MODULE trcbio_medusa !!====================================================================== !! *** MODULE trcbio *** !! TOP : MEDUSA !!====================================================================== !! History : !! - ! 1999-07 (M. Levy) original code !! - ! 2000-12 (E. Kestenare) assign parameters to name individual tracers !! - ! 2001-03 (M. Levy) LNO3 + dia2d !! 2.0 ! 2007-12 (C. Deltel, G. Madec) F90 !! - ! 2008-08 (K. Popova) adaptation for MEDUSA !! - ! 2008-11 (A. Yool) continuing adaptation for MEDUSA !! - ! 2010-03 (A. Yool) updated for branch inclusion !! - ! 2011-08 (A. Yool) updated for ROAM (see below) !! - ! 2013-03 (A. Yool) updated for iMARNET !! - ! 2013-05 (A. Yool) updated for v3.5 !! - ! 2014-08 (A. Yool, J. Palm) Add DMS module for UKESM1 model !! - ! 2015-06 (A. Yool) Update to include MOCSY !! - ! 2015-07 (A. Yool) Update for rolling averages !! - ! 2015-10 (J. Palm) Update for diag outputs through iom_use !! - ! 2016-11 (A. Yool) Updated diags for CMIP6 !! - ! 2017-05 (A. Yool) Added extra DMS calculation !!---------------------------------------------------------------------- !! #if defined key_roam !!---------------------------------------------------------------------- !! Updates for the ROAM project include: !! - addition of DIC, alkalinity, detrital carbon and oxygen tracers !! - addition of air-sea fluxes of CO2 and oxygen !! - periodic (monthly) calculation of full 3D carbonate chemistry !! - detrital C:N ratio now free to evolve dynamically !! - benthic storage pools !! !! Opportunity also taken to add functionality: !! - switch for Liebig Law (= most-limiting) nutrient uptake !! - switch for accelerated seafloor detritus remineralisation !! - switch for fast -> slow detritus transfer at seafloor !! - switch for ballast vs. Martin vs. Henson fast detritus remin. !! - per GMD referee remarks, xfdfrac3 introduced for grazed PDS !!---------------------------------------------------------------------- #endif !! #if defined key_mocsy !!---------------------------------------------------------------------- !! Updates with the addition of MOCSY include: !! - option to use PML or MOCSY carbonate chemistry (the latter is !! preferred) !! - central calculation of gas transfer velocity, f_kw660; previously !! this was done separately for CO2 and O2 with predictable results !! - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux !! calculations and to those for O2 air-sea flux !! - extra diagnostics included for MOCSY !!---------------------------------------------------------------------- #endif !! #if defined key_medusa !!---------------------------------------------------------------------- !! MEDUSA bio-model !!---------------------------------------------------------------------- !! trc_bio_medusa : !!---------------------------------------------------------------------- USE oce_trc USE trc USE sms_medusa USE lbclnk USE prtctl_trc ! Print control for debugging USE trcsed_medusa USE sbc_oce ! surface forcing USE sbcrnf ! surface boundary condition: runoff variables USE in_out_manager ! I/O manager # if defined key_iomput USE iom USE trcnam_medusa ! JPALM 13-11-2015 -- if iom_use for diag !!USE trc_nam_iom_medusa ! JPALM 13-11-2015 -- if iom_use for diag # endif # if defined key_roam USE gastransfer # if defined key_mocsy USE mocsy_wrapper # else USE trcco2_medusa # endif USE trcoxy_medusa !! Jpalm (08/08/2014) USE trcdms_medusa # endif !! AXY (18/01/12): brought in for benthic timestepping USE trcnam_trp ! AXY (24/05/2013) USE trdmxl_trc USE trdtrc_oce ! AXY (24/05/2013) !! AXY (30/01/14): necessary to find NaNs on HECTOR USE, INTRINSIC :: ieee_arithmetic !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm USE sbc_oce, ONLY: lk_oasis USE oce, ONLY: CO2Flux_out_cpl, DMS_out_cpl, PCO2a_in_cpl, chloro_out_cpl IMPLICIT NONE PRIVATE PUBLIC trc_bio_medusa ! called in ??? !!* Substitution # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trc_bio_medusa( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE trc_bio *** !! !! ** Purpose : compute the now trend due to biogeochemical processes !! and add it to the general trend of passive tracers equations !! !! ** Method : each now biological flux is calculated in function of now !! concentrations of tracers. !! depending on the tracer, these fluxes are sources or sinks. !! the total of the sources and sinks for each tracer !! is added to the general trend. !! !! tra = tra + zf...tra - zftra... !! | | !! | | !! source sink !! !! IF 'key_trc_diabio' defined , the biogeochemical trends !! for passive tracers are saved for futher diagnostics. !!--------------------------------------------------------------------- !! !! !!---------------------------------------------------------------------- !! Variable conventions !!---------------------------------------------------------------------- !! !! names: z*** - state variable !! f*** - function (or temporary variable used in part of a function) !! x*** - parameter !! b*** - right-hand part (sources and sinks) !! i*** - integer variable (usually used in yes/no flags) !! !! time (integer timestep) INTEGER, INTENT( in ) :: kt !! !! spatial array indices INTEGER :: ji,jj,jk,jn !! !! AXY (27/07/10): add in indices for depth horizons (for sinking flux !! and seafloor iron inputs) !! INTEGER :: i0100, i0200, i0500, i1000, i1100 !! !! model state variables REAL(wp) :: zchn,zchd,zphn,zphd,zpds,zzmi REAL(wp) :: zzme,zdet,zdtc,zdin,zsil,zfer REAL(wp) :: zage # if defined key_roam REAL(wp) :: zdic, zalk, zoxy REAL(wp) :: ztmp, zsal # endif # if defined key_mocsy REAL(wp) :: zpho # endif !! !! integrated source and sink terms REAL(wp) :: b0 !! AXY (23/08/13): changed from individual variables for each flux to !! an array that holds all fluxes REAL(wp), DIMENSION(jp_medusa) :: btra !! !! primary production and chl related quantities REAL(wp) :: fthetan,faln,fchn1,fchn,fjln,fprn,frn REAL(wp) :: fthetad,fald,fchd1,fchd,fjld,fprd,frd !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range) REAL(wp) :: fjlim_pn, fjlim_pd !! AXY (03/02/11): add in Liebig terms REAL(wp) :: fpnlim, fpdlim !! AXY (16/07/09): add in Eppley curve functionality REAL(wp) :: loc_T,fun_T,xvpnT,xvpdT INTEGER :: ieppley !! AXY (16/05/11): per Katya's prompting, add in new T-dependence !! for phytoplankton growth only (i.e. no change !! for remineralisation) REAL(wp) :: fun_Q10 !! AXY (01/03/10): add in mixed layer PP diagnostics REAL(wp), DIMENSION(jpi,jpj) :: fprn_ml,fprd_ml !! !! nutrient limiting factors REAL(wp) :: fnln,ffln !! N and Fe REAL(wp) :: fnld,ffld,fsld,fsld2 !! N, Fe and Si !! !! silicon cycle REAL(wp) :: fsin,fnsi,fsin1,fnsi1,fnsi2,fprds,fsdiss !! !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme REAL(wp) :: ffetop,ffebot,ffescav REAL(wp) :: xLgF, xFeT, xFeF, xFeL !! state variables for iron-ligand system REAL(wp), DIMENSION(jpi,jpj) :: xFree !! state variables for iron-ligand system REAL(wp) :: xb_coef_tmp, xb2M4ac !! iron-ligand parameters REAL(wp) :: xmaxFeF,fdeltaFe !! max Fe' parameters !! !! local parameters for Moore et al. (2004) alternative scavenging scheme REAL(wp) :: fbase_scav,fscal_sink,fscal_part,fscal_scav !! !! local parameters for Moore et al. (2008) alternative scavenging scheme REAL(wp) :: fscal_csink,fscal_sisink,fscal_casink !! !! local parameters for Galbraith et al. (2010) alternative scavenging scheme REAL(wp) :: xCscav1, xCscav2, xk_org, xORGscav !! organic portion of scavenging REAL(wp) :: xk_inorg, xINORGscav !! inorganic portion of scavenging !! !! microzooplankton grazing REAL(wp) :: fmi1,fmi,fgmipn,fgmid,fgmidc REAL(wp) :: finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp !! !! mesozooplankton grazing REAL(wp) :: fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc REAL(wp) :: finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp !! !! mortality/Remineralisation (defunct parameter "fz" removed) REAL(wp) :: fdpn,fdpd,fdpds,fdzmi,fdzme,fdd # if defined key_roam REAL(wp) :: fddc # endif REAL(wp) :: fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2 REAL(wp) :: fslown, fslowc REAL(wp), DIMENSION(jpi,jpj) :: fslownflux, fslowcflux REAL(wp) :: fregen,fregensi REAL(wp), DIMENSION(jpi,jpj) :: fregenfast,fregenfastsi # if defined key_roam REAL(wp) :: fregenc REAL(wp), DIMENSION(jpi,jpj) :: fregenfastc # endif !! !! particle flux REAL(WP) :: fthk,fdep,fdep1,fdep2,flat,fcaco3 REAL(WP) :: ftempn,ftempsi,ftempfe,ftempc,ftempca REAL(wp) :: freminn,freminsi,freminfe,freminc,freminca REAL(wp), DIMENSION(jpi,jpj) :: ffastn,ffastsi,ffastfe,ffastc,ffastca REAL(wp) :: fleftn,fleftsi,fleftfe,fleftc,fleftca REAL(wp) :: fheren,fheresi,fherefe,fherec,fhereca REAL(wp) :: fprotf REAL(wp), DIMENSION(jpi,jpj) :: fsedn,fsedsi,fsedfe,fsedc,fsedca REAL(wp), DIMENSION(jpi,jpj) :: fccd REAL(wp) :: fccd_dep !! AXY (28/11/16): fix mbathy bug INTEGER :: jmbathy !! !! AXY (06/07/11): alternative fast detritus schemes REAL(wp) :: fb_val, fl_sst !! !! AXY (08/07/11): fate of fast detritus reaching the seafloor REAL(wp) :: ffast2slown,ffast2slowfe,ffast2slowc !! !! conservation law REAL(wp) :: fnit0,fsil0,ffer0 # if defined key_roam REAL(wp) :: fcar0,falk0,foxy0 # endif !! !! temporary variables REAL(wp) :: fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8,fq9 !! !! water column nutrient and flux integrals REAL(wp), DIMENSION(jpi,jpj) :: ftot_n,ftot_si,ftot_fe REAL(wp), DIMENSION(jpi,jpj) :: fflx_n,fflx_si,fflx_fe REAL(wp), DIMENSION(jpi,jpj) :: fifd_n,fifd_si,fifd_fe REAL(wp), DIMENSION(jpi,jpj) :: fofd_n,fofd_si,fofd_fe # if defined key_roam REAL(wp), DIMENSION(jpi,jpj) :: ftot_c,ftot_a,ftot_o2 REAL(wp), DIMENSION(jpi,jpj) :: fflx_c,fflx_a,fflx_o2 REAL(wp), DIMENSION(jpi,jpj) :: fifd_c,fifd_a,fifd_o2 REAL(wp), DIMENSION(jpi,jpj) :: fofd_c,fofd_a,fofd_o2 # endif !! !! zooplankton grazing integrals REAL(wp), DIMENSION(jpi,jpj) :: fzmi_i,fzmi_o,fzme_i,fzme_o !! !! limitation term temporary variables REAL(wp), DIMENSION(jpi,jpj) :: ftot_pn,ftot_pd REAL(wp), DIMENSION(jpi,jpj) :: ftot_zmi,ftot_zme,ftot_det,ftot_dtc !! use ballast scheme (1) or simple exponential scheme (0; a conservation test) INTEGER :: iball !! use biological fluxes (1) or not (0) INTEGER :: ibio_switch !! !! diagnose fluxes (should only be used in 1D runs) INTEGER :: idf, idfval !! !! nitrogen and silicon production and consumption REAL(wp) :: fn_prod, fn_cons, fs_prod, fs_cons REAL(wp), DIMENSION(jpi,jpj) :: fnit_prod, fnit_cons, fsil_prod, fsil_cons # if defined key_roam !! !! flags to help with calculating the position of the CCD INTEGER, DIMENSION(jpi,jpj) :: i2_omcal,i2_omarg !! !! ROAM air-sea flux and diagnostic parameters REAL(wp) :: f_wind !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have) REAL(wp) :: f_xco2a REAL(wp) :: f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux REAL(wp) :: f_TDIC, f_TALK, f_dcf, f_henry REAL(wp) :: f_uwind, f_vwind, f_pp0 REAL(wp) :: f_kw660, f_o2flux, f_o2sat, f_o2sat3 REAL(wp), DIMENSION(jpi,jpj) :: f_omcal, f_omarg !! !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen REAL(wp) :: f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm REAL(wp) :: f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2 !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s REAL, PARAMETER :: weight_CO2_mol = 44.0095 !! g / mol REAL, PARAMETER :: secs_in_day = 86400.0 !! s / d REAL, PARAMETER :: CO2flux_conv = (1.e-6 * weight_CO2_mol) / secs_in_day !! INTEGER :: iters REAL(wp) :: f_year INTEGER :: i_year INTEGER :: iyr1, iyr2 !! !! carbon, alkalinity production and consumption REAL(wp) :: fc_prod, fc_cons, fa_prod, fa_cons REAL(wp), DIMENSION(jpi,jpj) :: fcomm_resp REAL(wp), DIMENSION(jpi,jpj) :: fcar_prod, fcar_cons !! !! oxygen production and consumption (and non-consumption) REAL(wp) :: fo2_prod, fo2_cons, fo2_ncons, fo2_ccons REAL(wp), DIMENSION(jpi,jpj) :: foxy_prod, foxy_cons, foxy_anox !! Jpalm (11-08-2014) !! add DMS in MEDUSA for UKESM1 model REAL(wp) :: dms_surf !! AXY (13/03/15): add in other DMS calculations REAL(wp) :: dms_andr, dms_simo, dms_aran, dms_hall, dms_andm, dms_nlim, dms_wtkn # endif !! !! benthic fluxes INTEGER :: ibenthic REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe, f_sbenin_c REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca REAL(wp) :: zfact !! !! benthic fluxes of CaCO3 that shouldn't happen because of lysocline REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca !! !! riverine fluxes REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk !! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface REAL(wp) :: f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk !! !! Jpalm -- 11-10-2015 -- adapt diag to iom_use !! 2D var for diagnostics. REAL(wp), POINTER, DIMENSION(:,: ) :: fprn2d, fdpn2d, fprd2d, fdpd2d, fprds2d, fsdiss2d, fgmipn2d REAL(wp), POINTER, DIMENSION(:,: ) :: fgmid2d, fdzmi2d, fgmepn2d, fgmepd2d, fgmezmi2d, fgmed2d REAL(wp), POINTER, DIMENSION(:,: ) :: fdzme2d, fslown2d, fdd2d, ffetop2d, ffebot2d, ffescav2d REAL(wp), POINTER, DIMENSION(:,: ) :: fjln2d, fnln2d, ffln2d, fjld2d, fnld2d, ffld2d, fsld2d2 REAL(wp), POINTER, DIMENSION(:,: ) :: fsld2d, fregen2d, fregensi2d, ftempn2d, ftempsi2d, ftempfe2d REAL(wp), POINTER, DIMENSION(:,: ) :: ftempc2d, ftempca2d, freminn2d, freminsi2d, freminfe2d REAL(wp), POINTER, DIMENSION(:,: ) :: freminc2d, freminca2d REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d # if defined key_roam REAL(wp), POINTER, DIMENSION(:,: ) :: ffastca2d, rivn2d, rivsi2d, rivc2d, rivalk2d, fslowc2d REAL(wp), POINTER, DIMENSION(:,: ) :: fdpn22d, fdpd22d, fdzmi22d, fdzme22d, zimesn2d, zimesd2d REAL(wp), POINTER, DIMENSION(:,: ) :: zimesc2d, zimesdc2d, ziexcr2d, ziresp2d, zigrow2d, zemesn2d REAL(wp), POINTER, DIMENSION(:,: ) :: zemesd2d, zemesc2d, zemesdc2d, zeexcr2d, zeresp2d, zegrow2d REAL(wp), POINTER, DIMENSION(:,: ) :: mdetc2d, gmidc2d, gmedc2d, f_pco2a2d, f_pco2w2d, f_co2flux2d REAL(wp), POINTER, DIMENSION(:,: ) :: f_TDIC2d, f_TALK2d, f_kw6602d, f_pp02d, f_o2flux2d, f_o2sat2d REAL(wp), POINTER, DIMENSION(:,: ) :: dms_andr2d, dms_simo2d, dms_aran2d, dms_hall2d, dms_andm2d, dms_surf2d REAL(wp), POINTER, DIMENSION(:,: ) :: iben_n2d, iben_fe2d, iben_c2d, iben_si2d, iben_ca2d, oben_n2d REAL(wp), POINTER, DIMENSION(:,: ) :: oben_fe2d, oben_c2d, oben_si2d, oben_ca2d, sfr_ocal2d REAL(wp), POINTER, DIMENSION(:,: ) :: sfr_oarg2d, lyso_ca2d !! AXY (23/11/16): extra MOCSY diagnostics REAL(wp), POINTER, DIMENSION(:,: ) :: f_xco2a_2d, f_fco2w_2d, f_fco2a_2d REAL(wp), POINTER, DIMENSION(:,: ) :: f_ocnrhosw_2d, f_ocnschco2_2d, f_ocnkwco2_2d REAL(wp), POINTER, DIMENSION(:,: ) :: f_ocnk0_2d, f_co2starair_2d, f_ocndpco2_2d # endif !! !! 3D var for diagnostics. REAL(wp), POINTER, DIMENSION(:,:,:) :: tpp3d, detflux3d, remin3dn !! # if defined key_roam !! AXY (04/11/16) !! 2D var for new CMIP6 diagnostics (behind a key_roam ifdef for simplicity) REAL(wp), POINTER, DIMENSION(:,: ) :: fgco2, intdissic, intdissin, intdissisi, inttalk, o2min, zo2min REAL(wp), POINTER, DIMENSION(:,: ) :: fbddtalk, fbddtdic, fbddtdife, fbddtdin, fbddtdisi !! !! 3D var for new CMIP6 diagnostics REAL(wp), POINTER, DIMENSION(:,:,:) :: tppd3 REAL(wp), POINTER, DIMENSION(:,:,:) :: bddtalk3, bddtdic3, bddtdife3, bddtdin3, bddtdisi3 REAL(wp), POINTER, DIMENSION(:,:,:) :: fd_nit3, fd_sil3, fd_car3, fd_cal3 REAL(wp), POINTER, DIMENSION(:,:,:) :: co33, co3satarag3, co3satcalc3, dcalc3 REAL(wp), POINTER, DIMENSION(:,:,:) :: expc3, expn3 REAL(wp), POINTER, DIMENSION(:,:,:) :: fediss3, fescav3 REAL(wp), POINTER, DIMENSION(:,:,:) :: migrazp3, migrazd3, megrazp3, megrazd3, megrazz3 REAL(wp), POINTER, DIMENSION(:,:,:) :: o2sat3, pbsi3, pcal3, remoc3 REAL(wp), POINTER, DIMENSION(:,:,:) :: pnlimj3, pnlimn3, pnlimfe3, pdlimj3, pdlimn3, pdlimfe3, pdlimsi3 # endif !!--------------------------------------------------------------------- # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined' CALL flush(numout) # endif !! AXY (20/11/14): alter this to report on first MEDUSA call !! IF( kt == nit000 ) THEN IF( kt == nittrc000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model' IF(lwp) WRITE(numout,*) ' ~~~~~~~' IF(lwp) WRITE(numout,*) ' kt =',kt ENDIF !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes ibenthic = 1 !! not sure what this is for; it's not used anywhere; commenting out !! fbodn(:,:) = 0.e0 !! IF( ln_diatrc ) THEN !! blank 2D diagnostic array trc2d(:,:,:) = 0.e0 !! !! blank 3D diagnostic array trc3d(:,:,:,:) = 0.e0 ENDIF !!---------------------------------------------------------------------- !! b0 is present for debugging purposes; using b0 = 0 sets the tendency !! terms of all biological equations to 0. !!---------------------------------------------------------------------- !! !! AXY (03/09/14): probably not the smartest move ever, but it'll fit !! the bill for now; another item on the things-to-sort- !! out-in-the-future list ... # if defined key_kill_medusa b0 = 0. # else b0 = 1. # endif !!---------------------------------------------------------------------- !! fast detritus ballast scheme (0 = no; 1 = yes) !! alternative to ballast scheme is same scheme but with no ballast !! protection (not dissimilar to Martin et al., 1987) !!---------------------------------------------------------------------- !! iball = 1 !!---------------------------------------------------------------------- !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output !! these should *only* be used in 1D since they give comprehensive !! output for ecological functions in the model; primarily used in !! debugging !!---------------------------------------------------------------------- !! idf = 0 !! !! timer mechanism if (kt/120*120.eq.kt) then idfval = 1 else idfval = 0 endif !!---------------------------------------------------------------------- !! blank fast-sinking detritus 2D fields !!---------------------------------------------------------------------- !! ffastn(:,:) = 0.0 !! organic nitrogen ffastsi(:,:) = 0.0 !! biogenic silicon ffastfe(:,:) = 0.0 !! organic iron ffastc(:,:) = 0.0 !! organic carbon ffastca(:,:) = 0.0 !! biogenic calcium carbonate !! fsedn(:,:) = 0.0 !! Seafloor flux of N fsedsi(:,:) = 0.0 !! Seafloor flux of Si fsedfe(:,:) = 0.0 !! Seafloor flux of Fe fsedc(:,:) = 0.0 !! Seafloor flux of C fsedca(:,:) = 0.0 !! Seafloor flux of CaCO3 !! fregenfast(:,:) = 0.0 !! integrated N regeneration (fast detritus) fregenfastsi(:,:) = 0.0 !! integrated Si regeneration (fast detritus) # if defined key_roam fregenfastc(:,:) = 0.0 !! integrated C regeneration (fast detritus) # endif !! fccd(:,:) = 0.0 !! last depth level before CCD !!---------------------------------------------------------------------- !! blank nutrient/flux inventories !!---------------------------------------------------------------------- !! fflx_n(:,:) = 0.0 !! nitrogen flux total fflx_si(:,:) = 0.0 !! silicon flux total fflx_fe(:,:) = 0.0 !! iron flux total fifd_n(:,:) = 0.0 !! nitrogen fast detritus production fifd_si(:,:) = 0.0 !! silicon fast detritus production fifd_fe(:,:) = 0.0 !! iron fast detritus production fofd_n(:,:) = 0.0 !! nitrogen fast detritus remineralisation fofd_si(:,:) = 0.0 !! silicon fast detritus remineralisation fofd_fe(:,:) = 0.0 !! iron fast detritus remineralisation # if defined key_roam fflx_c(:,:) = 0.0 !! carbon flux total fflx_a(:,:) = 0.0 !! alkalinity flux total fflx_o2(:,:) = 0.0 !! oxygen flux total ftot_c(:,:) = 0.0 !! carbon inventory ftot_a(:,:) = 0.0 !! alkalinity inventory ftot_o2(:,:) = 0.0 !! oxygen inventory fifd_c(:,:) = 0.0 !! carbon fast detritus production fifd_a(:,:) = 0.0 !! alkalinity fast detritus production fifd_o2(:,:) = 0.0 !! oxygen fast detritus production fofd_c(:,:) = 0.0 !! carbon fast detritus remineralisation fofd_a(:,:) = 0.0 !! alkalinity fast detritus remineralisation fofd_o2(:,:) = 0.0 !! oxygen fast detritus remineralisation !! fnit_prod(:,:) = 0.0 !! (organic) nitrogen production fnit_cons(:,:) = 0.0 !! (organic) nitrogen consumption fsil_prod(:,:) = 0.0 !! (inorganic) silicon production fsil_cons(:,:) = 0.0 !! (inorganic) silicon consumption fcar_prod(:,:) = 0.0 !! (organic) carbon production fcar_cons(:,:) = 0.0 !! (organic) carbon consumption !! foxy_prod(:,:) = 0.0 !! oxygen production foxy_cons(:,:) = 0.0 !! oxygen consumption foxy_anox(:,:) = 0.0 !! unrealised oxygen consumption !! # endif ftot_n(:,:) = 0.0 !! N inventory ftot_si(:,:) = 0.0 !! Si inventory ftot_fe(:,:) = 0.0 !! Fe inventory ftot_pn(:,:) = 0.0 !! integrated non-diatom phytoplankton ftot_pd(:,:) = 0.0 !! integrated diatom phytoplankton ftot_zmi(:,:) = 0.0 !! integrated microzooplankton ftot_zme(:,:) = 0.0 !! integrated mesozooplankton ftot_det(:,:) = 0.0 !! integrated slow detritus, nitrogen ftot_dtc(:,:) = 0.0 !! integrated slow detritus, carbon !! fzmi_i(:,:) = 0.0 !! material grazed by microzooplankton fzmi_o(:,:) = 0.0 !! ... sum of fate of this material fzme_i(:,:) = 0.0 !! material grazed by mesozooplankton fzme_o(:,:) = 0.0 !! ... sum of fate of this material !! f_sbenin_n(:,:) = 0.0 !! slow detritus N -> benthic pool f_sbenin_fe(:,:) = 0.0 !! slow detritus Fe -> benthic pool f_sbenin_c(:,:) = 0.0 !! slow detritus C -> benthic pool f_fbenin_n(:,:) = 0.0 !! fast detritus N -> benthic pool f_fbenin_fe(:,:) = 0.0 !! fast detritus Fe -> benthic pool f_fbenin_si(:,:) = 0.0 !! fast detritus Si -> benthic pool f_fbenin_c(:,:) = 0.0 !! fast detritus C -> benthic pool f_fbenin_ca(:,:) = 0.0 !! fast detritus Ca -> benthic pool !! f_benout_n(:,:) = 0.0 !! benthic N pool -> dissolved f_benout_fe(:,:) = 0.0 !! benthic Fe pool -> dissolved f_benout_si(:,:) = 0.0 !! benthic Si pool -> dissolved f_benout_c(:,:) = 0.0 !! benthic C pool -> dissolved f_benout_ca(:,:) = 0.0 !! benthic Ca pool -> dissolved !! f_benout_lyso_ca(:,:) = 0.0 !! benthic Ca pool -> dissolved (when it shouldn't!) !! f_runoff(:,:) = 0.0 !! riverine runoff f_riv_n(:,:) = 0.0 !! riverine N input f_riv_si(:,:) = 0.0 !! riverine Si input f_riv_c(:,:) = 0.0 !! riverine C input f_riv_alk(:,:) = 0.0 !! riverine alk input !! !! Jpalm -- 06-03-2017 -- Forgotten var to init f_omarg(:,:) = 0.0 !! f_omcal(:,:) = 0.0 xFree(:,:) = 0.0 !! state variables for iron-ligand system fcomm_resp(:,:) = 0.0 fprn_ml(:,:) = 0.0 !! mixed layer PP diagnostics fprd_ml(:,:) = 0.0 !! mixed layer PP diagnostics !! fslownflux(:,:) = 0.0 fslowcflux(:,:) = 0.0 !! !! allocate and initiate 2D diag !! ----------------------------- !! Juju :: add kt condition !! IF ( lk_iomput .AND. .NOT. ln_diatrc ) THEN !! if ( kt == nittrc000 ) CALL trc_nam_iom_medusa !! initialise iom_use test !! CALL wrk_alloc( jpi, jpj, zw2d ) zw2d(:,:) = 0.0 !! IF ( med_diag%PRN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fprn2d ) fprn2d(:,:) = 0.0 !! ENDIF IF ( med_diag%MPN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdpn2d ) fdpn2d(:,:) = 0.0 !! ENDIF IF ( med_diag%PRD%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fprd2d ) fprd2d(:,:) = 0.0 !! ENDIF IF( med_diag%MPD%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdpd2d ) fdpd2d(:,:) = 0.0 !! ENDIF IF( med_diag%OPAL%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fprds2d ) fprds2d(:,:) = 0.0 !! ENDIF IF( med_diag%OPALDISS%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fsdiss2d ) fsdiss2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMIPn%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgmipn2d ) fgmipn2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMID%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgmid2d ) fgmid2d(:,:) = 0.0 !! ENDIF IF( med_diag%MZMI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdzmi2d ) fdzmi2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMEPN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgmepn2d ) fgmepn2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMEPD%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgmepd2d ) fgmepd2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMEZMI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgmezmi2d ) fgmezmi2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMED%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgmed2d ) fgmed2d(:,:) = 0.0 !! ENDIF IF( med_diag%MZME%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdzme2d ) fdzme2d(:,:) = 0.0 !! ENDIF IF( med_diag%DETN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fslown2d ) fslown2d(:,:) = 0.0 !! ENDIF IF( med_diag%MDET%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdd2d ) fdd2d(:,:) = 0.0 !! ENDIF IF( med_diag%AEOLIAN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ffetop2d ) ffetop2d(:,:) = 0.0 !! ENDIF IF( med_diag%BENTHIC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ffebot2d ) ffebot2d(:,:) = 0.0 !! ENDIF IF( med_diag%SCAVENGE%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ffescav2d ) ffescav2d(:,:) = 0.0 !! ENDIF IF( med_diag%PN_JLIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fjln2d ) fjln2d(:,:) = 0.0 !! ENDIF IF( med_diag%PN_NLIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fnln2d ) fnln2d(:,:) = 0.0 !! ENDIF IF( med_diag%PN_FELIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ffln2d ) ffln2d(:,:) = 0.0 !! ENDIF IF( med_diag%PD_JLIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fjld2d ) fjld2d(:,:) = 0.0 !! ENDIF IF( med_diag%PD_NLIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fnld2d ) fnld2d(:,:) = 0.0 !! ENDIF IF( med_diag%PD_FELIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ffld2d ) ffld2d(:,:) = 0.0 !! ENDIF IF( med_diag%PD_SILIM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fsld2d2 ) fsld2d2(:,:) = 0.0 !! ENDIF IF( med_diag%PDSILIM2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fsld2d ) fsld2d(:,:) = 0.0 !! ENDIF !! !! skip SDT_XXXX diagnostics here !! IF( med_diag%TOTREG_N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fregen2d ) fregen2d(:,:) = 0.0 !! ENDIF IF( med_diag%TOTRG_SI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fregensi2d ) fregensi2d(:,:) = 0.0 !! ENDIF !! !! skip REG_XXXX diagnostics here !! IF( med_diag%FASTN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ftempn2d ) ftempn2d(:,:) = 0.0 !! ENDIF IF( med_diag%FASTSI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ftempsi2d ) ftempsi2d(:,:) = 0.0 !! ENDIF IF( med_diag%FASTFE%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ftempfe2d ) ftempfe2d(:,:) = 0.0 !! ENDIF IF( med_diag%FASTC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ftempc2d ) ftempc2d(:,:) = 0.0 !! ENDIF IF( med_diag%FASTCA%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ftempca2d ) ftempca2d(:,:) = 0.0 !! ENDIF !! !! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here !! IF( med_diag%REMINN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, freminn2d ) freminn2d(:,:) = 0.0 !! ENDIF IF( med_diag%REMINSI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, freminsi2d ) freminsi2d(:,:) = 0.0 !! ENDIF IF( med_diag%REMINFE%dgsave ) THEN CALL wrk_alloc( jpi, jpj, freminfe2d ) freminfe2d(:,:) = 0.0 !! ENDIF IF( med_diag%REMINC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, freminc2d ) freminc2d(:,:) = 0.0 !! ENDIF IF( med_diag%REMINCA%dgsave ) THEN CALL wrk_alloc( jpi, jpj, freminca2d ) freminca2d(:,:) = 0.0 !! ENDIF # if defined key_roam !! !! skip SEAFLRXX, MED_XXXX, INTFLX_XX, INT_XX, ML_XXX, OCAL_XXX, FE_XXXX, MED_XZE, WIND diagnostics here !! IF( med_diag%RR_0100%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ffastca2d ) ffastca2d(:,:) = 0.0 !! ENDIF IF( med_diag%ATM_PCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_pco2a2d ) f_pco2a2d(:,:) = 0.0 !! ENDIF !! !! skip OCN_PH diagnostic here !! IF( med_diag%OCN_PCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_pco2w2d ) f_pco2w2d(:,:) = 0.0 !! ENDIF !! !! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here !! IF( med_diag%CO2FLUX%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_co2flux2d ) f_co2flux2d(:,:) = 0.0 !! ENDIF !! !! skip OM_XXX diagnostics here !! IF( med_diag%TCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_TDIC2d ) f_TDIC2d(:,:) = 0.0 !! ENDIF IF( med_diag%TALK%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_TALK2d ) f_TALK2d(:,:) = 0.0 !! ENDIF IF( med_diag%KW660%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_kw6602d ) f_kw6602d(:,:) = 0.0 !! ENDIF IF( med_diag%ATM_PP0%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_pp02d ) f_pp02d(:,:) = 0.0 !! ENDIF IF( med_diag%O2FLUX%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_o2flux2d ) f_o2flux2d(:,:) = 0.0 !! ENDIF IF( med_diag%O2SAT%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_o2sat2d ) f_o2sat2d(:,:) = 0.0 !! ENDIF !! !! skip XXX_CCD diagnostics here !! IF( med_diag%SFR_OCAL%dgsave ) THEN CALL wrk_alloc( jpi, jpj, sfr_ocal2d ) sfr_ocal2d(:,:) = 0.0 !! ENDIF IF( med_diag%SFR_OARG%dgsave ) THEN CALL wrk_alloc( jpi, jpj, sfr_oarg2d ) sfr_oarg2d(:,:) = 0.0 !! ENDIF !! !! skip XX_PROD, XX_CONS, O2_ANOX, RR_XXXX diagnostics here !! IF( med_diag%IBEN_N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, iben_n2d ) iben_n2d(:,:) = 0.0 !! ENDIF IF( med_diag%IBEN_FE%dgsave ) THEN CALL wrk_alloc( jpi, jpj, iben_fe2d ) iben_fe2d(:,:) = 0.0 !! ENDIF IF( med_diag%IBEN_C%dgsave ) THEN CALL wrk_alloc( jpi, jpj, iben_c2d ) iben_c2d(:,:) = 0.0 !! ENDIF IF( med_diag%IBEN_SI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, iben_si2d ) iben_si2d(:,:) = 0.0 !! ENDIF IF( med_diag%IBEN_CA%dgsave ) THEN CALL wrk_alloc( jpi, jpj, iben_ca2d ) iben_ca2d(:,:) = 0.0 !! ENDIF IF( med_diag%OBEN_N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, oben_n2d ) oben_n2d(:,:) = 0.0 !! ENDIF IF( med_diag%OBEN_FE%dgsave ) THEN CALL wrk_alloc( jpi, jpj, oben_fe2d ) oben_fe2d(:,:) = 0.0 !! ENDIF IF( med_diag%OBEN_C%dgsave ) THEN CALL wrk_alloc( jpi, jpj, oben_c2d ) oben_c2d(:,:) = 0.0 !! ENDIF IF( med_diag%OBEN_SI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, oben_si2d ) oben_si2d(:,:) = 0.0 !! ENDIF IF( med_diag%OBEN_CA%dgsave ) THEN CALL wrk_alloc( jpi, jpj, oben_ca2d ) oben_ca2d(:,:) = 0.0 !! ENDIF !! !! skip BEN_XX diagnostics here !! IF( med_diag%RIV_N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, rivn2d ) rivn2d(:,:) = 0.0 !! ENDIF IF( med_diag%RIV_SI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, rivsi2d ) rivsi2d(:,:) = 0.0 !! ENDIF IF( med_diag%RIV_C%dgsave ) THEN CALL wrk_alloc( jpi, jpj, rivc2d ) rivc2d(:,:) = 0.0 !! ENDIF IF( med_diag%RIV_ALK%dgsave ) THEN CALL wrk_alloc( jpi, jpj, rivalk2d ) rivalk2d(:,:) = 0.0 !! ENDIF IF( med_diag%DETC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fslowc2d ) fslowc2d(:,:) = 0.0 !! ENDIF !! !! skip SDC_XXXX, INVTXXX diagnostics here !! IF( med_diag%LYSO_CA%dgsave ) THEN CALL wrk_alloc( jpi, jpj, lyso_ca2d ) lyso_ca2d(:,:) = 0.0 !! ENDIF !! !! skip COM_RESP diagnostic here !! IF( med_diag%PN_LLOSS%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdpn22d ) fdpn22d(:,:) = 0.0 !! ENDIF IF( med_diag%PD_LLOSS%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdpd22d ) fdpd22d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_LLOSS%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdzmi22d ) fdzmi22d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_LLOSS%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fdzme22d ) fdzme22d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_MES_N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zimesn2d ) zimesn2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_MES_D%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zimesd2d ) zimesd2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_MES_C%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zimesc2d ) zimesc2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_MESDC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zimesdc2d ) zimesdc2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_EXCR%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ziexcr2d ) ziexcr2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_RESP%dgsave ) THEN CALL wrk_alloc( jpi, jpj, ziresp2d ) ziresp2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZI_GROW%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zigrow2d ) zigrow2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_MES_N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zemesn2d ) zemesn2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_MES_D%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zemesd2d ) zemesd2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_MES_C%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zemesc2d ) zemesc2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_MESDC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zemesdc2d ) zemesdc2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_EXCR%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zeexcr2d ) zeexcr2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_RESP%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zeresp2d ) zeresp2d(:,:) = 0.0 !! ENDIF IF( med_diag%ZE_GROW%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zegrow2d ) zegrow2d(:,:) = 0.0 !! ENDIF IF( med_diag%MDETC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, mdetc2d ) mdetc2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMIDC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, gmidc2d ) gmidc2d(:,:) = 0.0 !! ENDIF IF( med_diag%GMEDC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, gmedc2d ) gmedc2d(:,:) = 0.0 !! ENDIF !! !! skip INT_XXX diagnostics here !! IF (jdms .eq. 1) THEN IF( med_diag%DMS_SURF%dgsave ) THEN CALL wrk_alloc( jpi, jpj, dms_surf2d ) dms_surf2d(:,:) = 0.0 !! ENDIF IF( med_diag%DMS_ANDR%dgsave ) THEN CALL wrk_alloc( jpi, jpj, dms_andr2d ) dms_andr2d(:,:) = 0.0 !! ENDIF IF( med_diag%DMS_SIMO%dgsave ) THEN CALL wrk_alloc( jpi, jpj, dms_simo2d ) dms_simo2d(:,:) = 0.0 !! ENDIF IF( med_diag%DMS_ARAN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, dms_aran2d ) dms_aran2d(:,:) = 0.0 !! ENDIF IF( med_diag%DMS_HALL%dgsave ) THEN CALL wrk_alloc( jpi, jpj, dms_hall2d ) dms_hall2d(:,:) = 0.0 !! ENDIF IF( med_diag%DMS_ANDM%dgsave ) THEN CALL wrk_alloc( jpi, jpj, dms_andm2d ) dms_andm2d(:,:) = 0.0 !! ENDIF ENDIF !! !! AXY (24/11/16): extra MOCSY diagnostics, 2D IF( med_diag%ATM_XCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_xco2a_2d ) f_xco2a_2d(:,:) = 0.0 !! ENDIF IF( med_diag%OCN_FCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_fco2w_2d ) f_fco2w_2d(:,:) = 0.0 !! ENDIF IF( med_diag%ATM_FCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_fco2a_2d ) f_fco2a_2d(:,:) = 0.0 !! ENDIF IF( med_diag%OCN_RHOSW%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d ) f_ocnrhosw_2d(:,:) = 0.0 !! ENDIF IF( med_diag%OCN_SCHCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d ) f_ocnschco2_2d(:,:) = 0.0 !! ENDIF IF( med_diag%OCN_KWCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d ) f_ocnkwco2_2d(:,:) = 0.0 !! ENDIF IF( med_diag%OCN_K0%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_ocnk0_2d ) f_ocnk0_2d(:,:) = 0.0 !! ENDIF IF( med_diag%CO2STARAIR%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_co2starair_2d ) f_co2starair_2d(:,:) = 0.0 !! ENDIF IF( med_diag%OCN_DPCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d ) f_ocndpco2_2d(:,:) = 0.0 !! ENDIF # endif IF( med_diag%TPP3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, tpp3d ) tpp3d(:,:,:) = 0.0 !! ENDIF IF( med_diag%DETFLUX3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, detflux3d ) detflux3d(:,:,:) = 0.0 !! ENDIF IF( med_diag%REMIN3N%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, remin3dn ) remin3dn(:,:,:) = 0.0 !! ENDIF !! !! AXY (10/11/16): CMIP6 diagnostics, 2D !! JPALM -- 17-11-16 -- put fgco2 alloc out of diag request !! needed for coupling/passed through restart !! IF( med_diag%FGCO2%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fgco2 ) fgco2(:,:) = 0.0 !! !! ENDIF IF( med_diag%INTDISSIC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, intdissic ) intdissic(:,:) = 0.0 !! ENDIF IF( med_diag%INTDISSIN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, intdissin ) intdissin(:,:) = 0.0 !! ENDIF IF( med_diag%INTDISSISI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, intdissisi ) intdissisi(:,:) = 0.0 !! ENDIF IF( med_diag%INTTALK%dgsave ) THEN CALL wrk_alloc( jpi, jpj, inttalk ) inttalk(:,:) = 0.0 !! ENDIF IF( med_diag%O2min%dgsave ) THEN CALL wrk_alloc( jpi, jpj, o2min ) o2min(:,:) = 1.e3 !! set to high value as we're looking for min(o2) ENDIF IF( med_diag%ZO2min%dgsave ) THEN CALL wrk_alloc( jpi, jpj, zo2min ) zo2min(:,:) = 0.0 !! ENDIF IF( med_diag%FBDDTALK%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fbddtalk ) fbddtalk(:,:) = 0.0 !! ENDIF IF( med_diag%FBDDTDIC%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fbddtdic ) fbddtdic(:,:) = 0.0 !! ENDIF IF( med_diag%FBDDTDIFE%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fbddtdife ) fbddtdife(:,:) = 0.0 !! ENDIF IF( med_diag%FBDDTDIN%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fbddtdin ) fbddtdin(:,:) = 0.0 !! ENDIF IF( med_diag%FBDDTDISI%dgsave ) THEN CALL wrk_alloc( jpi, jpj, fbddtdisi ) fbddtdisi(:,:) = 0.0 !! ENDIF !! !! AXY (10/11/16): CMIP6 diagnostics, 3D IF( med_diag%TPPD3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, tppd3 ) tppd3(:,:,:) = 0.0 !! ENDIF IF( med_diag%BDDTALK3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, bddtalk3 ) bddtalk3(:,:,:) = 0.0 !! ENDIF IF( med_diag%BDDTDIC3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, bddtdic3 ) bddtdic3(:,:,:) = 0.0 !! ENDIF IF( med_diag%BDDTDIFE3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, bddtdife3 ) bddtdife3(:,:,:) = 0.0 !! ENDIF IF( med_diag%BDDTDIN3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, bddtdin3 ) bddtdin3(:,:,:) = 0.0 !! ENDIF IF( med_diag%BDDTDISI3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, bddtdisi3 ) bddtdisi3(:,:,:) = 0.0 !! ENDIF IF( med_diag%FD_NIT3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, fd_nit3 ) fd_nit3(:,:,:) = 0.0 !! ENDIF IF( med_diag%FD_SIL3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, fd_sil3 ) fd_sil3(:,:,:) = 0.0 !! ENDIF IF( med_diag%FD_CAR3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, fd_car3 ) fd_car3(:,:,:) = 0.0 !! ENDIF IF( med_diag%FD_CAL3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, fd_cal3 ) fd_cal3(:,:,:) = 0.0 !! ENDIF IF( med_diag%DCALC3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, dcalc3 ) dcalc3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%EXPC3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, expc3 ) expc3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%EXPN3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, expn3 ) expn3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%FEDISS3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, fediss3 ) fediss3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%FESCAV3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, fescav3 ) fescav3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%MIGRAZP3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, migrazp3 ) migrazp3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%MIGRAZD3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, migrazd3 ) migrazd3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%MEGRAZP3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, megrazp3 ) megrazp3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%MEGRAZD3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, megrazd3 ) megrazd3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%MEGRAZZ3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, megrazz3 ) megrazz3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%O2SAT3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, o2sat3 ) o2sat3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PBSI3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pbsi3 ) pbsi3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PCAL3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pcal3 ) pcal3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%REMOC3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, remoc3 ) remoc3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PNLIMJ3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pnlimj3 ) pnlimj3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PNLIMN3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pnlimn3 ) pnlimn3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PNLIMFE3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pnlimfe3 ) pnlimfe3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PDLIMJ3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pdlimj3 ) pdlimj3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PDLIMN3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pdlimn3 ) pdlimn3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PDLIMFE3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pdlimfe3 ) pdlimfe3(:,:,: ) = 0.0 !! ENDIF IF( med_diag%PDLIMSI3%dgsave ) THEN CALL wrk_alloc( jpi, jpj, jpk, pdlimsi3 ) pdlimsi3(:,:,: ) = 0.0 !! ENDIF ENDIF !! lk_iomput !! # if defined key_axy_nancheck DO jn = 1,jptra !! fq0 = MINVAL(trn(:,:,:,jn)) !! fq1 = MAXVAL(trn(:,:,:,jn)) fq2 = SUM(trn(:,:,:,jn)) !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', & !! & kt, jn, fq0, fq1, fq2 !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR !! and has been replaced here with a specialist routine !! if (fq2 /= fq2 ) then if ( ieee_is_nan( fq2 ) ) then !! there's a NaN here if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:' DO jk = 1,jpk DO jj = 1,jpj DO ji = 1,jpi !! AXY (30/01/14): "isnan" problem on HECTOR !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', & & tmask(ji,jj,jk), ji, jj, jk, jn endif enddo enddo enddo CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' ) endif ENDDO CALL flush(numout) # endif # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked' CALL flush(numout) # endif # if defined key_roam !!---------------------------------------------------------------------- !! calculate atmospheric pCO2 !!---------------------------------------------------------------------- !! !! what's atmospheric pCO2 doing? (data start in 1859) iyr1 = nyear - 1859 + 1 iyr2 = iyr1 + 1 if (iyr1 .le. 1) then !! before 1860 f_xco2a = hist_pco2(1) elseif (iyr2 .ge. 242) then !! after 2099 f_xco2a = hist_pco2(242) else !! just right fq0 = hist_pco2(iyr1) fq1 = hist_pco2(iyr2) fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0) !! AXY (14/06/12): tweaked to make more sense (and be correct) # if defined key_bs_axy_yrlen fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 !! bugfix: for 360d year with HadGEM2-ES forcing # else fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 !! original use of 365 days (not accounting for leap year or 360d year) # endif fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) f_xco2a = fq4 endif # if defined key_axy_pi_co2 !! f_xco2a = 284.725 !! CMIP5 pre-industrial pCO2 f_xco2a = 284.317 !! CMIP6 pre-industrial pCO2 # endif !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', nyear !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day =', real(nsec_day) !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year) !! AXY (29/01/14): remove surplus diagnostics !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0 =', fq0 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1 =', fq1 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2 =', fq2 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3 =', fq3 IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a # endif # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry' IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000 CALL flush(numout) # endif # if defined key_roam !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every !! month (this is hardwired as 960 timesteps but should !! be calculated and done properly !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN !! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN !!============================= !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call : !! we don't want to call on the first time-step of all run submission, !! but only on the very first time-step, and then every month !! So we call on nittrc000 if not restarted run, !! else if one month after last call. !! assume one month is 30d --> 3600*24*30 : 2592000s !! try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt) !! ++ need to pass carb-chem output var through restarts If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN !!---------------------------------------------------------------------- !! Calculate the carbonate chemistry for the whole ocean on the first !! simulation timestep and every month subsequently; the resulting 3D !! field of omega calcite is used to determine the depth of the CCD !!---------------------------------------------------------------------- !! IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt CALL flush(numout) !! blank flags i2_omcal(:,:) = 0 i2_omarg(:,:) = 0 !! loop over 3D space DO jk = 1,jpk DO jj = 2,jpjm1 DO ji = 2,jpim1 !! OPEN wet point IF..THEN loop if (tmask(ji,jj,jk).eq.1) then IF (lk_oasis) THEN f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling ENDIF !! do carbonate chemistry !! fdep2 = fsdept(ji,jj,jk) !! set up level midpoint !! AXY (28/11/16): local seafloor depth !! previously mbathy(ji,jj) - 1, now mbathy(ji,jj) jmbathy = mbathy(ji,jj) !! !! set up required state variables zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity ztmp = tsn(ji,jj,jk,jp_tem) !! temperature zsal = tsn(ji,jj,jk,jp_sal) !! salinity # if defined key_mocsy zsil = max(0.,trn(ji,jj,jk,jpsil)) !! silicic acid zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield # endif !! !! AXY (28/02/14): check input fields if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 3D, ', & tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & ji, ',', jj, ',', jk, ') at time', kt IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', & tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) ztmp = tsb(ji,jj,jk,jp_tem) !! temperature endif if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 3D, ', & tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & ji, ',', jj, ',', jk, ') at time', kt endif !! !! blank input variables not used at this stage (they relate to air-sea flux) f_kw660 = 1.0 f_pp0 = 1.0 !! !! calculate carbonate chemistry at grid cell midpoint # if defined key_mocsy !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate !! chemistry package CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs f_pp0, fdep2, gphit(ji,jj), f_kw660, f_xco2a, 1, & ! inputs f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs f_co2starair, f_co2flux, f_dpco2 ) ! outputs !! f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 -> ueq / kg f_dcf = f_rhosw # else !! AXY (22/06/15): use old PML carbonate chemistry package (the !! MEDUSA-2 default) CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, f_kw660, & ! inputs f_xco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters) ! outputs !! !! AXY (28/02/14): check output fields if (iters .eq. 25) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: 3D ITERS WARNING, ', & iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt endif # endif !! !! store 3D outputs f3_pH(ji,jj,jk) = f_ph f3_h2co3(ji,jj,jk) = f_h2co3 f3_hco3(ji,jj,jk) = f_hco3 f3_co3(ji,jj,jk) = f_co3 f3_omcal(ji,jj,jk) = f_omcal(ji,jj) f3_omarg(ji,jj,jk) = f_omarg(ji,jj) !! !! CCD calculation: calcite if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then if (jk .eq. 1) then f2_ccd_cal(ji,jj) = fdep2 else fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj) fq1 = f3_omcal(ji,jj,jk-1) - 1.0 fq2 = fq1 / (fq0 + tiny(fq0)) fq3 = fdep2 - fsdept(ji,jj,jk-1) fq4 = fq2 * fq3 f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4 endif i2_omcal(ji,jj) = 1 endif if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then !! reached seafloor and still no dissolution; set to seafloor (W-point) f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1) i2_omcal(ji,jj) = 1 endif !! !! CCD calculation: aragonite if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then if (jk .eq. 1) then f2_ccd_arg(ji,jj) = fdep2 else fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj) fq1 = f3_omarg(ji,jj,jk-1) - 1.0 fq2 = fq1 / (fq0 + tiny(fq0)) fq3 = fdep2 - fsdept(ji,jj,jk-1) fq4 = fq2 * fq3 f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4 endif i2_omarg(ji,jj) = 1 endif if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then !! reached seafloor and still no dissolution; set to seafloor (W-point) f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1) i2_omarg(ji,jj) = 1 endif endif ENDDO ENDDO ENDDO ENDIF # endif # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations' CALL flush(numout) # endif !!---------------------------------------------------------------------- !! MEDUSA has unified equation through the water column !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers) !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1 !!---------------------------------------------------------------------- !! !! NOTE: the ordering of the loops below differs from that of some other !! models; looping over the vertical dimension is the outermost loop and !! this complicates some calculations (e.g. storage of vertical fluxes !! that can otherwise be done via a singular variable require 2D fields !! here); however, these issues are relatively easily resolved, but the !! loops CANNOT be reordered without potentially causing code efficiency !! problems (e.g. array indexing means that reordering the loops would !! require skipping between widely-spaced memory location; potentially !! outside those immediately cached) !! !! OPEN vertical loop DO jk = 1,jpk !! OPEN horizontal loops DO jj = 2,jpjm1 DO ji = 2,jpim1 !! OPEN wet point IF..THEN loop if (tmask(ji,jj,jk).eq.1) then !!====================================================================== !! SETUP LOCAL GRID CELL !!====================================================================== !! !!--------------------------------------------------------------------- !! Some notes on grid vertical structure !! - fsdepw(ji,jj,jk) is the depth of the upper surface of level jk !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of level jk !! - fse3t(ji,jj,jk) is the thickness of level jk !!--------------------------------------------------------------------- !! !! AXY (11/12/08): set up level thickness fthk = fse3t(ji,jj,jk) !! AXY (25/02/10): set up level depth (top of level) fdep = fsdepw(ji,jj,jk) !! AXY (01/03/10): set up level depth (bottom of level) fdep1 = fdep + fthk !! AXY (28/11/16): local seafloor depth !! previously mbathy(ji,jj) - 1, now mbathy(ji,jj) jmbathy = mbathy(ji,jj) !! !! set up model tracers !! negative values of state variables are not allowed to !! contribute to the calculated fluxes zchn = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll zchd = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll zphn = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms zphd = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms zpds = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon !! AXY (28/01/10): probably need to take account of chl/biomass connection if (zchn.eq.0.) zphn = 0. if (zchd.eq.0.) zphd = 0. if (zphn.eq.0.) zchn = 0. if (zphd.eq.0.) zchd = 0. !! AXY (23/01/14): duh - why did I forget diatom silicon? if (zpds.eq.0.) zphd = 0. if (zphd.eq.0.) zpds = 0. zzmi = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton zzme = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton zdet = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen zdin = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen zsil = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid zfer = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron" # if defined key_roam zdtc = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity zoxy = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen # if defined key_axy_carbchem && defined key_mocsy zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield # endif !! !! also need physical parameters for gas exchange calculations ztmp = tsn(ji,jj,jk,jp_tem) zsal = tsn(ji,jj,jk,jp_sal) !! !! AXY (28/02/14): check input fields if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & ji, ',', jj, ',', jk, ') at time', kt IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', & tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) ztmp = tsb(ji,jj,jk,jp_tem) !! temperature endif if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & ji, ',', jj, ',', jk, ') at time', kt endif # else zdtc = zdet * xthetad !! implicit detrital carbon # endif # if defined key_debug_medusa if (idf.eq.1) then !! AXY (15/01/10) if (trn(ji,jj,jk,jpdin).lt.0.) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', trn(ji,jj,jk,jpdin) IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', ji, jj, jk, kt endif if (trn(ji,jj,jk,jpsil).lt.0.) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', trn(ji,jj,jk,jpsil) IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', ji, jj, jk, kt endif # if defined key_roam if (trn(ji,jj,jk,jpdic).lt.0.) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', trn(ji,jj,jk,jpdic) IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', ji, jj, jk, kt endif if (trn(ji,jj,jk,jpalk).lt.0.) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', trn(ji,jj,jk,jpalk) IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', ji, jj, jk, kt endif if (trn(ji,jj,jk,jpoxy).lt.0.) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', trn(ji,jj,jk,jpoxy) IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', ji, jj, jk, kt endif # endif endif # endif # if defined key_debug_medusa !! report state variable values if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fthk(',jk,') = ', fthk IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer # if defined key_roam IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy # endif endif # endif # if defined key_debug_medusa if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'dust = ', dust(ji,jj) endif # endif !! sum tracers for inventory checks IF( lk_iomput ) THEN IF ( med_diag%INVTN%dgsave ) THEN ftot_n(ji,jj) = ftot_n(ji,jj) + & (fthk * ( zphn + zphd + zzmi + zzme + zdet + zdin ) ) ENDIF IF ( med_diag%INVTSI%dgsave ) THEN ftot_si(ji,jj) = ftot_si(ji,jj) + & (fthk * ( zpds + zsil ) ) ENDIF IF ( med_diag%INVTFE%dgsave ) THEN ftot_fe(ji,jj) = ftot_fe(ji,jj) + & (fthk * ( xrfn * ( zphn + zphd + zzmi + zzme + zdet ) + zfer ) ) ENDIF # if defined key_roam IF ( med_diag%INVTC%dgsave ) THEN ftot_c(ji,jj) = ftot_c(ji,jj) + & (fthk * ( (xthetapn * zphn) + (xthetapd * zphd) + & (xthetazmi * zzmi) + (xthetazme * zzme) + zdtc + & zdic ) ) ENDIF IF ( med_diag%INVTALK%dgsave ) THEN ftot_a(ji,jj) = ftot_a(ji,jj) + (fthk * ( zalk ) ) ENDIF IF ( med_diag%INVTO2%dgsave ) THEN ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fthk * ( zoxy ) ) ENDIF !! !! AXY (10/11/16): CMIP6 diagnostics IF ( med_diag%INTDISSIC%dgsave ) THEN intdissic(ji,jj) = intdissic(ji,jj) + (fthk * zdic) ENDIF IF ( med_diag%INTDISSIN%dgsave ) THEN intdissin(ji,jj) = intdissin(ji,jj) + (fthk * zdin) ENDIF IF ( med_diag%INTDISSISI%dgsave ) THEN intdissisi(ji,jj) = intdissisi(ji,jj) + (fthk * zsil) ENDIF IF ( med_diag%INTTALK%dgsave ) THEN inttalk(ji,jj) = inttalk(ji,jj) + (fthk * zalk) ENDIF IF ( med_diag%O2min%dgsave ) THEN if ( zoxy < o2min(ji,jj) ) then o2min(ji,jj) = zoxy IF ( med_diag%ZO2min%dgsave ) THEN zo2min(ji,jj) = (fdep + fdep1) / 2. !! layer midpoint ENDIF endif ENDIF # endif ENDIF CALL flush(numout) !!====================================================================== !! LOCAL GRID CELL CALCULATIONS !!====================================================================== !! # if defined key_roam if ( jk .eq. 1 ) then !!---------------------------------------------------------------------- !! Air-sea gas exchange !!---------------------------------------------------------------------- !! !! AXY (17/07/14): zwind_i and zwind_j do not exist in this !! version of NEMO because it does not include !! the SBC changes that our local version has !! for accessing the HadGEM2 forcing; they !! could be added, but an alternative approach !! is to make use of wndm from oce_trc.F90 !! which is wind speed at 10m (which is what !! is required here; this may need to be !! revisited when MEDUSA properly interacts !! with UKESM1 physics !! f_wind = wndm(ji,jj) IF (lk_oasis) THEN f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling ENDIF !! !! AXY (23/06/15): as part of an effort to update the carbonate chemistry !! in MEDUSA, the gas transfer velocity used in the carbon !! and oxygen cycles has been harmonised and is calculated !! by the same function here; this harmonisation includes !! changes to the PML carbonate chemistry scheme so that !! it too makes use of the same gas transfer velocity; the !! preferred parameterisation of this is Wanninkhof (2014), !! option 7 !! # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' CALL flush(numout) # endif CALL gas_transfer( f_wind, 1, 7, & ! inputs f_kw660 ) ! outputs # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' CALL flush(numout) # endif !! !! air pressure (atm); ultimately this will use air pressure at the base !! of the UKESM1 atmosphere !! f_pp0 = 1.0 !! !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) !! IF(lwp) WRITE(numout,*) ' MEDUSA f_wind =', f_wind !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj) !! # if defined key_axy_carbchem # if defined key_mocsy !! !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate !! chemistry package; note that depth is set to !! zero in this call CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs f_pp0, 0.0, gphit(ji,jj), f_kw660, f_xco2a, 1, & ! inputs f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs f_co2starair, f_co2flux, f_dpco2 ) ! outputs !! f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 -> ueq / kg f_dcf = f_rhosw # else iters = 0 !! !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_xco2a, & ! inputs f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters ) ! outputs !! !! AXY (09/01/14): removed iteration and NaN checks; these have !! been moved to trc_co2_medusa together with a !! fudge that amends erroneous values (this is !! intended to be a temporary fudge!); the !! output warnings are retained here so that !! failure position can be determined if (iters .eq. 25) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt endif # endif # else !! AXY (18/04/13): switch off carbonate chemistry calculations; provide !! quasi-sensible alternatives f_ph = 8.1 f_pco2w = f_xco2a f_h2co3 = 0.005 * zdic f_hco3 = 0.865 * zdic f_co3 = 0.130 * zdic f_omcal(ji,jj) = 4. f_omarg(ji,jj) = 2. f_co2flux = 0. f_TDIC = zdic f_TALK = zalk f_dcf = 1.026 f_henry = 1. !! AXY (23/06/15): add in some extra MOCSY diagnostics f_fco2w = f_xco2a f_BetaD = 1. f_rhosw = 1.026 f_opres = 0. f_insitut = ztmp f_pco2atm = f_xco2a f_fco2atm = f_xco2a f_schmidtco2 = 660. f_kwco2 = 0. f_K0 = 0. f_co2starair = f_xco2a f_dpco2 = 0. # endif !! !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux * 86400. / fthk !! !! oxygen (O2); OCMIP-2 code !! AXY (23/06/15): amend input list for oxygen to account for common gas !! transfer velocity !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs !! f_kw660, f_o2flux, f_o2sat ) ! outputs CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy, & ! inputs f_kwo2, f_o2flux, f_o2sat ) ! outputs !! !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness f_o2flux = (1. - fr_i(ji,jj)) * f_o2flux * 86400. / fthk !! !! Jpalm (08-2014) !! DMS surface concentration calculation !! initialy added for UKESM1 model. !! using MET-OFFICE subroutine. !! DMS module only needs Chl concentration and MLD !! to get an aproximate value of DMS concentration. !! air-sea fluxes are calculated by atmospheric chemitry model !! from atm and oc-surface concentrations. !! !! AXY (13/03/15): this is amended to calculate all of the DMS !! estimates examined during UKESM1 (see comments !! in trcdms_medusa.F90) !! !! AXY (25/05/17): amended to additionally pass DIN limitation as well as [DIN]; !! accounts for differences in nutrient half-saturations; changes !! also made in trc_dms_medusa; this permits an additional DMS !! calculation while retaining the existing Anderson one !! IF (jdms .eq. 1) THEN !! !! calculate weighted half-saturation for DIN uptake dms_wtkn = ((zphn * xnln) + (zphd * xnld)) / (zphn + zphd) !! !! feed in correct inputs if (jdms_input .eq. 0) then !! use instantaneous inputs dms_nlim = zdin / (zdin + dms_wtkn) !! CALL trc_dms_medusa( zchn, zchd, & ! inputs hmld(ji,jj), qsr(ji,jj), & ! inputs zdin, dms_nlim, & ! inputs dms_andr, dms_simo, dms_aran, dms_hall, dms_andm ) ! outputs else !! use diel-average inputs dms_nlim = zn_dms_din(ji,jj) / (zn_dms_din(ji,jj) + dms_wtkn) !! CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), & ! inputs zn_dms_din(ji,jj), dms_nlim, & ! inputs dms_andr, dms_simo, dms_aran, dms_hall, dms_andm ) ! outputs endif !! !! assign correct output to variable passed to atmosphere if (jdms_model .eq. 1) then dms_surf = dms_andr elseif (jdms_model .eq. 2) then dms_surf = dms_simo elseif (jdms_model .eq. 3) then dms_surf = dms_aran elseif (jdms_model .eq. 4) then dms_surf = dms_hall elseif (jdms_model .eq. 5) then dms_surf = dms_andm endif !! !! 2D diag through iom_use IF( lk_iomput ) THEN IF( med_diag%DMS_SURF%dgsave ) THEN dms_surf2d(ji,jj) = dms_surf ENDIF IF( med_diag%DMS_ANDR%dgsave ) THEN dms_andr2d(ji,jj) = dms_andr ENDIF IF( med_diag%DMS_SIMO%dgsave ) THEN dms_simo2d(ji,jj) = dms_simo ENDIF IF( med_diag%DMS_ARAN%dgsave ) THEN dms_aran2d(ji,jj) = dms_aran ENDIF IF( med_diag%DMS_HALL%dgsave ) THEN dms_hall2d(ji,jj) = dms_hall ENDIF IF( med_diag%DMS_ANDM%dgsave ) THEN dms_andm2d(ji,jj) = dms_andm ENDIF # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms' CALL flush(numout) # endif ENDIF !! End iom ENDIF !! End DMS Loop !! !! store 2D outputs !! !! JPALM -- 17-11-16 -- put fgco2 out of diag request !! is needed for coupling; pass through restart !! IF( med_diag%FGCO2%dgsave ) THEN !! convert from mol/m2/day to kg/m2/s fgco2(ji,jj) = f_co2flux * fthk * CO2flux_conv !! mmol-C/m3/d -> kg-CO2/m2/s !! ENDIF IF ( lk_iomput ) THEN IF( med_diag%ATM_PCO2%dgsave ) THEN f_pco2a2d(ji,jj) = f_pco2atm ENDIF IF( med_diag%OCN_PCO2%dgsave ) THEN f_pco2w2d(ji,jj) = f_pco2w ENDIF IF( med_diag%CO2FLUX%dgsave ) THEN f_co2flux2d(ji,jj) = f_co2flux * fthk !! mmol/m3/d -> mmol/m2/d ENDIF IF( med_diag%TCO2%dgsave ) THEN f_TDIC2d(ji,jj) = f_TDIC ENDIF IF( med_diag%TALK%dgsave ) THEN f_TALK2d(ji,jj) = f_TALK ENDIF IF( med_diag%KW660%dgsave ) THEN f_kw6602d(ji,jj) = f_kw660 ENDIF IF( med_diag%ATM_PP0%dgsave ) THEN f_pp02d(ji,jj) = f_pp0 ENDIF IF( med_diag%O2FLUX%dgsave ) THEN f_o2flux2d(ji,jj) = f_o2flux ENDIF IF( med_diag%O2SAT%dgsave ) THEN f_o2sat2d(ji,jj) = f_o2sat ENDIF !! AXY (24/11/16): add in extra MOCSY diagnostics IF( med_diag%ATM_XCO2%dgsave ) THEN f_xco2a_2d(ji,jj) = f_xco2a ENDIF IF( med_diag%OCN_FCO2%dgsave ) THEN f_fco2w_2d(ji,jj) = f_fco2w ENDIF IF( med_diag%ATM_FCO2%dgsave ) THEN f_fco2a_2d(ji,jj) = f_fco2atm ENDIF IF( med_diag%OCN_RHOSW%dgsave ) THEN f_ocnrhosw_2d(ji,jj) = f_rhosw ENDIF IF( med_diag%OCN_SCHCO2%dgsave ) THEN f_ocnschco2_2d(ji,jj) = f_schmidtco2 ENDIF IF( med_diag%OCN_KWCO2%dgsave ) THEN f_ocnkwco2_2d(ji,jj) = f_kwco2 ENDIF IF( med_diag%OCN_K0%dgsave ) THEN f_ocnk0_2d(ji,jj) = f_K0 ENDIF IF( med_diag%CO2STARAIR%dgsave ) THEN f_co2starair_2d(ji,jj) = f_co2starair ENDIF IF( med_diag%OCN_DPCO2%dgsave ) THEN f_ocndpco2_2d(ji,jj) = f_dpco2 ENDIF ENDIF !! endif !! End jk = 1 loop within ROAM key !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic IF ( med_diag%O2SAT3%dgsave ) THEN call oxy_sato( ztmp, zsal, f_o2sat3 ) o2sat3(ji, jj, jk) = f_o2sat3 ENDIF # endif if ( jk .eq. 1 ) then !!---------------------------------------------------------------------- !! River inputs !!---------------------------------------------------------------------- !! !! runoff comes in as kg / m2 / s !! used and written out as m3 / m2 / d (= m / d) !! where 1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d !! !! AXY (17/07/14): the compiler doesn't like this line for some reason; !! as MEDUSA doesn't even use runoff for riverine inputs, !! a temporary solution is to switch off runoff entirely !! here; again, this change is one of several that will !! need revisiting once MEDUSA has bedded down in UKESM1; !! particularly so if the land scheme provides information !! concerning nutrient fluxes !! !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. f_runoff(ji,jj) = 0.0 !! !! nutrients are added via rivers to the model in one of two ways: !! 1. via river concentration; i.e. the average nutrient concentration !! of a river water is described by a spatial file, and this is !! multiplied by runoff to give a nutrient flux !! 2. via direct river flux; i.e. the average nutrient flux due to !! rivers is described by a spatial file, and this is simply applied !! as a direct nutrient flux (i.e. it does not relate or respond to !! model runoff) !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and !! alkalinity are derived from continent-scale DIC estimates (Huang et al., !! 2012) and some Arctic river alkalinity estimates (Katya?) !! !! as of 19/07/12, riverine nutrients can now be spread vertically across !! several grid cells rather than just poured into the surface box; this !! block of code is still executed, however, to set up the total amounts !! of nutrient entering via rivers !! !! nitrogen if (jriver_n .eq. 1) then !! river concentration specified; use runoff to calculate input f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) elseif (jriver_n .eq. 2) then !! river flux specified; independent of runoff f_riv_n(ji,jj) = riv_n(ji,jj) endif !! !! silicon if (jriver_si .eq. 1) then !! river concentration specified; use runoff to calculate input f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) elseif (jriver_si .eq. 2) then !! river flux specified; independent of runoff f_riv_si(ji,jj) = riv_si(ji,jj) endif !! !! carbon if (jriver_c .eq. 1) then !! river concentration specified; use runoff to calculate input f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) elseif (jriver_c .eq. 2) then !! river flux specified; independent of runoff f_riv_c(ji,jj) = riv_c(ji,jj) endif !! !! alkalinity if (jriver_alk .eq. 1) then !! river concentration specified; use runoff to calculate input f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) elseif (jriver_alk .eq. 2) then !! river flux specified; independent of runoff f_riv_alk(ji,jj) = riv_alk(ji,jj) endif endif !!---------------------------------------------------------------------- !! Chlorophyll calculations !!---------------------------------------------------------------------- !! !! non-diatoms if (zphn.GT.rsmall) then fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn))) faln = xaln * fthetan else fthetan = 0. faln = 0. endif !! !! diatoms if (zphd.GT.rsmall) then fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd))) fald = xald * fthetad else fthetad = 0. fald = 0. endif # if defined key_debug_medusa !! report biological calculations if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'faln(',jk,') = ', faln IF (lwp) write (numout,*) 'fald(',jk,') = ', fald endif # endif !!---------------------------------------------------------------------- !! Phytoplankton light limitation !!---------------------------------------------------------------------- !! !! It is assumed xpar is the depth-averaged (vertical layer) PAR !! Light limitation (check self-shading) in W/m2 !! !! Note that there is no temperature dependence in phytoplankton !! growth rate or any other function. !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid !! NaNs in case of Phy==0. !! !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and non-diat: !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012 !! !! AXY (16/07/09) !! temperature for new Eppley style phytoplankton growth loc_T = tsn(ji,jj,jk,jp_tem) fun_T = 1.066**(1.0 * loc_T) !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for !phytoplankton !! growth; remin. unaffected fun_Q10 = jq10**((loc_T - 0.0) / 10.0) if (jphy.eq.1) then xvpnT = xvpn * fun_T xvpdT = xvpd * fun_T elseif (jphy.eq.2) then xvpnT = xvpn * fun_Q10 xvpdT = xvpd * fun_Q10 else xvpnT = xvpn xvpdT = xvpd endif !! !! non-diatoms fchn1 = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk)) if (fchn1.GT.rsmall) then fchn = xvpnT / (sqrt(fchn1) + tiny(fchn1)) else fchn = 0. endif fjln = fchn * faln * xpar(ji,jj,jk) !! non-diatom J term fjlim_pn = fjln / xvpnT !! !! diatoms fchd1 = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk)) if (fchd1.GT.rsmall) then fchd = xvpdT / (sqrt(fchd1) + tiny(fchd1)) else fchd = 0. endif fjld = fchd * fald * xpar(ji,jj,jk) !! diatom J term fjlim_pd = fjld / xvpdT # if defined key_debug_medusa !! report phytoplankton light limitation if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld endif # endif !!---------------------------------------------------------------------- !! Phytoplankton nutrient limitation !!---------------------------------------------------------------------- !! !! non-diatoms (N, Fe) fnln = zdin / (zdin + xnln) !! non-diatom Qn term ffln = zfer / (zfer + xfln) !! non-diatom Qf term !! !! diatoms (N, Si, Fe) fnld = zdin / (zdin + xnld) !! diatom Qn term fsld = zsil / (zsil + xsld) !! diatom Qs term ffld = zfer / (zfer + xfld) !! diatom Qf term # if defined key_debug_medusa !! report phytoplankton nutrient limitation if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld IF (lwp) write (numout,*) 'ffln(',jk,') = ', ffln IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld endif # endif !!---------------------------------------------------------------------- !! Primary production (non-diatoms) !! (note: still needs multiplying by phytoplankton concentration) !!---------------------------------------------------------------------- !! if (jliebig .eq. 0) then !! multiplicative nutrient limitation fpnlim = fnln * ffln elseif (jliebig .eq. 1) then !! Liebig Law (= most limiting) nutrient limitation fpnlim = min(fnln, ffln) endif fprn = fjln * fpnlim !!---------------------------------------------------------------------- !! Primary production (diatoms) !! (note: still needs multiplying by phytoplankton concentration) !! !! production here is split between nitrogen production and that of !! silicon; depending upon the "intracellular" ratio of Si:N, model !! diatoms will uptake nitrogen/silicon differentially; this borrows !! from the diatom model of Mongin et al. (2006) !!---------------------------------------------------------------------- !! if (jliebig .eq. 0) then !! multiplicative nutrient limitation fpdlim = fnld * ffld elseif (jliebig .eq. 1) then !! Liebig Law (= most limiting) nutrient limitation fpdlim = min(fnld, ffld) endif !! if (zphd.GT.rsmall .AND. zpds.GT.rsmall) then !! "intracellular" elemental ratios ! fsin = zpds / (zphd + tiny(zphd)) ! fnsi = zphd / (zpds + tiny(zpds)) fsin = 0.0 IF( zphd .GT. rsmall) fsin = zpds / zphd fnsi = 0.0 IF( zpds .GT. rsmall) fnsi = zphd / zpds !! AXY (23/02/10): these next variables derive from Mongin et al. (2003) fsin1 = 3.0 * xsin0 !! = 0.6 fnsi1 = 1.0 / fsin1 !! = 1.667 fnsi2 = 1.0 / xsin0 !! = 5.0 !! !! conditionalities based on ratios !! nitrogen (and iron and carbon) if (fsin.le.xsin0) then fprd = 0.0 fsld2 = 0.0 elseif (fsin.lt.fsin1) then fprd = xuif * ((fsin - xsin0) / (fsin + tiny(fsin))) * (fjld * fpdlim) fsld2 = xuif * ((fsin - xsin0) / (fsin + tiny(fsin))) elseif (fsin.ge.fsin1) then fprd = (fjld * fpdlim) fsld2 = 1.0 endif !! !! silicon if (fsin.lt.fnsi1) then fprds = (fjld * fsld) elseif (fsin.lt.fnsi2) then fprds = xuif * ((fnsi - xnsi0) / (fnsi + tiny(fnsi))) * (fjld * fsld) else fprds = 0.0 endif else fsin = 0.0 fnsi = 0.0 fprd = 0.0 fsld2 = 0.0 fprds = 0.0 endif # if defined key_debug_medusa !! report phytoplankton growth (including diatom silicon submodel) if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fsin(',jk,') = ', fsin IF (lwp) write (numout,*) 'fnsi(',jk,') = ', fnsi IF (lwp) write (numout,*) 'fsld2(',jk,') = ', fsld2 IF (lwp) write (numout,*) 'fprn(',jk,') = ', fprn IF (lwp) write (numout,*) 'fprd(',jk,') = ', fprd IF (lwp) write (numout,*) 'fprds(',jk,') = ', fprds endif # endif !!---------------------------------------------------------------------- !! Mixed layer primary production !! this block calculates the amount of primary production that occurs !! within the upper mixed layer; this allows the separate diagnosis !! of "sub-surface" primary production; it does assume that short- !! term variability in mixed layer depth doesn't mess with things !! though !!---------------------------------------------------------------------- !! if (fdep1.le.hmld(ji,jj)) then !! this level is entirely in the mixed layer fq0 = 1.0 elseif (fdep.ge.hmld(ji,jj)) then !! this level is entirely below the mixed layer fq0 = 0.0 else !! this level straddles the mixed layer fq0 = (hmld(ji,jj) - fdep) / fthk endif !! fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn * zphn * fthk * fq0) fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd * zphd * fthk * fq0) !!---------------------------------------------------------------------- !! Vertical Integral -- !!---------------------------------------------------------------------- ftot_pn(ji,jj) = ftot_pn(ji,jj) + (zphn * fthk) !! vertical integral non-diatom phytoplankton ftot_pd(ji,jj) = ftot_pd(ji,jj) + (zphd * fthk) !! vertical integral diatom phytoplankton ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi * fthk) !! vertical integral microzooplankton ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme * fthk) !! vertical integral mesozooplankton ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet * fthk) !! vertical integral slow detritus, nitrogen ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc * fthk) !! vertical integral slow detritus, carbon !!---------------------------------------------------------------------- !! More chlorophyll calculations !!---------------------------------------------------------------------- !! !! frn = (xthetam / fthetan) * (fprn / (fthetan * xpar(ji,jj,jk))) !! frd = (xthetam / fthetad) * (fprd / (fthetad * xpar(ji,jj,jk))) frn = (xthetam * fchn * fnln * ffln ) / (fthetan + tiny(fthetan)) !! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid !! limitation, is used in the following line to regulate chlorophyll !! growth in a manner that is inconsistent with its use in the regulation !! of biomass growth; the Mongin term term used in growth is more complex !! than the simple multiplicative function used below !! frd = (xthetam * fchd * fnld * ffld * fsld) / (fthetad + tiny(fthetad)) !! AXY (12/05/09): this replacement line uses the new variable, fsld2, to !! regulate chlorophyll growth frd = (xthetamd * fchd * fnld * ffld * fsld2) / (fthetad + tiny(fthetad)) # if defined key_debug_medusa !! report chlorophyll calculations if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad IF (lwp) write (numout,*) 'frn(',jk,') = ', frn IF (lwp) write (numout,*) 'frd(',jk,') = ', frd endif # endif !!---------------------------------------------------------------------- !! Zooplankton Grazing !! this code supplements the base grazing model with one that !! considers the C:N ratio of grazed food and balances this against !! the requirements of zooplankton growth; this model is derived !! from that of Anderson & Pondaven (2003) !! !! the current version of the code assumes a fixed C:N ratio for !! detritus (in contrast to Anderson & Pondaven, 2003), though the !! full equations are retained for future extension !!---------------------------------------------------------------------- !! !!---------------------------------------------------------------------- !! Microzooplankton first !!---------------------------------------------------------------------- !! fmi1 = (xkmi * xkmi) + (xpmipn * zphn * zphn) + (xpmid * zdet * zdet) fmi = xgmi * zzmi / fmi1 fgmipn = fmi * xpmipn * zphn * zphn !! grazing on non-diatoms fgmid = fmi * xpmid * zdet * zdet !! grazing on detrital nitrogen # if defined key_roam fgmidc = rsmall !acc IF ( zdet .GT. rsmall ) fgmidc = (zdtc / (zdet + tiny(zdet))) * fgmid !! grazing on detrital carbon # else !! AXY (26/11/08): implicit detrital carbon change fgmidc = xthetad * fgmid !! grazing on detrital carbon # endif !! !! which translates to these incoming N and C fluxes finmi = (1.0 - xphi) * (fgmipn + fgmid) ficmi = (1.0 - xphi) * ((xthetapn * fgmipn) + fgmidc) !! !! the ideal food C:N ratio for microzooplankton !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 fstarmi = (xbetan * xthetazmi) / (xbetac * xkc) !! !! process these to determine proportioning of grazed N and C !! (since there is no explicit consideration of respiration, !! only growth and excretion are calculated here) fmith = (ficmi / (finmi + tiny(finmi))) if (fmith.ge.fstarmi) then fmigrow = xbetan * finmi fmiexcr = 0.0 else fmigrow = (xbetac * xkc * ficmi) / xthetazmi fmiexcr = ficmi * ((xbetan / (fmith + tiny(fmith))) - ((xbetac * xkc) / xthetazmi)) endif # if defined key_roam fmiresp = (xbetac * ficmi) - (xthetazmi * fmigrow) # endif # if defined key_debug_medusa !! report microzooplankton grazing if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fmi1(',jk,') = ', fmi1 IF (lwp) write (numout,*) 'fmi(',jk,') = ', fmi IF (lwp) write (numout,*) 'fgmipn(',jk,') = ', fgmipn IF (lwp) write (numout,*) 'fgmid(',jk,') = ', fgmid IF (lwp) write (numout,*) 'fgmidc(',jk,') = ', fgmidc IF (lwp) write (numout,*) 'finmi(',jk,') = ', finmi IF (lwp) write (numout,*) 'ficmi(',jk,') = ', ficmi IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi IF (lwp) write (numout,*) 'fmith(',jk,') = ', fmith IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr # if defined key_roam IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp # endif endif # endif !!---------------------------------------------------------------------- !! Mesozooplankton second !!---------------------------------------------------------------------- !! fme1 = (xkme * xkme) + (xpmepn * zphn * zphn) + (xpmepd * zphd * zphd) + & (xpmezmi * zzmi * zzmi) + (xpmed * zdet * zdet) fme = xgme * zzme / fme1 fgmepn = fme * xpmepn * zphn * zphn !! grazing on non-diatoms fgmepd = fme * xpmepd * zphd * zphd !! grazing on diatoms fgmepds = fsin * fgmepd !! grazing on diatom silicon fgmezmi = fme * xpmezmi * zzmi * zzmi !! grazing on microzooplankton fgmed = fme * xpmed * zdet * zdet !! grazing on detrital nitrogen # if defined key_roam fgmedc = rsmall !acc IF ( zdet .GT. rsmall ) fgmedc = (zdtc / (zdet + tiny(zdet))) * fgmed !! grazing on detrital carbon # else !! AXY (26/11/08): implicit detrital carbon change fgmedc = xthetad * fgmed !! grazing on detrital carbon # endif !! !! which translates to these incoming N and C fluxes finme = (1.0 - xphi) * (fgmepn + fgmepd + fgmezmi + fgmed) ficme = (1.0 - xphi) * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & (xthetazmi * fgmezmi) + fgmedc) !! !! the ideal food C:N ratio for mesozooplankton !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 fstarme = (xbetan * xthetazme) / (xbetac * xkc) !! !! process these to determine proportioning of grazed N and C !! (since there is no explicit consideration of respiration, !! only growth and excretion are calculated here) fmeth = (ficme / (finme + tiny(finme))) if (fmeth.ge.fstarme) then fmegrow = xbetan * finme fmeexcr = 0.0 else fmegrow = (xbetac * xkc * ficme) / xthetazme fmeexcr = ficme * ((xbetan / (fmeth + tiny(fmeth))) - ((xbetac * xkc) / xthetazme)) endif # if defined key_roam fmeresp = (xbetac * ficme) - (xthetazme * fmegrow) # endif # if defined key_debug_medusa !! report mesozooplankton grazing if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fme1(',jk,') = ', fme1 IF (lwp) write (numout,*) 'fme(',jk,') = ', fme IF (lwp) write (numout,*) 'fgmepn(',jk,') = ', fgmepn IF (lwp) write (numout,*) 'fgmepd(',jk,') = ', fgmepd IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi IF (lwp) write (numout,*) 'fgmed(',jk,') = ', fgmed IF (lwp) write (numout,*) 'fgmedc(',jk,') = ', fgmedc IF (lwp) write (numout,*) 'finme(',jk,') = ', finme IF (lwp) write (numout,*) 'ficme(',jk,') = ', ficme IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme IF (lwp) write (numout,*) 'fmeth(',jk,') = ', fmeth IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr # if defined key_roam IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp # endif endif # endif fzmi_i(ji,jj) = fzmi_i(ji,jj) + fthk * ( & fgmipn + fgmid ) fzmi_o(ji,jj) = fzmi_o(ji,jj) + fthk * ( & fmigrow + (xphi * (fgmipn + fgmid)) + fmiexcr + ((1.0 - xbetan) * finmi) ) fzme_i(ji,jj) = fzme_i(ji,jj) + fthk * ( & fgmepn + fgmepd + fgmezmi + fgmed ) fzme_o(ji,jj) = fzme_o(ji,jj) + fthk * ( & fmegrow + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + fmeexcr + ((1.0 - xbetan) * finme) ) !!---------------------------------------------------------------------- !! Plankton metabolic losses !! Linear loss processes assumed to be metabolic in origin !!---------------------------------------------------------------------- !! fdpn2 = xmetapn * zphn fdpd2 = xmetapd * zphd fdpds2 = xmetapd * zpds fdzmi2 = xmetazmi * zzmi fdzme2 = xmetazme * zzme !!---------------------------------------------------------------------- !! Plankton mortality losses !! EKP (26/02/09): phytoplankton hyperbolic mortality term introduced !! to improve performance in gyres !!---------------------------------------------------------------------- !! !! non-diatom phytoplankton if (jmpn.eq.1) fdpn = xmpn * zphn !! linear if (jmpn.eq.2) fdpn = xmpn * zphn * zphn !! quadratic if (jmpn.eq.3) fdpn = xmpn * zphn * & !! hyperbolic (zphn / (xkphn + zphn)) if (jmpn.eq.4) fdpn = xmpn * zphn * & !! sigmoid ((zphn * zphn) / (xkphn + (zphn * zphn))) !! !! diatom phytoplankton if (jmpd.eq.1) fdpd = xmpd * zphd !! linear if (jmpd.eq.2) fdpd = xmpd * zphd * zphd !! quadratic if (jmpd.eq.3) fdpd = xmpd * zphd * & !! hyperbolic (zphd / (xkphd + zphd)) if (jmpd.eq.4) fdpd = xmpd * zphd * & !! sigmoid ((zphd * zphd) / (xkphd + (zphd * zphd))) fdpds = fdpd * fsin !! !! microzooplankton if (jmzmi.eq.1) fdzmi = xmzmi * zzmi !! linear if (jmzmi.eq.2) fdzmi = xmzmi * zzmi * zzmi !! quadratic if (jmzmi.eq.3) fdzmi = xmzmi * zzmi * & !! hyperbolic (zzmi / (xkzmi + zzmi)) if (jmzmi.eq.4) fdzmi = xmzmi * zzmi * & !! sigmoid ((zzmi * zzmi) / (xkzmi + (zzmi * zzmi))) !! !! mesozooplankton if (jmzme.eq.1) fdzme = xmzme * zzme !! linear if (jmzme.eq.2) fdzme = xmzme * zzme * zzme !! quadratic if (jmzme.eq.3) fdzme = xmzme * zzme * & !! hyperbolic (zzme / (xkzme + zzme)) if (jmzme.eq.4) fdzme = xmzme * zzme * & !! sigmoid ((zzme * zzme) / (xkzme + (zzme * zzme))) !!---------------------------------------------------------------------- !! Detritus remineralisation !! Constant or temperature-dependent !!---------------------------------------------------------------------- !! if (jmd.eq.1) then !! temperature-dependent fdd = xmd * fun_T * zdet # if defined key_roam fddc = xmdc * fun_T * zdtc # endif elseif (jmd.eq.2) then !! AXY (16/05/13): add in Q10-based parameterisation (def in nmlst) !! temperature-dependent fdd = xmd * fun_Q10 * zdet # if defined key_roam fddc = xmdc * fun_Q10 * zdtc # endif else !! temperature-independent fdd = xmd * zdet # if defined key_roam fddc = xmdc * zdtc # endif endif !! !! AXY (22/07/09): accelerate detrital remineralisation in the bottom box if ((jk.eq.jmbathy) .and. jsfd.eq.1) then fdd = 1.0 * zdet # if defined key_roam fddc = 1.0 * zdtc # endif endif # if defined key_debug_medusa !! report plankton mortality and remineralisation if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fdpn2(',jk,') = ', fdpn2 IF (lwp) write (numout,*) 'fdpd2(',jk,') = ', fdpd2 IF (lwp) write (numout,*) 'fdpds2(',jk,')= ', fdpds2 IF (lwp) write (numout,*) 'fdzmi2(',jk,')= ', fdzmi2 IF (lwp) write (numout,*) 'fdzme2(',jk,')= ', fdzme2 IF (lwp) write (numout,*) 'fdpn(',jk,') = ', fdpn IF (lwp) write (numout,*) 'fdpd(',jk,') = ', fdpd IF (lwp) write (numout,*) 'fdpds(',jk,') = ', fdpds IF (lwp) write (numout,*) 'fdzmi(',jk,') = ', fdzmi IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme IF (lwp) write (numout,*) 'fdd(',jk,') = ', fdd # if defined key_roam IF (lwp) write (numout,*) 'fddc(',jk,') = ', fddc # endif endif # endif !!---------------------------------------------------------------------- !! Detritus addition to benthos !! If activated, slow detritus in the bottom box will enter the !! benthic pool !!---------------------------------------------------------------------- !! if ((jk.eq.jmbathy) .and. jorgben.eq.1) then !! this is the BOTTOM OCEAN BOX -> into the benthic pool! !! f_sbenin_n(ji,jj) = (zdet * vsed * 86400.) f_sbenin_fe(ji,jj) = (zdet * vsed * 86400. * xrfn) # if defined key_roam f_sbenin_c(ji,jj) = (zdtc * vsed * 86400.) # else f_sbenin_c(ji,jj) = (zdet * vsed * 86400. * xthetad) # endif endif !!---------------------------------------------------------------------- !! Iron chemistry and fractionation !! following the Parekh et al. (2004) scheme adopted by the Met. !! Office, Medusa models total iron but considers "free" and !! ligand-bound forms for the purposes of scavenging (only "free" !! iron can be scavenged !!---------------------------------------------------------------------- !! !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) xFeT = zfer * 1.e3 !! !! calculate fractionation (based on Diat-HadOCC; in turn based on Parekh et al., 2004) xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0 xb2M4ac = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0) !! !! "free" ligand concentration xLgF = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL !! !! ligand-bound iron concentration xFeL = xLgT - xLgF !! !! "free" iron concentration (and convert to mmol Fe / m3) xFeF = (xFeT - xFeL) * 1.e-3 xFree(ji,jj)= xFeF / (zfer + tiny(zfer)) !! !! scavenging of iron (multiple schemes); I'm only really happy with the !! first one at the moment - the others involve assumptions (sometimes !! guessed at by me) that are potentially questionable !! if (jiron.eq.1) then !!---------------------------------------------------------------------- !! Scheme 1: Dutkiewicz et al. (2005) !! This scheme includes a single scavenging term based solely on a !! fixed rate and the availablility of "free" iron !!---------------------------------------------------------------------- !! ffescav = xk_sc_Fe * xFeF ! = mmol/m3/d !! !!---------------------------------------------------------------------- !! !! Mick's code contains a further (optional) implicit "scavenging" of !! iron that sets an upper bound on "free" iron concentration, and !! essentially caps the concentration of total iron as xFeL + "free" !! iron; since the former is constrained by a fixed total ligand !! concentration (= 1.0 umol/m3), and the latter isn't allowed above !! this upper bound, total iron is constrained to a maximum of ... !! !! xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3 !! !! In Mick's code, the actual value of total iron is reset to this !! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't !! our favoured approach to tracer updating here (not least because !! of the leapfrog), so here the amount scavenged is augmented by an !! additional amount that serves to drag total iron back towards that !! expected from this limitation on iron concentration ... !! xmaxFeF = min((xFeF * 1.e3), 0.3) ! = umol/m3 !! !! Here, the difference between current total Fe and (FeL + Fe') is !! calculated and added to the scavenging flux already calculated !! above ... !! fdeltaFe = (xFeT - (xFeL + xmaxFeF)) * 1.e-3 ! = mmol/m3 !! !! This assumes that the "excess" iron is dissipated with a time- !! scale of 1 day; seems reasonable to me ... (famous last words) !! ffescav = ffescav + fdeltaFe ! = mmol/m3/d !! # if defined key_deep_fe_fix !! AXY (17/01/13) !! stop scavenging for iron concentrations below 0.5 umol / m3 !! at depths greater than 1000 m; this aims to end MEDUSA's !! continual loss of iron at depth without impacting things !! at the surface too much; the justification for this is that !! it appears to be what Mick Follows et al. do in their work !! (as evidenced by the iron initial condition they supplied !! me with); to be honest, it looks like Follow et al. do this !! at shallower depths than 1000 m, but I'll stick with this !! for now; I suspect that this seemingly arbitrary approach !! effectively "parameterises" the particle-based scavenging !! rates that other models use (i.e. at depth there are no !! sinking particles, so scavenging stops); it might be fun !! justifying this in a paper though! !! if ((fdep.gt.1000.) .and. (xFeT.lt.0.5)) then ffescav = 0. endif # endif !! elseif (jiron.eq.2) then !!---------------------------------------------------------------------- !! Scheme 2: Moore et al. (2004) !! This scheme includes a single scavenging term that accounts for !! both suspended and sinking particles in the water column; this !! term scavenges total iron rather than "free" iron !!---------------------------------------------------------------------- !! !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) xFeT = zfer * 1.e3 !! !! this has a base scavenging rate (12% / y) which is modified by local !! particle concentration and sinking flux (and dust - but I'm ignoring !! that here for now) and which is accelerated when Fe concentration gets !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as concentrations !! below 0.4 nM (= 0.4 umol/m3 = 0.0004 mmol/m3) !! !! base scavenging rate (0.12 / y) fbase_scav = 0.12 / 365.25 !! !! calculate sinking particle part of scaling factor !! this takes local fast sinking carbon (mmol C / m2 / d) and !! gets it into nmol C / cm3 / s ("rdt" below is the number of seconds in !! a model timestep) !! !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.) !! !! ... actually, re-reading Moore et al.'s equations, it looks like he uses !! his sinking flux directly, without scaling it by time-step or anything, !! so I'll copy this here ... !! fscal_sink = ffastc(ji,jj) * 1.e2 !! !! calculate particle part of scaling factor !! this totals up the carbon in suspended particles (Pn, Pd, Zmi, Zme, D), !! which comes out in mmol C / m3 (= nmol C / cm3), and then multiplies it !! by a magic factor, 0.002, to get it into nmol C / cm2 / s !! fscal_part = ((xthetapn * zphn) + (xthetapd * zphd) + (xthetazmi * zzmi) + & (xthetazme * zzme) + (xthetad * zdet)) * 0.002 !! !! calculate scaling factor for base scavenging rate !! this uses the (now correctly scaled) sinking flux and standing !! particle concentration, divides through by some sort of reference !! value (= 0.0066 nmol C / cm2 / s) and then uses this, or not if its !! too high, to rescale the base scavenging rate !! fscal_scav = fbase_scav * min(((fscal_sink + fscal_part) / 0.0066), 4.0) !! !! the resulting scavenging rate is then scaled further according to the !! local iron concentration (i.e. diminished in low iron regions; enhanced !! in high iron regions; less alone in intermediate iron regions) !! if (xFeT.lt.0.4) then !! !! low iron region !! fscal_scav = fscal_scav * (xFeT / 0.4) !! elseif (xFeT.gt.0.6) then !! !! high iron region !! fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4)) !! else !! !! intermediate iron region: do nothing !! endif !! !! apply the calculated scavenging rate ... !! ffescav = fscal_scav * zfer !! elseif (jiron.eq.3) then !!---------------------------------------------------------------------- !! Scheme 3: Moore et al. (2008) !! This scheme includes a single scavenging term that accounts for !! sinking particles in the water column, and includes organic C, !! biogenic opal, calcium carbonate and dust in this (though the !! latter is ignored here until I work out what units the incoming !! "dust" flux is in); this term scavenges total iron rather than !! "free" iron !!---------------------------------------------------------------------- !! !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) xFeT = zfer * 1.e3 !! !! this has a base scavenging rate which is modified by local !! particle sinking flux (including dust - but I'm ignoring that !! here for now) and which is accelerated when Fe concentration !! is > 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as !! concentrations < 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3) !! !! base scavenging rate (Fe_b in paper; units may be wrong there) fbase_scav = 0.00384 ! (ng)^-1 cm !! !! calculate sinking particle part of scaling factor; this converts !! mmol / m2 / d fluxes of organic carbon, silicon and calcium !! carbonate into ng / cm2 / s fluxes; it is assumed here that the !! mass conversions simply consider the mass of the main element !! (C, Si and Ca) and *not* the mass of the molecules that they are !! part of; Moore et al. (2008) is unclear on the conversion that !! should be used !! !! milli -> nano; mol -> gram; /m2 -> /cm2; /d -> /s fscal_csink = (ffastc(ji,jj) * 1.e6 * xmassc * 1.e-4 / 86400.) ! ng C / cm2 / s fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi * 1.e-4 / 86400.) ! ng Si / cm2 / s fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca * 1.e-4 / 86400.) ! ng Ca / cm2 / s !! !! sum up these sinking fluxes and convert to ng / cm by dividing !! through by a sinking rate of 100 m / d = 1.157 cm / s fscal_sink = ((fscal_csink * 6.) + fscal_sisink + fscal_casink) / & (100. * 1.e3 / 86400) ! ng / cm !! !! now calculate the scavenging rate based upon the base rate and !! this particle flux scaling; according to the published units, !! the result actually has *no* units, but as it must be expressed !! per unit time for it to make any sense, I'm assuming a missing !! "per second" fscal_scav = fbase_scav * fscal_sink ! / s !! !! the resulting scavenging rate is then scaled further according to the !! local iron concentration (i.e. diminished in low iron regions; enhanced !! in high iron regions; less alone in intermediate iron regions) !! if (xFeT.lt.0.5) then !! !! low iron region (0.5 instead of the 0.4 in Moore et al., 2004) !! fscal_scav = fscal_scav * (xFeT / 0.5) !! elseif (xFeT.gt.0.6) then !! !! high iron region (functional form different in Moore et al., 2004) !! fscal_scav = fscal_scav + ((xFeT - 0.6) * 0.00904) !! else !! !! intermediate iron region: do nothing !! endif !! !! apply the calculated scavenging rate ... !! ffescav = fscal_scav * zfer !! elseif (jiron.eq.4) then !!---------------------------------------------------------------------- !! Scheme 4: Galbraith et al. (2010) !! This scheme includes two scavenging terms, one for organic, !! particle-based scavenging, and another for inorganic scavenging; !! both terms scavenge "free" iron only !!---------------------------------------------------------------------- !! !! Galbraith et al. (2010) present a more straightforward outline of !! the scheme in Parekh et al. (2005) ... !! !! sinking particulate carbon available for scavenging !! this assumes a sinking rate of 100 m / d (Moore & Braucher, 2008), xCscav1 = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3 !! !! scale by Honeyman et al. (1981) exponent coefficient !! multiply by 1.e-3 to express C flux in g C rather than mg C xCscav2 = (xCscav1 * 1.e-3)**0.58 !! !! multiply by Galbraith et al. (2010) scavenging rate xk_org = 0.5 ! ((g C m/3)^-1) / d xORGscav = xk_org * xCscav2 * xFeF !! !! Galbraith et al. (2010) also include an inorganic bit ... !! !! this occurs at a fixed rate, again based on the availability of !! "free" iron !! !! k_inorg = 1000 d**-1 nmol Fe**-0.5 kg**-0.5 !! !! to implement this here, scale xFeF by 1026 to put in units of !! umol/m3 which approximately equal nmol/kg !! xk_inorg = 1000. ! ((nmol Fe / kg)^1.5) xINORGscav = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e-3 !! !! sum these two terms together ffescav = xORGscav + xINORGscav else !!---------------------------------------------------------------------- !! No Scheme: you coward! !! This scheme puts its head in the sand and eskews any decision about !! how iron is removed from the ocean; prepare to get deluged in iron !! you fool! !!---------------------------------------------------------------------- ffescav = 0. endif !!---------------------------------------------------------------------- !! Other iron cycle processes !!---------------------------------------------------------------------- !! !! aeolian iron deposition if (jk.eq.1) then !! zirondep is in mmol-Fe / m2 / day !! ffetop is in mmol-dissolved-Fe / m3 / day ffetop = zirondep(ji,jj) * xfe_sol / fthk else ffetop = 0.0 endif !! !! seafloor iron addition !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then if ((jk.eq.jmbathy).AND.jk.le.i0500) then !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value !! but apply it everywhere !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 1-37) ffebot = (xfe_sed / fthk) else ffebot = 0.0 endif !! AXY (16/12/09): remove iron addition/removal processes !! For the purposes of the quarter degree run, the iron cycle is being pegged to the !! initial condition supplied by Mick Follows via restoration with a 30 day period; !! iron addition at the seafloor is still permitted with the idea that this extra !! iron will be removed by the restoration away from the source !! ffescav = 0.0 !! ffetop = 0.0 !! ffebot = 0.0 # if defined key_debug_medusa !! report miscellaneous calculations if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'xfe_sol = ', xfe_sol IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass IF (lwp) write (numout,*) 'ffetop(',jk,') = ', ffetop IF (lwp) write (numout,*) 'ffebot(',jk,') = ', ffebot IF (lwp) write (numout,*) 'xFree(',jk,') = ', xFree(ji,jj) IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav endif # endif !!---------------------------------------------------------------------- !! Miscellaneous !!---------------------------------------------------------------------- !! !! diatom frustule dissolution fsdiss = xsdiss * zpds # if defined key_debug_medusa !! report miscellaneous calculations if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fsdiss(',jk,') = ', fsdiss endif # endif !!---------------------------------------------------------------------- !! Slow detritus creation !!---------------------------------------------------------------------- !! this variable integrates the creation of slow sinking detritus !! to allow the split between fast and slow detritus to be !! diagnosed fslown = fdpn + fdzmi + ((1.0 - xfdfrac1) * fdpd) + & ((1.0 - xfdfrac2) * fdzme) + ((1.0 - xbetan) * (finmi + finme)) !! !! this variable records the slow detrital sinking flux at this !! particular depth; it is used in the output of this flux at !! standard depths in the diagnostic outputs; needs to be !! adjusted from per second to per day because of parameter vsed fslownflux(ji,jj) = zdet * vsed * 86400. # if defined key_roam !! !! and the same for detrital carbon fslowc = (xthetapn * fdpn) + (xthetazmi * fdzmi) + & (xthetapd * (1.0 - xfdfrac1) * fdpd) + & (xthetazme * (1.0 - xfdfrac2) * fdzme) + & ((1.0 - xbetac) * (ficmi + ficme)) !! !! this variable records the slow detrital sinking flux at this !! particular depth; it is used in the output of this flux at !! standard depths in the diagnostic outputs; needs to be !! adjusted from per second to per day because of parameter vsed fslowcflux(ji,jj) = zdtc * vsed * 86400. # endif !!---------------------------------------------------------------------- !! Nutrient regeneration !! this variable integrates total nitrogen regeneration down the !! watercolumn; its value is stored and output as a 2D diagnostic; !! the corresponding dissolution flux of silicon (from sources !! other than fast detritus) is also integrated; note that, !! confusingly, the linear loss terms from plankton compartments !! are labelled as fdX2 when one might have expected fdX or fdX1 !!---------------------------------------------------------------------- !! !! nitrogen fregen = (( (xphi * (fgmipn + fgmid)) + & ! messy feeding (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + & ! messy feeding fmiexcr + fmeexcr + fdd + & ! excretion + D remin. fdpn2 + fdpd2 + fdzmi2 + fdzme2) * fthk) ! linear mortality !! !! silicon fregensi = (( fsdiss + ((1.0 - xfdfrac1) * fdpds) + & ! dissolution + non-lin. mortality ((1.0 - xfdfrac3) * fgmepds) + & ! egestion by zooplankton fdpds2) * fthk) ! linear mortality # if defined key_roam !! !! carbon fregenc = (( (xphi * ((xthetapn * fgmipn) + fgmidc)) + & ! messy feeding (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & ! messy feeding (xthetazmi * fgmezmi) + fgmedc)) + & ! messy feeding fmiresp + fmeresp + fddc + & ! respiration + D remin. (xthetapn * fdpn2) + (xthetapd * fdpd2) + & ! linear mortality (xthetazmi * fdzmi2) + (xthetazme * fdzme2)) * fthk) ! linear mortality # endif !!---------------------------------------------------------------------- !! Fast-sinking detritus terms !! "local" variables declared so that conservation can be checked; !! the calculated terms are added to the fast-sinking flux later on !! only after the flux entering this level has experienced some !! remineralisation !! note: these fluxes need to be scaled by the level thickness !!---------------------------------------------------------------------- !! !! nitrogen: diatom and mesozooplankton mortality ftempn = b0 * ((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) !! !! silicon: diatom mortality and grazed diatoms ftempsi = b0 * ((xfdfrac1 * fdpds) + (xfdfrac3 * fgmepds)) !! !! iron: diatom and mesozooplankton mortality ftempfe = b0 * (((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) * xrfn) !! !! carbon: diatom and mesozooplankton mortality ftempc = b0 * ((xfdfrac1 * xthetapd * fdpd) + & (xfdfrac2 * xthetazme * fdzme)) !! # if defined key_roam if (jrratio.eq.0) then !! CaCO3: latitudinally-based fraction of total primary production !! absolute latitude of current grid cell flat = abs(gphit(ji,jj)) !! 0.10 at equator; 0.02 at pole fcaco3 = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0)) elseif (jrratio.eq.1) then !! CaCO3: Ridgwell et al. (2007) submodel, version 1 !! this uses SURFACE omega calcite to regulate rain ratio if (f_omcal(ji,jj).ge.1.0) then fq1 = (f_omcal(ji,jj) - 1.0)**0.81 else fq1 = 0. endif fcaco3 = xridg_r0 * fq1 elseif (jrratio.eq.2) then !! CaCO3: Ridgwell et al. (2007) submodel, version 2 !! this uses FULL 3D omega calcite to regulate rain ratio if (f3_omcal(ji,jj,jk).ge.1.0) then fq1 = (f3_omcal(ji,jj,jk) - 1.0)**0.81 else fq1 = 0. endif fcaco3 = xridg_r0 * fq1 endif # else !! CaCO3: latitudinally-based fraction of total primary production !! absolute latitude of current grid cell flat = abs(gphit(ji,jj)) !! 0.10 at equator; 0.02 at pole fcaco3 = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0)) # endif !! AXY (09/03/09): convert CaCO3 production from function of !! primary production into a function of fast-sinking material; !! technically, this is what Dunne et al. (2007) do anyway; they !! convert total primary production estimated from surface !! chlorophyll to an export flux for which they apply conversion !! factors to estimate the various elemental fractions (Si, Ca) ftempca = ftempc * fcaco3 # if defined key_debug_medusa !! integrate total fast detritus production if (idf.eq.1) then fifd_n(ji,jj) = fifd_n(ji,jj) + (ftempn * fthk) fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi * fthk) fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe * fthk) # if defined key_roam fifd_c(ji,jj) = fifd_c(ji,jj) + (ftempc * fthk) # endif endif !! report quantities of fast-sinking detritus for each component if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'fdpd(',jk,') = ', fdpd IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme IF (lwp) write (numout,*) 'ftempn(',jk,') = ', ftempn IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe IF (lwp) write (numout,*) 'ftempc(',jk,') = ', ftempc IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca IF (lwp) write (numout,*) 'flat(',jk,') = ', flat IF (lwp) write (numout,*) 'fcaco3(',jk,') = ', fcaco3 endif # endif !!---------------------------------------------------------------------- !! This version of MEDUSA offers a choice of three methods for !! handling the remineralisation of fast detritus. All three !! do so in broadly the same way: !! !! 1. Fast detritus is stored as a 2D array [ ffastX ] !! 2. Fast detritus is added level-by-level [ ftempX ] !! 3. Fast detritus is not remineralised in the top box [ freminX ] !! 4. Remaining fast detritus is remineralised in the bottom [ fsedX ] !! box !! !! The three remineralisation methods are: !! !! 1. Ballast model (i.e. that published in Yool et al., 2011) !! (1b. Ballast-sans-ballast model) !! 2. Martin et al. (1987) !! 3. Henson et al. (2011) !! !! The first of these couples C, N and Fe remineralisation to !! the remineralisation of particulate Si and CaCO3, but the !! latter two treat remineralisation of C, N, Fe, Si and CaCO3 !! completely separately. At present a switch within the code !! regulates which submodel is used, but this should be moved !! to the namelist file. !! !! The ballast-sans-ballast submodel is an original development !! feature of MEDUSA in which the ballast submodel's general !! framework and parameterisation is used, but in which there !! is no protection of organic material afforded by ballasting !! minerals. While similar, it is not the same as the Martin !! et al. (1987) submodel. !! !! Since the three submodels behave the same in terms of !! accumulating sinking material and remineralising it all at !! the seafloor, these portions of the code below are common to !! all three. !!---------------------------------------------------------------------- if (jexport.eq.1) then !!====================================================================== !! BALLAST SUBMODEL !!====================================================================== !! !!---------------------------------------------------------------------- !! Fast-sinking detritus fluxes, pt. 1: REMINERALISATION !! aside from explicitly modelled, slow-sinking detritus, the !! model includes an implicit representation of detrital !! particles that sink too quickly to be modelled with !! explicit state variables; this sinking flux is instead !! instantaneously remineralised down the water column using !! the version of Armstrong et al. (2002)'s ballast model !! used by Dunne et al. (2007); the version of this model !! here considers silicon and calcium carbonate ballast !! minerals; this section of the code redistributes the fast !! sinking material generated locally down the water column; !! this differs from Dunne et al. (2007) in that fast sinking !! material is distributed at *every* level below that it is !! generated, rather than at every level below some fixed !! depth; this scheme is also different in that sinking material !! generated in one level is aggregated with that generated by !! shallower levels; this should make the ballast model more !! self-consistent (famous last words) !!---------------------------------------------------------------------- !! if (jk.eq.1) then !! this is the SURFACE OCEAN BOX (no remineralisation) !! freminc = 0.0 freminn = 0.0 freminfe = 0.0 freminsi = 0.0 freminca = 0.0 elseif (jk.le.jmbathy) then !! this is an OCEAN BOX (remineralise some material) !! !! set up CCD depth to be used depending on user choice if (jocalccd.eq.0) then !! use default CCD field fccd_dep = ocal_ccd(ji,jj) elseif (jocalccd.eq.1) then !! use calculated CCD field fccd_dep = f2_ccd_cal(ji,jj) endif !! !! === organic carbon === fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol) if (iball.eq.1) then fq1 = (fq0 * xmassc) !! how much it weighs (mass) fq2 = (ffastca(ji,jj) * xmassca) !! how much CaCO3 enters this box (mass) fq3 = (ffastsi(ji,jj) * xmasssi) !! how much opal enters this box (mass) fq4 = (fq2 * xprotca) + (fq3 * xprotsi) !! total protected organic C (mass) !! this next term is calculated for C but used for N and Fe as well !! it needs to be protected in case ALL C is protected if (fq4.lt.fq1) then fprotf = (fq4 / (fq1 + tiny(fq1))) !! protected fraction of total organic C (non-dim) else fprotf = 1.0 !! all organic C is protected (non-dim) endif fq5 = (1.0 - fprotf) !! unprotected fraction of total organic C (non-dim) fq6 = (fq0 * fq5) !! how much organic C is unprotected (mol) fq7 = (fq6 * exp(-(fthk / xfastc))) !! how much unprotected C leaves this box (mol) fq8 = (fq7 + (fq0 * fprotf)) !! how much total C leaves this box (mol) freminc = (fq0 - fq8) / fthk !! C remineralisation in this box (mol) ffastc(ji,jj) = fq8 # if defined key_debug_medusa !! report in/out/remin fluxes of carbon for this level if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'totalC(',jk,') = ', fq1 IF (lwp) write (numout,*) 'prtctC(',jk,') = ', fq4 IF (lwp) write (numout,*) 'fprotf(',jk,') = ', fprotf IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'IN C(',jk,') = ', fq0 IF (lwp) write (numout,*) 'LOST C(',jk,') = ', freminc * fthk IF (lwp) write (numout,*) 'OUT C(',jk,') = ', fq8 IF (lwp) write (numout,*) 'NEW C(',jk,') = ', ftempc * fthk endif # endif else fq1 = fq0 * exp(-(fthk / xfastc)) !! how much organic C leaves this box (mol) freminc = (fq0 - fq1) / fthk !! C remineralisation in this box (mol) ffastc(ji,jj) = fq1 endif !! !! === organic nitrogen === fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol) if (iball.eq.1) then fq5 = (1.0 - fprotf) !! unprotected fraction of total organic N (non-dim) fq6 = (fq0 * fq5) !! how much organic N is unprotected (mol) fq7 = (fq6 * exp(-(fthk / xfastc))) !! how much unprotected N leaves this box (mol) fq8 = (fq7 + (fq0 * fprotf)) !! how much total N leaves this box (mol) freminn = (fq0 - fq8) / fthk !! N remineralisation in this box (mol) ffastn(ji,jj) = fq8 # if defined key_debug_medusa !! report in/out/remin fluxes of carbon for this level if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'totalN(',jk,') = ', fq1 IF (lwp) write (numout,*) 'prtctN(',jk,') = ', fq4 IF (lwp) write (numout,*) 'fprotf(',jk,') = ', fprotf IF (lwp) write (numout,*) '------------------------------' if (freminn < 0.0) then IF (lwp) write (numout,*) '** FREMIN ERROR **' endif IF (lwp) write (numout,*) 'IN N(',jk,') = ', fq0 IF (lwp) write (numout,*) 'LOST N(',jk,') = ', freminn * fthk IF (lwp) write (numout,*) 'OUT N(',jk,') = ', fq8 IF (lwp) write (numout,*) 'NEW N(',jk,') = ', ftempn * fthk endif # endif else fq1 = fq0 * exp(-(fthk / xfastc)) !! how much organic N leaves this box (mol) freminn = (fq0 - fq1) / fthk !! N remineralisation in this box (mol) ffastn(ji,jj) = fq1 endif !! !! === organic iron === fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol) if (iball.eq.1) then fq5 = (1.0 - fprotf) !! unprotected fraction of total organic Fe (non-dim) fq6 = (fq0 * fq5) !! how much organic Fe is unprotected (mol) fq7 = (fq6 * exp(-(fthk / xfastc))) !! how much unprotected Fe leaves this box (mol) fq8 = (fq7 + (fq0 * fprotf)) !! how much total Fe leaves this box (mol) freminfe = (fq0 - fq8) / fthk !! Fe remineralisation in this box (mol) ffastfe(ji,jj) = fq8 else fq1 = fq0 * exp(-(fthk / xfastc)) !! how much total Fe leaves this box (mol) freminfe = (fq0 - fq1) / fthk !! Fe remineralisation in this box (mol) ffastfe(ji,jj) = fq1 endif !! !! === biogenic silicon === fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol) fq1 = fq0 * exp(-(fthk / xfastsi)) !! how much opal leaves this box (mol) freminsi = (fq0 - fq1) / fthk !! Si remineralisation in this box (mol) ffastsi(ji,jj) = fq1 !! !! === biogenic calcium carbonate === fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol) if (fdep.le.fccd_dep) then !! whole grid cell above CCD fq1 = fq0 !! above lysocline, no Ca dissolves (mol) freminca = 0.0 !! above lysocline, no Ca dissolves (mol) fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#) elseif (fdep.ge.fccd_dep) then !! whole grid cell below CCD fq1 = fq0 * exp(-(fthk / xfastca)) !! how much CaCO3 leaves this box (mol) freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol) else !! partial grid cell below CCD fq2 = fdep1 - fccd_dep !! amount of grid cell below CCD (m) fq1 = fq0 * exp(-(fq2 / xfastca)) !! how much CaCO3 leaves this box (mol) freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol) endif ffastca(ji,jj) = fq1 else !! this is BELOW THE LAST OCEAN BOX (do nothing) freminc = 0.0 freminn = 0.0 freminfe = 0.0 freminsi = 0.0 freminca = 0.0 endif elseif (jexport.eq.2.or.jexport.eq.3) then if (jexport.eq.2) then !!====================================================================== !! MARTIN ET AL. (1987) SUBMODEL !!====================================================================== !! !!---------------------------------------------------------------------- !! This submodel uses the classic Martin et al. (1987) curve !! to determine the attenuation of fast-sinking detritus down !! the water column. All three organic elements, C, N and Fe, !! are handled identically, and their quantities in sinking !! particles attenuate according to a power relationship !! governed by parameter "b". This is assigned a canonical !! value of -0.858. Biogenic opal and calcium carbonate are !! attentuated using the same function as in the ballast !! submodel !!---------------------------------------------------------------------- !! fb_val = -0.858 elseif (jexport.eq.3) then !!====================================================================== !! HENSON ET AL. (2011) SUBMODEL !!====================================================================== !! !!---------------------------------------------------------------------- !! This submodel reconfigures the Martin et al. (1987) curve by !! allowing the "b" value to vary geographically. Its value is !! set, following Henson et al. (2011), as a function of local !! sea surface temperature: !! b = -1.06 + (0.024 * SST) !! This means that remineralisation length scales are longer in !! warm, tropical areas and shorter in cold, polar areas. This !! does seem back-to-front (i.e. one would expect GREATER !! remineralisation in warmer waters), but is an outcome of !! analysis of sediment trap data, and it may reflect details !! of ecosystem structure that pertain to particle production !! rather than simply Q10. !!---------------------------------------------------------------------- !! fl_sst = tsn(ji,jj,1,jp_tem) fb_val = -1.06 + (0.024 * fl_sst) endif !! if (jk.eq.1) then !! this is the SURFACE OCEAN BOX (no remineralisation) !! freminc = 0.0 freminn = 0.0 freminfe = 0.0 freminsi = 0.0 freminca = 0.0 elseif (jk.le.jmbathy) then !! this is an OCEAN BOX (remineralise some material) !! !! === organic carbon === fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol) fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic C leaves this box (mol) freminc = (fq0 - fq1) / fthk !! C remineralisation in this box (mol) ffastc(ji,jj) = fq1 !! !! === organic nitrogen === fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol) fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic N leaves this box (mol) freminn = (fq0 - fq1) / fthk !! N remineralisation in this box (mol) ffastn(ji,jj) = fq1 !! !! === organic iron === fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol) fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic Fe leaves this box (mol) freminfe = (fq0 - fq1) / fthk !! Fe remineralisation in this box (mol) ffastfe(ji,jj) = fq1 !! !! === biogenic silicon === fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol) fq1 = fq0 * exp(-(fthk / xfastsi)) !! how much opal leaves this box (mol) freminsi = (fq0 - fq1) / fthk !! Si remineralisation in this box (mol) ffastsi(ji,jj) = fq1 !! !! === biogenic calcium carbonate === fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol) if (fdep.le.ocal_ccd(ji,jj)) then !! whole grid cell above CCD fq1 = fq0 !! above lysocline, no Ca dissolves (mol) freminca = 0.0 !! above lysocline, no Ca dissolves (mol) fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#) elseif (fdep.ge.ocal_ccd(ji,jj)) then !! whole grid cell below CCD fq1 = fq0 * exp(-(fthk / xfastca)) !! how much CaCO3 leaves this box (mol) freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol) else !! partial grid cell below CCD fq2 = fdep1 - ocal_ccd(ji,jj) !! amount of grid cell below CCD (m) fq1 = fq0 * exp(-(fq2 / xfastca)) !! how much CaCO3 leaves this box (mol) freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol) endif ffastca(ji,jj) = fq1 else !! this is BELOW THE LAST OCEAN BOX (do nothing) freminc = 0.0 freminn = 0.0 freminfe = 0.0 freminsi = 0.0 freminca = 0.0 endif endif !!---------------------------------------------------------------------- !! Fast-sinking detritus fluxes, pt. 2: UPDATE FAST FLUXES !! here locally calculated additions to the fast-sinking flux are added !! to the total fast-sinking flux; this is done here such that material !! produced in a particular layer is only remineralised below this !! layer !!---------------------------------------------------------------------- !! !! add sinking material generated in this layer to running totals !! !! === organic carbon === (diatom and mesozooplankton mortality) ffastc(ji,jj) = ffastc(ji,jj) + (ftempc * fthk) !! !! === organic nitrogen === (diatom and mesozooplankton mortality) ffastn(ji,jj) = ffastn(ji,jj) + (ftempn * fthk) !! !! === organic iron === (diatom and mesozooplankton mortality) ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe * fthk) !! !! === biogenic silicon === (diatom mortality and grazed diatoms) ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi * fthk) !! !! === biogenic calcium carbonate === (latitudinally-based fraction of total primary production) ffastca(ji,jj) = ffastca(ji,jj) + (ftempca * fthk) !!---------------------------------------------------------------------- !! Fast-sinking detritus fluxes, pt. 3: SEAFLOOR !! remineralise all remaining fast-sinking detritus to dissolved !! nutrients; the sedimentation fluxes calculated here allow the !! separation of what's remineralised sinking through the final !! ocean box from that which is added to the final box by the !! remineralisation of material that reaches the seafloor (i.e. !! the model assumes that *all* material that hits the seafloor !! is remineralised and that none is permanently buried; hey, !! this is a giant GCM model that can't be run for long enough !! to deal with burial fluxes!) !! !! in a change to this process, in part so that MEDUSA behaves !! a little more like ERSEM et al., fast-sinking detritus (N, Fe !! and C) is converted to slow sinking detritus at the seafloor !! instead of being remineralised; the rationale is that in !! shallower shelf regions (... that are not fully mixed!) this !! allows the detrital material to return slowly to dissolved !! nutrient rather than instantaneously as now; the alternative !! would be to explicitly handle seafloor organic material - a !! headache I don't wish to experience at this point; note that !! fast-sinking Si and Ca detritus is just remineralised as !! per usual !! !! AXY (13/01/12) !! in a further change to this process, again so that MEDUSA is !! a little more like ERSEM et al., material that reaches the !! seafloor can now be added to sediment pools and stored for !! slow release; there are new 2D arrays for organic nitrogen, !! iron and carbon and inorganic silicon and carbon that allow !! fast and slow detritus that reaches the seafloor to be held !! and released back to the water column more slowly; these arrays !! are transferred via the tracer restart files between repeat !! submissions of the model !!---------------------------------------------------------------------- !! ffast2slowc = 0.0 ffast2slown = 0.0 ffast2slowfe = 0.0 !! if (jk.eq.jmbathy) then !! this is the BOTTOM OCEAN BOX (remineralise everything) !! !! AXY (17/01/12): tweaked to include benthos pools !! !! === organic carbon === if (jfdfate.eq.0 .and. jorgben.eq.0) then freminc = freminc + (ffastc(ji,jj) / fthk) !! C remineralisation in this box (mol/m3) elseif (jfdfate.eq.1 .and. jorgben.eq.0) then ffast2slowc = ffastc(ji,jj) / fthk !! fast C -> slow C (mol/m3) fslowc = fslowc + ffast2slowc elseif (jfdfate.eq.0 .and. jorgben.eq.1) then f_fbenin_c(ji,jj) = ffastc(ji,jj) !! fast C -> benthic C (mol/m2) endif fsedc(ji,jj) = ffastc(ji,jj) !! record seafloor C (mol/m2) ffastc(ji,jj) = 0.0 !! !! === organic nitrogen === if (jfdfate.eq.0 .and. jorgben.eq.0) then freminn = freminn + (ffastn(ji,jj) / fthk) !! N remineralisation in this box (mol/m3) elseif (jfdfate.eq.1 .and. jorgben.eq.0) then ffast2slown = ffastn(ji,jj) / fthk !! fast N -> slow N (mol/m3) fslown = fslown + ffast2slown elseif (jfdfate.eq.0 .and. jorgben.eq.1) then f_fbenin_n(ji,jj) = ffastn(ji,jj) !! fast N -> benthic N (mol/m2) endif fsedn(ji,jj) = ffastn(ji,jj) !! record seafloor N (mol/m2) ffastn(ji,jj) = 0.0 !! !! === organic iron === if (jfdfate.eq.0 .and. jorgben.eq.0) then freminfe = freminfe + (ffastfe(ji,jj) / fthk) !! Fe remineralisation in this box (mol/m3) elseif (jfdfate.eq.1 .and. jorgben.eq.0) then ffast2slowfe = ffastn(ji,jj) / fthk !! fast Fe -> slow Fe (mol/m3) elseif (jfdfate.eq.0 .and. jorgben.eq.1) then f_fbenin_fe(ji,jj) = ffastfe(ji,jj) !! fast Fe -> benthic Fe (mol/m2) endif fsedfe(ji,jj) = ffastfe(ji,jj) !! record seafloor Fe (mol/m2) ffastfe(ji,jj) = 0.0 !! !! === biogenic silicon === if (jinorgben.eq.0) then freminsi = freminsi + (ffastsi(ji,jj) / fthk) !! Si remineralisation in this box (mol/m3) elseif (jinorgben.eq.1) then f_fbenin_si(ji,jj) = ffastsi(ji,jj) !! fast Si -> benthic Si (mol/m2) endif fsedsi(ji,jj) = ffastsi(ji,jj) !! record seafloor Si (mol/m2) ffastsi(ji,jj) = 0.0 !! !! === biogenic calcium carbonate === if (jinorgben.eq.0) then freminca = freminca + (ffastca(ji,jj) / fthk) !! Ca remineralisation in this box (mol/m3) elseif (jinorgben.eq.1) then f_fbenin_ca(ji,jj) = ffastca(ji,jj) !! fast Ca -> benthic Ca (mol/m2) endif fsedca(ji,jj) = ffastca(ji,jj) !! record seafloor Ca (mol/m2) ffastca(ji,jj) = 0.0 endif # if defined key_debug_medusa if (idf.eq.1) then !!---------------------------------------------------------------------- !! Integrate total fast detritus remineralisation !!---------------------------------------------------------------------- !! fofd_n(ji,jj) = fofd_n(ji,jj) + (freminn * fthk) fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi * fthk) fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe * fthk) # if defined key_roam fofd_c(ji,jj) = fofd_c(ji,jj) + (freminc * fthk) # endif endif # endif !!---------------------------------------------------------------------- !! Sort out remineralisation tally of fast-sinking detritus !!---------------------------------------------------------------------- !! !! update fast-sinking regeneration arrays fregenfast(ji,jj) = fregenfast(ji,jj) + (freminn * fthk) fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi * fthk) # if defined key_roam fregenfastc(ji,jj) = fregenfastc(ji,jj) + (freminc * fthk) # endif !!---------------------------------------------------------------------- !! Benthic remineralisation fluxes !!---------------------------------------------------------------------- !! if (jk.eq.jmbathy) then !! !! organic components if (jorgben.eq.1) then f_benout_n(ji,jj) = xsedn * zn_sed_n(ji,jj) f_benout_fe(ji,jj) = xsedfe * zn_sed_fe(ji,jj) f_benout_c(ji,jj) = xsedc * zn_sed_c(ji,jj) endif !! !! inorganic components if (jinorgben.eq.1) then f_benout_si(ji,jj) = xsedsi * zn_sed_si(ji,jj) f_benout_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj) !! !! account for CaCO3 that dissolves when it shouldn't if ( fdep .le. fccd_dep ) then f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj) endif endif endif CALL flush(numout) !!====================================================================== !! LOCAL GRID CELL TRENDS !!====================================================================== !! !!---------------------------------------------------------------------- !! Determination of trends !!---------------------------------------------------------------------- !! !!---------------------------------------------------------------------- !! chlorophyll btra(jpchn) = b0 * ( & + ((frn * fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2) * (fthetan / xxi) ) btra(jpchd) = b0 * ( & + ((frd * fprd * zphd) - fgmepd - fdpd - fdpd2) * (fthetad / xxi) ) !! !!---------------------------------------------------------------------- !! phytoplankton btra(jpphn) = b0 * ( & + (fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2 ) btra(jpphd) = b0 * ( & + (fprd * zphd) - fgmepd - fdpd - fdpd2 ) btra(jppds) = b0 * ( & + (fprds * zpds) - fgmepds - fdpds - fsdiss - fdpds2 ) !! !!---------------------------------------------------------------------- !! zooplankton btra(jpzmi) = b0 * ( & + fmigrow - fgmezmi - fdzmi - fdzmi2 ) btra(jpzme) = b0 * ( & + fmegrow - fdzme - fdzme2 ) !! !!---------------------------------------------------------------------- !! detritus btra(jpdet) = b0 * ( & + fdpn + ((1.0 - xfdfrac1) * fdpd) & ! mort. losses + fdzmi + ((1.0 - xfdfrac2) * fdzme) & ! mort. losses + ((1.0 - xbetan) * (finmi + finme)) & ! assim. inefficiency - fgmid - fgmed - fdd & ! grazing and remin. + ffast2slown ) ! seafloor fast->slow !! !!---------------------------------------------------------------------- !! dissolved inorganic nitrogen nutrient fn_cons = 0.0 & - (fprn * zphn) - (fprd * zphd) ! primary production fn_prod = 0.0 & + (xphi * (fgmipn + fgmid)) & ! messy feeding remin. + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) & ! messy feeding remin. + fmiexcr + fmeexcr + fdd + freminn & ! excretion and remin. + fdpn2 + fdpd2 + fdzmi2 + fdzme2 ! metab. losses !! !! riverine flux if ( jriver_n .gt. 0 ) then f_riv_loc_n = f_riv_n(ji,jj) * friver_dep(jk,jmbathy) / fthk fn_prod = fn_prod + f_riv_loc_n endif !! !! benthic remineralisation if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then fn_prod = fn_prod + (f_benout_n(ji,jj) / fthk) endif !! btra(jpdin) = b0 * ( & fn_prod + fn_cons ) !! fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved nitrogen fn_cons ) ) fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fthk * ( & ! production of dissolved nitrogen fn_prod ) ) !! !!---------------------------------------------------------------------- !! dissolved silicic acid nutrient fs_cons = 0.0 & - (fprds * zpds) ! opal production fs_prod = 0.0 & + fsdiss & ! opal dissolution + ((1.0 - xfdfrac1) * fdpds) & ! mort. loss + ((1.0 - xfdfrac3) * fgmepds) & ! egestion of grazed Si + freminsi + fdpds2 ! fast diss. and metab. losses !! !! riverine flux if ( jriver_si .gt. 0 ) then f_riv_loc_si = f_riv_si(ji,jj) * friver_dep(jk,jmbathy) / fthk fs_prod = fs_prod + f_riv_loc_si endif !! !! benthic remineralisation if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then fs_prod = fs_prod + (f_benout_si(ji,jj) / fthk) endif !! btra(jpsil) = b0 * ( & fs_prod + fs_cons ) !! fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved silicon fs_cons ) ) fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fthk * ( & ! production of dissolved silicon fs_prod ) ) !! !!---------------------------------------------------------------------- !! dissolved "iron" nutrient btra(jpfer) = b0 * ( & + (xrfn * btra(jpdin)) + ffetop + ffebot - ffescav ) # if defined key_roam !! !!---------------------------------------------------------------------- !! AXY (26/11/08): implicit detrital carbon change btra(jpdtc) = b0 * ( & + (xthetapn * fdpn) + ((1.0 - xfdfrac1) * (xthetapd * fdpd)) & ! mort. losses + (xthetazmi * fdzmi) + ((1.0 - xfdfrac2) * (xthetazme * fdzme)) & ! mort. losses + ((1.0 - xbetac) * (ficmi + ficme)) & ! assim. inefficiency - fgmidc - fgmedc - fddc & ! grazing and remin. + ffast2slowc ) ! seafloor fast->slow !! !!---------------------------------------------------------------------- !! dissolved inorganic carbon fc_cons = 0.0 & - (xthetapn * fprn * zphn) - (xthetapd * fprd * zphd) ! primary production fc_prod = 0.0 & + (xthetapn * xphi * fgmipn) + (xphi * fgmidc) & ! messy feeding remin + (xthetapn * xphi * fgmepn) + (xthetapd * xphi * fgmepd) & ! messy feeding remin + (xthetazmi * xphi * fgmezmi) + (xphi * fgmedc) & ! messy feeding remin + fmiresp + fmeresp + fddc + freminc + (xthetapn * fdpn2) & ! resp., remin., losses + (xthetapd * fdpd2) + (xthetazmi * fdzmi2) & ! losses + (xthetazme * fdzme2) ! losses !! !! riverine flux if ( jriver_c .gt. 0 ) then f_riv_loc_c = f_riv_c(ji,jj) * friver_dep(jk,jmbathy) / fthk fc_prod = fc_prod + f_riv_loc_c endif !! !! benthic remineralisation if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then fc_prod = fc_prod + (f_benout_c(ji,jj) / fthk) endif if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then fc_prod = fc_prod + (f_benout_ca(ji,jj) / fthk) endif !! !! community respiration (does not include CaCO3 terms - obviously!) fcomm_resp(ji,jj) = fcomm_resp(ji,jj) + fc_prod !! !! CaCO3 fc_prod = fc_prod - ftempca + freminca !! !! riverine flux if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then fc_prod = fc_prod + f_riv_c(ji,jj) endif !! btra(jpdic) = b0 * ( & fc_prod + fc_cons ) !! fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved carbon fc_cons ) ) fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fthk * ( & ! production of dissolved carbon fc_prod ) ) !! !!---------------------------------------------------------------------- !! alkalinity fa_prod = 0.0 & + (2.0 * freminca) ! CaCO3 dissolution fa_cons = 0.0 & - (2.0 * ftempca) ! CaCO3 production !! !! riverine flux if ( jriver_alk .gt. 0 ) then f_riv_loc_alk = f_riv_alk(ji,jj) * friver_dep(jk,jmbathy) / fthk fa_prod = fa_prod + f_riv_loc_alk endif !! !! benthic remineralisation if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fthk) endif !! btra(jpalk) = b0 * ( & fa_prod + fa_cons ) !! !!---------------------------------------------------------------------- !! oxygen (has protection at low O2 concentrations; OCMIP-2 style) fo2_prod = 0.0 & + (xthetanit * fprn * zphn) & ! Pn primary production, N + (xthetanit * fprd * zphd) & ! Pd primary production, N + (xthetarem * xthetapn * fprn * zphn) & ! Pn primary production, C + (xthetarem * xthetapd * fprd * zphd) ! Pd primary production, C fo2_ncons = 0.0 & - (xthetanit * xphi * fgmipn) & ! Pn messy feeding remin., N - (xthetanit * xphi * fgmid) & ! D messy feeding remin., N - (xthetanit * xphi * fgmepn) & ! Pn messy feeding remin., N - (xthetanit * xphi * fgmepd) & ! Pd messy feeding remin., N - (xthetanit * xphi * fgmezmi) & ! Zi messy feeding remin., N - (xthetanit * xphi * fgmed) & ! D messy feeding remin., N - (xthetanit * fmiexcr) & ! microzoo excretion, N - (xthetanit * fmeexcr) & ! mesozoo excretion, N - (xthetanit * fdd) & ! slow detritus remin., N - (xthetanit * freminn) & ! fast detritus remin., N - (xthetanit * fdpn2) & ! Pn losses, N - (xthetanit * fdpd2) & ! Pd losses, N - (xthetanit * fdzmi2) & ! Zmi losses, N - (xthetanit * fdzme2) ! Zme losses, N !! !! benthic remineralisation if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then fo2_ncons = fo2_ncons - (xthetanit * f_benout_n(ji,jj) / fthk) endif fo2_ccons = 0.0 & - (xthetarem * xthetapn * xphi * fgmipn) & ! Pn messy feeding remin., C - (xthetarem * xphi * fgmidc) & ! D messy feeding remin., C - (xthetarem * xthetapn * xphi * fgmepn) & ! Pn messy feeding remin., C - (xthetarem * xthetapd * xphi * fgmepd) & ! Pd messy feeding remin., C - (xthetarem * xthetazmi * xphi * fgmezmi) & ! Zi messy feeding remin., C - (xthetarem * xphi * fgmedc) & ! D messy feeding remin., C - (xthetarem * fmiresp) & ! microzoo respiration, C - (xthetarem * fmeresp) & ! mesozoo respiration, C - (xthetarem * fddc) & ! slow detritus remin., C - (xthetarem * freminc) & ! fast detritus remin., C - (xthetarem * xthetapn * fdpn2) & ! Pn losses, C - (xthetarem * xthetapd * fdpd2) & ! Pd losses, C - (xthetarem * xthetazmi * fdzmi2) & ! Zmi losses, C - (xthetarem * xthetazme * fdzme2) ! Zme losses, C !! !! benthic remineralisation if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then fo2_ccons = fo2_ccons - (xthetarem * f_benout_c(ji,jj) / fthk) endif fo2_cons = fo2_ncons + fo2_ccons !! !! is this a suboxic zone? if (zoxy.lt.xo2min) then ! deficient O2; production fluxes only btra(jpoxy) = b0 * ( & fo2_prod ) foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod ) foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fthk * fo2_cons ) else ! sufficient O2; production + consumption fluxes btra(jpoxy) = b0 * ( & fo2_prod + fo2_cons ) foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod ) foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fthk * fo2_cons ) endif !! !! air-sea fluxes (if this is the surface box) if (jk.eq.1) then !! !! CO2 flux btra(jpdic) = btra(jpdic) + (b0 * f_co2flux) !! !! O2 flux (mol/m3/s -> mmol/m3/d) btra(jpoxy) = btra(jpoxy) + (b0 * f_o2flux) endif # endif # if defined key_debug_medusa !! report state variable fluxes (not including fast-sinking detritus) if (idf.eq.1.AND.idfval.eq.1) then IF (lwp) write (numout,*) '------------------------------' IF (lwp) write (numout,*) 'btra(jpchn)(',jk,') = ', btra(jpchn) IF (lwp) write (numout,*) 'btra(jpchd)(',jk,') = ', btra(jpchd) IF (lwp) write (numout,*) 'btra(jpphn)(',jk,') = ', btra(jpphn) IF (lwp) write (numout,*) 'btra(jpphd)(',jk,') = ', btra(jpphd) IF (lwp) write (numout,*) 'btra(jppds)(',jk,') = ', btra(jppds) IF (lwp) write (numout,*) 'btra(jpzmi)(',jk,') = ', btra(jpzmi) IF (lwp) write (numout,*) 'btra(jpzme)(',jk,') = ', btra(jpzme) IF (lwp) write (numout,*) 'btra(jpdet)(',jk,') = ', btra(jpdet) IF (lwp) write (numout,*) 'btra(jpdin)(',jk,') = ', btra(jpdin) IF (lwp) write (numout,*) 'btra(jpsil)(',jk,') = ', btra(jpsil) IF (lwp) write (numout,*) 'btra(jpfer)(',jk,') = ', btra(jpfer) # if defined key_roam IF (lwp) write (numout,*) 'btra(jpdtc)(',jk,') = ', btra(jpdtc) IF (lwp) write (numout,*) 'btra(jpdic)(',jk,') = ', btra(jpdic) IF (lwp) write (numout,*) 'btra(jpalk)(',jk,') = ', btra(jpalk) IF (lwp) write (numout,*) 'btra(jpoxy)(',jk,') = ', btra(jpoxy) # endif endif # endif !!---------------------------------------------------------------------- !! Integrate calculated fluxes for mass balance !!---------------------------------------------------------------------- !! !! === nitrogen === fflx_n(ji,jj) = fflx_n(ji,jj) + & fthk * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) ) !! === silicon === fflx_si(ji,jj) = fflx_si(ji,jj) + & fthk * ( btra(jppds) + btra(jpsil) ) !! === iron === fflx_fe(ji,jj) = fflx_fe(ji,jj) + & fthk * ( ( xrfn * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet)) ) + btra(jpfer) ) # if defined key_roam !! === carbon === fflx_c(ji,jj) = fflx_c(ji,jj) + & fthk * ( (xthetapn * btra(jpphn)) + (xthetapd * btra(jpphd)) + & (xthetazmi * btra(jpzmi)) + (xthetazme * btra(jpzme)) + btra(jpdtc) + btra(jpdic) ) !! === alkalinity === fflx_a(ji,jj) = fflx_a(ji,jj) + & fthk * ( btra(jpalk) ) !! === oxygen === fflx_o2(ji,jj) = fflx_o2(ji,jj) + & fthk * ( btra(jpoxy) ) # endif !!---------------------------------------------------------------------- !! Apply calculated tracer fluxes !!---------------------------------------------------------------------- !! !! units: [unit of tracer] per second (fluxes are calculated above per day) !! ibio_switch = 1 # if defined key_gulf_finland !! AXY (17/05/13): fudge in a Gulf of Finland correction; uses longitude- !! latitude range to establish if this is a Gulf of Finland !! grid cell; if so, then BGC fluxes are ignored (though !! still calculated); for reference, this is meant to be a !! temporary fix to see if all of my problems can be done !! away with if I switch off BGC fluxes in the Gulf of !! Finland, which currently appears the source of trouble if ( glamt(ji,jj).gt.24.7 .and. glamt(ji,jj).lt.27.8 .and. & & gphit(ji,jj).gt.59.2 .and. gphit(ji,jj).lt.60.2 ) then ibio_switch = 0 endif # endif if (ibio_switch.eq.1) then tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(jpchn) / 86400.) tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(jpchd) / 86400.) tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(jpphn) / 86400.) tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(jpphd) / 86400.) tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(jppds) / 86400.) tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(jpzmi) / 86400.) tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(jpzme) / 86400.) tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(jpdet) / 86400.) tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(jpdin) / 86400.) tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(jpsil) / 86400.) tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(jpfer) / 86400.) # if defined key_roam tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(jpdtc) / 86400.) tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(jpdic) / 86400.) tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(jpalk) / 86400.) tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(jpoxy) / 86400.) # endif endif !! AXY (18/11/16): CMIP6 diagnostics IF( med_diag%FBDDTALK%dgsave ) THEN fbddtalk(ji,jj) = fbddtalk(ji,jj) + (btra(jpalk) * fthk) ENDIF IF( med_diag%FBDDTDIC%dgsave ) THEN fbddtdic(ji,jj) = fbddtdic(ji,jj) + (btra(jpdic) * fthk) ENDIF IF( med_diag%FBDDTDIFE%dgsave ) THEN fbddtdife(ji,jj) = fbddtdife(ji,jj) + (btra(jpfer) * fthk) ENDIF IF( med_diag%FBDDTDIN%dgsave ) THEN fbddtdin(ji,jj) = fbddtdin(ji,jj) + (btra(jpdin) * fthk) ENDIF IF( med_diag%FBDDTDISI%dgsave ) THEN fbddtdisi(ji,jj) = fbddtdisi(ji,jj) + (btra(jpsil) * fthk) ENDIF !! IF( med_diag%BDDTALK3%dgsave ) THEN bddtalk3(ji,jj,jk) = btra(jpalk) ENDIF IF( med_diag%BDDTDIC3%dgsave ) THEN bddtdic3(ji,jj,jk) = btra(jpdic) ENDIF IF( med_diag%BDDTDIFE3%dgsave ) THEN bddtdife3(ji,jj,jk) = btra(jpfer) ENDIF IF( med_diag%BDDTDIN3%dgsave ) THEN bddtdin3(ji,jj,jk) = btra(jpdin) ENDIF IF( med_diag%BDDTDISI3%dgsave ) THEN bddtdisi3(ji,jj,jk) = btra(jpsil) ENDIF # if defined key_debug_medusa IF (lwp) write (numout,*) '------' IF (lwp) write (numout,*) 'trc_bio_medusa: end all calculations' IF (lwp) write (numout,*) 'trc_bio_medusa: now outputs' CALL flush(numout) # endif # if defined key_axy_nancheck !!---------------------------------------------------------------------- !! Check calculated tracer fluxes !!---------------------------------------------------------------------- !! DO jn = 1,jptra fq0 = btra(jn) !! AXY (30/01/14): "isnan" problem on HECTOR !! if (fq0 /= fq0 ) then if ( ieee_is_nan( fq0 ) ) then !! there's a NaN here if (lwp) write(numout,*) 'NAN detected in btra(', ji, ',', & & jj, ',', jk, ',', jn, ') at time', kt CALL ctl_stop( 'trcbio_medusa, NAN in btra field' ) endif ENDDO DO jn = 1,jptra fq0 = tra(ji,jj,jk,jn) !! AXY (30/01/14): "isnan" problem on HECTOR !! if (fq0 /= fq0 ) then if ( ieee_is_nan( fq0 ) ) then !! there's a NaN here if (lwp) write(numout,*) 'NAN detected in tra(', ji, ',', & & jj, ',', jk, ',', jn, ') at time', kt CALL ctl_stop( 'trcbio_medusa, NAN in tra field' ) endif ENDDO CALL flush(numout) # endif !!---------------------------------------------------------------------- !! Check model conservation !! these terms merely sum up the tendency terms of the relevant !! state variables, which should sum to zero; the iron cycle is !! complicated by fluxes that add (aeolian deposition and seafloor !! remineralisation) and remove (scavenging) dissolved iron from !! the model (i.e. the sum of iron fluxes is unlikely to be zero) !!---------------------------------------------------------------------- !! !! fnit0 = btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) ! + ftempn !! fsil0 = btra(jppds) + btra(jpsil) ! + ftempsi !! ffer0 = (xrfn * fnit0) + btra(jpfer) # if defined key_roam !! fcar0 = 0. !! falk0 = 0. !! foxy0 = 0. # endif !! !! if (kt/240*240.eq.kt) then !! if (ji.eq.2.and.jj.eq.2.and.jk.eq.1) then !! IF (lwp) write (*,*) '*******!MEDUSA Conservation!*******',kt # if defined key_roam !! IF (lwp) write (*,*) fnit0,fsil0,ffer0,fcar0,falk0,foxy0 # else !! IF (lwp) write (*,*) fnit0,fsil0,ffer0 # endif !! endif !! endif # if defined key_trc_diabio !!====================================================================== !! LOCAL GRID CELL DIAGNOSTICS !!====================================================================== !! !!---------------------------------------------------------------------- !! Full diagnostics key_trc_diabio: !! LOBSTER and PISCES support full diagnistics option key_trc_diabio !! which gives an option of FULL output of biological sourses and sinks. !! I cannot see any reason for doing this. If needed, it can be done !! as shown below. !!---------------------------------------------------------------------- !! IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio' !! trbio(ji,jj,jk, 1) = fprn # endif IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN !!---------------------------------------------------------------------- !! Add in XML diagnostics stuff !!---------------------------------------------------------------------- !! !! ** 2D diagnostics # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk loop' CALL flush(numout) # endif IF ( med_diag%PRN%dgsave ) THEN fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn * zphn * fthk) ENDIF IF ( med_diag%MPN%dgsave ) THEN fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn * fthk) ENDIF IF ( med_diag%PRD%dgsave ) THEN fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd * zphd * fthk) ENDIF IF( med_diag%MPD%dgsave ) THEN fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd * fthk) ENDIF ! IF( med_diag%DSED%dgsave ) THEN ! CALL iom_put( "DSED" , ftot_n ) ! ENDIF IF( med_diag%OPAL%dgsave ) THEN fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds * zpds * fthk) ENDIF IF( med_diag%OPALDISS%dgsave ) THEN fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss * fthk) ENDIF IF( med_diag%GMIPn%dgsave ) THEN fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn * fthk) ENDIF IF( med_diag%GMID%dgsave ) THEN fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid * fthk) ENDIF IF( med_diag%MZMI%dgsave ) THEN fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi * fthk) ENDIF IF( med_diag%GMEPN%dgsave ) THEN fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn * fthk) ENDIF IF( med_diag%GMEPD%dgsave ) THEN fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd * fthk) ENDIF IF( med_diag%GMEZMI%dgsave ) THEN fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi * fthk) ENDIF IF( med_diag%GMED%dgsave ) THEN fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed * fthk) ENDIF IF( med_diag%MZME%dgsave ) THEN fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme * fthk) ENDIF ! IF( med_diag%DEXP%dgsave ) THEN ! CALL iom_put( "DEXP" , ftot_n ) ! ENDIF IF( med_diag%DETN%dgsave ) THEN fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown * fthk) ENDIF IF( med_diag%MDET%dgsave ) THEN fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd * fthk) ENDIF IF( med_diag%AEOLIAN%dgsave ) THEN ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop * fthk) ENDIF IF( med_diag%BENTHIC%dgsave ) THEN ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot * fthk) ENDIF IF( med_diag%SCAVENGE%dgsave ) THEN ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav * fthk) ENDIF IF( med_diag%PN_JLIM%dgsave ) THEN ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln * zphn * fthk) fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn * zphn * fthk) ENDIF IF( med_diag%PN_NLIM%dgsave ) THEN fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln * zphn * fthk) ENDIF IF( med_diag%PN_FELIM%dgsave ) THEN ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln * zphn * fthk) ENDIF IF( med_diag%PD_JLIM%dgsave ) THEN ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld * zphd * fthk) fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd * zphd * fthk) ENDIF IF( med_diag%PD_NLIM%dgsave ) THEN fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld * zphd * fthk) ENDIF IF( med_diag%PD_FELIM%dgsave ) THEN ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld * zphd * fthk) ENDIF IF( med_diag%PD_SILIM%dgsave ) THEN fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2 * zphd * fthk) ENDIF IF( med_diag%PDSILIM2%dgsave ) THEN fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld * zphd * fthk) ENDIF !! IF( med_diag%TOTREG_N%dgsave ) THEN fregen2d(ji,jj) = fregen2d(ji,jj) + fregen ENDIF IF( med_diag%TOTRG_SI%dgsave ) THEN fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi ENDIF !! IF( med_diag%FASTN%dgsave ) THEN ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn * fthk) ENDIF IF( med_diag%FASTSI%dgsave ) THEN ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi * fthk) ENDIF IF( med_diag%FASTFE%dgsave ) THEN ftempfe2d(ji,jj) =ftempfe2d(ji,jj) + (ftempfe * fthk) ENDIF IF( med_diag%FASTC%dgsave ) THEN ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc * fthk) ENDIF IF( med_diag%FASTCA%dgsave ) THEN ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca * fthk) ENDIF !! IF( med_diag%REMINN%dgsave ) THEN freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn * fthk) ENDIF IF( med_diag%REMINSI%dgsave ) THEN freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi * fthk) ENDIF IF( med_diag%REMINFE%dgsave ) THEN freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe * fthk) ENDIF IF( med_diag%REMINC%dgsave ) THEN freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc * fthk) ENDIF IF( med_diag%REMINCA%dgsave ) THEN freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca * fthk) ENDIF !! # if defined key_roam !! !! AXY (09/11/16): CMIP6 diagnostics IF( med_diag%FD_NIT3%dgsave ) THEN fd_nit3(ji,jj,jk) = ffastn(ji,jj) ENDIF IF( med_diag%FD_SIL3%dgsave ) THEN fd_sil3(ji,jj,jk) = ffastsi(ji,jj) ENDIF IF( med_diag%FD_CAR3%dgsave ) THEN fd_car3(ji,jj,jk) = ffastc(ji,jj) ENDIF IF( med_diag%FD_CAL3%dgsave ) THEN fd_cal3(ji,jj,jk) = ffastca(ji,jj) ENDIF !! IF (jk.eq.i0100) THEN IF( med_diag%RR_0100%dgsave ) THEN ffastca2d(ji,jj) = & ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) ENDIF ELSE IF (jk.eq.i0500) THEN IF( med_diag%RR_0500%dgsave ) THEN ffastca2d(ji,jj) = & ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) ENDIF ELSE IF (jk.eq.i1000) THEN IF( med_diag%RR_1000%dgsave ) THEN ffastca2d(ji,jj) = & ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) ENDIF ELSE IF (jk.eq.jmbathy) THEN IF( med_diag%IBEN_N%dgsave ) THEN iben_n2d(ji,jj) = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) ENDIF IF( med_diag%IBEN_FE%dgsave ) THEN iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj) ENDIF IF( med_diag%IBEN_C%dgsave ) THEN iben_c2d(ji,jj) = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) ENDIF IF( med_diag%IBEN_SI%dgsave ) THEN iben_si2d(ji,jj) = f_fbenin_si(ji,jj) ENDIF IF( med_diag%IBEN_CA%dgsave ) THEN iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj) ENDIF IF( med_diag%OBEN_N%dgsave ) THEN oben_n2d(ji,jj) = f_benout_n(ji,jj) ENDIF IF( med_diag%OBEN_FE%dgsave ) THEN oben_fe2d(ji,jj) = f_benout_fe(ji,jj) ENDIF IF( med_diag%OBEN_C%dgsave ) THEN oben_c2d(ji,jj) = f_benout_c(ji,jj) ENDIF IF( med_diag%OBEN_SI%dgsave ) THEN oben_si2d(ji,jj) = f_benout_si(ji,jj) ENDIF IF( med_diag%OBEN_CA%dgsave ) THEN oben_ca2d(ji,jj) = f_benout_ca(ji,jj) ENDIF IF( med_diag%SFR_OCAL%dgsave ) THEN sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk) ENDIF IF( med_diag%SFR_OARG%dgsave ) THEN sfr_oarg2d(ji,jj) = f3_omarg(ji,jj,jk) ENDIF IF( med_diag%LYSO_CA%dgsave ) THEN lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj) ENDIF ENDIF !! end bathy-1 diags !! IF( med_diag%RIV_N%dgsave ) THEN rivn2d(ji,jj) = rivn2d(ji,jj) + (f_riv_loc_n * fthk) ENDIF IF( med_diag%RIV_SI%dgsave ) THEN rivsi2d(ji,jj) = rivsi2d(ji,jj) + (f_riv_loc_si * fthk) ENDIF IF( med_diag%RIV_C%dgsave ) THEN rivc2d(ji,jj) = rivc2d(ji,jj) + (f_riv_loc_c * fthk) ENDIF IF( med_diag%RIV_ALK%dgsave ) THEN rivalk2d(ji,jj) = rivalk2d(ji,jj) + (f_riv_loc_alk * fthk) ENDIF IF( med_diag%DETC%dgsave ) THEN fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc * fthk) ENDIF !! !! !! IF( med_diag%PN_LLOSS%dgsave ) THEN fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2 * fthk) ENDIF IF( med_diag%PD_LLOSS%dgsave ) THEN fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2 * fthk) ENDIF IF( med_diag%ZI_LLOSS%dgsave ) THEN fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2 * fthk) ENDIF IF( med_diag%ZE_LLOSS%dgsave ) THEN fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2 * fthk) ENDIF IF( med_diag%ZI_MES_N%dgsave ) THEN zimesn2d(ji,jj) = zimesn2d(ji,jj) + & (xphi * (fgmipn + fgmid) * fthk) ENDIF IF( med_diag%ZI_MES_D%dgsave ) THEN zimesd2d(ji,jj) = zimesd2d(ji,jj) + & ((1. - xbetan) * finmi * fthk) ENDIF IF( med_diag%ZI_MES_C%dgsave ) THEN zimesc2d(ji,jj) = zimesc2d(ji,jj) + & (xphi * ((xthetapn * fgmipn) + fgmidc) * fthk) ENDIF IF( med_diag%ZI_MESDC%dgsave ) THEN zimesdc2d(ji,jj) = zimesdc2d(ji,jj) + & ((1. - xbetac) * ficmi * fthk) ENDIF IF( med_diag%ZI_EXCR%dgsave ) THEN ziexcr2d(ji,jj) = ziexcr2d(ji,jj) + (fmiexcr * fthk) ENDIF IF( med_diag%ZI_RESP%dgsave ) THEN ziresp2d(ji,jj) = ziresp2d(ji,jj) + (fmiresp * fthk) ENDIF IF( med_diag%ZI_GROW%dgsave ) THEN zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow * fthk) ENDIF IF( med_diag%ZE_MES_N%dgsave ) THEN zemesn2d(ji,jj) = zemesn2d(ji,jj) + & (xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk) ENDIF IF( med_diag%ZE_MES_D%dgsave ) THEN zemesd2d(ji,jj) = zemesd2d(ji,jj) + & ((1. - xbetan) * finme * fthk) ENDIF IF( med_diag%ZE_MES_C%dgsave ) THEN zemesc2d(ji,jj) = zemesc2d(ji,jj) + & (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & (xthetazmi * fgmezmi) + fgmedc) * fthk) ENDIF IF( med_diag%ZE_MESDC%dgsave ) THEN zemesdc2d(ji,jj) = zemesdc2d(ji,jj) + & ((1. - xbetac) * ficme * fthk) ENDIF IF( med_diag%ZE_EXCR%dgsave ) THEN zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr * fthk) ENDIF IF( med_diag%ZE_RESP%dgsave ) THEN zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp * fthk) ENDIF IF( med_diag%ZE_GROW%dgsave ) THEN zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow * fthk) ENDIF IF( med_diag%MDETC%dgsave ) THEN mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc * fthk) ENDIF IF( med_diag%GMIDC%dgsave ) THEN gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc * fthk) ENDIF IF( med_diag%GMEDC%dgsave ) THEN gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc * fthk) ENDIF !! # endif !! !! ** 3D diagnostics IF( med_diag%TPP3%dgsave ) THEN tpp3d(ji,jj,jk) = (fprn * zphn) + (fprd * zphd) !CALL iom_put( "TPP3" , tpp3d ) ENDIF IF( med_diag%TPPD3%dgsave ) THEN tppd3(ji,jj,jk) = (fprd * zphd) ENDIF IF( med_diag%REMIN3N%dgsave ) THEN remin3dn(ji,jj,jk) = fregen + (freminn * fthk) !! remineralisation !CALL iom_put( "REMIN3N" , remin3dn ) ENDIF !! IF( med_diag%PH3%dgsave ) THEN !! CALL iom_put( "PH3" , f3_pH ) !! ENDIF !! IF( med_diag%OM_CAL3%dgsave ) THEN !! CALL iom_put( "OM_CAL3" , f3_omcal ) !! ENDIF !! !! AXY (09/11/16): CMIP6 diagnostics IF ( med_diag%DCALC3%dgsave ) THEN dcalc3(ji,jj,jk) = freminca ENDIF IF ( med_diag%FEDISS3%dgsave ) THEN fediss3(ji,jj,jk) = ffetop ENDIF IF ( med_diag%FESCAV3%dgsave ) THEN fescav3(ji,jj,jk) = ffescav ENDIF IF ( med_diag%MIGRAZP3%dgsave ) THEN migrazp3(ji,jj,jk) = fgmipn * xthetapn ENDIF IF ( med_diag%MIGRAZD3%dgsave ) THEN migrazd3(ji,jj,jk) = fgmidc ENDIF IF ( med_diag%MEGRAZP3%dgsave ) THEN megrazp3(ji,jj,jk) = (fgmepn * xthetapn) + (fgmepd * xthetapd) ENDIF IF ( med_diag%MEGRAZD3%dgsave ) THEN megrazd3(ji,jj,jk) = fgmedc ENDIF IF ( med_diag%MEGRAZZ3%dgsave ) THEN megrazz3(ji,jj,jk) = (fgmezmi * xthetazmi) ENDIF IF ( med_diag%PBSI3%dgsave ) THEN pbsi3(ji,jj,jk) = (fprds * zpds) ENDIF IF ( med_diag%PCAL3%dgsave ) THEN pcal3(ji,jj,jk) = ftempca ENDIF IF ( med_diag%REMOC3%dgsave ) THEN remoc3(ji,jj,jk) = freminc ENDIF IF ( med_diag%PNLIMJ3%dgsave ) THEN ! pnlimj3(ji,jj,jk) = fjln pnlimj3(ji,jj,jk) = fjlim_pn ENDIF IF ( med_diag%PNLIMN3%dgsave ) THEN pnlimn3(ji,jj,jk) = fnln ENDIF IF ( med_diag%PNLIMFE3%dgsave ) THEN pnlimfe3(ji,jj,jk) = ffln ENDIF IF ( med_diag%PDLIMJ3%dgsave ) THEN ! pdlimj3(ji,jj,jk) = fjld pdlimj3(ji,jj,jk) = fjlim_pd ENDIF IF ( med_diag%PDLIMN3%dgsave ) THEN pdlimn3(ji,jj,jk) = fnld ENDIF IF ( med_diag%PDLIMFE3%dgsave ) THEN pdlimfe3(ji,jj,jk) = ffld ENDIF IF ( med_diag%PDLIMSI3%dgsave ) THEN pdlimsi3(ji,jj,jk) = fsld2 ENDIF !! !! ** Without using iom_use ELSE IF( ln_diatrc ) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk ln_diatrc' CALL flush(numout) # endif !!---------------------------------------------------------------------- !! Prepare 2D diagnostics !!---------------------------------------------------------------------- !! !! if ((kt / 240*240).eq.kt) then !! IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt !! endif trc2d(ji,jj,1) = ftot_n(ji,jj) !! nitrogen inventory trc2d(ji,jj,2) = ftot_si(ji,jj) !! silicon inventory trc2d(ji,jj,3) = ftot_fe(ji,jj) !! iron inventory trc2d(ji,jj,4) = trc2d(ji,jj,4) + (fprn * zphn * fthk) !! non-diatom production trc2d(ji,jj,5) = trc2d(ji,jj,5) + (fdpn * fthk) !! non-diatom non-grazing losses trc2d(ji,jj,6) = trc2d(ji,jj,6) + (fprd * zphd * fthk) !! diatom production trc2d(ji,jj,7) = trc2d(ji,jj,7) + (fdpd * fthk) !! diatom non-grazing losses !! diagnostic field 8 is (ostensibly) supplied by trcsed.F trc2d(ji,jj,9) = trc2d(ji,jj,9) + (fprds * zpds * fthk) !! diatom silicon production trc2d(ji,jj,10) = trc2d(ji,jj,10) + (fsdiss * fthk) !! diatom silicon dissolution trc2d(ji,jj,11) = trc2d(ji,jj,11) + (fgmipn * fthk) !! microzoo grazing on non-diatoms trc2d(ji,jj,12) = trc2d(ji,jj,12) + (fgmid * fthk) !! microzoo grazing on detrital nitrogen trc2d(ji,jj,13) = trc2d(ji,jj,13) + (fdzmi * fthk) !! microzoo non-grazing losses trc2d(ji,jj,14) = trc2d(ji,jj,14) + (fgmepn * fthk) !! mesozoo grazing on non-diatoms trc2d(ji,jj,15) = trc2d(ji,jj,15) + (fgmepd * fthk) !! mesozoo grazing on diatoms trc2d(ji,jj,16) = trc2d(ji,jj,16) + (fgmezmi * fthk) !! mesozoo grazing on microzoo trc2d(ji,jj,17) = trc2d(ji,jj,17) + (fgmed * fthk) !! mesozoo grazing on detrital nitrogen trc2d(ji,jj,18) = trc2d(ji,jj,18) + (fdzme * fthk) !! mesozoo non-grazing losses !! diagnostic field 19 is (ostensibly) supplied by trcexp.F trc2d(ji,jj,20) = trc2d(ji,jj,20) + (fslown * fthk) !! slow sinking detritus N production trc2d(ji,jj,21) = trc2d(ji,jj,21) + (fdd * fthk) !! detrital remineralisation trc2d(ji,jj,22) = trc2d(ji,jj,22) + (ffetop * fthk) !! aeolian iron addition trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot * fthk) !! seafloor iron addition trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav * fthk) !! "free" iron scavenging trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn * zphn * fthk) !! non-diatom J limitation term trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln * zphn * fthk) !! non-diatom N limitation term trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln * zphn * fthk) !! non-diatom Fe limitation term trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd * zphd * fthk) !! diatom J limitation term trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld * zphd * fthk) !! diatom N limitation term trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld * zphd * fthk) !! diatom Fe limitation term trc2d(ji,jj,31) = trc2d(ji,jj,31) + (fsld2 * zphd * fthk) !! diatom Si limitation term trc2d(ji,jj,32) = trc2d(ji,jj,32) + (fsld * zphd * fthk) !! diatom Si uptake limitation term if (jk.eq.i0100) trc2d(ji,jj,33) = fslownflux(ji,jj) !! slow detritus flux at 100 m if (jk.eq.i0200) trc2d(ji,jj,34) = fslownflux(ji,jj) !! slow detritus flux at 200 m if (jk.eq.i0500) trc2d(ji,jj,35) = fslownflux(ji,jj) !! slow detritus flux at 500 m if (jk.eq.i1000) trc2d(ji,jj,36) = fslownflux(ji,jj) !! slow detritus flux at 1000 m trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen !! non-fast N full column regeneration trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi !! non-fast Si full column regeneration if (jk.eq.i0100) trc2d(ji,jj,39) = trc2d(ji,jj,37) !! non-fast N regeneration to 100 m if (jk.eq.i0200) trc2d(ji,jj,40) = trc2d(ji,jj,37) !! non-fast N regeneration to 200 m if (jk.eq.i0500) trc2d(ji,jj,41) = trc2d(ji,jj,37) !! non-fast N regeneration to 500 m if (jk.eq.i1000) trc2d(ji,jj,42) = trc2d(ji,jj,37) !! non-fast N regeneration to 1000 m trc2d(ji,jj,43) = trc2d(ji,jj,43) + (ftempn * fthk) !! fast sinking detritus N production trc2d(ji,jj,44) = trc2d(ji,jj,44) + (ftempsi * fthk) !! fast sinking detritus Si production trc2d(ji,jj,45) = trc2d(ji,jj,45) + (ftempfe * fthk) !! fast sinking detritus Fe production trc2d(ji,jj,46) = trc2d(ji,jj,46) + (ftempc * fthk) !! fast sinking detritus C production trc2d(ji,jj,47) = trc2d(ji,jj,47) + (ftempca * fthk) !! fast sinking detritus CaCO3 production if (jk.eq.i0100) trc2d(ji,jj,48) = ffastn(ji,jj) !! fast detritus N flux at 100 m if (jk.eq.i0200) trc2d(ji,jj,49) = ffastn(ji,jj) !! fast detritus N flux at 200 m if (jk.eq.i0500) trc2d(ji,jj,50) = ffastn(ji,jj) !! fast detritus N flux at 500 m if (jk.eq.i1000) trc2d(ji,jj,51) = ffastn(ji,jj) !! fast detritus N flux at 1000 m if (jk.eq.i0100) trc2d(ji,jj,52) = fregenfast(ji,jj) !! N regeneration to 100 m if (jk.eq.i0200) trc2d(ji,jj,53) = fregenfast(ji,jj) !! N regeneration to 200 m if (jk.eq.i0500) trc2d(ji,jj,54) = fregenfast(ji,jj) !! N regeneration to 500 m if (jk.eq.i1000) trc2d(ji,jj,55) = fregenfast(ji,jj) !! N regeneration to 1000 m if (jk.eq.i0100) trc2d(ji,jj,56) = ffastsi(ji,jj) !! fast detritus Si flux at 100 m if (jk.eq.i0200) trc2d(ji,jj,57) = ffastsi(ji,jj) !! fast detritus Si flux at 200 m if (jk.eq.i0500) trc2d(ji,jj,58) = ffastsi(ji,jj) !! fast detritus Si flux at 500 m if (jk.eq.i1000) trc2d(ji,jj,59) = ffastsi(ji,jj) !! fast detritus Si flux at 1000 m if (jk.eq.i0100) trc2d(ji,jj,60) = fregenfastsi(ji,jj) !! Si regeneration to 100 m if (jk.eq.i0200) trc2d(ji,jj,61) = fregenfastsi(ji,jj) !! Si regeneration to 200 m if (jk.eq.i0500) trc2d(ji,jj,62) = fregenfastsi(ji,jj) !! Si regeneration to 500 m if (jk.eq.i1000) trc2d(ji,jj,63) = fregenfastsi(ji,jj) !! Si regeneration to 1000 m trc2d(ji,jj,64) = trc2d(ji,jj,64) + (freminn * fthk) !! sum of fast-sinking N fluxes trc2d(ji,jj,65) = trc2d(ji,jj,65) + (freminsi * fthk) !! sum of fast-sinking Si fluxes trc2d(ji,jj,66) = trc2d(ji,jj,66) + (freminfe * fthk) !! sum of fast-sinking Fe fluxes trc2d(ji,jj,67) = trc2d(ji,jj,67) + (freminc * fthk) !! sum of fast-sinking C fluxes trc2d(ji,jj,68) = trc2d(ji,jj,68) + (freminca * fthk) !! sum of fast-sinking Ca fluxes if (jk.eq.jmbathy) then trc2d(ji,jj,69) = fsedn(ji,jj) !! N sedimentation flux trc2d(ji,jj,70) = fsedsi(ji,jj) !! Si sedimentation flux trc2d(ji,jj,71) = fsedfe(ji,jj) !! Fe sedimentation flux trc2d(ji,jj,72) = fsedc(ji,jj) !! C sedimentation flux trc2d(ji,jj,73) = fsedca(ji,jj) !! Ca sedimentation flux endif if (jk.eq.1) trc2d(ji,jj,74) = qsr(ji,jj) if (jk.eq.1) trc2d(ji,jj,75) = xpar(ji,jj,jk) !! if (jk.eq.1) trc2d(ji,jj,75) = real(iters) !! diagnostic fields 76 to 80 calculated below trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj) !! mixed layer non-diatom production trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj) !! mixed layer diatom production # if defined key_gulf_finland if (jk.eq.1) trc2d(ji,jj,83) = real(ibio_switch) !! Gulf of Finland check # else trc2d(ji,jj,83) = ocal_ccd(ji,jj) !! calcite CCD depth # endif trc2d(ji,jj,84) = fccd(ji,jj) !! last model level above calcite CCD depth if (jk.eq.1) trc2d(ji,jj,85) = xFree(ji,jj) !! surface "free" iron if (jk.eq.i0200) trc2d(ji,jj,86) = xFree(ji,jj) !! "free" iron at 100 m if (jk.eq.i0200) trc2d(ji,jj,87) = xFree(ji,jj) !! "free" iron at 200 m if (jk.eq.i0500) trc2d(ji,jj,88) = xFree(ji,jj) !! "free" iron at 500 m if (jk.eq.i1000) trc2d(ji,jj,89) = xFree(ji,jj) !! "free" iron at 1000 m !! AXY (27/06/12): extract "euphotic depth" if (jk.eq.1) trc2d(ji,jj,90) = xze(ji,jj) !! # if defined key_roam !! ROAM provisionally has access to a further 20 2D diagnostics if (jk .eq. 1) then trc2d(ji,jj,91) = trc2d(ji,jj,91) + f_wind !! surface wind trc2d(ji,jj,92) = trc2d(ji,jj,92) + f_pco2atm !! atmospheric pCO2 trc2d(ji,jj,93) = trc2d(ji,jj,93) + f_ph !! ocean pH trc2d(ji,jj,94) = trc2d(ji,jj,94) + f_pco2w !! ocean pCO2 trc2d(ji,jj,95) = trc2d(ji,jj,95) + f_h2co3 !! ocean H2CO3 conc. trc2d(ji,jj,96) = trc2d(ji,jj,96) + f_hco3 !! ocean HCO3 conc. trc2d(ji,jj,97) = trc2d(ji,jj,97) + f_co3 !! ocean CO3 conc. trc2d(ji,jj,98) = trc2d(ji,jj,98) + f_co2flux !! air-sea CO2 flux trc2d(ji,jj,99) = trc2d(ji,jj,99) + f_omcal(ji,jj) !! ocean omega calcite trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj) !! ocean omega aragonite trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC !! ocean TDIC trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK !! ocean TALK trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660 !! surface kw660 trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0 !! surface pressure trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux !! air-sea O2 flux trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat !! ocean O2 saturation trc2d(ji,jj,107) = f2_ccd_cal(ji,jj) !! depth calcite CCD trc2d(ji,jj,108) = f2_ccd_arg(ji,jj) !! depth aragonite CCD endif if (jk .eq. jmbathy) then trc2d(ji,jj,109) = f3_omcal(ji,jj,jk) !! seafloor omega calcite trc2d(ji,jj,110) = f3_omarg(ji,jj,jk) !! seafloor omega aragonite endif !! diagnostic fields 111 to 117 calculated below if (jk.eq.i0100) trc2d(ji,jj,118) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) !! rain ratio at 100 m if (jk.eq.i0500) trc2d(ji,jj,119) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) !! rain ratio at 500 m if (jk.eq.i1000) trc2d(ji,jj,120) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) !! rain ratio at 1000 m !! AXY (18/01/12): benthic flux diagnostics if (jk.eq.jmbathy) then trc2d(ji,jj,121) = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj) trc2d(ji,jj,123) = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) trc2d(ji,jj,124) = f_fbenin_si(ji,jj) trc2d(ji,jj,125) = f_fbenin_ca(ji,jj) trc2d(ji,jj,126) = f_benout_n(ji,jj) trc2d(ji,jj,127) = f_benout_fe(ji,jj) trc2d(ji,jj,128) = f_benout_c(ji,jj) trc2d(ji,jj,129) = f_benout_si(ji,jj) trc2d(ji,jj,130) = f_benout_ca(ji,jj) endif !! diagnostics fields 131 to 135 calculated below trc2d(ji,jj,136) = f_runoff(ji,jj) !! AXY (19/07/12): amended to allow for riverine nutrient addition below surface trc2d(ji,jj,137) = trc2d(ji,jj,137) + (f_riv_loc_n * fthk) trc2d(ji,jj,138) = trc2d(ji,jj,138) + (f_riv_loc_si * fthk) trc2d(ji,jj,139) = trc2d(ji,jj,139) + (f_riv_loc_c * fthk) trc2d(ji,jj,140) = trc2d(ji,jj,140) + (f_riv_loc_alk * fthk) trc2d(ji,jj,141) = trc2d(ji,jj,141) + (fslowc * fthk) !! slow sinking detritus C production if (jk.eq.i0100) trc2d(ji,jj,142) = fslowcflux(ji,jj) !! slow detritus flux at 100 m if (jk.eq.i0200) trc2d(ji,jj,143) = fslowcflux(ji,jj) !! slow detritus flux at 200 m if (jk.eq.i0500) trc2d(ji,jj,144) = fslowcflux(ji,jj) !! slow detritus flux at 500 m if (jk.eq.i1000) trc2d(ji,jj,145) = fslowcflux(ji,jj) !! slow detritus flux at 1000 m trc2d(ji,jj,146) = trc2d(ji,jj,146) + ftot_c(ji,jj) !! carbon inventory trc2d(ji,jj,147) = trc2d(ji,jj,147) + ftot_a(ji,jj) !! alkalinity inventory trc2d(ji,jj,148) = trc2d(ji,jj,148) + ftot_o2(ji,jj) !! oxygen inventory if (jk.eq.jmbathy) then trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj) endif trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fthk !! community respiration !! !! AXY (14/02/14): a Valentines Day gift to BASIN - a shedload of new !! diagnostics that they'll most likely never need! !! (actually, as with all such gifts, I'm giving them !! some things I'd like myself!) !! !! ---------------------------------------------------------------------- !! linear losses !! non-diatom trc2d(ji,jj,151) = trc2d(ji,jj,151) + (fdpn2 * fthk) !! diatom trc2d(ji,jj,152) = trc2d(ji,jj,152) + (fdpd2 * fthk) !! microzooplankton trc2d(ji,jj,153) = trc2d(ji,jj,153) + (fdzmi2 * fthk) !! mesozooplankton trc2d(ji,jj,154) = trc2d(ji,jj,154) + (fdzme2 * fthk) !! ---------------------------------------------------------------------- !! microzooplankton grazing !! microzooplankton messy -> N trc2d(ji,jj,155) = trc2d(ji,jj,155) + (xphi * (fgmipn + fgmid) * fthk) !! microzooplankton messy -> D trc2d(ji,jj,156) = trc2d(ji,jj,156) + ((1. - xbetan) * finmi * fthk) !! microzooplankton messy -> DIC trc2d(ji,jj,157) = trc2d(ji,jj,157) + (xphi * ((xthetapn * fgmipn) + fgmidc) * fthk) !! microzooplankton messy -> Dc trc2d(ji,jj,158) = trc2d(ji,jj,158) + ((1. - xbetac) * ficmi * fthk) !! microzooplankton excretion trc2d(ji,jj,159) = trc2d(ji,jj,159) + (fmiexcr * fthk) !! microzooplankton respiration trc2d(ji,jj,160) = trc2d(ji,jj,160) + (fmiresp * fthk) !! microzooplankton growth trc2d(ji,jj,161) = trc2d(ji,jj,161) + (fmigrow * fthk) !! ---------------------------------------------------------------------- !! mesozooplankton grazing !! mesozooplankton messy -> N trc2d(ji,jj,162) = trc2d(ji,jj,162) + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk) !! mesozooplankton messy -> D trc2d(ji,jj,163) = trc2d(ji,jj,163) + ((1. - xbetan) * finme * fthk) !! mesozooplankton messy -> DIC trc2d(ji,jj,164) = trc2d(ji,jj,164) + (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & & (xthetazmi * fgmezmi) + fgmedc) * fthk) !! mesozooplankton messy -> Dc trc2d(ji,jj,165) = trc2d(ji,jj,165) + ((1. - xbetac) * ficme * fthk) !! mesozooplankton excretion trc2d(ji,jj,166) = trc2d(ji,jj,166) + (fmeexcr * fthk) !! mesozooplankton respiration trc2d(ji,jj,167) = trc2d(ji,jj,167) + (fmeresp * fthk) !! mesozooplankton growth trc2d(ji,jj,168) = trc2d(ji,jj,168) + (fmegrow * fthk) !! ---------------------------------------------------------------------- !! miscellaneous trc2d(ji,jj,169) = trc2d(ji,jj,169) + (fddc * fthk) !! detrital C remineralisation trc2d(ji,jj,170) = trc2d(ji,jj,170) + (fgmidc * fthk) !! microzoo grazing on detrital carbon trc2d(ji,jj,171) = trc2d(ji,jj,171) + (fgmedc * fthk) !! mesozoo grazing on detrital carbon !! !! ---------------------------------------------------------------------- !! !! AXY (23/10/14): extract primary production related surface fields to !! deal with diel cycle issues; hijacking BASIN 150m !! diagnostics to do so (see commented out diagnostics !! below this section) !! !! extract fields at surface !! if (jk .eq. 1) then !! trc2d(ji,jj,172) = zchn !! Pn chlorophyll !! trc2d(ji,jj,173) = zphn !! Pn biomass !! trc2d(ji,jj,174) = fjln !! Pn J-term !! trc2d(ji,jj,175) = (fprn * zphn) !! Pn PP !! trc2d(ji,jj,176) = zchd !! Pd chlorophyll !! trc2d(ji,jj,177) = zphd !! Pd biomass !! trc2d(ji,jj,178) = fjld !! Pd J-term !! trc2d(ji,jj,179) = xpar(ji,jj,jk) !! Pd PP !! trc2d(ji,jj,180) = loc_T !! local temperature !! endif !! !! !! !! extract fields at 50m (actually 44-50m) !! if (jk .eq. 18) then !! trc2d(ji,jj,181) = zchn !! Pn chlorophyll !! trc2d(ji,jj,182) = zphn !! Pn biomass !! trc2d(ji,jj,183) = fjln !! Pn J-term !! trc2d(ji,jj,184) = (fprn * zphn) !! Pn PP !! trc2d(ji,jj,185) = zchd !! Pd chlorophyll !! trc2d(ji,jj,186) = zphd !! Pd biomass !! trc2d(ji,jj,187) = fjld !! Pd J-term !! trc2d(ji,jj,188) = xpar(ji,jj,jk) !! Pd PP !! trc2d(ji,jj,189) = loc_T !! local temperature !! endif !! !! !! !! extract fields at 100m !! if (jk .eq. i0100) then !! trc2d(ji,jj,190) = zchn !! Pn chlorophyll !! trc2d(ji,jj,191) = zphn !! Pn biomass !! trc2d(ji,jj,192) = fjln !! Pn J-term !! trc2d(ji,jj,193) = (fprn * zphn) !! Pn PP !! trc2d(ji,jj,194) = zchd !! Pd chlorophyll !! trc2d(ji,jj,195) = zphd !! Pd biomass !! trc2d(ji,jj,196) = fjld !! Pd J-term !! trc2d(ji,jj,197) = xpar(ji,jj,jk) !! Pd PP !! trc2d(ji,jj,198) = loc_T !! local temperature !! endif !! !! extract relevant BASIN fields at 150m if (jk .eq. i0150) then trc2d(ji,jj,172) = trc2d(ji,jj,4) !! Pn PP trc2d(ji,jj,173) = trc2d(ji,jj,151) !! Pn linear loss trc2d(ji,jj,174) = trc2d(ji,jj,5) !! Pn non-linear loss trc2d(ji,jj,175) = trc2d(ji,jj,11) !! Pn grazing to Zmi trc2d(ji,jj,176) = trc2d(ji,jj,14) !! Pn grazing to Zme trc2d(ji,jj,177) = trc2d(ji,jj,6) !! Pd PP trc2d(ji,jj,178) = trc2d(ji,jj,152) !! Pd linear loss trc2d(ji,jj,179) = trc2d(ji,jj,7) !! Pd non-linear loss trc2d(ji,jj,180) = trc2d(ji,jj,15) !! Pd grazing to Zme trc2d(ji,jj,181) = trc2d(ji,jj,12) !! Zmi grazing on D trc2d(ji,jj,182) = trc2d(ji,jj,170) !! Zmi grazing on Dc trc2d(ji,jj,183) = trc2d(ji,jj,155) !! Zmi messy feeding loss to N trc2d(ji,jj,184) = trc2d(ji,jj,156) !! Zmi messy feeding loss to D trc2d(ji,jj,185) = trc2d(ji,jj,157) !! Zmi messy feeding loss to DIC trc2d(ji,jj,186) = trc2d(ji,jj,158) !! Zmi messy feeding loss to Dc trc2d(ji,jj,187) = trc2d(ji,jj,159) !! Zmi excretion trc2d(ji,jj,188) = trc2d(ji,jj,160) !! Zmi respiration trc2d(ji,jj,189) = trc2d(ji,jj,161) !! Zmi growth trc2d(ji,jj,190) = trc2d(ji,jj,153) !! Zmi linear loss trc2d(ji,jj,191) = trc2d(ji,jj,13) !! Zmi non-linear loss trc2d(ji,jj,192) = trc2d(ji,jj,16) !! Zmi grazing to Zme trc2d(ji,jj,193) = trc2d(ji,jj,17) !! Zme grazing on D trc2d(ji,jj,194) = trc2d(ji,jj,171) !! Zme grazing on Dc trc2d(ji,jj,195) = trc2d(ji,jj,162) !! Zme messy feeding loss to N trc2d(ji,jj,196) = trc2d(ji,jj,163) !! Zme messy feeding loss to D trc2d(ji,jj,197) = trc2d(ji,jj,164) !! Zme messy feeding loss to DIC trc2d(ji,jj,198) = trc2d(ji,jj,165) !! Zme messy feeding loss to Dc trc2d(ji,jj,199) = trc2d(ji,jj,166) !! Zme excretion trc2d(ji,jj,200) = trc2d(ji,jj,167) !! Zme respiration trc2d(ji,jj,201) = trc2d(ji,jj,168) !! Zme growth trc2d(ji,jj,202) = trc2d(ji,jj,154) !! Zme linear loss trc2d(ji,jj,203) = trc2d(ji,jj,18) !! Zme non-linear loss trc2d(ji,jj,204) = trc2d(ji,jj,20) !! Slow detritus production, N trc2d(ji,jj,205) = trc2d(ji,jj,21) !! Slow detritus remineralisation, N trc2d(ji,jj,206) = trc2d(ji,jj,141) !! Slow detritus production, C trc2d(ji,jj,207) = trc2d(ji,jj,169) !! Slow detritus remineralisation, C trc2d(ji,jj,208) = trc2d(ji,jj,43) !! Fast detritus production, N trc2d(ji,jj,209) = trc2d(ji,jj,21) !! Fast detritus remineralisation, N trc2d(ji,jj,210) = trc2d(ji,jj,64) !! Fast detritus production, C trc2d(ji,jj,211) = trc2d(ji,jj,67) !! Fast detritus remineralisation, C trc2d(ji,jj,212) = trc2d(ji,jj,150) !! Community respiration trc2d(ji,jj,213) = fslownflux(ji,jj) !! Slow detritus N flux at 150 m trc2d(ji,jj,214) = fslowcflux(ji,jj) !! Slow detritus C flux at 150 m trc2d(ji,jj,215) = ffastn(ji,jj) !! Fast detritus N flux at 150 m trc2d(ji,jj,216) = ffastc(ji,jj) !! Fast detritus C flux at 150 m endif !! !! Jpalm (11-08-2014) !! Add UKESM1 diagnoatics !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ if ((jk .eq. 1) .and.( jdms.eq.1)) then trc2d(ji,jj,221) = dms_surf !! DMS surface concentration !! AXY (13/03/15): add in other DMS estimates trc2d(ji,jj,222) = dms_andr !! DMS surface concentration trc2d(ji,jj,223) = dms_simo !! DMS surface concentration trc2d(ji,jj,224) = dms_aran !! DMS surface concentration trc2d(ji,jj,225) = dms_hall !! DMS surface concentration endif # endif !! other possible future diagnostics include: !! - integrated tracer values (esp. biological) !! - mixed layer tracer values !! - sub-surface chlorophyll maxima (plus depth) !! - different mixed layer depth criteria (T, sigma, var. sigma) !!---------------------------------------------------------------------- !! Prepare 3D diagnostics !!---------------------------------------------------------------------- !! trc3d(ji,jj,jk,1) = ((fprn + fprd) * zphn) !! primary production trc3d(ji,jj,jk,2) = fslownflux(ji,jj) + ffastn(ji,jj) !! detrital flux trc3d(ji,jj,jk,3) = fregen + (freminn * fthk) !! remineralisation # if defined key_roam trc3d(ji,jj,jk,4) = f3_pH(ji,jj,jk) !! pH trc3d(ji,jj,jk,5) = f3_omcal(ji,jj,jk) !! omega calcite # else trc3d(ji,jj,jk,4) = ffastsi(ji,jj) !! fast Si flux # endif ENDIF ! end of ln_diatrc option !! CLOSE wet point IF..THEN loop endif !! CLOSE horizontal loops ENDDO ENDDO !! IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN !! first - 2D diag implemented !! on every K level !!----------------------------------------- !! -- !!second - 2d specific k level diags !! !!----------------------------------------- IF (jk.eq.1) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 1' CALL flush(numout) # endif !! JPALM -- 02-06-2017 -- !! add Chl surf coupling !! no need to output, just pass to cpl var IF (lk_oasis) THEN zn_chl_srf(:,:) = (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6 !! surf Chl in Kg-chl/m3 as needed for cpl chloro_out_cpl(:,:) = zn_chl_srf(:,:) !! Coupling Chl END IF IF( med_diag%MED_QSR%dgsave ) THEN CALL iom_put( "MED_QSR" , qsr ) ! ENDIF IF( med_diag%MED_XPAR%dgsave ) THEN CALL iom_put( "MED_XPAR" , xpar(:,:,jk) ) ! ENDIF IF( med_diag%OCAL_CCD%dgsave ) THEN CALL iom_put( "OCAL_CCD" , ocal_ccd ) ! ENDIF IF( med_diag%FE_0000%dgsave ) THEN CALL iom_put( "FE_0000" , xFree ) ! ENDIF IF( med_diag%MED_XZE%dgsave ) THEN CALL iom_put( "MED_XZE" , xze ) ! ENDIF # if defined key_roam IF( med_diag%WIND%dgsave ) THEN CALL iom_put( "WIND" , wndm ) ENDIF IF( med_diag%ATM_PCO2%dgsave ) THEN CALL iom_put( "ATM_PCO2" , f_pco2a2d ) CALL wrk_dealloc( jpi, jpj, f_pco2a2d ) ENDIF IF( med_diag%OCN_PH%dgsave ) THEN zw2d(:,:) = f3_pH(:,:,jk) CALL iom_put( "OCN_PH" , zw2d ) ENDIF IF( med_diag%OCN_PCO2%dgsave ) THEN CALL iom_put( "OCN_PCO2" , f_pco2w2d ) CALL wrk_dealloc( jpi, jpj, f_pco2w2d ) ENDIF IF( med_diag%OCNH2CO3%dgsave ) THEN zw2d(:,:) = f3_h2co3(:,:,jk) CALL iom_put( "OCNH2CO3" , zw2d ) ENDIF IF( med_diag%OCN_HCO3%dgsave ) THEN zw2d(:,:) = f3_hco3(:,:,jk) CALL iom_put( "OCN_HCO3" , zw2d ) ENDIF IF( med_diag%OCN_CO3%dgsave ) THEN zw2d(:,:) = f3_co3(:,:,jk) CALL iom_put( "OCN_CO3" , zw2d ) ENDIF IF( med_diag%CO2FLUX%dgsave ) THEN CALL iom_put( "CO2FLUX" , f_co2flux2d ) CALL wrk_dealloc( jpi, jpj, f_co2flux2d ) ENDIF !! !! AXY (10/11/16): repeat CO2 flux diagnostic in UKMO/CMIP6 units; this !! both outputs the CO2 flux in specified units and !! sends the resulting field to the coupler !! JPALM (17/11/16): put CO2 flux (fgco2) alloc/unalloc/pass to zn !! out of diag list request CALL lbc_lnk( fgco2(:,:),'T',1. ) IF( med_diag%FGCO2%dgsave ) THEN CALL iom_put( "FGCO2" , fgco2 ) ENDIF !! JPALM (17/11/16): should mv this fgco2 part !! out of lk_iomput loop zb_co2_flx = zn_co2_flx zn_co2_flx = fgco2 IF (lk_oasis) THEN CO2Flux_out_cpl = zn_co2_flx ENDIF CALL wrk_dealloc( jpi, jpj, fgco2 ) !! --- IF( med_diag%OM_CAL%dgsave ) THEN CALL iom_put( "OM_CAL" , f_omcal ) ENDIF IF( med_diag%OM_ARG%dgsave ) THEN CALL iom_put( "OM_ARG" , f_omarg ) ENDIF IF( med_diag%TCO2%dgsave ) THEN CALL iom_put( "TCO2" , f_TDIC2d ) CALL wrk_dealloc( jpi, jpj, f_TDIC2d ) ENDIF IF( med_diag%TALK%dgsave ) THEN CALL iom_put( "TALK" , f_TALK2d ) CALL wrk_dealloc( jpi, jpj, f_TALK2d ) ENDIF IF( med_diag%KW660%dgsave ) THEN CALL iom_put( "KW660" , f_kw6602d ) CALL wrk_dealloc( jpi, jpj, f_kw6602d ) ENDIF IF( med_diag%ATM_PP0%dgsave ) THEN CALL iom_put( "ATM_PP0" , f_pp02d ) CALL wrk_dealloc( jpi, jpj, f_pp02d ) ENDIF IF( med_diag%O2FLUX%dgsave ) THEN CALL iom_put( "O2FLUX" , f_o2flux2d ) CALL wrk_dealloc( jpi, jpj, f_o2flux2d ) ENDIF IF( med_diag%O2SAT%dgsave ) THEN CALL iom_put( "O2SAT" , f_o2sat2d ) CALL wrk_dealloc( jpi, jpj, f_o2sat2d ) ENDIF IF( med_diag%CAL_CCD%dgsave ) THEN CALL iom_put( "CAL_CCD" , f2_ccd_cal ) ENDIF IF( med_diag%ARG_CCD%dgsave ) THEN CALL iom_put( "ARG_CCD" , f2_ccd_arg ) ENDIF IF (jdms .eq. 1) THEN IF( med_diag%DMS_SURF%dgsave ) THEN CALL lbc_lnk(dms_surf2d(:,:),'T',1. ) CALL iom_put( "DMS_SURF" , dms_surf2d ) zb_dms_srf = zn_dms_srf zn_dms_srf = dms_surf2d IF (lk_oasis) THEN DMS_out_cpl = zn_dms_srf ENDIF CALL wrk_dealloc( jpi, jpj, dms_surf2d ) ENDIF IF( med_diag%DMS_ANDR%dgsave ) THEN CALL iom_put( "DMS_ANDR" , dms_andr2d ) CALL wrk_dealloc( jpi, jpj, dms_andr2d ) ENDIF IF( med_diag%DMS_SIMO%dgsave ) THEN CALL iom_put( "DMS_SIMO" , dms_simo2d ) CALL wrk_dealloc( jpi, jpj, dms_simo2d ) ENDIF IF( med_diag%DMS_ARAN%dgsave ) THEN CALL iom_put( "DMS_ARAN" , dms_aran2d ) CALL wrk_dealloc( jpi, jpj, dms_aran2d ) ENDIF IF( med_diag%DMS_HALL%dgsave ) THEN CALL iom_put( "DMS_HALL" , dms_hall2d ) CALL wrk_dealloc( jpi, jpj, dms_hall2d ) ENDIF IF( med_diag%DMS_ANDM%dgsave ) THEN CALL iom_put( "DMS_ANDM" , dms_andm2d ) CALL wrk_dealloc( jpi, jpj, dms_andm2d ) ENDIF ENDIF !! AXY (24/11/16): extra MOCSY diagnostics IF( med_diag%ATM_XCO2%dgsave ) THEN CALL iom_put( "ATM_XCO2" , f_xco2a_2d ) CALL wrk_dealloc( jpi, jpj, f_xco2a_2d ) ENDIF IF( med_diag%OCN_FCO2%dgsave ) THEN CALL iom_put( "OCN_FCO2" , f_fco2w_2d ) CALL wrk_dealloc( jpi, jpj, f_fco2w_2d ) ENDIF IF( med_diag%ATM_FCO2%dgsave ) THEN CALL iom_put( "ATM_FCO2" , f_fco2a_2d ) CALL wrk_dealloc( jpi, jpj, f_fco2a_2d ) ENDIF IF( med_diag%OCN_RHOSW%dgsave ) THEN CALL iom_put( "OCN_RHOSW" , f_ocnrhosw_2d ) CALL wrk_dealloc( jpi, jpj, f_ocnrhosw_2d ) ENDIF IF( med_diag%OCN_SCHCO2%dgsave ) THEN CALL iom_put( "OCN_SCHCO2" , f_ocnschco2_2d ) CALL wrk_dealloc( jpi, jpj, f_ocnschco2_2d ) ENDIF IF( med_diag%OCN_KWCO2%dgsave ) THEN CALL iom_put( "OCN_KWCO2" , f_ocnkwco2_2d ) CALL wrk_dealloc( jpi, jpj, f_ocnkwco2_2d ) ENDIF IF( med_diag%OCN_K0%dgsave ) THEN CALL iom_put( "OCN_K0" , f_ocnk0_2d ) CALL wrk_dealloc( jpi, jpj, f_ocnk0_2d ) ENDIF IF( med_diag%CO2STARAIR%dgsave ) THEN CALL iom_put( "CO2STARAIR" , f_co2starair_2d ) CALL wrk_dealloc( jpi, jpj, f_co2starair_2d ) ENDIF IF( med_diag%OCN_DPCO2%dgsave ) THEN CALL iom_put( "OCN_DPCO2" , f_ocndpco2_2d ) CALL wrk_dealloc( jpi, jpj, f_ocndpco2_2d ) ENDIF # endif ELSE IF (jk.eq.i0100) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 100' CALL flush(numout) # endif IF( med_diag%SDT__100%dgsave ) THEN zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDT__100" , zw2d ) ENDIF IF( med_diag%REG__100%dgsave ) THEN CALL iom_put( "REG__100" , fregen2d ) ENDIF IF( med_diag%FDT__100%dgsave ) THEN CALL iom_put( "FDT__100" , ffastn ) ENDIF IF( med_diag%RG__100F%dgsave ) THEN CALL iom_put( "RG__100F" , fregenfast ) ENDIF IF( med_diag%FDS__100%dgsave ) THEN CALL iom_put( "FDS__100" , ffastsi ) ENDIF IF( med_diag%RGS_100F%dgsave ) THEN CALL iom_put( "RGS_100F" , fregenfastsi ) ENDIF IF( med_diag%FE_0100%dgsave ) THEN CALL iom_put( "FE_0100" , xFree ) ENDIF # if defined key_roam IF( med_diag%RR_0100%dgsave ) THEN CALL iom_put( "RR_0100" , ffastca2d ) ENDIF IF( med_diag%SDC__100%dgsave ) THEN zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDC__100" , zw2d ) ENDIF IF( med_diag%epC100%dgsave ) THEN zw2d(:,:) = (fslowcflux + ffastc) * tmask(:,:,jk) CALL iom_put( "epC100" , zw2d ) ENDIF IF( med_diag%epCALC100%dgsave ) THEN CALL iom_put( "epCALC100" , ffastca ) ENDIF IF( med_diag%epN100%dgsave ) THEN zw2d(:,:) = (fslownflux + ffastn) * tmask(:,:,jk) CALL iom_put( "epN100" , zw2d ) ENDIF IF( med_diag%epSI100%dgsave ) THEN CALL iom_put( "epSI100" , ffastsi ) ENDIF ELSE IF (jk.eq.i0150) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 150' CALL flush(numout) # endif # endif ELSE IF (jk.eq.i0200) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 200' CALL flush(numout) # endif IF( med_diag%SDT__200%dgsave ) THEN zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDT__200" , zw2d ) ENDIF IF( med_diag%REG__200%dgsave ) THEN CALL iom_put( "REG__200" , fregen2d ) ENDIF IF( med_diag%FDT__200%dgsave ) THEN CALL iom_put( "FDT__200" , ffastn ) ENDIF IF( med_diag%RG__200F%dgsave ) THEN CALL iom_put( "RG__200F" , fregenfast ) ENDIF IF( med_diag%FDS__200%dgsave ) THEN CALL iom_put( "FDS__200" , ffastsi ) ENDIF IF( med_diag%RGS_200F%dgsave ) THEN CALL iom_put( "RGS_200F" , fregenfastsi ) ENDIF IF( med_diag%FE_0200%dgsave ) THEN CALL iom_put( "FE_0200" , xFree ) ENDIF # if defined key_roam IF( med_diag%SDC__200%dgsave ) THEN zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDC__200" , zw2d ) ENDIF # endif ELSE IF (jk.eq.i0500) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 500' CALL flush(numout) # endif IF( med_diag%SDT__500%dgsave ) THEN zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDT__500" , zw2d ) ENDIF IF( med_diag%REG__500%dgsave ) THEN CALL iom_put( "REG__500" , fregen2d ) ENDIF IF( med_diag%FDT__500%dgsave ) THEN CALL iom_put( "FDT__500" , ffastn ) ENDIF IF( med_diag%RG__500F%dgsave ) THEN CALL iom_put( "RG__500F" , fregenfast ) ENDIF IF( med_diag%FDS__500%dgsave ) THEN CALL iom_put( "FDS__500" , ffastsi ) ENDIF IF( med_diag%RGS_500F%dgsave ) THEN CALL iom_put( "RGS_500F" , fregenfastsi ) ENDIF IF( med_diag%FE_0500%dgsave ) THEN CALL iom_put( "FE_0500" , xFree ) ENDIF # if defined key_roam IF( med_diag%RR_0500%dgsave ) THEN CALL iom_put( "RR_0500" , ffastca2d ) ENDIF IF( med_diag%SDC__500%dgsave ) THEN zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDC__500" , zw2d ) ENDIF # endif ELSE IF (jk.eq.i1000) THEN # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 1000' CALL flush(numout) # endif IF( med_diag%SDT_1000%dgsave ) THEN zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDT_1000" , zw2d ) ENDIF IF( med_diag%REG_1000%dgsave ) THEN CALL iom_put( "REG_1000" , fregen2d ) ENDIF IF( med_diag%FDT_1000%dgsave ) THEN CALL iom_put( "FDT_1000" , ffastn ) ENDIF IF( med_diag%RG_1000F%dgsave ) THEN CALL iom_put( "RG_1000F" , fregenfast ) ENDIF IF( med_diag%FDS_1000%dgsave ) THEN CALL iom_put( "FDS_1000" , ffastsi ) ENDIF IF( med_diag%RGS1000F%dgsave ) THEN CALL iom_put( "RGS1000F" , fregenfastsi ) ENDIF IF( med_diag%FE_1000%dgsave ) THEN CALL iom_put( "FE_1000" , xFree ) ENDIF # if defined key_roam IF( med_diag%RR_1000%dgsave ) THEN CALL iom_put( "RR_1000" , ffastca2d ) CALL wrk_dealloc( jpi, jpj, ffastca2d ) ENDIF IF( med_diag%SDC_1000%dgsave ) THEN zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk) CALL iom_put( "SDC_1000" , zw2d ) ENDIF # endif ENDIF !! to do on every k loop : IF( med_diag%DETFLUX3%dgsave ) THEN detflux3d(:,:,jk) = (fslownflux(:,:) + ffastn(:,:)) * tmask(:,:,jk) !! detrital flux !CALL iom_put( "DETFLUX3" , ftot_n ) ENDIF # if defined key_roam IF( med_diag%EXPC3%dgsave ) THEN expc3(:,:,jk) = (fslowcflux(:,:) + ffastc(:,:)) * tmask(:,:,jk) ENDIF IF( med_diag%EXPN3%dgsave ) THEN expn3(:,:,jk) = (fslownflux(:,:) + ffastn(:,:)) * tmask(:,:,jk) ENDIF # endif ENDIF !! CLOSE vertical loop ENDDO !!---------------------------------------------------------------------- !! Process benthic in/out fluxes !! These can be handled outside of the 3D calculations since the !! benthic pools (and fluxes) are 2D in nature; this code is !! (shamelessly) borrowed from corresponding code in the LOBSTER !! model !!---------------------------------------------------------------------- !! !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt if (jorgben.eq.1) then za_sed_n(:,:) = zn_sed_n(:,:) + & & ( f_sbenin_n(:,:) + f_fbenin_n(:,:) - f_benout_n(:,:) ) * (rdt / 86400.) zn_sed_n(:,:) = za_sed_n(:,:) !! za_sed_fe(:,:) = zn_sed_fe(:,:) + & & ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * (rdt / 86400.) zn_sed_fe(:,:) = za_sed_fe(:,:) !! za_sed_c(:,:) = zn_sed_c(:,:) + & & ( f_sbenin_c(:,:) + f_fbenin_c(:,:) - f_benout_c(:,:) ) * (rdt / 86400.) zn_sed_c(:,:) = za_sed_c(:,:) endif if (jinorgben.eq.1) then za_sed_si(:,:) = zn_sed_si(:,:) + & & ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * (rdt / 86400.) zn_sed_si(:,:) = za_sed_si(:,:) !! za_sed_ca(:,:) = zn_sed_ca(:,:) + & & ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * (rdt / 86400.) zn_sed_ca(:,:) = za_sed_ca(:,:) endif IF( ln_diatrc ) THEN DO jj = 2,jpjm1 DO ji = 2,jpim1 trc2d(ji,jj,131) = za_sed_n(ji,jj) trc2d(ji,jj,132) = za_sed_fe(ji,jj) trc2d(ji,jj,133) = za_sed_c(ji,jj) trc2d(ji,jj,134) = za_sed_si(ji,jj) trc2d(ji,jj,135) = za_sed_ca(ji,jj) ENDDO ENDDO !! AXY (07/07/15): temporary hijacking # if defined key_roam !! trc2d(:,:,126) = zn_dms_chn(:,:) !! trc2d(:,:,127) = zn_dms_chd(:,:) !! trc2d(:,:,128) = zn_dms_mld(:,:) !! trc2d(:,:,129) = zn_dms_qsr(:,:) !! trc2d(:,:,130) = zn_dms_din(:,:) # endif ENDIF !! if (ibenthic.eq.2) then !! The code below (in this if ... then ... endif loop) is !! effectively commented out because it does not work as !! anticipated; it can be deleted at a later date if (jorgben.eq.1) then za_sed_n(:,:) = ( f_sbenin_n(:,:) + f_fbenin_n(:,:) - f_benout_n(:,:) ) * rdt za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt za_sed_c(:,:) = ( f_sbenin_c(:,:) + f_fbenin_c(:,:) - f_benout_c(:,:) ) * rdt endif if (jinorgben.eq.1) then za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt endif !! !! Leap-frog scheme - only in explicit case, otherwise the time stepping !! is already being done in trczdf !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN !! zfact = 2. * rdttra(jk) * FLOAT( ndttrc ) !! IF( neuler == 0 .AND. kt == nittrc000 ) zfact = rdttra(jk) * FLOAT(ndttrc) !! if (jorgben.eq.1) then !! za_sed_n(:,:) = zb_sed_n(:,:) + ( zfact * za_sed_n(:,:) ) !! za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) ) !! za_sed_c(:,:) = zb_sed_c(:,:) + ( zfact * za_sed_c(:,:) ) !! endif !! if (jinorgben.eq.1) then !! za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) ) !! za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) ) !! endif !! ENDIF !! !! Time filter and swap of arrays IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN ! centred or tvd scheme IF( neuler == 0 .AND. kt == nittrc000 ) THEN if (jorgben.eq.1) then zb_sed_n(:,:) = zn_sed_n(:,:) zn_sed_n(:,:) = za_sed_n(:,:) za_sed_n(:,:) = 0.0 !! zb_sed_fe(:,:) = zn_sed_fe(:,:) zn_sed_fe(:,:) = za_sed_fe(:,:) za_sed_fe(:,:) = 0.0 !! zb_sed_c(:,:) = zn_sed_c(:,:) zn_sed_c(:,:) = za_sed_c(:,:) za_sed_c(:,:) = 0.0 endif if (jinorgben.eq.1) then zb_sed_si(:,:) = zn_sed_si(:,:) zn_sed_si(:,:) = za_sed_si(:,:) za_sed_si(:,:) = 0.0 !! zb_sed_ca(:,:) = zn_sed_ca(:,:) zn_sed_ca(:,:) = za_sed_ca(:,:) za_sed_ca(:,:) = 0.0 endif ELSE if (jorgben.eq.1) then zb_sed_n(:,:) = (atfp * ( zb_sed_n(:,:) + za_sed_n(:,:) )) + (atfp1 * zn_sed_n(:,:) ) zn_sed_n(:,:) = za_sed_n(:,:) za_sed_n(:,:) = 0.0 !! zb_sed_fe(:,:) = (atfp * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:)) zn_sed_fe(:,:) = za_sed_fe(:,:) za_sed_fe(:,:) = 0.0 !! zb_sed_c(:,:) = (atfp * ( zb_sed_c(:,:) + za_sed_c(:,:) )) + (atfp1 * zn_sed_c(:,:) ) zn_sed_c(:,:) = za_sed_c(:,:) za_sed_c(:,:) = 0.0 endif if (jinorgben.eq.1) then zb_sed_si(:,:) = (atfp * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:)) zn_sed_si(:,:) = za_sed_si(:,:) za_sed_si(:,:) = 0.0 !! zb_sed_ca(:,:) = (atfp * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:)) zn_sed_ca(:,:) = za_sed_ca(:,:) za_sed_ca(:,:) = 0.0 endif ENDIF ELSE ! case of smolar scheme or muscl if (jorgben.eq.1) then zb_sed_n(:,:) = za_sed_n(:,:) zn_sed_n(:,:) = za_sed_n(:,:) za_sed_n(:,:) = 0.0 !! zb_sed_fe(:,:) = za_sed_fe(:,:) zn_sed_fe(:,:) = za_sed_fe(:,:) za_sed_fe(:,:) = 0.0 !! zb_sed_c(:,:) = za_sed_c(:,:) zn_sed_c(:,:) = za_sed_c(:,:) za_sed_c(:,:) = 0.0 endif if (jinorgben.eq.1) then zb_sed_si(:,:) = za_sed_si(:,:) zn_sed_si(:,:) = za_sed_si(:,:) za_sed_si(:,:) = 0.0 !! zb_sed_ca(:,:) = za_sed_ca(:,:) zn_sed_ca(:,:) = za_sed_ca(:,:) za_sed_ca(:,:) = 0.0 endif ENDIF endif IF( ln_diatrc ) THEN !!---------------------------------------------------------------------- !! Output several accumulated diagnostics !! - biomass-average phytoplankton limitation terms !! - integrated tendency terms !!---------------------------------------------------------------------- !! DO jj = 2,jpjm1 DO ji = 2,jpim1 !! non-diatom phytoplankton limitations trc2d(ji,jj,25) = trc2d(ji,jj,25) / MAX(ftot_pn(ji,jj), rsmall) trc2d(ji,jj,26) = trc2d(ji,jj,26) / MAX(ftot_pn(ji,jj), rsmall) trc2d(ji,jj,27) = trc2d(ji,jj,27) / MAX(ftot_pn(ji,jj), rsmall) !! diatom phytoplankton limitations trc2d(ji,jj,28) = trc2d(ji,jj,28) / MAX(ftot_pd(ji,jj), rsmall) trc2d(ji,jj,29) = trc2d(ji,jj,29) / MAX(ftot_pd(ji,jj), rsmall) trc2d(ji,jj,30) = trc2d(ji,jj,30) / MAX(ftot_pd(ji,jj), rsmall) trc2d(ji,jj,31) = trc2d(ji,jj,31) / MAX(ftot_pd(ji,jj), rsmall) trc2d(ji,jj,32) = trc2d(ji,jj,32) / MAX(ftot_pd(ji,jj), rsmall) !! tendency terms trc2d(ji,jj,76) = fflx_n(ji,jj) trc2d(ji,jj,77) = fflx_si(ji,jj) trc2d(ji,jj,78) = fflx_fe(ji,jj) !! integrated biomass trc2d(ji,jj,79) = ftot_pn(ji,jj) !! integrated non-diatom phytoplankton trc2d(ji,jj,80) = ftot_pd(ji,jj) !! integrated diatom phytoplankton trc2d(ji,jj,217) = ftot_zmi(ji,jj) !! Integrated microzooplankton trc2d(ji,jj,218) = ftot_zme(ji,jj) !! Integrated mesozooplankton trc2d(ji,jj,219) = ftot_det(ji,jj) !! Integrated slow detritus, nitrogen trc2d(ji,jj,220) = ftot_dtc(ji,jj) !! Integrated slow detritus, carbon # if defined key_roam !! the balance of nitrogen production/consumption trc2d(ji,jj,111) = fnit_prod(ji,jj) !! integrated nitrogen production trc2d(ji,jj,112) = fnit_cons(ji,jj) !! integrated nitrogen consumption !! the balance of carbon production/consumption trc2d(ji,jj,113) = fcar_prod(ji,jj) !! integrated carbon production trc2d(ji,jj,114) = fcar_cons(ji,jj) !! integrated carbon consumption !! the balance of oxygen production/consumption trc2d(ji,jj,115) = foxy_prod(ji,jj) !! integrated oxygen production trc2d(ji,jj,116) = foxy_cons(ji,jj) !! integrated oxygen consumption trc2d(ji,jj,117) = foxy_anox(ji,jj) !! integrated unrealised oxygen consumption # endif ENDDO ENDDO # if defined key_roam # if defined key_axy_nancheck !!---------------------------------------------------------------------- !! Check for NaNs in diagnostic outputs !!---------------------------------------------------------------------- !! !! 2D diagnostics DO jn = 1,150 fq0 = SUM(trc2d(:,:,jn)) !! AXY (30/01/14): "isnan" problem on HECTOR !! if (fq0 /= fq0 ) then if ( ieee_is_nan( fq0 ) ) then !! there's a NaN here if (lwp) write(numout,*) 'NAN detected in 2D diagnostic field', jn, 'at time', kt, 'at position:' DO jj = 1,jpj DO ji = 1,jpi if ( ieee_is_nan( trc2d(ji,jj,jn) ) ) then if (lwp) write (numout,'(a,3i6)') 'NAN-CHECK', & & ji, jj, jn endif ENDDO ENDDO CALL ctl_stop( 'trcbio_medusa, NAN in 2D diagnostic field' ) endif ENDDO !! !! 3D diagnostics DO jn = 1,5 fq0 = SUM(trc3d(:,:,:,jn)) !! AXY (30/01/14): "isnan" problem on HECTOR !! if (fq0 /= fq0 ) then if ( ieee_is_nan( fq0 ) ) then !! there's a NaN here if (lwp) write(numout,*) 'NAN detected in 3D diagnostic field', jn, 'at time', kt, 'at position:' DO jk = 1,jpk DO jj = 1,jpj DO ji = 1,jpi if ( ieee_is_nan( trc3d(ji,jj,jk,jn) ) ) then if (lwp) write (numout,'(a,4i6)') 'NAN-CHECK', & & ji, jj, jk, jn endif ENDDO ENDDO ENDDO CALL ctl_stop( 'trcbio_medusa, NAN in 3D diagnostic field' ) endif ENDDO CALL flush(numout) # endif # endif !!---------------------------------------------------------------------- !! Don't know what this does; belongs to someone else ... !!---------------------------------------------------------------------- !! !! Lateral boundary conditions on trc2d DO jn=1,jp_medusa_2d CALL lbc_lnk(trc2d(:,:,jn),'T',1. ) ENDDO !! Lateral boundary conditions on trc3d DO jn=1,jp_medusa_3d CALL lbc_lnk(trc3d(:,:,1,jn),'T',1. ) ENDDO # if defined key_axy_nodiag !!---------------------------------------------------------------------- !! Blank diagnostics as a NaN-trap !!---------------------------------------------------------------------- !! !! blank 2D diagnostic array trc2d(:,:,:) = 0.e0 !! !! blank 3D diagnostic array trc3d(:,:,:,:) = 0.e0 # endif !!---------------------------------------------------------------------- !! Add in XML diagnostics stuff !!---------------------------------------------------------------------- !! !! ** 2D diagnostics DO jn=1,jp_medusa_2d CALL iom_put(TRIM(ctrc2d(jn)), trc2d(:,:,jn)) END DO !! AXY (17/02/14): don't think I need this if I modify the above for all diagnostics !! # if defined key_roam !! DO jn=91,jp_medusa_2d !! CALL iom_put(TRIM(ctrc2d(jn)), trc2d(:,:,jn)) !! END DO !! # endif !! !! ** 3D diagnostics DO jn=1,jp_medusa_3d CALL iom_put(TRIM(ctrc3d(jn)), trc3d(:,:,:,jn)) END DO !! AXY (17/02/14): don't think I need this if I modify the above for all diagnostics !! # if defined key_roam !! CALL iom_put(TRIM(ctrc3d(5)), trc3d(:,:,:,5)) !! # endif ELSE IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN !!!---------------------------------------------------------------------- !! Add very last diag calculations !!!---------------------------------------------------------------------- DO jj = 2,jpjm1 DO ji = 2,jpim1 !! IF( med_diag%PN_JLIM%dgsave ) THEN fjln2d(ji,jj) = fjln2d(ji,jj) / MAX(ftot_pn(ji,jj), rsmall) ENDIF IF( med_diag%PN_NLIM%dgsave ) THEN fnln2d(ji,jj) = fnln2d(ji,jj) / MAX(ftot_pn(ji,jj), rsmall) ENDIF IF( med_diag%PN_FELIM%dgsave ) THEN ffln2d(ji,jj) = ffln2d(ji,jj) / MAX(ftot_pn(ji,jj), rsmall) ENDIF IF( med_diag%PD_JLIM%dgsave ) THEN fjld2d(ji,jj) = fjld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall) ENDIF IF( med_diag%PD_NLIM%dgsave ) THEN fnld2d(ji,jj) = fnld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall) ENDIF IF( med_diag%PD_FELIM%dgsave ) THEN ffld2d(ji,jj) = ffld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall) ENDIF IF( med_diag%PD_SILIM%dgsave ) THEN fsld2d2(ji,jj) = fsld2d2(ji,jj) / MAX(ftot_pd(ji,jj), rsmall) ENDIF IF( med_diag%PDSILIM2%dgsave ) THEN fsld2d(ji,jj) = fsld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall) ENDIF ENDDO ENDDO !!---------------------------------------------------------------------- !! Add in XML diagnostics stuff !!---------------------------------------------------------------------- !! !! ** 2D diagnostics # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: export all diag.' CALL flush(numout) # endif IF ( med_diag%INVTN%dgsave ) THEN CALL iom_put( "INVTN" , ftot_n ) ENDIF IF ( med_diag%INVTSI%dgsave ) THEN CALL iom_put( "INVTSI" , ftot_si ) ENDIF IF ( med_diag%INVTFE%dgsave ) THEN CALL iom_put( "INVTFE" , ftot_fe ) ENDIF IF ( med_diag%ML_PRN%dgsave ) THEN CALL iom_put( "ML_PRN" , fprn_ml ) ENDIF IF ( med_diag%ML_PRD%dgsave ) THEN CALL iom_put( "ML_PRD" , fprd_ml ) ENDIF IF ( med_diag%OCAL_LVL%dgsave ) THEN CALL iom_put( "OCAL_LVL" , fccd ) ENDIF IF ( med_diag%PN_JLIM%dgsave ) THEN CALL iom_put( "PN_JLIM" , fjln2d ) CALL wrk_dealloc( jpi, jpj, fjln2d ) ENDIF IF ( med_diag%PN_NLIM%dgsave ) THEN CALL iom_put( "PN_NLIM" , fnln2d ) CALL wrk_dealloc( jpi, jpj, fnln2d ) ENDIF IF ( med_diag%PN_FELIM%dgsave ) THEN CALL iom_put( "PN_FELIM" , ffln2d ) CALL wrk_dealloc( jpi, jpj, ffln2d ) ENDIF IF ( med_diag%PD_JLIM%dgsave ) THEN CALL iom_put( "PD_JLIM" , fjld2d ) CALL wrk_dealloc( jpi, jpj, fjld2d ) ENDIF IF ( med_diag%PD_NLIM%dgsave ) THEN CALL iom_put( "PD_NLIM" , fnld2d ) CALL wrk_dealloc( jpi, jpj, fnld2d ) ENDIF IF ( med_diag%PD_FELIM%dgsave ) THEN CALL iom_put( "PD_FELIM" , ffld2d ) CALL wrk_dealloc( jpi, jpj, ffld2d ) ENDIF IF ( med_diag%PD_SILIM%dgsave ) THEN CALL iom_put( "PD_SILIM" , fsld2d2 ) CALL wrk_dealloc( jpi, jpj, fsld2d2 ) ENDIF IF ( med_diag%PDSILIM2%dgsave ) THEN CALL iom_put( "PDSILIM2" , fsld2d ) CALL wrk_dealloc( jpi, jpj, fsld2d ) ENDIF IF ( med_diag%INTFLX_N%dgsave ) THEN CALL iom_put( "INTFLX_N" , fflx_n ) ENDIF IF ( med_diag%INTFLX_SI%dgsave ) THEN CALL iom_put( "INTFLX_SI" , fflx_si ) ENDIF IF ( med_diag%INTFLX_FE%dgsave ) THEN CALL iom_put( "INTFLX_FE" , fflx_fe ) ENDIF IF ( med_diag%INT_PN%dgsave ) THEN CALL iom_put( "INT_PN" , ftot_pn ) ENDIF IF ( med_diag%INT_PD%dgsave ) THEN CALL iom_put( "INT_PD" , ftot_pd ) ENDIF IF ( med_diag%INT_ZMI%dgsave ) THEN CALL iom_put( "INT_ZMI" , ftot_zmi ) ENDIF IF ( med_diag%INT_ZME%dgsave ) THEN CALL iom_put( "INT_ZME" , ftot_zme ) ENDIF IF ( med_diag%INT_DET%dgsave ) THEN CALL iom_put( "INT_DET" , ftot_det ) ENDIF IF ( med_diag%INT_DTC%dgsave ) THEN CALL iom_put( "INT_DTC" , ftot_dtc ) ENDIF IF ( med_diag%BEN_N%dgsave ) THEN CALL iom_put( "BEN_N" , za_sed_n ) ENDIF IF ( med_diag%BEN_FE%dgsave ) THEN CALL iom_put( "BEN_FE" , za_sed_fe ) ENDIF IF ( med_diag%BEN_C%dgsave ) THEN CALL iom_put( "BEN_C" , za_sed_c ) ENDIF IF ( med_diag%BEN_SI%dgsave ) THEN CALL iom_put( "BEN_SI" , za_sed_si ) ENDIF IF ( med_diag%BEN_CA%dgsave ) THEN CALL iom_put( "BEN_CA" , za_sed_ca ) ENDIF IF ( med_diag%RUNOFF%dgsave ) THEN CALL iom_put( "RUNOFF" , f_runoff ) ENDIF # if defined key_roam IF ( med_diag%N_PROD%dgsave ) THEN CALL iom_put( "N_PROD" , fnit_prod ) ENDIF IF ( med_diag%N_CONS%dgsave ) THEN CALL iom_put( "N_CONS" , fnit_cons ) ENDIF IF ( med_diag%C_PROD%dgsave ) THEN CALL iom_put( "C_PROD" , fcar_prod ) ENDIF IF ( med_diag%C_CONS%dgsave ) THEN CALL iom_put( "C_CONS" , fcar_cons ) ENDIF IF ( med_diag%O2_PROD%dgsave ) THEN CALL iom_put( "O2_PROD" , foxy_prod ) ENDIF IF ( med_diag%O2_CONS%dgsave ) THEN CALL iom_put( "O2_CONS" , foxy_cons ) ENDIF IF ( med_diag%O2_ANOX%dgsave ) THEN CALL iom_put( "O2_ANOX" , foxy_anox ) ENDIF IF ( med_diag%INVTC%dgsave ) THEN CALL iom_put( "INVTC" , ftot_c ) ENDIF IF ( med_diag%INVTALK%dgsave ) THEN CALL iom_put( "INVTALK" , ftot_a ) ENDIF IF ( med_diag%INVTO2%dgsave ) THEN CALL iom_put( "INVTO2" , ftot_o2 ) ENDIF IF ( med_diag%COM_RESP%dgsave ) THEN CALL iom_put( "COM_RESP" , fcomm_resp ) ENDIF # endif !! !! diagnostic filled in the i-j-k main loop !!-------------------------------------------- IF ( med_diag%PRN%dgsave ) THEN CALL iom_put( "PRN" , fprn2d ) CALL wrk_dealloc( jpi, jpj, fprn2d ) ENDIF IF ( med_diag%MPN%dgsave ) THEN CALL iom_put( "MPN" ,fdpn2d ) CALL wrk_dealloc( jpi, jpj, fdpn2d ) ENDIF IF ( med_diag%PRD%dgsave ) THEN CALL iom_put( "PRD" ,fprd2d ) CALL wrk_dealloc( jpi, jpj, fprd2d ) ENDIF IF( med_diag%MPD%dgsave ) THEN CALL iom_put( "MPD" , fdpd2d ) CALL wrk_dealloc( jpi, jpj, fdpd2d ) ENDIF ! IF( med_diag%DSED%dgsave ) THEN ! CALL iom_put( "DSED" , ftot_n ) ! ENDIF IF( med_diag%OPAL%dgsave ) THEN CALL iom_put( "OPAL" , fprds2d ) CALL wrk_dealloc( jpi, jpj, fprds2d ) ENDIF IF( med_diag%OPALDISS%dgsave ) THEN CALL iom_put( "OPALDISS" , fsdiss2d ) CALL wrk_dealloc( jpi, jpj, fsdiss2d ) ENDIF IF( med_diag%GMIPn%dgsave ) THEN CALL iom_put( "GMIPn" , fgmipn2d ) CALL wrk_dealloc( jpi, jpj, fgmipn2d ) ENDIF IF( med_diag%GMID%dgsave ) THEN CALL iom_put( "GMID" , fgmid2d ) CALL wrk_dealloc( jpi, jpj, fgmid2d ) ENDIF IF( med_diag%MZMI%dgsave ) THEN CALL iom_put( "MZMI" , fdzmi2d ) CALL wrk_dealloc( jpi, jpj, fdzmi2d ) ENDIF IF( med_diag%GMEPN%dgsave ) THEN CALL iom_put( "GMEPN" , fgmepn2d ) CALL wrk_dealloc( jpi, jpj, fgmepn2d ) ENDIF IF( med_diag%GMEPD%dgsave ) THEN CALL iom_put( "GMEPD" , fgmepd2d ) CALL wrk_dealloc( jpi, jpj, fgmepd2d ) ENDIF IF( med_diag%GMEZMI%dgsave ) THEN CALL iom_put( "GMEZMI" , fgmezmi2d ) CALL wrk_dealloc( jpi, jpj, fgmezmi2d ) ENDIF IF( med_diag%GMED%dgsave ) THEN CALL iom_put( "GMED" , fgmed2d ) CALL wrk_dealloc( jpi, jpj, fgmed2d ) ENDIF IF( med_diag%MZME%dgsave ) THEN CALL iom_put( "MZME" , fdzme2d ) CALL wrk_dealloc( jpi, jpj, fdzme2d ) ENDIF ! IF( med_diag%DEXP%dgsave ) THEN ! CALL iom_put( "DEXP" , ftot_n ) ! ENDIF IF( med_diag%DETN%dgsave ) THEN CALL iom_put( "DETN" , fslown2d ) CALL wrk_dealloc( jpi, jpj, fslown2d ) ENDIF IF( med_diag%MDET%dgsave ) THEN CALL iom_put( "MDET" , fdd2d ) CALL wrk_dealloc( jpi, jpj, fdd2d ) ENDIF IF( med_diag%AEOLIAN%dgsave ) THEN CALL iom_put( "AEOLIAN" , ffetop2d ) CALL wrk_dealloc( jpi, jpj, ffetop2d ) ENDIF IF( med_diag%BENTHIC%dgsave ) THEN CALL iom_put( "BENTHIC" , ffebot2d ) CALL wrk_dealloc( jpi, jpj, ffebot2d ) ENDIF IF( med_diag%SCAVENGE%dgsave ) THEN CALL iom_put( "SCAVENGE" , ffescav2d ) CALL wrk_dealloc( jpi, jpj, ffescav2d ) ENDIF !! IF( med_diag%TOTREG_N%dgsave ) THEN CALL iom_put( "TOTREG_N" , fregen2d ) CALL wrk_dealloc( jpi, jpj, fregen2d ) ENDIF IF( med_diag%TOTRG_SI%dgsave ) THEN CALL iom_put( "TOTRG_SI" , fregensi2d ) CALL wrk_dealloc( jpi, jpj, fregensi2d ) ENDIF !! IF( med_diag%FASTN%dgsave ) THEN CALL iom_put( "FASTN" , ftempn2d ) CALL wrk_dealloc( jpi, jpj, ftempn2d ) ENDIF IF( med_diag%FASTSI%dgsave ) THEN CALL iom_put( "FASTSI" , ftempsi2d ) CALL wrk_dealloc( jpi, jpj, ftempsi2d ) ENDIF IF( med_diag%FASTFE%dgsave ) THEN CALL iom_put( "FASTFE" , ftempfe2d ) CALL wrk_dealloc( jpi, jpj, ftempfe2d ) ENDIF IF( med_diag%FASTC%dgsave ) THEN CALL iom_put( "FASTC" , ftempc2d ) CALL wrk_dealloc( jpi, jpj, ftempc2d ) ENDIF IF( med_diag%FASTCA%dgsave ) THEN CALL iom_put( "FASTCA" , ftempca2d ) CALL wrk_dealloc( jpi, jpj, ftempca2d ) ENDIF !! IF( med_diag%REMINN%dgsave ) THEN CALL iom_put( "REMINN" , freminn2d ) CALL wrk_dealloc( jpi, jpj, freminn2d ) ENDIF IF( med_diag%REMINSI%dgsave ) THEN CALL iom_put( "REMINSI" , freminsi2d ) CALL wrk_dealloc( jpi, jpj, freminsi2d ) ENDIF IF( med_diag%REMINFE%dgsave ) THEN CALL iom_put( "REMINFE" , freminfe2d ) CALL wrk_dealloc( jpi, jpj, freminfe2d ) ENDIF IF( med_diag%REMINC%dgsave ) THEN CALL iom_put( "REMINC" , freminc2d ) CALL wrk_dealloc( jpi, jpj, freminc2d ) ENDIF IF( med_diag%REMINCA%dgsave ) THEN CALL iom_put( "REMINCA" , freminca2d ) CALL wrk_dealloc( jpi, jpj, freminca2d ) ENDIF IF( med_diag%SEAFLRN%dgsave ) THEN CALL iom_put( "SEAFLRN" , fsedn ) ENDIF IF( med_diag%SEAFLRSI%dgsave ) THEN CALL iom_put( "SEAFLRSI" , fsedsi ) ENDIF IF( med_diag%SEAFLRFE%dgsave ) THEN CALL iom_put( "SEAFLRFE" , fsedfe ) ENDIF IF( med_diag%SEAFLRC%dgsave ) THEN CALL iom_put( "SEAFLRC" , fsedc ) ENDIF IF( med_diag%SEAFLRCA%dgsave ) THEN CALL iom_put( "SEAFLRCA" , fsedca ) ENDIF !! # if defined key_roam !! IF( med_diag%RIV_N%dgsave ) THEN CALL iom_put( "RIV_N" , rivn2d ) CALL wrk_dealloc( jpi, jpj, rivn2d ) ENDIF IF( med_diag%RIV_SI%dgsave ) THEN CALL iom_put( "RIV_SI" , rivsi2d ) CALL wrk_dealloc( jpi, jpj, rivsi2d ) ENDIF IF( med_diag%RIV_C%dgsave ) THEN CALL iom_put( "RIV_C" , rivc2d ) CALL wrk_dealloc( jpi, jpj, rivc2d ) ENDIF IF( med_diag%RIV_ALK%dgsave ) THEN CALL iom_put( "RIV_ALK" , rivalk2d ) CALL wrk_dealloc( jpi, jpj, rivalk2d ) ENDIF IF( med_diag%DETC%dgsave ) THEN CALL iom_put( "DETC" , fslowc2d ) CALL wrk_dealloc( jpi, jpj, fslowc2d ) ENDIF !! IF( med_diag%PN_LLOSS%dgsave ) THEN CALL iom_put( "PN_LLOSS" , fdpn22d ) CALL wrk_dealloc( jpi, jpj, fdpn22d ) ENDIF IF( med_diag%PD_LLOSS%dgsave ) THEN CALL iom_put( "PD_LLOSS" , fdpd22d ) CALL wrk_dealloc( jpi, jpj, fdpd22d ) ENDIF IF( med_diag%ZI_LLOSS%dgsave ) THEN CALL iom_put( "ZI_LLOSS" , fdzmi22d ) CALL wrk_dealloc( jpi, jpj, fdzmi22d ) ENDIF IF( med_diag%ZE_LLOSS%dgsave ) THEN CALL iom_put( "ZE_LLOSS" , fdzme22d ) CALL wrk_dealloc( jpi, jpj, fdzme22d ) ENDIF IF( med_diag%ZI_MES_N%dgsave ) THEN CALL iom_put( "ZI_MES_N" , zimesn2d ) CALL wrk_dealloc( jpi, jpj, zimesn2d ) ENDIF IF( med_diag%ZI_MES_D%dgsave ) THEN CALL iom_put( "ZI_MES_D" , zimesd2d ) CALL wrk_dealloc( jpi, jpj, zimesd2d ) ENDIF IF( med_diag%ZI_MES_C%dgsave ) THEN CALL iom_put( "ZI_MES_C" , zimesc2d ) CALL wrk_dealloc( jpi, jpj, zimesc2d ) ENDIF IF( med_diag%ZI_MESDC%dgsave ) THEN CALL iom_put( "ZI_MESDC" ,zimesdc2d ) CALL wrk_dealloc( jpi, jpj, zimesdc2d ) ENDIF IF( med_diag%ZI_EXCR%dgsave ) THEN CALL iom_put( "ZI_EXCR" , ziexcr2d ) CALL wrk_dealloc( jpi, jpj, ziexcr2d ) ENDIF IF( med_diag%ZI_RESP%dgsave ) THEN CALL iom_put( "ZI_RESP" , ziresp2d ) CALL wrk_dealloc( jpi, jpj, ziresp2d ) ENDIF IF( med_diag%ZI_GROW%dgsave ) THEN CALL iom_put( "ZI_GROW" , zigrow2d ) CALL wrk_dealloc( jpi, jpj, zigrow2d ) ENDIF IF( med_diag%ZE_MES_N%dgsave ) THEN CALL iom_put( "ZE_MES_N" , zemesn2d ) CALL wrk_dealloc( jpi, jpj, zemesn2d ) ENDIF IF( med_diag%ZE_MES_D%dgsave ) THEN CALL iom_put( "ZE_MES_D" , zemesd2d ) CALL wrk_dealloc( jpi, jpj, zemesd2d ) ENDIF IF( med_diag%ZE_MES_C%dgsave ) THEN CALL iom_put( "ZE_MES_C" , zemesc2d ) CALL wrk_dealloc( jpi, jpj, zemesc2d ) ENDIF IF( med_diag%ZE_MESDC%dgsave ) THEN CALL iom_put( "ZE_MESDC" , zemesdc2d ) CALL wrk_dealloc( jpi, jpj, zemesdc2d ) ENDIF IF( med_diag%ZE_EXCR%dgsave ) THEN CALL iom_put( "ZE_EXCR" , zeexcr2d ) CALL wrk_dealloc( jpi, jpj, zeexcr2d ) ENDIF IF( med_diag%ZE_RESP%dgsave ) THEN CALL iom_put( "ZE_RESP" , zeresp2d ) CALL wrk_dealloc( jpi, jpj, zeresp2d ) ENDIF IF( med_diag%ZE_GROW%dgsave ) THEN CALL iom_put( "ZE_GROW" , zegrow2d ) CALL wrk_dealloc( jpi, jpj, zegrow2d ) ENDIF IF( med_diag%MDETC%dgsave ) THEN CALL iom_put( "MDETC" , mdetc2d ) CALL wrk_dealloc( jpi, jpj, mdetc2d ) ENDIF IF( med_diag%GMIDC%dgsave ) THEN CALL iom_put( "GMIDC" , gmidc2d ) CALL wrk_dealloc( jpi, jpj, gmidc2d ) ENDIF IF( med_diag%GMEDC%dgsave ) THEN CALL iom_put( "GMEDC" , gmedc2d ) CALL wrk_dealloc( jpi, jpj, gmedc2d ) ENDIF IF( med_diag%IBEN_N%dgsave ) THEN CALL iom_put( "IBEN_N" , iben_n2d ) CALL wrk_dealloc( jpi, jpj, iben_n2d ) ENDIF IF( med_diag%IBEN_FE%dgsave ) THEN CALL iom_put( "IBEN_FE" , iben_fe2d ) CALL wrk_dealloc( jpi, jpj, iben_fe2d ) ENDIF IF( med_diag%IBEN_C%dgsave ) THEN CALL iom_put( "IBEN_C" , iben_c2d ) CALL wrk_dealloc( jpi, jpj, iben_c2d ) ENDIF IF( med_diag%IBEN_SI%dgsave ) THEN CALL iom_put( "IBEN_SI" , iben_si2d ) CALL wrk_dealloc( jpi, jpj, iben_si2d ) ENDIF IF( med_diag%IBEN_CA%dgsave ) THEN CALL iom_put( "IBEN_CA" , iben_ca2d ) CALL wrk_dealloc( jpi, jpj, iben_ca2d ) ENDIF IF( med_diag%OBEN_N%dgsave ) THEN CALL iom_put( "OBEN_N" , oben_n2d ) CALL wrk_dealloc( jpi, jpj, oben_n2d ) ENDIF IF( med_diag%OBEN_FE%dgsave ) THEN CALL iom_put( "OBEN_FE" , oben_fe2d ) CALL wrk_dealloc( jpi, jpj, oben_fe2d ) ENDIF IF( med_diag%OBEN_C%dgsave ) THEN CALL iom_put( "OBEN_C" , oben_c2d ) CALL wrk_dealloc( jpi, jpj, oben_c2d ) ENDIF IF( med_diag%OBEN_SI%dgsave ) THEN CALL iom_put( "OBEN_SI" , oben_si2d ) CALL wrk_dealloc( jpi, jpj, oben_si2d ) ENDIF IF( med_diag%OBEN_CA%dgsave ) THEN CALL iom_put( "OBEN_CA" , oben_ca2d ) CALL wrk_dealloc( jpi, jpj, oben_ca2d ) ENDIF IF( med_diag%SFR_OCAL%dgsave ) THEN CALL iom_put( "SFR_OCAL" , sfr_ocal2d ) CALL wrk_dealloc( jpi, jpj, sfr_ocal2d ) ENDIF IF( med_diag%SFR_OARG%dgsave ) THEN CALL iom_put( "SFR_OARG" , sfr_oarg2d ) CALL wrk_dealloc( jpi, jpj, sfr_oarg2d ) ENDIF IF( med_diag%LYSO_CA%dgsave ) THEN CALL iom_put( "LYSO_CA" , lyso_ca2d ) CALL wrk_dealloc( jpi, jpj, lyso_ca2d ) ENDIF # endif !! !! ** 3D diagnostics IF( med_diag%TPP3%dgsave ) THEN CALL iom_put( "TPP3" , tpp3d ) CALL wrk_dealloc( jpi, jpj, jpk, tpp3d ) ENDIF IF( med_diag%DETFLUX3%dgsave ) THEN CALL iom_put( "DETFLUX3" , detflux3d ) CALL wrk_dealloc( jpi, jpj, jpk, detflux3d ) ENDIF IF( med_diag%REMIN3N%dgsave ) THEN CALL iom_put( "REMIN3N" , remin3dn ) CALL wrk_dealloc( jpi, jpj, jpk, remin3dn ) ENDIF # if defined key_roam IF( med_diag%PH3%dgsave ) THEN CALL iom_put( "PH3" , f3_pH ) ENDIF IF( med_diag%OM_CAL3%dgsave ) THEN CALL iom_put( "OM_CAL3" , f3_omcal ) ENDIF !! !! AXY (09/11/16): 2D CMIP6 diagnostics IF( med_diag%INTDISSIC%dgsave ) THEN CALL iom_put( "INTDISSIC" , intdissic ) CALL wrk_dealloc( jpi, jpj, intdissic ) ENDIF IF( med_diag%INTDISSIN%dgsave ) THEN CALL iom_put( "INTDISSIN" , intdissin ) CALL wrk_dealloc( jpi, jpj, intdissin ) ENDIF IF( med_diag%INTDISSISI%dgsave ) THEN CALL iom_put( "INTDISSISI" , intdissisi ) CALL wrk_dealloc( jpi, jpj, intdissisi ) ENDIF IF( med_diag%INTTALK%dgsave ) THEN CALL iom_put( "INTTALK" , inttalk ) CALL wrk_dealloc( jpi, jpj, inttalk ) ENDIF IF( med_diag%O2min%dgsave ) THEN CALL iom_put( "O2min" , o2min ) CALL wrk_dealloc( jpi, jpj, o2min ) ENDIF IF( med_diag%ZO2min%dgsave ) THEN CALL iom_put( "ZO2min" , zo2min ) CALL wrk_dealloc( jpi, jpj, zo2min ) ENDIF IF( med_diag%FBDDTALK%dgsave ) THEN CALL iom_put( "FBDDTALK" , fbddtalk ) CALL wrk_dealloc( jpi, jpj, fbddtalk ) ENDIF IF( med_diag%FBDDTDIC%dgsave ) THEN CALL iom_put( "FBDDTDIC" , fbddtdic ) CALL wrk_dealloc( jpi, jpj, fbddtdic ) ENDIF IF( med_diag%FBDDTDIFE%dgsave ) THEN CALL iom_put( "FBDDTDIFE" , fbddtdife ) CALL wrk_dealloc( jpi, jpj, fbddtdife ) ENDIF IF( med_diag%FBDDTDIN%dgsave ) THEN CALL iom_put( "FBDDTDIN" , fbddtdin ) CALL wrk_dealloc( jpi, jpj, fbddtdin ) ENDIF IF( med_diag%FBDDTDISI%dgsave ) THEN CALL iom_put( "FBDDTDISI" , fbddtdisi ) CALL wrk_dealloc( jpi, jpj, fbddtdisi ) ENDIF !! !! AXY (09/11/16): 3D CMIP6 diagnostics IF( med_diag%TPPD3%dgsave ) THEN CALL iom_put( "TPPD3" , tppd3 ) CALL wrk_dealloc( jpi, jpj, jpk, tppd3 ) ENDIF IF( med_diag%BDDTALK3%dgsave ) THEN CALL iom_put( "BDDTALK3" , bddtalk3 ) CALL wrk_dealloc( jpi, jpj, jpk, bddtalk3 ) ENDIF IF( med_diag%BDDTDIC3%dgsave ) THEN CALL iom_put( "BDDTDIC3" , bddtdic3 ) CALL wrk_dealloc( jpi, jpj, jpk, bddtdic3 ) ENDIF IF( med_diag%BDDTDIFE3%dgsave ) THEN CALL iom_put( "BDDTDIFE3" , bddtdife3 ) CALL wrk_dealloc( jpi, jpj, jpk, bddtdife3 ) ENDIF IF( med_diag%BDDTDIN3%dgsave ) THEN CALL iom_put( "BDDTDIN3" , bddtdin3 ) CALL wrk_dealloc( jpi, jpj, jpk, bddtdin3 ) ENDIF IF( med_diag%BDDTDISI3%dgsave ) THEN CALL iom_put( "BDDTDISI3" , bddtdisi3 ) CALL wrk_dealloc( jpi, jpj, jpk, bddtdisi3 ) ENDIF IF( med_diag%FD_NIT3%dgsave ) THEN CALL iom_put( "FD_NIT3" , fd_nit3 ) CALL wrk_dealloc( jpi, jpj, jpk, fd_nit3 ) ENDIF IF( med_diag%FD_SIL3%dgsave ) THEN CALL iom_put( "FD_SIL3" , fd_sil3 ) CALL wrk_dealloc( jpi, jpj, jpk, fd_sil3 ) ENDIF IF( med_diag%FD_CAL3%dgsave ) THEN CALL iom_put( "FD_CAL3" , fd_cal3 ) CALL wrk_dealloc( jpi, jpj, jpk, fd_cal3 ) ENDIF IF( med_diag%FD_CAR3%dgsave ) THEN CALL iom_put( "FD_CAR3" , fd_car3 ) CALL wrk_dealloc( jpi, jpj, jpk, fd_car3 ) ENDIF IF( med_diag%CO33%dgsave ) THEN CALL iom_put( "CO33" , f3_co3 ) ENDIF IF( med_diag%CO3SATARAG3%dgsave ) THEN CALL iom_put( "CO3SATARAG3" , f3_omarg ) ENDIF IF( med_diag%CO3SATCALC3%dgsave ) THEN CALL iom_put( "CO3SATCALC3" , f3_omcal ) ENDIF IF( med_diag%EXPC3%dgsave ) THEN CALL iom_put( "EXPC3" , expc3 ) CALL wrk_dealloc( jpi, jpj, jpk, expc3 ) ENDIF IF( med_diag%EXPN3%dgsave ) THEN CALL iom_put( "EXPN3" , expn3 ) CALL wrk_dealloc( jpi, jpj, jpk, expn3 ) ENDIF IF( med_diag%DCALC3%dgsave ) THEN CALL iom_put( "DCALC3" , dcalc3 ) CALL wrk_dealloc( jpi, jpj, jpk, dcalc3 ) ENDIF IF( med_diag%FEDISS3%dgsave ) THEN CALL iom_put( "FEDISS3" , fediss3 ) CALL wrk_dealloc( jpi, jpj, jpk, fediss3 ) ENDIF IF( med_diag%FESCAV3%dgsave ) THEN CALL iom_put( "FESCAV3" , fescav3 ) CALL wrk_dealloc( jpi, jpj, jpk, fescav3 ) ENDIF IF( med_diag%MIGRAZP3%dgsave ) THEN CALL iom_put( "MIGRAZP3" , migrazp3 ) CALL wrk_dealloc( jpi, jpj, jpk, migrazp3 ) ENDIF IF( med_diag%MIGRAZD3%dgsave ) THEN CALL iom_put( "MIGRAZD3" , migrazd3 ) CALL wrk_dealloc( jpi, jpj, jpk, migrazd3 ) ENDIF IF( med_diag%MEGRAZP3%dgsave ) THEN CALL iom_put( "MEGRAZP3" , megrazp3 ) CALL wrk_dealloc( jpi, jpj, jpk, megrazp3 ) ENDIF IF( med_diag%MEGRAZD3%dgsave ) THEN CALL iom_put( "MEGRAZD3" , megrazd3 ) CALL wrk_dealloc( jpi, jpj, jpk, megrazd3 ) ENDIF IF( med_diag%MEGRAZZ3%dgsave ) THEN CALL iom_put( "MEGRAZZ3" , megrazz3 ) CALL wrk_dealloc( jpi, jpj, jpk, megrazz3 ) ENDIF IF( med_diag%O2SAT3%dgsave ) THEN CALL iom_put( "O2SAT3" , o2sat3 ) CALL wrk_dealloc( jpi, jpj, jpk, o2sat3 ) ENDIF IF( med_diag%PBSI3%dgsave ) THEN CALL iom_put( "PBSI3" , pbsi3 ) CALL wrk_dealloc( jpi, jpj, jpk, pbsi3 ) ENDIF IF( med_diag%PCAL3%dgsave ) THEN CALL iom_put( "PCAL3" , pcal3 ) CALL wrk_dealloc( jpi, jpj, jpk, pcal3 ) ENDIF IF( med_diag%REMOC3%dgsave ) THEN CALL iom_put( "REMOC3" , remoc3 ) CALL wrk_dealloc( jpi, jpj, jpk, remoc3 ) ENDIF IF( med_diag%PNLIMJ3%dgsave ) THEN CALL iom_put( "PNLIMJ3" , pnlimj3 ) CALL wrk_dealloc( jpi, jpj, jpk, pnlimj3 ) ENDIF IF( med_diag%PNLIMN3%dgsave ) THEN CALL iom_put( "PNLIMN3" , pnlimn3 ) CALL wrk_dealloc( jpi, jpj, jpk, pnlimn3 ) ENDIF IF( med_diag%PNLIMFE3%dgsave ) THEN CALL iom_put( "PNLIMFE3" , pnlimfe3 ) CALL wrk_dealloc( jpi, jpj, jpk, pnlimfe3 ) ENDIF IF( med_diag%PDLIMJ3%dgsave ) THEN CALL iom_put( "PDLIMJ3" , pdlimj3 ) CALL wrk_dealloc( jpi, jpj, jpk, pdlimj3 ) ENDIF IF( med_diag%PDLIMN3%dgsave ) THEN CALL iom_put( "PDLIMN3" , pdlimn3 ) CALL wrk_dealloc( jpi, jpj, jpk, pdlimn3 ) ENDIF IF( med_diag%PDLIMFE3%dgsave ) THEN CALL iom_put( "PDLIMFE3" , pdlimfe3 ) CALL wrk_dealloc( jpi, jpj, jpk, pdlimfe3 ) ENDIF IF( med_diag%PDLIMSI3%dgsave ) THEN CALL iom_put( "PDLIMSI3" , pdlimsi3 ) CALL wrk_dealloc( jpi, jpj, jpk, pdlimsi3 ) ENDIF # endif CALL wrk_dealloc( jpi, jpj, zw2d ) ENDIF ! end of ln_diatrc option # if defined key_trc_diabio !! Lateral boundary conditions on trcbio DO jn=1,jp_medusa_trd CALL lbc_lnk(trbio(:,:,1,jn),'T',1. ) ENDDO # endif # if defined key_debug_medusa IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt CALL flush(numout) # endif END SUBROUTINE trc_bio_medusa #else !!====================================================================== !! Dummy module : No MEDUSA bio-model !!====================================================================== CONTAINS SUBROUTINE trc_bio_medusa( kt ) ! Empty routine INTEGER, INTENT( in ) :: kt WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt END SUBROUTINE trc_bio_medusa #endif !!====================================================================== END MODULE trcbio_medusa