! ************************************************************************************************* ! biogem_main.f90 ! BioGeM ! MAIN ! ************************************************************************************************* ! *** SETUP BioGeM *** SUBROUTINE setup_biogem(dum_t0,dum_filestring,dum_ans,dum_nyear, & & dum_maxi,dum_maxj,dum_maxk,dum_maxl,dum_imax,dum_jmax,dum_kmax,dum_lmax, & & dum_nsteps,dum_npstp,dum_iwstp,dum_itstp,dum_ianav,dum_pi,dum_saln0,dum_rhoair,dum_cd,dum_ds,dum_dphi, & & dum_usc,dum_rsc,dum_tsc,dum_dsc,dum_fsc,dum_gsc,dum_rh0sc,dum_rhosc,dum_cpsc,dum_solconst, & & dum_ips,dum_ipf,dum_ias,dum_iaf,dum_jsf, & & dum_k1,dum_tau,dum_dz,dum_dza, & & dum_c,dum_cv,dum_s,dum_sv, & & dum_ts, & & dum_ts1, & & dum_sfcatm1,dum_sfxatm1, & & dum_sfcocn1,dum_sfcsed1, & & dum_sfxocn1,dum_sfxsed1) USE biogem_lib USE biogem_data IMPLICIT NONE ! dummy arguments real,intent(in)::dum_t0 CHARACTER(LEN=*),INTENT(in)::dum_filestring CHARACTER(len=1),intent(in)::dum_ans INTEGER,INTENT(inout)::dum_nyear INTEGER,INTENT(in)::dum_maxi,dum_maxj,dum_maxk,dum_maxl INTEGER,INTENT(inout)::dum_imax,dum_jmax,dum_kmax,dum_lmax INTEGER,INTENT(inout)::dum_nsteps,dum_npstp,dum_iwstp,dum_itstp,dum_ianav REAL,INTENT(in)::dum_pi,dum_saln0,dum_rhoair,dum_cd,dum_ds,dum_dphi real,INTENT(in)::dum_usc,dum_rsc,dum_tsc,dum_dsc,dum_fsc,dum_gsc,dum_rh0sc,dum_rhosc,dum_cpsc,dum_solconst INTEGER,INTENT(in),DIMENSION(n_maxj)::dum_ips,dum_ipf,dum_ias,dum_iaf INTEGER,INTENT(in)::dum_jsf integer,DIMENSION(n_maxi,n_maxj),INTENT(in)::dum_k1 REAL,DIMENSION(2,n_maxi,n_maxj),INTENT(in)::dum_tau REAL,DIMENSION(n_maxk),INTENT(in)::dum_dz,dum_dza REAL,DIMENSION(0:n_maxj),INTENT(in)::dum_c,dum_cv,dum_s,dum_sv REAL,DIMENSION(n_maxl,n_maxi,n_maxj,n_maxk),INTENT(inout)::dum_ts,dum_ts1 REAL,DIMENSION(0:n_atm,n_maxi,n_maxj),INTENT(inout)::dum_sfcatm1 REAL,DIMENSION(0:n_atm,n_maxi,n_maxj),INTENT(inout)::dum_sfxatm1 REAL,DIMENSION(0:n_ocn,n_maxi,n_maxj),INTENT(inout)::dum_sfcocn1 REAL,DIMENSION(0:n_sed,n_maxi,n_maxj),INTENT(inout)::dum_sfcsed1 REAL,DIMENSION(0:n_sed,n_maxi,n_maxj),INTENT(inout)::dum_sfxocn1 REAL,DIMENSION(0:n_sed,n_maxi,n_maxj),INTENT(inout)::dum_sfxsed1 ! local variables real::loc_t ! ******************************************************* ! *** Synchronize common cgoldstein-biogem parameters *** ! ******************************************************* ! check primary array dimensions IF ((n_maxi /= dum_maxi) .OR. (n_maxj /= dum_maxj) .OR. (n_maxk /= dum_maxk) .OR. (n_maxl /= dum_maxl)) THEN CALL sub_report_error( & & 'biogem_main','setup_biogem','ocean grid dimensions do not match between '// & & 'BioGeM (n_maxl,n_maxi,n_maxj,n_maxk in file biogem_lib.f90; dir genie-biogem) and '// & & 'GOLDSTEIn (maxl,maxi,maxj,maxk in file var.cmn; dir genie-cgoldstein); '// & & 'make respective definitions consistent and re-compile model', & & 'STOPPING', & & (/const_real_null/),.TRUE. & & ) END IF ! move 'goin' file-read one line on if not a continuing run ! (GOLDSTEIn reads a different number of lines depending on whether it is a continuing run or not) IF ((dum_ans == 'n') .or. (dum_ans == 'N')) then read(unit=5,fmt=*) end if ! load ('goin') run-time options CALL sub_load_goin_biogem() ! set GeM time ! NOTE: modify 'par_misc_t_start' according to the run-time accumulated in any requested restart, ! so that the time that BioGeM starts with is the same as the actual start time requested in goin_biogem ! (BioGeM calculates its internal time as equal to par_misc_t_start + elapsed GOLDSTEIn time, ! BUT, elapsed GOLDSTEIn time will include any accumulated restart time, ! hence accumulated restart time is subtracted from par_misc_t_start) ! NOTE: I'm sure that this doesn't make any sense at all ... :( if (opt_misc(iopt_misc_t_timescale_BP)) then par_misc_t_end = par_misc_t_start - par_misc_t_runtime if ((dum_ans == 'c') .or. (dum_ans == 'C')) then par_misc_t_start = par_misc_t_start + (dum_tsc*dum_t0)/conv_yr_s end if else par_misc_t_end = par_misc_t_start + par_misc_t_runtime if ((dum_ans == 'c') .or. (dum_ans == 'C')) then par_misc_t_start = par_misc_t_start - (dum_tsc*dum_t0)/conv_yr_s end if end if ! over-write goldstein time-step control if requested IF (opt_misc(iopt_misc_t_timescale_BioGeM)) THEN dum_nsteps = int(real(dum_nyear)*par_misc_t_runtime) dum_iwstp = dum_nsteps - 1 dum_npstp = dum_nsteps + 1 dum_ianav = dum_nsteps + 1 END IF ! *** copy GOLDSTEIn parameters *** ! NOTE: this is so that keep GOLDSTEIn parameters are made 'global' to BioGeM ! (they are all defined in the GOLDSTEIn var.cmn common block, ! and cannot be seen unless explicitly passed and copied) ! grid dimensions n_imax = dum_imax n_jmax = dum_jmax n_kmax = dum_kmax n_lmax = dum_lmax ! NOTE: instead of using value of GODLSTEIn take default value from ! (since it can be reset later if too small) n_lmax = n_maxl ! copy time stepping and time control parameters goldstein_nsteps = dum_nsteps goldstein_npstp = dum_npstp goldstein_iwstp = dum_iwstp goldstein_itstp = dum_itstp goldstein_t0 = dum_t0 ! copy dimensional scale factors for ocean goldstein_usc = dum_usc goldstein_rsc = dum_rsc goldstein_tsc = dum_tsc goldstein_dsc = dum_dsc goldstein_fsc = dum_fsc goldstein_gsc = dum_gsc goldstein_rh0sc = dum_rh0sc goldstein_rhosc = dum_rhosc goldstein_cpsc = dum_cpsc ! copy miscellaneous constants ! NOTE: use identical value of pi to avoid possible GOLDETSIn/BioGeM mis-match in the calculation of ocean areas and volumes goldstein_pi = dum_pi goldstein_saln0 = dum_saln0 goldstein_rhoair = dum_rhoair goldstein_cd = dum_cd goldstein_ds = dum_ds goldstein_dphi = dum_dphi goldstein_solconst = dum_solconst ! copy ocean positions goldstein_jsf = dum_jsf goldstein_ips(:) = dum_ips(:) goldstein_ipf(:) = dum_ipf(:) goldstein_ias(:) = dum_ias(:) goldstein_iaf(:) = dum_iaf(:) ! copy ocean bottom index grid goldstein_k1(:,:) = dum_k1(:,:) ! copy wind stress field goldstein_tau(:,:,:) = dum_tau(:,:,:) ! miscellaneous goldstein_dz(:) = dum_dz(:) goldstein_dza(:) = dum_dza(:) goldstein_c(:) = dum_c(:) goldstein_cv(:) = dum_cv(:) goldstein_s(:) = dum_s(:) goldstein_sv(:) = dum_sv(:) ! *** define ts->ocn offset and copy array information *** ! define ofset between GOLDSTEIn tracer units and BioGeM ! => temperature is degrees C in GOLDSTEIn, but K in BioGeM ! => salinity is as a relative deviation (o/oo) (from saln0) in GOLDSTEIn, but absolute (o/oo) in BioGeM tstoocn_offset(:) = 0.0 tstoocn_offset(1) = const_zeroC tstoocn_offset(2) = goldstein_saln0 ! *** set-up biogem-sedgem and biogem-atchem interfacing *** ! initialize interface arrays dum_sfcatm1(:,:,:) = 0.0 dum_sfxatm1(:,:,:) = 0.0 dum_sfcocn1(:,:,:) = 0.0 dum_sfcocn1(:,:,:) = 0.0 dum_sfxocn1(:,:,:) = 0.0 dum_sfxsed1(:,:,:) = 0.0 ! define relationships between tracers CALL sub_def_tracerrelationships() ! ************************* ! *** initialize BioGeM *** ! ************************* ! set initial (BioGeM) model time loc_t = ABS(par_misc_t_end - par_misc_t_start) - (dum_tsc*dum_t0)/conv_yr_s ! initialize arrays par_force_flux_weather(:) = 0.0 CALL sub_init_int_timeslice() CALL sub_init_int_timeseries() CALL sub_init_force() CALL sub_init_bio() IF (opt_misc(iopt_misc_audit)) CALL sub_init_audit() ! initialize ocean and ocean-atmosphere interfacing grids and 'physics' CALL sub_init_phys_ocn() CALL sub_init_phys_ocnatm() ! load default values for ocean, atmosphere, and sediment tracers and and initialize tracer arrays and options CALL sub_init_tracer_ocn() call sub_biogem_copy_tstoocn(dum_ts) CALL sub_init_tracer_ocnatm() CALL sub_init_tracer_ocnsed() ! load primary BioGeM run-time options CALL sub_load_biogem_config() ! load bio options (biological scheme and details of remineralization) call sub_load_biogem_bio_config() ! set meta-options and verify self-consistency of chosen parameters call sub_check_par() ! initialize default 'Redfield' ratios call sub_init_bio_Redfield ! initialize carbon cycle sub-systems (if selected) IF (opt_select(iopt_select_carbchem)) CALL sub_init_carb() ! initialize the time scale save string CALL sub_init_char() ! initialize basic data saving CALL sub_init_data_save() ! open units for runtime file I/O CALL sub_init_data_save_runtime(loc_t) ! initialize system forcings CALL sub_init_force_restore_atm() CALL sub_init_force_flux_atm() CALL sub_init_force_restore_ocn() CALL sub_init_force_flux_ocn() !!$ CALL sub_init_force_restore_sed() CALL sub_init_force_flux_sed() ! *** set the value of in GOLDSTEIn *** ! make copy of initial GOLDSTEIn lmax goldstein_lmax = dum_lmax ! re-set GOLDSTEIn lmax dum_lmax = n_lmax ! *** load restart information *** IF ((dum_ans == 'c') .or. (dum_ans == 'C')) then call sub_load_biogem_restart(trim(dum_filestring)) end if ! *** initialize tracer auditing *** ! carry out initial tracer inventory audit IF (opt_misc(iopt_misc_audit)) THEN audit_ocn_init(:) = fun_calc_ocn_tot() audit_ocn_old(:) = audit_ocn_init(:) end if ! *** write ocean tracer data to GOLDSTEIn arrays *** call sub_biogem_copy_ocntots(dum_ts,dum_ts1) ! *** enable BioGeM *** par_misc_t_go = .TRUE. END SUBROUTINE setup_biogem ! *** TIMESTEP BioGeM *** SUBROUTINE tstep_biogem(dum_t, & & dum_t0, & & dum_dt, & & dum_ts, & & dum_ts1, & & dum_varice, & & dum_cost, & & dum_solfor, & & dum_u, & & dum_sfcatm1,dum_sfxatm1, & & dum_sfcocn1,dum_sfcsed1, & & dum_sfxocn1,dum_sfxsed1) USE biogem_lib USE biogem_box USE biogem_data IMPLICIT NONE ! dummy arguments REAL,INTENT(in)::dum_t,dum_t0,dum_dt REAL,DIMENSION(n_maxl,n_maxi,n_maxj,n_maxk),INTENT(inout)::dum_ts REAL,DIMENSION(n_maxl,n_maxi,n_maxj,n_maxk),INTENT(inout)::dum_ts1 REAL,DIMENSION(n_maxi,n_maxj),INTENT(in)::dum_varice REAL,DIMENSION(n_maxi,n_maxj),INTENT(in)::dum_cost REAL,DIMENSION(n_maxj),INTENT(in)::dum_solfor REAL,DIMENSION(3,n_maxi,n_maxj,n_maxk),INTENT(in)::dum_u REAL,DIMENSION(0:n_atm,n_maxi,n_maxj),INTENT(in)::dum_sfcatm1 REAL,DIMENSION(0:n_atm,n_maxi,n_maxj),INTENT(inout)::dum_sfxatm1 REAL,DIMENSION(0:n_ocn,n_maxi,n_maxj),INTENT(inout)::dum_sfcocn1 REAL,DIMENSION(0:n_sed,n_maxi,n_maxj),INTENT(in)::dum_sfcsed1 REAL,DIMENSION(0:n_sed,n_maxi,n_maxj),INTENT(in)::dum_sfxocn1 REAL,DIMENSION(0:n_sed,n_maxi,n_maxj),INTENT(inout)::dum_sfxsed1 ! local variables INTEGER::i,j,k,l,io,ia,is integer::loc_k,loc_k1 integer::loc_k_min,loc_k_max real::loc_min,loc_max real::loc_t,loc_dts,loc_dtyr,loc_yr ! local time and time step BLAH actual year logical,DIMENSION(0:n_ocn)::locio_mask REAL,DIMENSION(0:n_ocn,n_maxi,n_maxj,n_maxk)::locijk_ocn REAL,DIMENSION(0:n_ocn,n_maxi,n_maxj,n_maxk)::locijk_focn REAL,DIMENSION(0:n_atm,n_maxi,n_maxj)::locij_fatm REAL,DIMENSION(0:n_ocn,n_maxk)::loc_docn_restore REAL,DIMENSION(0:n_atm)::loc_datm_restore !!$ REAL,DIMENSION(0:n_sed)::loc_dsed_restore REAL,DIMENSION(0:n_ocn)::loc_force_flux_weather real::loc_ocn_tot_M,loc_ocn_tot_A,loc_ocnatm_tot_A,loc_ocn_mean_S REAL,DIMENSION(0:n_ocn)::loc_ocn_tot_OLD,loc_ocn_tot_NEW REAL,DIMENSION(0:n_ocn)::loc_force_restore_ocn_tmod ! relaxation modifier REAL,DIMENSION(0:n_atm)::loc_force_restore_atm_tmod ! relaxation modifier !!$ REAL,DIMENSION(0:n_sed)::loc_force_restore_sed_tmod ! relaxation modifier REAL,DIMENSION(0:n_atm,n_maxi,n_maxj)::locij_focnatm ! local ocn->atm flux (atm tracer currency), units (mol yr-1) REAL,DIMENSION(0:n_sed,n_maxi,n_maxj)::locij_focnsed ! local ocn->sed change (sed tracer currency), units (mol) REAL,DIMENSION(0:n_ocn,n_maxi,n_maxj)::locij_fsedocn ! local sed->ocean change (ocn tracer currency), units (mol) REAL,DIMENSION(0:n_ocn)::loc_ocnsed_audit ! temporary ocn->sed auditing array (ocn tracer currency), units (mol) REAL,DIMENSION(0:n_maxj,0:n_maxk)::loc_opsi,loc_opsia,loc_opsip,loc_zpsi real::loc_opsi_scale ! *** INITIALIZE LOCAL ARRAYS *** ! NOTE: only bother initializing for selected tracers DO ia=1,n_atm IF (ocnatm_select(ia)) THEN locij_fatm(:,:,:) = 0.0 locij_focnatm(:,:,:) = 0.0 end if end do DO io=1,n_ocn IF (ocn_select(io)) THEN locijk_ocn(io,:,:,:) = 0.0 locijk_focn(io,:,:,:) = 0.0 locij_fsedocn(io,:,:) = 0.0 end if end do DO is=1,n_sed IF (ocnsed_select(is)) THEN locij_focnsed(is,:,:) = 0.0 end if end do ! *** CALCULATE GEM TIME *** ! update gemchemical model time ! NOTE: counted down in years ! calculate time step length ! NOTE: convert between time units of BioGeM (years) and GOLDSTEIn (use scale factor to convert to seconds) ! NOTE: subtract time par_misc_t_err (yrs) off local time to account for finite numerical inaccurary ! in making a comparison between loc_t with a prescribed point in time loc_t = ABS(par_misc_t_end - par_misc_t_start) - (goldstein_tsc*dum_t)/conv_yr_s - par_misc_t_err loc_dts = goldstein_tsc*dum_dt loc_dtyr = loc_dts/conv_yr_s IF (opt_misc(iopt_misc_t_timescale_BP)) THEN loc_yr = loc_t + par_misc_t_end ELSE loc_yr = par_misc_t_end - loc_t END IF ! *** UPDATE PHYSICS *** DO i=1,n_imax DO j=1,n_jmax loc_k1 = goldstein_k1(i,j) IF (n_kmax >= loc_k1) THEN phys_ocnatm(ipoa_solfor,i,j) = dum_solfor(j) phys_ocn(ipo_dcost,i,j,n_kmax) = dum_cost(i,j) - phys_ocn(ipo_cost,i,j,n_kmax) phys_ocn(ipo_cost,i,j,n_kmax) = dum_cost(i,j) phys_ocnatm(ipoa_seaice,i,j) = dum_varice(i,j) do k=loc_k1,n_kmax phys_ocn(ipo_rho,i,j,k) = fun_calc_rho(ocn(io_T,i,j,k),ocn(io_S,i,j,k)) end do end IF end DO end DO ! *** CALCULATE LOCAL CONSTANTS *** ! fractional relaxation DO io=1,n_ocn IF (ocn_select(io)) THEN loc_force_restore_ocn_tmod(io) = 1.0 - EXP(-loc_dtyr/force_restore_ocn_tconst(io)) end if end do DO ia=1,n_atm IF (ocnatm_select(ia)) THEN loc_force_restore_atm_tmod(ia) = 1.0 - EXP(-loc_dtyr/force_restore_atm_tconst(ia)) end if end do !!$ DO is=1,n_sed !!$ IF (ocnsed_select(is)) THEN !!$ loc_force_restore_sed_tmod(is) = 1.0 - EXP(-loc_dtyr/force_restore_sed_tconst(is)) !!$ end if !!$ end do ! total ocean mass loc_ocn_tot_M = sum(phys_ocn(ipo_M,:,:,:)) ! total global surface area loc_ocnatm_tot_A = sum(phys_ocnatm(ipoa_A,:,:)) ! ocean surface area (ice-free) loc_ocn_tot_A = sum((1.0 - phys_ocnatm(ipoa_seaice,:,:))*phys_ocn(ipo_A,:,:,n_kmax)) ! local global fixed flux forcing ! NOTE: units of are (mol yr-1) ! NOTE: convert to dissolved tracer currency and into units of (mol yr-1) loc_force_flux_weather(:) = matmul(conv_sed_ocn(:,:),par_force_flux_weather(:)) ! calculate local opsi conversion constant loc_opsi_scale = goldstein_dsc*goldstein_usc*goldstein_rsc*1.0E-6 ! *** START-UP REPORTING *** if ((dum_t - dum_t0 - conv_yr_s*par_misc_t_err/goldstein_tsc) <= dum_dt) then print*,' ' print*,'*** BioGeM run-time daignostics ***' CALL sub_calc_psi(dum_u,loc_opsi,loc_opsia,loc_opsip,loc_zpsi) par_misc_t_echo_header = .TRUE. call sub_echo_runtime(loc_yr,loc_ocn_tot_M,loc_ocn_tot_A,loc_ocnatm_tot_A,loc_opsi_scale,loc_opsia(:,:),dum_sfcatm1(:,:,:)) print*,' ' END IF ! ****************************** ! *** UPDATE BIOGEOCHEMISTRY *** ! ****************************** ! main BioGeM loop - all updating of ocean-atmosphere-sediment biogeochemistry occurs within this ! NOTE: update conditional on end of biogeochemical simulation not having been reached ! NOTE: code to excecute every iterations (i.e., at every GOLDSTEIn time step) unless otherwise stated ! ocean biogeochemistry is updated on a grid point by grid point ((i,j) loop) basis ! NOTE: tracer loops follow the GOLDSTEIn notation, with indices converted to their BioGeM equivalents ! NOTE: should there be tracer provision in GOLDSTEIn not used in BioGeM, ! i.e, 'spare' array indices arrising when the number of tracers selected in BioGeM is less than ! the correspsonding compiled dimensions in GOLDSTEIn, ! GOLDSTEIn indices are mapped onto the '0' index in BioGeM (essentially a null tracer) to avoid array overflow ! (and thus memory overwriting) IF (par_misc_t_go) THEN ! *** UPDATE FORCING TIME SERIES DATA *** ! recalculate time-varying restoring and flux forcings of the system ! ATMOSPHERE (applied at the ocean-atmospere interface) DO ia=1,n_atm IF (ocnatm_select(ia)) THEN IF (force_restore_atm_select(ia)) THEN CALL sub_update_force_restore_atm(loc_t,ia) END IF IF (force_flux_atm_select(ia)) THEN CALL sub_update_force_flux_atm(loc_t,ia) END IF END IF END DO ! OCEAN DO io=1,n_ocn IF (ocn_select(io)) THEN IF (force_restore_ocn_select(io)) THEN CALL sub_update_force_restore_ocn(loc_t,io) END IF IF (force_flux_ocn_select(io)) THEN CALL sub_update_force_flux_ocn(loc_t,io) END IF end if END DO ! SEDIMENTS (applied at the ocean surface) DO is=1,n_sed IF (ocnsed_select(is)) THEN !!$ IF (force_restore_sed_select(is)) THEN !!$ CALL sub_update_force_restore_sed(loc_t,is) !!$ END IF IF (force_flux_sed_select(is)) THEN CALL sub_update_force_flux_sed(loc_t,is) END IF END IF END DO ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! *** (i,j) GRID PT LOOP START *** ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DO i=1,n_imax DO j=1,n_jmax ! *** INITIALIZE LOOP VARIABLES *** ! set local depth loop limit loc_k1 = goldstein_k1(i,j) ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! *** WET GRID PT CONDITIONALITY START *** ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF (n_kmax >= loc_k1) THEN ! *** UPDATE OCEAN AQUEOUS CARBONATE SYSTEM *** IF (opt_select(iopt_select_carbchem)) THEN ! re-calculate carbonate constants ! NOTE: only re-calculate is T or S has changed significantly ! (currently the threshold for both is set to 0.1) ! NOTE: saving of cumulative update frequency is called from end_biogem (file; 'gemdiag_carb_n.res') if ( & & (ABS(carb_TSn(1,i,j,n_kmax) - ocn(io_T,i,j,n_kmax)) > 0.1) & & .OR. & & (ABS(carb_TSn(2,i,j,n_kmax) - ocn(io_S,i,j,n_kmax)) > 0.1) & & ) then CALL sub_calc_carbconst( & & phys_ocn(ipo_Dmid,i,j,n_kmax), & & ocn(io_T,i,j,n_kmax), & & ocn(io_S,i,j,n_kmax), & & carbconst(:,i,j,n_kmax) & & ) carb_TSn(1,i,j,n_kmax) = ocn(io_T,i,j,n_kmax) carb_TSn(2,i,j,n_kmax) = ocn(io_S,i,j,n_kmax) carb_TSn(3,i,j,n_kmax) = carb_TSn(3,i,j,n_kmax) + 1.0 end if ! re-estimate Ca and borate concentrations from salinity (if not selected and therefore explicitly treated) IF (.NOT. ocn_select(io_Ca)) ocn(io_Ca,i,j,n_kmax) = fun_calc_Ca(ocn(io_S,i,j,n_kmax)) IF (.NOT. ocn_select(io_B)) ocn(io_B,i,j,n_kmax) = fun_calc_B(ocn(io_S,i,j,n_kmax)) ! re-calculate surface ocean carbonate chemistry CALL sub_calc_carb( & & ocn(io_CO2,i,j,n_kmax), & & ocn(io_ALK,i,j,n_kmax), & & ocn(io_PO4,i,j,n_kmax), & & ocn(io_Si,i,j,n_kmax), & & ocn(io_B,i,j,n_kmax), & & ocn(io_Ca,i,j,n_kmax), & & carbconst(:,i,j,n_kmax), & & carb(:,i,j,n_kmax) & & ) END IF ! *** CALCULATE RESTORING BOUNDARY CONDITIONS *** ! ATMOSPHERE ! calculate tracer changes necessary to meet any imposed atmospheric restoring boundary conditions ! NOTE: do not apply restoring forcing when the two indices of the forcing array are identical, ! as this typically indicates that the current time is outside of the interval ! defined in the forcing time series input file ! -> this allows forcings to be 'turned' off and on in a time-depentent manner (although only once) ! NOTE: divide restoring atmospheric flux by to produce a value in units of mol a-1, ! (since in updating the atmosphere subsequently, the flux is multiplied by the timestep ! to get the required increment) ! NOTE: use a crude conversion factor for partial pressure (atm) -> mol, BUT ! take into account that only 1/(imax x jmax) of the atmosphere is being forced under each grid point ! NOTE: catch missing restoring data (tracer values < 0.0) => force restorign flux to zero DO ia=1,n_atm IF (ocnatm_select(ia)) THEN IF (force_restore_atm_select(ia) .AND. (force_restore_atm_sig_i(ia,1) /= force_restore_atm_sig_i(ia,2))) THEN if (force_restore_atm(ia,i,j) < 0.0) force_restore_atm(ia,i,j) = dum_sfcatm1(ia,i,j) loc_datm_restore(ia) = (force_restore_atm(ia,i,j) - dum_sfcatm1(ia,i,j))*loc_force_restore_atm_tmod(ia) locij_fatm(ia,i,j) = (1.0/real(n_imax*n_jmax))*conv_atm_mol*loc_datm_restore(ia)/loc_dtyr END IF end if END DO ! OCEAN ! hack to overwrite surface ocean [O2] restoring field and equilibrate with atmosphere if requested IF (opt_misc(iopt_misc_O2_equil)) THEN force_restore_ocn(io_O2,i,j,n_kmax) = carbconst(icc_QO2,i,j,n_kmax)*dum_sfcatm1(ia_pO2,i,j) end if ! calculate tracer changes necessary to meet any imposed oceanic restoring boundary conditions ! NOTE: do not take into account fractional ice-free area ! NOTE: catch missing restoring data (tracer values < 0.0) => force restorign flux to zero DO io=1,n_ocn IF (ocn_select(io)) THEN IF (force_restore_ocn_select(io) .AND. (force_restore_ocn_sig_i(io,1) /= force_restore_ocn_sig_i(io,2))) THEN DO k=force_restore_ocn_k1(io,i,j),n_kmax if (force_restore_ocn(io,i,j,k) < 0.0) force_restore_ocn(io,i,j,k) = ocn(io,i,j,k) loc_docn_restore(io,k) = (force_restore_ocn(io,i,j,k) - ocn(io,i,j,k))*loc_force_restore_ocn_tmod(io) locijk_focn(io,i,j,k) = loc_docn_restore(io,k)*phys_ocn(ipo_M,i,j,k)/loc_dtyr END DO END IF end if END DO !!$ ! SEDIMENTS !!$ DO is=1,n_sed !!$ IF (ocnsed_select(is)) THEN !!$ IF (force_restore_sed_select(is) .AND. (force_restore_sed_sig_i(is,1) /= force_restore_sed_sig_i(is,2))) THEN !!$ if (force_restore_sed(is,i,j) < 0.0) force_restore_sed(is,i,j) = dum_sfcsed1(is,i,j) !!$ loc_dsed_restore(is) = (force_restore_sed(is,i,j) - dum_sfcsed1(is,i,j))*loc_force_restore_sed_tmod(is) !!$ locij_fsed(is,i,j) = (1.0/real(n_imax*n_jmax))*conv_sed_mol*loc_dsed_restore(is)/loc_dtyr !!$ END IF !!$ end if !!$ END DO ! *** CH4/CLATHRATES *** ! deal with clathrate destabalization and CH4 bubbling ! ********************* ! *** *** ! ********************* ! *** WATER COLUMN REMINERALIZATION *** ! NOTE: remineralize particulate material left over from PREVIOUS time step, ! i.e., do not immediately remineralize material newly created in the current time step ! (which is why the subroutine is called before ) ! NOTE: units of (mol kg-1) call sub_calc_bio_remin(i,j,loc_dtyr) ! *** SURFACE OCEAN BIOLOGICAL PRODUCTIVITY *** ! NOTE: units of (mol kg-1) call sub_calc_bio_uptake(i,j,loc_dtyr,loc_docn_restore(:,n_kmax),locijk_focn(:,i,j,n_kmax)) ! *** OCEAN-ATMOPSHERE EXCHANGE FLUXES *** ! calculate ocean-atmosphere exchange (restricted to wet grid points) ! NOTE: a positive value represents net ocean to atmosphere transfer ! NOTE: units of (mol yr-1) ! ******************************************* ! *** *** ! ******************************************* ! \/\/\/ SPECIAL CASE FOR CO2 AND C ISOTOPES \/\/\/ ! calculate surface ocean CO2 piston velocity and air-sea exchange CALL sub_calc_pv_CO2(i,j) ! calculate CO2 air-sea exchange locij_focnatm(ia_pCO2,i,j) = fun_calc_ocnatm_flux_CO2(i,j,dum_sfcatm1(ia_pCO2,i,j)) ! /\/\/\ SPECIAL CASE FOR CO2 AND C ISOTOPES /\/\/\ ! set update of ocean and atmosphere reservoirs ! >>> GENERIC ALGORITHM BUT CONTAINS USER-DEFINABLE INFORMATION IN ARRAY DO ia=1,n_atm IF (ocnatm_select(ia) .AND. (ocnatm_type(ia) /= 0)) THEN locij_fatm(ia,i,j) = locij_fatm(ia,i,j) + locij_focnatm(ia,i,j) locio_mask(:) = .TRUE. do io = maxval(maxloc(abs(conv_atm_ocn(:,ia)),locio_mask(:)))-1 if (io == 0) then exit else locijk_focn(io,i,j,n_kmax) = locijk_focn(io,i,j,n_kmax) - conv_atm_ocn(io,ia)*locij_focnatm(ia,i,j) locio_mask(io) = .FALSE. end if end do end IF end DO ! *** EXTERNAL FLUX FORCING *** ! calculate value of applied flux forcings ! NOTE: in units of (mol yr-1) ! ATMOSPHERE (applied at the ocean-atmospere interface) DO ia=1,n_atm IF (ocnatm_select(ia)) THEN IF (force_flux_atm_select(ia)) THEN locij_fatm(ia,i,j) = locij_fatm(ia,i,j) + force_flux_atm(ia,i,j) END IF end if END DO ! OCEAN DO io=1,n_ocn IF (ocn_select(io)) THEN IF (force_flux_ocn_select(io)) THEN DO k=loc_k1,n_kmax locijk_focn(io,i,j,k) = locijk_focn(io,i,j,k) + force_flux_ocn(io,i,j,k) END DO END IF ! add prescribed total system flux forcing (aka 'weathering') ! NOTE: units of (mol yr-1) ! NOTE: distribute total uniformly throughout ocean according to fractional ocean cell mass locijk_focn(io,i,j,loc_k1:n_kmax) = locijk_focn(io,i,j,loc_k1:n_kmax) + & & loc_force_flux_weather(io)*phys_ocn(ipo_M,i,j,loc_k1:n_kmax)/loc_ocn_tot_M end if END DO ! SEDIMENTS (applied at the ocean surface) ! NOTE: addition is made directly to particulate sedimentary tracer array (scaled by time-step and cell mass) DO is=1,n_sed IF (ocnsed_select(is)) THEN IF (force_flux_sed_select(is)) THEN bio_part(is,i,j,n_kmax) = bio_part(is,i,j,n_kmax) + & & loc_dtyr*phys_ocn(ipo_rM,i,j,n_kmax)*force_flux_sed(is,i,j) END IF end if END DO ! *** SEDIMENT DISSOLUTION INPUT *** ! modify remineralization array according to dissolution input from sediments ! NOTE: in units of (mol m-2 s-1) - needs to be converted to (mol per time step) ! NOTE: if the model is configured as a 'closed' system (biogem_congif.par), ! set dissolution flux equal to rain flux and add to overlying ocean If (opt_misc(iopt_misc_sed_select)) then locij_fsedocn(:,i,j) = 0.0 DO is=1,n_sed IF (ocnsed_select(is)) THEN locio_mask(:) = .TRUE. do io = maxval(maxloc(abs(conv_sed_ocn(:,is)),locio_mask(:)))-1 if (io == 0) then exit else if (opt_misc(iopt_misc_sed_closedsystem)) then locij_fsedocn(io,i,j) = locij_fsedocn(io,i,j) + & & conv_sed_ocn(io,is)*bio_settle(is,i,j,loc_k1) else locij_fsedocn(io,i,j) = locij_fsedocn(io,i,j) + & & loc_dts*phys_ocn(ipo_A,i,j,loc_k1)*conv_sed_ocn(io,is)*dum_sfxocn1(is,i,j) end if locio_mask(io) = .FALSE. end if end do end IF end DO DO io=1,n_ocn IF (ocn_select(io)) THEN bio_remin(io,i,j,loc_k1) = bio_remin(io,i,j,loc_k1) + & & phys_ocn(ipo_rM,i,j,loc_k1)*locij_fsedocn(io,i,j) end if end do end if end if ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! *** WET GRID PT CONDITIONALITY END *** ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< END DO END DO ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! *** (i,j) GRID PT LOOP END *** ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! *** OCEAN TRACER UPDATE *** ! *** [SALINITY NORMALIZED SCHEME] *** ! NOTE: units of aqueous concentration (mol kg-1) ! NOTE: loop though only tracers actually selected (i.e., 3:lmax) ! calculate mean ocean salinity loc_ocn_mean_S = SUM(ocn(io_S,:,:,:)*phys_ocn(ipo_M,:,:,:))/loc_ocn_tot_M ! calculate unadjusted biogeochem tracer inventory DO l=3,n_lmax io = goldstein_l_io(l) loc_ocn_tot_OLD(io) = SUM(ocn(io,:,:,:)*phys_ocn(ipo_M,:,:,:)) end DO ! (1) salinity-adjust circulation-updated GOLDSTEIn tracer array ! NOTE: no offset (array: ) required for biogeochem-only tracers DO i=1,n_imax DO j=1,n_jmax DO k=goldstein_k1(i,j),n_kmax DO l=3,n_lmax io = goldstein_l_io(l) locijk_ocn(io,i,j,k) = dum_ts(l,i,j,k)*(ocn(io_S,i,j,k)/loc_ocn_mean_S) end DO end DO end do end DO ! calculate salinity-adjusted biogeochem tracer inventory ! NOTE: be sure to avoid potential divide-by-zero problems DO l=3,n_lmax io = goldstein_l_io(l) loc_ocn_tot_NEW(io) = SUM(locijk_ocn(io,:,:,:)*phys_ocn(ipo_M,:,:,:)) if (loc_ocn_tot_NEW(io) < const_real_nullsmall) loc_ocn_tot_NEW(io) = const_real_nullhigh end do ! (2) update T,S ! NOTE: includes hack for T,S boundary condition forcing in biogeochemical model only (i.e., not see by GOLDSTEIN) DO i=1,n_imax DO j=1,n_jmax DO l=1,2 io = goldstein_l_io(l) DO k=goldstein_k1(i,j),n_kmax locijk_ocn(io,i,j,k) = dum_ts(l,i,j,k) + tstoocn_offset(l) If (opt_force(iopt_force_GOLDSTEInTS)) then ocn(io,i,j,k) = locijk_ocn(io,i,j,k) + loc_dtyr*phys_ocn(ipo_rM,i,j,k)*locijk_focn(io,i,j,k) else ocn(io,i,j,k) = locijk_ocn(io,i,j,k) end if dum_ts(l,i,j,k) = ocn(io,i,j,k) - tstoocn_offset(l) dum_ts1(l,i,j,k) = dum_ts(l,i,j,k) end DO end DO end DO end DO ! re-calculate mean ocean salinity loc_ocn_mean_S = SUM(ocn(io_S,:,:,:)*phys_ocn(ipo_M,:,:,:))/loc_ocn_tot_M ! (3a) normalize salinity-adjusted tracer concentrations and update biogeochem tracers ! NOTE: this is to preserve mass-balance - the final net effect is a re-distribution of tracer concentrations ! according to the relative salinity deviation from the mean but normalized so that the overall scale factor is unity ! (3b) salinity normalize biogeochem tracers and copy to DO i=1,n_imax DO j=1,n_jmax DO k=goldstein_k1(i,j),n_kmax DO l=3,n_lmax io = goldstein_l_io(l) ocn(io,i,j,k) = locijk_ocn(io,i,j,k)*(loc_ocn_tot_OLD(io)/loc_ocn_tot_NEW(io)) + & & bio_remin(io,i,j,k) + loc_dtyr*phys_ocn(ipo_rM,i,j,k)*locijk_focn(io,i,j,k) dum_ts(l,i,j,k) = ocn(io,i,j,k)*(loc_ocn_mean_S/ocn(io_S,i,j,k)) dum_ts1(l,i,j,k) = dum_ts(l,i,j,k) end do end DO end do end DO ! *** INTERFACE ARRAY UPDATE *** DO i=1,n_imax DO j=1,n_jmax loc_k1 = goldstein_k1(i,j) IF (n_kmax >= loc_k1) THEN ! (1) set ocn->atm flux ! NOTE: convert units from (mol yr-1) to (mol m-2 s-1) ! NOTE: update external interface array with TOTAL flux to atmosphere ! (i.e., due to both air-sea exchange and any forcing of the atmosphere) DO ia=1,n_atm IF (ocnatm_select(ia)) THEN dum_sfxatm1(ia,i,j) = phys_ocnatm(ipoa_rA,i,j)*conv_s_yr*locij_fatm(ia,i,j) end if end do ! (2) set bottom-water tracers ! NOTE: estimate Ca and borate concentrations from salinity (if not selected and therefore not explicitly treated) do io=1,n_ocn IF (ocn_select(io)) THEN dum_sfcocn1(io,i,j) = ocn(io,i,j,loc_k1) end if end do IF (.NOT. ocn_select(io_Ca)) dum_sfcocn1(io_Ca,i,j) = fun_calc_Ca(dum_sfcocn1(io_S,i,j)) IF (.NOT. ocn_select(io_B)) dum_sfcocn1(io_B,i,j) = fun_calc_B(dum_sfcocn1(io_S,i,j)) ! (3) set ocn->sed flux ! NOTE: convert units from (mol per timestep) to (mol m-2 s-1) ! NOTE: if 'allow particulate flux to sediments' option in biogem_config is not selected, ! the value of the particulate flux to sediments local array is zero do is=1,n_sed IF (ocnsed_select(is)) THEN If (opt_misc(iopt_misc_sed_select)) then locij_focnsed(is,i,j) = bio_settle(is,i,j,loc_k1) else locij_focnsed(is,i,j) = 0.0 end if dum_sfxsed1(is,i,j) = phys_ocn(ipo_rA,i,j,loc_k1)*locij_focnsed(is,i,j)/loc_dts end if end do ! (4) set CaCO3 age tracer if (ocnsed_select(is_CaCO3_age)) then dum_sfxsed1(is_CaCO3_age,i,j) = loc_t*dum_sfxsed1(is_CaCO3,i,j) end if end if END DO END DO ! *** AUDIT FLUX UPDATE & TRACER PLAUSIBILITY CHECK *** IF (opt_misc(iopt_misc_audit)) THEN DO i=1,n_imax DO j=1,n_jmax loc_k1 = goldstein_k1(i,j) IF (n_kmax >= loc_k1) THEN ! (1) Are we living in Lah-lah Land? (the TellyTubby check) ! i.e., are the new ocean tracer values calculated within 'pluasible' limits; ! this is useful pre-screening to catch potential numerical instability problems and ! situations likely to trip up the implicit solution of aqueous carbonate chemistry ! NOTE: carry out check only if tracer auditing has been enabled DO io=3,n_ocn IF (ocn_select(io)) THEN IF (minval(ocn(io,i,j,loc_k1:n_kmax)) < par_misc_ocn_min(io)) then ! find minimum value loc_min = par_misc_ocn_min(io) do loc_k=loc_k1,n_kmax if (ocn(io,i,j,loc_k) < loc_min) then loc_min = ocn(io,i,j,loc_k) loc_k_min = loc_k end if end do ! report details CALL sub_report_error( & & 'biogem_main','tstep_biogem', & & 'Whoops; things not looking so hot - tracer values < par_misc_ocn_min(io) '// & & '[tracer #ID, , , , par_misc_ocn_min(io), value]', & & 'OH DEAR - I AM NOT FEELING SO WELL ANYMORE ...', & & (/real(io),real(i),real(j),real(loc_k_min),par_misc_ocn_min(io),ocn(io,i,j,loc_k_min)/), & & .false. & & ) ! set run termination flag par_misc_t_go = .FALSE. end if IF (maxval(ocn(io,i,j,loc_k1:n_kmax)) > par_misc_ocn_max(io)) then ! find maximum value loc_max = par_misc_ocn_max(io) do loc_k=loc_k1,n_kmax if (ocn(io,i,j,loc_k) > loc_max) then loc_max = ocn(io,i,j,loc_k) loc_k_max = loc_k end if end do ! report details CALL sub_report_error( & & 'biogem_main','tstep_biogem', & & 'Whoops; things not looking so hot - tracer values > par_misc_ocn_max(io) '// & & '[tracer #ID, , , , par_misc_ocn_max(io), value]', & & 'OH DEAR - I AM NOT FEELING SO WELL ANYMORE ...', & & (/real(io),real(i),real(j),real(loc_k_max),par_misc_ocn_max(io),ocn(io,i,j,loc_k_max)/), & & .false. & & ) ! set run termination flag par_misc_t_go = .FALSE. end if if (.NOT. par_misc_t_go) then ! set for full data save opt_data(:) = .TRUE. ! initialize time-slice CALL sub_init_int_timeslice() par_data_save_timeslice_i = 1 par_data_save_timeslice(par_data_save_timeslice_i) = loc_t + loc_dtyr end if end if end do ! (2) update net integrated fluxes into the ocean for tracer inventory auditing (if selected) ! NOTE: take into account any sediment tracer flux forcing; ! the is implemented by incrementing the array, and wont show up in the auditing otherwise ! >>> GENERIC ALGORITHM BUT CONTAINS USER-DEFINABLE INFORMATION IN ARRAY loc_ocnsed_audit(:) = -locij_fsedocn(:,i,j) DO is=1,n_sed IF (ocnsed_select(is)) THEN locio_mask(:) = .TRUE. do io = maxval(maxloc(abs(conv_sed_ocn(:,is)),locio_mask(:)))-1 if (io == 0) then exit else If (opt_misc(iopt_misc_sed_select)) then loc_ocnsed_audit(io) = loc_ocnsed_audit(io) + conv_sed_ocn(io,is)*locij_focnsed(is,i,j) else loc_ocnsed_audit(io) = 0.0 end if locijk_focn(io,i,j,n_kmax) = locijk_focn(io,i,j,n_kmax) + conv_sed_ocn(io,is)*force_flux_sed(is,i,j) locio_mask(io) = .FALSE. end if end do end IF END DO DO io=1,n_ocn IF (ocn_select(io)) THEN audit_ocn_delta(io) = audit_ocn_delta(io) + & & loc_dtyr*SUM(locijk_focn(io,i,j,loc_k1:n_kmax)) - loc_ocnsed_audit(io) end if END DO end if END DO END DO end if ! *** TIME-SLICE DATA UPDATE *** ! update time slice data ! NOTE: carried out only when the local (BioGeM) time falls between a selected time slice time plus integration time, ! and the time slice time itself ! NOTE: do not construct and save a time slice if == 0, ! i.e., no valid (in-range) time slices have been requested in the time slice input file 'biogem_save_timeslice.dat', ! or the final requested time slice has been made IF ( & & ((par_data_save_timeslice(par_data_save_timeslice_i) + par_data_save_timeslice_dt) > (loc_t + loc_dtyr)) & & .AND. & & (par_data_save_timeslice_i > 0) & & ) THEN ! update integrated time int_t_timeslice = int_t_timeslice + loc_dtyr int_t_timeslice_count = int_t_timeslice_count + 1 ! update whole-ocean carbonate equilibrium DO i=1,n_imax DO j=1,n_jmax DO k=goldstein_k1(i,j),n_kmax CALL sub_calc_carb( & & ocn(io_CO2,i,j,k), & & ocn(io_ALK,i,j,k), & & ocn(io_PO4,i,j,k), & & ocn(io_Si,i,j,k), & & ocn(io_B,i,j,k), & & ocn(io_Ca,i,j,k), & & carbconst(:,i,j,k), & & carb(:,i,j,k) & & ) end do end do end do ! update time slice data - ocean ! NOTE: do not time increment weight quantities such as or , ! because they represent discrete increments rather than a flux or weightable concentration value int_ocn_timeslice(:,:,:,:) = int_ocn_timeslice(:,:,:,:) + loc_dtyr*ocn(:,:,:,:) int_bio_part_timeslice(:,:,:,:) = int_bio_part_timeslice(:,:,:,:) + loc_dtyr*bio_part(:,:,:,:) int_bio_settle_timeslice(:,:,:,:) = int_bio_settle_timeslice(:,:,:,:) + bio_settle(:,:,:,:) int_bio_remin_timeslice(:,:,:,:) = int_bio_remin_timeslice(:,:,:,:) + bio_remin(:,:,:,:) int_phys_ocn_timeslice(:,:,:,:) = int_phys_ocn_timeslice(:,:,:,:) + loc_dtyr*phys_ocn(:,:,:,:) int_carb_timeslice(:,:,:,:) = int_carb_timeslice(:,:,:,:) + loc_dtyr*carb(:,:,:,:) int_carbconst_timeslice(:,:,:,:) = int_carbconst_timeslice(:,:,:,:) + loc_dtyr*carbconst(:,:,:,:) ! update time slice data - ocean-atmosphere interface int_sfcatm1_timeslice(:,:,:) = int_sfcatm1_timeslice(:,:,:) + loc_dtyr*dum_sfcatm1(:,:,:) int_focnatm_timeslice(:,:,:) = int_focnatm_timeslice(:,:,:) + loc_dtyr*locij_focnatm(:,:,:) int_phys_ocnatm_timeslice(:,:,:) = int_phys_ocnatm_timeslice(:,:,:) + loc_dtyr*phys_ocnatm(:,:,:) ! update time slice data - ocean-sediment interface int_sfcsed1_timeslice(:,:,:) = int_sfcsed1_timeslice(:,:,:) + loc_dtyr*dum_sfcsed1(:,:,:) int_focnsed_timeslice(:,:,:) = int_focnsed_timeslice(:,:,:) + locij_focnsed(:,:,:) int_fsedocn_timeslice(:,:,:) = int_fsedocn_timeslice(:,:,:) + locij_fsedocn(:,:,:) ! update time-slice data - GOLDSTEIn CALL sub_calc_psi(dum_u,loc_opsi,loc_opsia,loc_opsip,loc_zpsi) int_opsi_timeslice(:,:) = int_opsi_timeslice(:,:) + loc_dtyr*loc_opsi(:,:) int_opsia_timeslice(:,:) = int_opsia_timeslice(:,:) + loc_dtyr*loc_opsia(:,:) int_opsip_timeslice(:,:) = int_opsip_timeslice(:,:) + loc_dtyr*loc_opsip(:,:) int_zpsi_timeslice(:,:) = int_zpsi_timeslice(:,:) + loc_dtyr*loc_zpsi(:,:) int_u_timeslice(:,:,:,:) = int_u_timeslice(:,:,:,:) + loc_dtyr*dum_u(:,:,:,:) IF (loc_t < par_data_save_timeslice(par_data_save_timeslice_i)) THEN ! save time-slice data CALL sub_data_save_timeslice() If (opt_misc(iopt_misc_sed_select)) CALL sub_data_save_timeslice_sed() ! update time slice index par_data_save_timeslice_i = par_data_save_timeslice_i - 1 ! reset array values CALL sub_init_int_timeslice() END IF END IF ! *** RUN-TIME DATA UPDATE *** ! update signal-averaging ingetrated time int_t_sig = int_t_sig + loc_dtyr int_t_sig_count = int_t_sig_count + 1 ! update signal-averaging ingetrated tracers ! NOTE: this is performed in a seperate series of loops, rather than from which the main program loop ! (where computational efficiency would be greater) ! to avoid numerical errors where the number of integration time steps >> 1 and ! the model is compiled with single precision reals ! (i.e., problems with loss of decimal accuracy when adding small quantities to a large quantity - ! helped by summing first over the entire ocean or atmosphere, and then adding this to the running integrated total) ! NOTE: mass-weight ocean tracers, and surface area-weight atmospheric tracers ! NOTE: for sea surface (SS) tracer values weight by fractional ice-free cover IF (opt_data(iopt_data_save_sig_ocn)) THEN DO io=1,n_ocn IF (ocn_select(io)) THEN int_ocn_sig(io) = int_ocn_sig(io) + & & loc_dtyr*SUM(phys_ocn(ipo_M,:,:,:)*ocn(io,:,:,:))/loc_ocn_tot_M end if END DO end if IF (opt_data(iopt_data_save_sig_atm)) THEN DO ia=1,n_atm IF (ocnatm_select(ia)) THEN int_ocnatm_sig(ia) = int_ocnatm_sig(ia) + & & loc_dtyr*SUM(phys_ocnatm(ipoa_A,:,:)*dum_sfcatm1(ia,:,:))/loc_ocnatm_tot_A end if END DO end if IF (opt_data(iopt_data_save_sig_focnatm)) THEN DO ia=1,n_atm IF (ocnatm_select(ia)) THEN int_ocnatm_flux_sig(ia) = int_ocnatm_flux_sig(ia) + & & loc_dtyr*SUM(locij_focnatm(ia,:,:)) end if END DO end if IF (opt_data(iopt_data_save_sig_focnsed)) THEN DO is=1,n_sed IF (ocnsed_select(is)) THEN int_ocnsed_flux_sig(is) = int_ocnsed_flux_sig(is) + & & SUM(locij_focnsed(is,:,:)) end if END DO end if IF (opt_data(iopt_data_save_sig_ocnSS)) THEN DO io=1,n_ocn IF (ocn_select(io)) THEN int_ocnSS_sig(io) = int_ocnSS_sig(io) + & & loc_dtyr*SUM((1.0 - phys_ocnatm(ipoa_seaice,:,:))*phys_ocn(ipo_A,:,:,n_kmax)*ocn(io,:,:,n_kmax))/loc_ocn_tot_A end if END DO end if IF (opt_data(iopt_data_save_sig_misc)) THEN CALL sub_calc_psi(dum_u,loc_opsi,loc_opsia,loc_opsip,loc_zpsi) int_misc_seaice_sig = int_misc_seaice_sig + & & loc_dtyr*SUM(phys_ocn(ipo_A,:,:,n_kmax)*phys_ocnatm(ipoa_seaice,:,:)) int_misc_THCmin_sig = int_misc_THCmin_sig + & & loc_dtyr*minval(loc_opsia(:,:)) int_misc_THCmax_sig = int_misc_THCmax_sig + & & loc_dtyr*maxval(loc_opsia(:,:)) end if ! data save + audit + run-time echo ! NOTE: data save at mid-points of data save time series control point intervals ! NOTE: also carry out an updated tracer inventory audit at this time (if selected) IF (loc_t < par_data_save_sig(par_data_save_sig_i)) THEN ! run-time data echo-ing CALL sub_calc_psi(dum_u,loc_opsi,loc_opsia,loc_opsip,loc_zpsi) ! \/\/\/ UN-COMMENT TO PERIODICALLY RE-PRINT HEADER INFORMATION \/\/\/ !!$ if (mod(par_data_save_sig_i,50) == 0) par_misc_t_echo_header = .TRUE. ! /\/\/\ UN-COMMENT TO PERIODICALLY RE-PRINT HEADER INFORMATION /\/\/\ call sub_echo_runtime(loc_yr,loc_ocn_tot_M,loc_ocn_tot_A,loc_ocnatm_tot_A,loc_opsi_scale,loc_opsia(:,:),dum_sfcatm1(:,:,:)) IF (int_t_sig > 0.0) then CALL sub_data_save_runtime(loc_t + (int_t_sig - (par_data_save_sig(par_data_save_sig_i) - loc_t))/2.0) end if IF (par_data_save_sig_i > 1) par_data_save_sig_i = par_data_save_sig_i - 1 ! reset array values CALL sub_init_int_timeseries() ! audit IF (opt_misc(iopt_misc_audit)) THEN CALL sub_audit_update() END IF END IF END IF ! *** TEST FOR END-OF-RUN *** ! if local (BioGeM) model time has surppassed the set model end time, then set the flag updating system biogeochemistry to false IF (loc_t <= 0.0) THEN par_misc_t_go = .FALSE. END IF ! clean-up and terminate model if biogeochemistry 'go' is FALSE and BioGeM rather than goldstein model time-control is requested IF (.NOT. par_misc_t_go .AND. opt_misc(iopt_misc_t_timescale_BioGeM)) THEN ! end reporting CALL sub_calc_psi(dum_u,loc_opsi,loc_opsia,loc_opsip,loc_zpsi) par_misc_t_echo_header = .TRUE. call sub_echo_runtime(loc_yr,loc_ocn_tot_M,loc_ocn_tot_A,loc_ocnatm_tot_A,loc_opsi_scale,loc_opsia(:,:),dum_sfcatm1(:,:,:)) CALL end_biogem(dum_sfcatm1(:,:,:)) STOP END IF END SUBROUTINE tstep_biogem ! *** biogem DIAGNOSTICS *** SUBROUTINE diag_biogem(dum_t,dum_dt,dum_u,dum_opsi,dum_opsia,dum_opsip,dum_zpsi) USE biogem_lib USE biogem_box use biogem_data IMPLICIT NONE ! dummy arguments real,intent(in)::dum_t,dum_dt REAL,DIMENSION(3,n_maxi,n_maxj,n_maxk),INTENT(in)::dum_u REAL,DIMENSION(0:n_maxj,0:n_maxk),INTENT(in)::dum_opsi,dum_opsia,dum_opsip,dum_zpsi ! local variables real::loc_t,loc_dts,loc_dtyr ! *** CALCULATE LOCAL CONSTANTS *** ! calculate gemchemical model time loc_t = ABS(par_misc_t_end - par_misc_t_start) - (goldstein_tsc*dum_t)/conv_yr_s - par_misc_t_err loc_dts = goldstein_tsc*dum_dt loc_dtyr = loc_dts/conv_yr_s ! *** TIME-SLICE DATA UPDATE *** ! update time slice data IF ( & & ((par_data_save_timeslice(par_data_save_timeslice_i) + par_data_save_timeslice_dt) > loc_t) & & .AND. & & (par_data_save_timeslice_i > 0) & & ) THEN ! update time-slice data - GOLDSTEIn int_opsi_timeslice(:,:) = int_opsi_timeslice(:,:) + loc_dtyr*dum_opsi(:,:) int_opsia_timeslice(:,:) = int_opsia_timeslice(:,:) + loc_dtyr*dum_opsia(:,:) int_opsip_timeslice(:,:) = int_opsip_timeslice(:,:) + loc_dtyr*dum_opsip(:,:) int_zpsi_timeslice(:,:) = int_zpsi_timeslice(:,:) + loc_dtyr*dum_zpsi(:,:) int_u_timeslice(:,:,:,:) = int_u_timeslice(:,:,:,:) + loc_dtyr*dum_u(:,:,:,:) IF (loc_t <= (par_data_save_timeslice(par_data_save_timeslice_i) + loc_dtyr)) THEN ! save time-slice data ! CALL sub_data_save_timeslice_goldstein() ! reset array values ! CALL sub_init_int_timeslice_goldstein() END IF END IF ! ********************* ! *** *** ! ********************* end SUBROUTINE diag_biogem ! *** RESTART BioGeM (save data) *** SUBROUTINE rest_biogem(dum_filestring) USE biogem_lib IMPLICIT NONE ! dummy arguments CHARACTER(LEN=*),INTENT(in)::dum_filestring ! local variables CHARACTER(len=255)::loc_filename ! dump restart data ! NOTE: data is saved unformatted for minimal file size ! also means that arrays can be written directly to file without needing to loop thought data loc_filename = TRIM(string_results_dir)//trim(dum_filestring)//'.biogem' OPEN(unit=out,status='replace',file=loc_filename,form='unformatted',action='write') WRITE(unit=out) & & ocn(:,:,:,:), & & bio_part(:,:,:,:) close(unit=out) end SUBROUTINE rest_biogem ! *** END BioGeM *** SUBROUTINE end_biogem(dum_sfcsumatm1) USE biogem_lib USE biogem_data IMPLICIT NONE ! dummy arguments REAL,DIMENSION(0:n_atm,n_maxi,n_maxj),INTENT(in)::dum_sfcsumatm1 ! local variables integer::ia,io CHARACTER(len=255)::loc_filename,loc_filename_in,loc_filename_out real,DIMENSION(n_maxi,n_maxj)::loc_ij real,DIMENSION(n_maxi,n_maxj,n_maxk)::loc_ijk ! close all currently open units (used for runtime data save) CALL sub_end_data_save_runtime() ! save copy of parameter files - model configuration IF (opt_data(iopt_data_save_config)) THEN loc_filename_in = TRIM(string_biogem_dir)// & & 'biogem_config.par' loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_CFG_'//'biogem_config'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) end IF ! save copy of parameter files - tracer configuration IF (opt_data(iopt_data_save_config)) THEN loc_filename_in = TRIM(string_data_dir)// & & 'gem_config_atm.par' loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_CFG_'//'gem_config_atm'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) loc_filename_in = TRIM(string_data_dir)// & & 'gem_config_ocn.par' loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_CFG_'//'gem_config_ocn'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) loc_filename_in = TRIM(string_data_dir)// & & 'gem_config_sed.par' loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_CFG_'//'gem_config_sed'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) end IF ! save copy of forcing files IF (opt_data(iopt_data_save_config)) THEN DO ia=1,n_atm IF (ocnatm_select(ia)) THEN IF (force_restore_atm_select(ia)) THEN loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_restore_atm_'//TRIM(string_atm(ia))//'_I'//string_data_ext CALL sub_load_data_ij(loc_filename_in,n_maxi,n_maxj,loc_ij) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_restore_atm_'//TRIM(string_atm(ia))//'_I'//string_results_ext CALL sub_save_data_ij(loc_filename_out,n_maxi,n_maxj,loc_ij) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_restore_atm_'//TRIM(string_atm(ia))//'_II'//string_data_ext CALL sub_load_data_ij(loc_filename_in,n_maxi,n_maxj,loc_ij) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_restore_atm_'//TRIM(string_atm(ia))//'_II'//string_results_ext CALL sub_save_data_ij(loc_filename_out,n_maxi,n_maxj,loc_ij) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_restore_atm_'//TRIM(string_atm(ia))//'_sig'//string_data_ext loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_restore_atm_'//TRIM(string_atm(ia))//'_sig'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) end if IF (force_flux_atm_select(ia)) THEN loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_flux_atm_'//TRIM(string_atm(ia))//'_I'//string_data_ext CALL sub_load_data_ij(loc_filename_in,n_maxi,n_maxj,loc_ij) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_flux_atm_'//TRIM(string_atm(ia))//'_I'//string_results_ext CALL sub_save_data_ij(loc_filename_out,n_maxi,n_maxj,loc_ij) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_flux_atm_'//TRIM(string_atm(ia))//'_II'//string_data_ext CALL sub_load_data_ij(loc_filename_in,n_maxi,n_maxj,loc_ij) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_flux_atm_'//TRIM(string_atm(ia))//'_II'//string_results_ext CALL sub_save_data_ij(loc_filename_out,n_maxi,n_maxj,loc_ij) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_flux_atm_'//TRIM(string_atm(ia))//'_sig'//string_data_ext loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_flux_atm_'//TRIM(string_atm(ia))//'_sig'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) END IF end if END DO DO io=1,n_ocn IF (ocn_select(io)) THEN IF (force_restore_ocn_select(io)) THEN loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_restore_ocn_'//TRIM(string_ocn(io))//'_I'//string_data_ext CALL sub_load_data_ijk(loc_filename_in,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_restore_ocn_'//TRIM(string_ocn(io))//'_I'//string_results_ext CALL sub_save_data_ijk(loc_filename_out,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_restore_ocn_'//TRIM(string_ocn(io))//'_II'//string_data_ext CALL sub_load_data_ijk(loc_filename_in,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_restore_ocn_'//TRIM(string_ocn(io))//'_II'//string_results_ext CALL sub_save_data_ijk(loc_filename_out,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_restore_ocn_'//TRIM(string_ocn(io))//'_sig'//string_data_ext loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_restore_ocn_'//TRIM(string_ocn(io))//'_sig'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) END IF IF (force_flux_ocn_select(io)) THEN loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_flux_ocn_'//TRIM(string_ocn(io))//'_I'//string_data_ext CALL sub_load_data_ijk(loc_filename_in,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_flux_ocn_'//TRIM(string_ocn(io))//'_I'//string_results_ext CALL sub_save_data_ijk(loc_filename_out,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_flux_ocn_'//TRIM(string_ocn(io))//'_II'//string_data_ext CALL sub_load_data_ijk(loc_filename_in,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_flux_ocn_'//TRIM(string_ocn(io))//'_II'//string_results_ext CALL sub_save_data_ijk(loc_filename_out,n_maxi,n_maxj,n_maxk,loc_ijk) loc_filename_in = TRIM(string_data_dir)// & & 'biogem_force_flux_ocn_'//TRIM(string_ocn(io))//'_sig'//string_data_ext loc_filename_out = TRIM(string_results_dir)// & & TRIM(string_runid)//'_FORCE_flux_ocn_'//TRIM(string_ocn(io))//'_sig'//string_results_ext call sub_copy_ascii_file(trim(loc_filename_in),trim(loc_filename_out)) END IF end if END DO end IF print*,'! save copy of forcing files' ! save misc diagnostics IF (opt_data(iopt_data_save_misc)) THEN loc_filename = string_results_dir//trim(string_runid)//'_DIAG_totalcarbiterations'//string_results_ext CALL sub_save_data_ijk(loc_filename,n_maxi,n_maxj,n_maxk,carb_TSn(3,:,:,:)) end IF print*,'! save misc diagnostics' ! output audit diagnostics (final reporting of initial/final tracer inventories) IF (opt_misc(iopt_misc_audit) .AND. opt_misc(iopt_misc_audit_reportfinal)) THEN PRINT*,' ' PRINT*,'*** BioGeM tracer audit diagnostics ***' CALL sub_audit_update() CALL sub_data_audit_diagnostics() END IF print*,'! output audit diagnostics (final reporting of initial/final tracer inventories)' ! print indicative carbon cycle diagnostics (if ocean CO2 is selected) (just some random useful info) IF (opt_misc(iopt_misc_audit_reportfinal)) THEN IF (ocn_select(io_CO2)) THEN PRINT*,' ' PRINT*,'*** BioGeM carbon cycle diagnostics ***' PRINT*,'FINAL (area-weighted) mean surface pCO2 (uatm): ', & & conv_mol_umol*SUM(dum_sfcsumatm1(ia_pCO2,:,:)*phys_ocnatm(ipoa_A,:,:))/SUM(phys_ocnatm(ipoa_A,:,:)) IF (opt_misc(iopt_misc_audit)) then PRINT*,'INITIAL / FINAL ocean CO2 inventory (mol): ', & & audit_ocn_init(io_CO2), & & ' ', & & audit_ocn_new(io_CO2) end IF PRINT*,'***************************************' PRINT*,' ' ELSE PRINT*,'***************************************' PRINT*,' ' END IF end if END SUBROUTINE end_biogem ! *** EnKF HACK *** SUBROUTINE enkf_biogem() USE biogem_lib IMPLICIT NONE ! local variables real,DIMENSION(12)::loc_enkf ! *** LOAD ENKF PARAMETERS *** ! organic carbon biogeochem cycle read(unit=5,fmt=*) loc_enkf(1) read(unit=5,fmt=*) loc_enkf(2) read(unit=5,fmt=*) loc_enkf(11) read(unit=5,fmt=*) loc_enkf(12) read(unit=5,fmt=*) loc_enkf(3) read(unit=5,fmt=*) loc_enkf(4) read(unit=5,fmt=*) loc_enkf(5) ! inorganic carbon biogeochem cycle read(unit=5,fmt=*) loc_enkf(6) read(unit=5,fmt=*) loc_enkf(7) read(unit=5,fmt=*) loc_enkf(8) read(unit=5,fmt=*) loc_enkf(9) read(unit=5,fmt=*) loc_enkf(10) ! *** SET ENKF PARAMETERS *** ! organic carbon biogeochem cycle par_bio_k0_PO4 = loc_enkf(1) par_bio_c0_PO4 = loc_enkf(2) par_bio_red_DOCfrac = loc_enkf(11) par_bio_remin_DOClifetime = loc_enkf(12) par_bio_remin_POC_frac2 = loc_enkf(3) par_bio_remin_POC_eL1 = loc_enkf(4) par_bio_remin_POC_eL2 = loc_enkf(5) ! inorganic carbon biogeochem cycle par_bio_red_POC_CaCO3 = loc_enkf(6) par_bio_red_POC_CaCO3_pP = loc_enkf(7) par_bio_remin_CaCO3_frac2 = loc_enkf(8) par_bio_remin_CaCO3_eL1 = loc_enkf(9) par_bio_remin_CaCO3_eL2 = loc_enkf(10) ! DISPLAY print*,' ' print*,'***************************' print*,'*** biogem EnFK version ***' print*,'***************************' print*,' ' print*,'par_bio_k0_PO4 : ',par_bio_k0_PO4 print*,'par_bio_c0_PO4 : ',par_bio_c0_PO4 print*,'par_bio_red_DOCfrac : ',par_bio_red_DOCfrac print*,'par_bio_remin_DOClifetime : ',par_bio_remin_DOClifetime print*,'par_bio_remin_POC_frac2 : ',par_bio_remin_POC_frac2 print*,'par_bio_remin_POC_eL1 : ',par_bio_remin_POC_eL1 print*,'par_bio_remin_POC_eL2 : ',par_bio_remin_POC_eL2 print*,'par_bio_red_POC_CaCO3 : ',par_bio_red_POC_CaCO3 print*,'par_bio_red_POC_CaCO3_pP : ',par_bio_red_POC_CaCO3_pP print*,'par_bio_remin_CaCO3_frac2 : ',par_bio_remin_CaCO3_frac2 print*,'par_bio_remin_CaCO3_eL1 : ',par_bio_remin_CaCO3_eL1 print*,'par_bio_remin_CaCO3_eL2 : ',par_bio_remin_CaCO3_eL2 print*,' ' end SUBROUTINE enkf_biogem ! ******* ! *************** ! ******************** ! ********************** ! *** CODE GRAVEYARD *** ! ********************** ! ********************** ! ** ** ** ** ! ** ** ** ** ** ** ! ** ** ** ** ** ** ! ** * *** ** ** ! ** ** ** ** ****** ! ** ** ** ** ****** ! ** ** ** ** ****** ! ********************** ! ********************** ! ********************** ! ********************** !!$! *** copy and make available to GeM, variable values local to in GOLDSTEIn *** !!$SUBROUTINE sub_biogem_init_copyvar( & !!$ & dum_nsteps,dum_npstp,dum_iwstp,dum_itstp,dum_t0, & !!$ & dum_usc,dum_rsc,dum_dsc,dum_fsc,dum_gsc,dum_rh0sc,dum_rhosc,dum_cpsc,dum_tsc, & !!$ & dum_pi,dum_saln0,dum_rhoair,dum_cd,dum_ds,dum_dphi, & !!$ & dum_ips,dum_ipf,dum_ias,dum_iaf,dum_jsf, & !!$ & dum_k1,dum_tau, & !!$ & dum_dz,dum_dza, & !!$ & dum_c,dum_cv,dum_s,dum_sv & !!$ & ) !!$ USE biogem_lib !!$ IMPLICIT NONE !!$ ! dummy arguments !!$ INTEGER,INTENT(in)::dum_nsteps,dum_npstp,dum_iwstp,dum_itstp !!$ real,intent(in)::dum_t0 !!$ real,INTENT(in)::dum_usc,dum_rsc,dum_dsc,dum_fsc,dum_gsc,dum_rh0sc,dum_rhosc,dum_cpsc,dum_tsc !!$ REAL,INTENT(in)::dum_pi,dum_saln0,dum_rhoair,dum_cd,dum_ds,dum_dphi !!$ INTEGER,INTENT(in),DIMENSION(1:n_maxj)::dum_ips,dum_ipf,dum_ias,dum_iaf !!$ INTEGER,INTENT(in)::dum_jsf !!$ integer,DIMENSION(n_maxi,n_maxj),INTENT(in)::dum_k1 !!$ REAL,DIMENSION(2,n_maxi,n_maxj),INTENT(in)::dum_tau !!$ REAL,DIMENSION(n_maxk),INTENT(in)::dum_dz,dum_dza !!$ REAL,DIMENSION(0:n_maxj),INTENT(in)::dum_c,dum_cv,dum_s,dum_sv !!$ ! *** copy GOLDSTEIn parameters *** !!$ ! !!$ goldstein_nsteps = dum_nsteps !!$ goldstein_npstp = dum_npstp !!$ goldstein_iwstp = dum_iwstp !!$ goldstein_itstp = dum_itstp !!$ goldstein_t0 = dum_t0 !!$ ! copy dimensional scale factors for ocean !!$ goldstein_usc = dum_usc !!$ goldstein_rsc = dum_rsc !!$ goldstein_dsc = dum_dsc !!$ goldstein_fsc = dum_fsc !!$ goldstein_gsc = dum_gsc !!$ goldstein_rh0sc = dum_rh0sc !!$ goldstein_rhosc = dum_rhosc !!$ goldstein_cpsc = dum_cpsc !!$ goldstein_tsc = dum_tsc !!$ ! copy miscellaneous constants !!$ ! NOTE: use identical value of pi to avoid possible GOLDETSIn/BioGeM mis-match in the calculation of ocean areas and volumes !!$ goldstein_pi = dum_pi !!$ goldstein_saln0 = dum_saln0 !!$ goldstein_rhoair = dum_rhoair !!$ goldstein_cd = dum_cd !!$ goldstein_ds = dum_ds !!$ goldstein_dphi = dum_dphi !!$ ! copy ocean positions !!$ goldstein_jsf = dum_jsf !!$ goldstein_ips(:) = dum_ips(:) !!$ goldstein_ipf(:) = dum_ipf(:) !!$ goldstein_ias(:) = dum_ias(:) !!$ goldstein_iaf(:) = dum_iaf(:) !!$ ! copy ocean bottom index grid !!$ goldstein_k1(:,:) = dum_k1(:,:) !!$ ! copy wind stress field !!$ goldstein_tau(:,:,:) = dum_tau(:,:,:) !!$ ! miscellaneous !!$ goldstein_dz(:) = dum_dz(:) !!$ goldstein_dza(:) = dum_dza(:) !!$ goldstein_c(:) = dum_c(:) !!$ goldstein_cv(:) = dum_cv(:) !!$ goldstein_s(:) = dum_s(:) !!$ goldstein_sv(:) = dum_sv(:) !!$END SUBROUTINE sub_biogem_init_copyvar !!$! *** initialize BioGeM *** !!$SUBROUTINE sub_biogem_init(dum_t,dum_tsc) !!$ USE biogem_lib !!$ USE biogem_data !!$ USE gem_data !!$ IMPLICIT NONE !!$ ! dummy arguments !!$ REAL,INTENT(in)::dum_t,dum_tsc !!$ ! local variables !!$ INTEGER::l,m,io,ia !!$ real::loc_t !!$ ! initialize variables !!$ par_force_flux_weather(:) = 0.0 !!$ ! enable biogeochemistry to be updated !!$ par_misc_t_go = .TRUE. !!$ ! load BioGeM run-time options and set initial (BioGeM) model time !!$ CALL sub_load_biogem_config() !!$ loc_t = ABS(par_misc_t_end - par_misc_t_start) - (dum_tsc*dum_t)/conv_yr_s !!$ ! initialize common (i.e., those defined in biogem_lib.mod) BioGeM arrays !!$ CALL sub_init_misc() !!$ CALL sub_init_char() !!$ CALL sub_init_int_timeslice() !!$ CALL sub_init_int_timeseries() !!$ CALL sub_init_phys_ocn() !!$ CALL sub_init_phys_atm() !!$ IF (opt_misc(iopt_misc_audit)) CALL sub_init_audit() !!$ CALL sub_init_force() !!$ ! load default values for ocean, atmosphere, and sediment tracers !!$ CALL sub_init_tracer_ocn() !!$ CALL sub_init_tracer_atm() !!$ CALL sub_init_tracer_sed() !!$ CALL sub_init_bio() !!$ ! set derived tracer selection options !!$ CALL sub_init_select() !!$ ! initialize carbon cycle sub-systems (if selected) !!$ IF (opt_select(iopt_select_carbchem)) CALL sub_init_carb() !!$ ! initialize basic data saving !!$ CALL sub_init_data_save() !!$ ! open units for runtime file I/O !!$ CALL sub_init_data_save_runtime(loc_t) !!$ ! initialize forcings !!$ ! NOTE: tracer loops follow the GOLDSTEIn notation, with indices converted to their BioGeM equivalents !!$ ! 1a. set up forcing functions - atmosphere !!$ DO m=1,n_mmax !!$ ia = goldstein_m_ia(m) !!$ IF (force_restore_atm_select(ia)) CALL sub_init_force_restore_atm(ia) !!$ IF (force_flux_atm_select(ia)) CALL sub_init_force_flux_atm(ia) !!$ END DO !!$ ! 1b. set up forcing functions - ocean !!$ DO l=1,n_lmax !!$ io = goldstein_l_io(l) !!$ IF (force_restore_ocn_select(io)) CALL sub_init_force_restore_ocn(io) !!$ IF (force_flux_ocn_select(io)) CALL sub_init_force_flux_ocn(io) !!$ END DO !!$ ! 2a. initialize forcings - atmosphere !!$ DO m=1,n_mmax !!$ ia = goldstein_m_ia(m) !!$ IF (force_restore_atm_select(ia)) THEN !!$ CALL sub_update_force_restore_atm(loc_t,ia) !!$ end if !!$ IF (force_flux_atm_select(ia)) THEN !!$ CALL sub_update_force_flux_atm(loc_t,ia) !!$ END IF !!$ END DO !!$ ! 2b. initialize forcings - ocean !!$ DO l=1,n_lmax !!$ io = goldstein_l_io(l) !!$ IF (force_restore_ocn_select(io)) THEN !!$ CALL sub_update_force_restore_ocn(loc_t,io) !!$ END IF !!$ IF (force_flux_ocn_select(io)) THEN !!$ CALL sub_update_force_flux_ocn(loc_t,io) !!$ END IF !!$ END DO !!$ ! misc init !!$ ! pre-calculate particulate transformation ratios !!$ If (opt_bio(iopt_bio_remin_fixed)) then !!$ call sub_init_bio_remin_frac() !!$ end if !!$ ! save model topology data !!$ CALL sub_data_save_topology() !!$END SUBROUTINE sub_biogem_init !!$! *** set the value of in GOLDSTEIn *** !!$SUBROUTINE sub_biogem_set_lmax(dum_lmax) !!$ USE biogem_lib !!$ IMPLICIT NONE !!$ INTEGER,INTENT(inout)::dum_lmax !!$ goldstein_lmax = dum_lmax !!$ dum_lmax = n_lmax !!$END SUBROUTINE sub_biogem_set_lmax !!$! *** reset the value of in GOLDSTEIn *** !!$SUBROUTINE sub_biogem_reset_lmax(dum_lmax) !!$ USE biogem_lib !!$ IMPLICIT NONE !!$ INTEGER,INTENT(inout)::dum_lmax !!$ dum_lmax = goldstein_lmax !!$END SUBROUTINE sub_biogem_reset_lmax !!$! *** hack C-GOLDSTEIn CO2 *** !!$SUBROUTINE sub_biogem_co2_update(dum_co2) !!$ USE biogem_lib !!$ IMPLICIT NONE !!$ ! dummy arguments !!$ REAL,DIMENSION(1:n_maxi,1:n_maxj),INTENT(inout)::dum_co2 !!$ ! update co2 if requested !!$ IF (opt_force(iopt_force_GOLDSTEIn_CO2) .AND. opt_select(iopt_select_ocnatm_CO2)) THEN !!$ dum_co2(:,:) = conv_mol_umol*atm(ia_pCO2,:,:) !!$ END IF !!$END SUBROUTINE sub_biogem_co2_update !!$! *** INITIAL INVENTORY AUDIT *** !!$SUBROUTINE sub_biogem_audit_init() !!$ USE biogem_box !!$ IMPLICIT NONE !!$ ! carry out initial tracer inventory audit !!$ IF (opt_misc(iopt_misc_audit)) THEN !!$ CALL sub_audit_init() !!$ END IF !!$END SUBROUTINE sub_biogem_audit_init !!$! *** INVENTORY AUDIT UPDATE *** !!$SUBROUTINE sub_biogem_audit_update(dum_istep) !!$ USE biogem_box !!$ IMPLICIT NONE !!$ ! dummy arguments !!$ INTEGER,INTENT(in)::dum_istep !!$ ! carry out an updated tracer inventory audit !!$ ! NOTE: code to excecute every iterations !!$ IF (opt_misc(iopt_misc_audit) .AND. MOD(dum_istep,par_misc_auditstp) == 0) THEN !!$ CALL sub_audit_update() !!$ END IF !!$END SUBROUTINE sub_biogem_audit_update !!$! *** RESTART GeM (save data) *** !!$SUBROUTINE rest_gem(dum_filestring) !!$ USE biogem_lib !!$ IMPLICIT NONE !!$ ! dummy arguments !!$ CHARACTER(LEN=*),INTENT(in)::dum_filestring !!$ ! local variables !!$ CHARACTER(len=255)::loc_filename !!$ ! dump data !!$ loc_filename = TRIM(string_results_dir)//trim(dum_filestring)//'.gem' !!$ OPEN(unit=out,status='replace',file=loc_filename,form='unformatted',action='write') !!$ WRITE(unit=out) & !!$ & interf_ocnsed_ocn(:,:,:), & !!$ & interf_ocnsed_phys(:,:,:), & !!$ & interf_sedtop(:,:,:), & !!$ & interf_ocnsed_flux(:,:,:), & !!$ & interf_sedocn_flux(:,:,:), & !!$ & interf_ocnatm_flux(:,:,:) !!$ close(unit=out) !!$end SUBROUTINE rest_gem