subroutine tracer (joff, js, je, is, ie) !#if defined 1 !======================================================================= ! compute tracers at "tau+1" for rows js through je in the MW. ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW !======================================================================= implicit none character(120) :: fname, new_file_name integer istrt, iend, i, k, j, ip, kr, jq, n, jp, jrow, iou, js integer je, limit, joff, is, ie, kmx, m, kb, idiag, index integer it(10), iu(10), ib(10), ic(10), nfnpzd, mfnpzd, mxfnpzd integer id_time, id_xt, id_yt, id_zt, fe_jlo, fe_m, fe_k, fe_n parameter (fe_n = 14) logical inqvardef, exists real rctheta, declin, gl, impo, expo, npp, time real remi, excr, graz, morp, morpt, morz, temp, swr, dayfrac real graz_Det, graz_Z, avej, avej_D, gmax, no3P, po4P, po4_D real phin, dz, prca, dprca, nud, bct, tap, fo2, so2, ai, hi, hs real npp_D, graz_D, morp_D, no3flag, deni, nfix, felimit real t_i, t_j, dz_t2r, dz_tr, dz_wtr, dx_t2r, dx_tr, dy_t2r real dy_tr, adv_tx, adv_ty, adv_tz, adv_txiso, adv_tyiso real adv_tziso, diff_tx, diff_ty, diff_tz, zmax, cont, drho real drhom1, wt, ahbi_cstr, ahbi_csu_dyur, gamma, rrstd, fy, fyz real fe_dy, fe_conc, fe_x(fe_n), fe_y(fe_n), bctz, felimit_D real Paulmier_a, Paulmier_z, Paulmier_R0 real expofe, impofe, feorgads, remife, thetamax, deffe, fecol real thetachl, chl, chl_D, feprime, fesed, bfe, sgb real ofehydro !wyao 20170427 !{ wyao 20170713 ! wyao 20170713} !======================= include file "size.h" ========================= !----------------------------------------------------------------------- ! USER INPUT: !----------------------------------------------------------------------- ! imt = number of grid points in the longitudinal direction ! (calculated points are from 2 through imt-1. End points ! are boundaries) ! jmt = number of grid points (latitude rows) in the latitudinal ! direction (calculated points are from 2 through jmt-1. ! End points are boundaries) ! km = number of grid points in the vertical direction ! (calculated points are from 1 through km) ! nt = number of tracers (temperature, salinity, ...) ! nsrc = number of tracer with sources ! kpzd = depth for limited npzd model ! jmz = size for "unrotated" zonal averages ! jmzm1 = jmz minus one ! mnisle = maximum number of islands (unconnected land masses) ! maxipp = maximum number of all island perimeter points !----------------------------------------------------------------------- integer imt, jmt, km, nt, nsrc, kpzd, nat, jmz, jmzm1, mnisle integer maxipp, jmw, jsmw, jemw parameter (imt= 1, jmt= 1, km= 19) parameter (nt=2 $ +1 $ +1 $ +1 $ +4 $ +2 $ +2 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +2 $ +2 $ ) parameter (kpzd=km) parameter (nat=2 $, jmz=jmt, jmzm1=jmz-1) parameter (mnisle=50, maxipp=5000) parameter (jmw=jmt) !----------------------------------------------------------------------- ! set first and last calculated row within the MW. other rows ! are used as buffers !----------------------------------------------------------------------- ! jsmw = 1st calculated row within the MW ! jemw = last calculated row within the MW parameter (jsmw=1, jemw=1) ! Moses-Triffid land model ! POINTS = Maximum number of points in grid. ! STEPSM = Maximum number of timesteps in a day. ! klmax = maximum ocean depth levels over which the land model can exist integer POINTS, STEPSM, klmax parameter (POINTS=14300, STEPSM=24, klmax=0) ! NNVG = Number of non-vegetation surface types. ! NPFT = Number of plant functional types. ! NTYPE = Number of surface types. ! SOIL = Index of the surface type 'Soil' ! Land surface types : ! 1 - Broadleaf Tree ! 2 - Needleleaf Tree ! 3 - C3 Grass ! 4 - C4 Grass ! 5 - Shrub ! 6 - Soil integer NNVG, NPFT, NTYPE, SOIL parameter (NNVG=4, NPFT=5, NTYPE=6, SOIL=6) !======================= include file "param.h" ======================== ! main parameter file which sets ocean characteristics: ! nvar = number of prognostic variables ! lseg = maximum number of longitudinal stream function segments ! nlatpr = maximum number of latitudes for matrix printouts ! on diagnostic time steps ! nhreg = number of regions in the horizontal used for averaging ! tracers. ! nvreg = number of regions in the vertical used for term balance ! calculations. note "nvreg" is not used for tracer ! averages ! numreg = total number of regions ( = product of nhreg & nvreg) ! used for term balance calculations ! nvarbh = number of prognostic variables using biharmonic mixing ! ncrows = number of calculated rows within the MW. ! (the remaining rows are buffer rows). ! rstd = standard c14/c12 ratio integer lseg, nlatpr, nhreg, nvreg, numreg, nvar, nvarbh integer imtm1, kmm1, imtp1, imtm2, jmtp1, jmtm1, jmtm2, jscan integer kmp1, kmp2, imtkm, nwds, nkflds, nslab, ntmin2, ncrows real rstd parameter (lseg=5, nlatpr=10) parameter (nhreg=3, nvreg=1, numreg=nhreg*nvreg) parameter (nvar=nt+2) parameter (imtm1=imt-1, kmm1=km-1) parameter (imtp1=imt+1, imtm2=imt-2 &, jmtp1=jmt+1, jmtm1=jmt-1, jmtm2=jmt-2 &, jscan=jmtm2 &, kmp1=km+1, kmp2=km+2 &, imtkm=imt*km, nwds=imt*jmt, nkflds=2 &, nslab=imt*nvar*km, ntmin2=nt+1/nt) parameter (ncrows = jmw - 3 + jmw/jmt) !====================== include file "pconst.h" ======================== ! rules for parameter constants ! use prefix of "c" for whole real numbers (eg: c57 for 57.0) ! use "m" after prefix to designate negative values (minus sign) ! (eg: cm7 for -7.0) ! use prefix of "p" for non repeating fractions (eg: p5 for 0.5) ! use prefix of "r" for reciprocals (eg: r3 for 1/3.0) ! combine use of prefix above and "e" for scientific notation, with ! (eg: c5e4 for 5.0e4, c1em10 for 1.0e-10) real c0, c1, c2, c3, c4, c5, c7, c8, c14, c16, c360, p125, p25 real p5, p75, epsln, c24, c60, c1440, r24, r60, r1440, secday parameter (c0=0.0, c1=1.0, c2=2.0, c3=3.0, c4=4.0, c5=5.0, c7=7.0) parameter (c8=8.0) parameter (c14=14.0, c16=16.0, c360=360.0) parameter (p125=0.125, p25=0.25, p5=0.5, p75=0.75) parameter (epsln=1.0e-20) parameter (c24=24.0, c60=60.0, c1440=1440.0) parameter (r24=c1/c24, r60=c1/c60, r1440=c1/c1440) parameter (secday=c1/(c60*c1440)) !====================== include file "stdunits.h" ====================== ! stdin = unit number for standard input. ! stdout = unit number for standard output. ! stderr = unit number for standard error. integer stdin, stdout, stderr parameter (stdin = 5, stdout = 6, stderr = 6) !======================= include file "accel.h" ======================== ! depth dependent tracer timestep acceleration multipliers used to ! hasten the convergence to equilibrium of the deeper portions of ! ocean-climate models. ! accelerate abyssal processes by varying the length of the tracer ! timestep with depth. by using longer timesteps at depth, one can ! in effect reduce the heat capacity of the deeper levels and speed ! convergence to equilibrium. ! note: ! by applying this method, one is assuming that there is a single ! steady-state solution to the model being considered. ! also, since the diagnostic timestep calculations of "termbt" do ! not attempt to account for depth variant timestep lengths, the ! truncation error reported will increase, because it will include ! the tracer changes due to variations in "dtxcel". ! reference: ! Bryan, K., 1984: accelerating the convergence to equilibrium ! of ocean climate models, J. Phys. Oceanogr., 14, 666-673. ! ("dtxcel" here is the same as 1/gamma in the above reference) ! set "dtxcel" to 1.0 at the surface and for upper levels not ! to be accelerated ! set "dtxcel" to values greater than 1.0 at deeper levels to ! accelerate convergence if above requirements are met ! dtxcel = model level dependent tracer timestep multipliers ! dtxsqr = square root of "dtxcel" (used in computation of ! maximum slope constraint for isopycnal mixing) ! dztxcl = layer thickness divided by the timestep multiplier ! (needed for convection code) ! dzwxcl = multiplication factor relating to the vertical ! distance between ts points, scaled according ! to timestep multipliers for use in convection code real dtxcel, dtxsqr, dztxcl, dzwxcl common /accel/ dtxcel(km) common /accel/ dtxsqr(km) common /accel/ dztxcl(km), dzwxcl(km) !======================= include file "coord.h" ======================== ! model grid point coordinates ! grid definition: ! the model uses a staggered arakawa "b" grid which is setup and ! generated by the "grids.F" module. ! xt(i) = longitude of the ith "t" point in degrees. i=1..imt ! xu(i) = longitude of the ith "u,v" point in degrees. i=1..imt ! yt(j) = latitude of the jth "t" point in degrees. j=1..jmt ! yu(j) = latitude of the jth "u,v" point in degrees. j=1..jmt ! zt(k) = distance from surface down to centre of level k (in cm) ! (for depth of "t" and "u,v" grid points: k=1,km) ! zw(k) = distance from surface down to bottom of level k (in cm) ! (for depth of "t" and "u,v" grid points: k=1,km) ! dxtdeg = widths for "t" grid cells (degrees) ! dytdeg = heights for "t" grid cells (degrees) ! dxudeg = widths for "u" grid cells (degrees) ! dyudeg = heights for "u" grid cells (degrees) ! dzt(k) = vertical resolution of "t" and "u" grid cells (in cm) ! dzw(k) = vertical resolution of "w" grid cells (in cm) ! "i" increases in an eastward direction, "j" increases in a ! northward direction, and "k" increases downward. real xt, yt, xu, yu, zw, zt real dxtdeg, dytdeg, dzt real dxudeg, dyudeg, dzw common /coord/ xt(imt), yt(jmt), xu(imt), yu(jmt), zw(km), zt(km) common /coord/ dxtdeg(imt), dytdeg(jmt), dzt(km) common /coord/ dxudeg(imt), dyudeg(jmt), dzw(0:km) !====================== include file "cregin.h" ======================== ! variables used for computing regional tracer averages (see ! "reg1st.F" & sub "diagt1" in "tracer.F") and for computing term ! balances for tracer and momentum equations (see "clinic.F", ! "tracer.F" and "diag2.F") ! mskhr = mask field defining regions in the horizontal ! (eg: mskhr(i,j) = n indicates point (i,j) is in the ! "nth" horizontal region where n=1..nhreg) ! The "mskhr" masks are used in "diagt1" when ! computing volume weighted tracer averages and in ! "clinic.F", "tracer.F" and "diag2.F" when computing ! term balances for tracers and momentum. ! mskvr = mask field defining regions in the vertical ! (eg: mskvr(k) = m indicates all points under a horizontal ! mask at level "k" are in the "mth" vertical region ! where m=1..nvreg) ! The "mskvr" masks are used in "diag.F", but not ! in "diagt1", where tracer averages are calculated ! for each k-level. ! hregnm = horizontal region name ! vregnm = vertical region name ! volbt = total volume under a given horizontal region ! volbk = volume contained in a horizontal region at level "k" ! volgt = total ocean volume ! volgk = total ocean volume at level "k" ! areab = total ocean surface area for a given horizontal region ! areag = total ocean surface area ! volt = ocean volume within a particular horizontal & vertical ! region (on the "t" grid) for tracer term balances ! volt(0) represents the sum of all regions ! rvolt = 1/volt ( 0.0 if volt = 0.0) ! areat = horizontal ocean surface area corresponding to "volt" ! areat(0) represents the sum of all regions ! rareat = 1/areat ( 0.0 if areat = 0.0) ! volu = ocean volume within a particular horizontal & vertical ! region (on the "u" grid) for momentum term balances ! volu(0) represents the sum of all regions ! rvolu = 1/volu ( 0.0 if volu = 0.0) ! areau = horizontal ocean area corresponding to "volu" ! areau(0) represents the sum of all regions ! rareau = 1/areau ( 0.0 if areau = 0.0) ! llvreg = level limits for defining vertical regions in term ! balance calculations (not used in computing volume ! weighted tracer averages) ! (eg: llvreg(3,1) = 4... means that starting level for ! the third region in the vertical is 4. similarly, ! llvreg(3,2) = 6 means the ending level is 6 for that ! region. note regions should not overlap.) character(40) :: hregnm character(20) :: vregnm common /cregn_c/ hregnm(nhreg), vregnm(nvreg) integer mskhr, mskvr, llvreg common /cregn_i/ mskhr(imt,jmt), mskvr(km), llvreg(numreg,2) real volbk, volbt, volgk, volgt, areab, areag, volt, volu real areat, areau, rvolt, rvolu, rareat, rareau common /cregn_r/ volbk(nhreg,km), volbt(nhreg), volgk(km) common /cregn_r/ volgt, areab(nhreg), areag, volt(0:numreg) common /cregn_r/ volu(0:numreg), areat(0:numreg), areau(0:numreg) common /cregn_r/ rvolt(0:numreg), rvolu(0:numreg) common /cregn_r/ rareat(0:numreg), rareau(0:numreg) !======================= include file "csbc.h" ========================= ! surface boundary conditions (S.B.C.) ! numsbc = total number of surface boundary conditions ! indices and names set in UVic_ESCM.F: ! itaux = x component of wind stress (dynes cm-2) ! itauy = y component of wind stress (dynes cm-2) ! iws = surface wind speed (cm s-1) ! iaca = atmospheric coalbedo ! isca = surface coalbedo (%/100) ! ihflx = heat flux ! isflx = salt flux ! isst = ocean model SST (C) ! isss = ocean model SSS (PSU*0.001-0.035) ! iwa = surface wind angle (degrees) ! iro = surface runoff (g cm-2 s-1) ! iwxq = x component of wind for moisture advection (cm s-1) ! iwyq = y component of wind for moisture advection (cm s-1) ! iwxt = x component of wind for temperature advection (cm s-1) ! iwyt = y component of wind for temperature advection (cm s-1) ! iwxc = x component of wind for carbon advection (cm s-1) ! iwyc = y component of wind for carbon advection (cm s-1) ! ipsw = penetrating shortwave (cal cm-2 s-1) ! isu = x component of ocean first layer velocity (cm s-1) ! isv = y component of ocean first layer velocity (cm s-1) ! igu = x component of ocean second layer velocity (cm s-1) ! igv = y component of ocean second layer velocity (cm s-1) ! issdic = sea surface concentration of dic (umol cm-3) ! idicflx = sea surface flux of carbon (umol cm-2 s-1) ! issalk = sea surface concentration of alkalinity (umol cm-3) ! ialkflx = sea surface flux of alkalinity (umol cm-2 s-1) ! isso2 = sea surface concentration of oxygen (umol cm-3) ! io2flx = sea surface flux of oxygen (umol cm-2 s-1) ! isspo4 = sea surface concentration of phosphate (nmol cm-3) ! ipo4flx = sea surface flux of phosphate (nmol cm-2 s-1) !{ wyao 20170704 ! wyao 20170704} ! issphyt = sea surface concentration of phytoplankton (nmol P cm-3) ! iphytflx = sea surface flux of phytoplankton (nmol P cm-2 s-1) ! isszoop = sea surface concentration of zooplankton (nmol P cm-3) ! izoopflx = sea surface flux of zooplankton (nmol P cm-2 s-1) ! issdetr = sea surface concentration of detritus (nmol P cm-3) ! idetrflx = sea surface flux of detritus (nmol P cm-2 s-1) ! issdetrfe = sea surface concentration of particualte Fe (nmol Fe cm-3) ! idetrfeflx= sea surface flux of particulate Fe (nmol Fe cm-2 s-1) ! issno3 = sea surface concentration of nitrate (nmol P cm-3) ! ino3flx = sea surface flux of nitrate (nmol P cm-2 s-1) !{ wyao 20170704 ! wyao 20170704} ! issdfe = sea surface concentration of iron (nmol Fe cm-3) ! idfeflx = sea surface flux of iron (nmol Fe cm-2 s-1) ! idfeadep = sea surface flux of iron (nmol Fe cm-2 s-1) ! issdiaz = sea surface concentration of diazotraphs (nmol P cm-3) ! idiazflx = sea surface flux of diazotraphs (nmol P cm-2 s-1) ! issc14 = sea surface concentration of carbon 14 (umol cm-3) ! ic14flx = sea surface flux of carbon 14 (umol cm-2 s-1) ! isscfc11 = sea surface concentration of cfc11 (umol cm-3) ! icfc11flx = sea surface flux of cfc11 (umol cm-2 s-1) ! isscfc12 = sea surface concentration of cfc12 (umol cm-3) ! icfc12flx = sea surface flux of cfc12 (umol cm-2 s-1) ! iat = surface air temperature (C) ! irh = surface relative humidity (%/100) ! ipr = precipitation as rain (g cm-2 s-1) ! ips = precipitation as snow (g cm-2 s-1) ! iaws = averaged surface wind speed (cm s-1) ! iswr = surface shortwave radiation (g s-3) ! ilwr = surface longwave radiation (g s-3) ! isens = surface sensible heat flux (g s-3) ! ievap = surface evaporation (g cm-2 s-1) ! idtr = diurnal temperature range (C) ! inpp = carbon flux from net primary production (kg C m-2 s-1) ! isr = carbon flux from soil respiration (kg m-2 s-1) ! iburn = carbon flux from burning vegetation (kg m-2 s-1) ! ibtemp = ocean bottom temperature (C) ! ibsalt = ocean bottom temperature (PSU*0.001-0.035) ! ibdic = ocean bottom dissolved organic carbon (umol cm-3) ! ibdicfx = ocean bottom flux of carbon (umol cm-2 s-1) ! ibalk = ocean bottom alkalinity (umol cm-3) ! ibalkfx = ocean bottom flux of alkalinity (umol cm-2 s-1) ! ibo2 = ocean bottom oxygen (umol cm-3) ! ircal = rain rate of calcite (umol cm-2 s-1) ! irorg = rain rate of organic carbon (umol cm-2 s-1) ! mapsbc = surface boundary conditions names ! sbc = surface boundary condition data. ! gaost = global average ocean surface tracer ! socn = average ocean sea surface salinity ! ntspas = the number of ocean time steps per atmosphere segment ! ntspls = the number of ocean time steps per land segment ! ntspos = the number of ocean time steps per ocean segment ! nats = number of atm time steps since mixing time step ! nots = number of ocn time steps since mixing time step ! addflxa = logical flag for adding only even mode fluxes from atm ! addflxo = logical flag for adding only even mode fluxes from ocn ! dtatm = atmosphere time step (s) ! dtism = icesheet time step (s) ! dtismyr = icesheet time step (years) ! nismacc = icesheet acceleration time (years) ! dtlnd = land time step (s) ! dtocn = ocean time step (s) ! dtsed = sediment time step (s) ! dtsedyr = sediment time step (years) ! nsedacc = sediment acceleration time (years) ! atatm = atmosphere boundary condition averaging time (s) ! atlnd = land boundary condition averaging time (s) ! atocn = ocean boundary condition averaging time (s) ! atsed = sediment boundary condition averaging time (s) ! atism = icesheet boundary condition averaging time (s) ! dampts = time scale for damping surface tracers (days) to data ! dampdz = ocean level thickness for converting Newtonian damping ! to a surface flux ! land_map = map with indices for coupling to land arrays ! dtoih = total system heat lost minus heat gained ! dtoic = total system carbon lost minus carbon gained ! avgpertavg = averaging period for time averages ! avgtimtavg = averaging time for time averages ! avgpertsi = averaging period for time integrals ! avgtimtsi = averaging time for time integrals integer numsbc parameter (numsbc = 14 & + 4 & + 2 & + 2 & + 2 & + 2 & + 2 & + 3 & ) integer itaux, itauy, iws, iaca, isca, ihflx, isflx, isst, isss integer iwa, iro, iwxq, iwyq, iwxt, iwyt, iwxc, iwyc, ipsw, isu integer isv, igu, igv, issdic, idicflx, issalk, ialkflx, isso2 integer io2flx, isspo4, ipo4flx, issphyt, iphytflx, isszoop integer izoopflx, issdetr, idetrflx, issno3, ino3flx, issdiaz integer idiazflx, issc14, ic14flx, isscfc11, icfc11flx, isscfc12 integer icfc12flx, iat, irh, ipr, ips, iaws, iswr, ilwr, isens integer ievap, idtr, inpp, isr, iburn, ibtemp, ibsalt, ibdic integer ibdicfx, ibalk, ibalkfx, ibo2, ircal, irorg !{ wyao 20170704 integer issdfe ! wyao 20170704} integer idfeflx, issdetrfe, idetrfeflx, idfeadep common /csbc_i/ itaux, itauy, iws, iaca, isca, ihflx, isflx, isst common /csbc_i/ isss, iwa, iro, iwxq, iwyq, iwxt, iwyt, iwxc, iwyc common /csbc_i/ ipsw, isu, isv, igu, igv, issdic, idicflx, issalk common /csbc_i/ ialkflx, isso2, io2flx, isspo4, ipo4flx, issphyt common /csbc_i/ iphytflx, isszoop, izoopflx, issdetr, idetrflx common /csbc_i/ issno3, ino3flx, issdiaz, idiazflx, issc14 common /csbc_i/ ic14flx, isscfc11, icfc11flx, isscfc12, icfc12flx common /csbc_i/ iat, irh, ipr, ips, iaws, iswr, ilwr, isens, ievap common /csbc_i/ idtr, inpp, isr, iburn, ibtemp, ibsalt, ibdic common /csbc_i/ ibdicfx, ibalk, ibalkfx, ibo2, ircal, irorg !{ wyao 20170704 ! wyao 20170704} common /csbc_i/ issdfe, idfeflx, issdetrfe, idetrfeflx common /csbc_i/ idfeadep character(20) :: mapsbc common /csbc_c/ mapsbc(numsbc) real sbc common /csbc_r/ sbc(imt,jmt,numsbc) real gaost, socn common /csbc_r/ gaost(nt), socn integer ntspas, ntspls, ntspos, nats, nots, nismacc, nsedacc common /csbc_i/ ntspas, ntspls, ntspos, nats, nots, nismacc common /csbc_i/ nsedacc logical addflxa, addflxo common /csbc_l/ addflxa, addflxo real dtatm, dtism, dtismyr, dtlnd, dtocn, dtsed, dtsedyr, atatm real atism, atlnd, atocn, atsed common /csbc_r/ dtatm, dtism, dtismyr, dtlnd, dtocn, dtsed common /csbc_r/ dtsedyr, atatm, atism, atlnd, atocn, atsed real dampts, dampdz, subflux, subz common /csbc_r/ dampts(nt), dampdz(nt) integer ntlbc common /csbc_i/ ntlbc real dtoih, dtoic common /csbc_r/ dtoih, dtoic real avgpertavg, avgtimtavg, avgpertsi, avgtimtsi common /csbc_r/ avgpertavg, avgtimtavg, avgpertsi, avgtimtsi !======================= include file "emode.h" ======================== ! variables for external mode ! psi = stream function (,,1) is for tau; (,,2) is for tau-1 ! zu = vertically averaged forcing from momentum equations ! (,,1) is zonal and (,,2) is meridional component ! ztd = curl of "zu" for the stream function equation ! ptd = time change of stream function ! h = depth over "u" points ! hr = reciprocal depth over "u" points ! res = residual from elliptic solver ! map = land mass map distinguishing, ocean, land, and perimeters ! mxscan = max number of allowable scans for Poisson solvers ! mscan = actual number of scans taken by Poisson solvers ! sor = successive over-relaxation constant ! tolrsf = tolerance for stream function calculation. ! the solution is halted when it is within "tolrsf" ! of the "true" solution assuming geometric convergence. ! tolrsp = tolerance for surface pressure calculation ! the solution is halted when it is within "tolrsp" ! of the "true" solution assuming geometric convergence. ! tolrfs = tolerance for implicit free surface calculation ! the solution is halted when it is within "tolrfs" ! of the "true" solution assuming geometric convergence. ! esterr = estimated maximum error in elliptic solver assuming ! geometric convergence ! nisle = number of land masses ! nippts= number of land mass perimeter points ! iperm = "i" coordinate for the land mass perimeter point ! jperm = "j" coordinate for the land mass perimeter point ! iofs = offset for indexing into the land mass perimeter points ! imask = controls whether calculations get done on perimeters ! set mask for land mass perimeters on which to perform calculations ! imask(-n) = .false. [no equations ever on dry land mass n] ! imask(0) = .true. [equations at all mid ocean points] ! imask(n) = .true./.false [controls whether there will be ! equations on the ocean perimeter of ! land mass n] ! note: land mass 1 is the northwest-most land mass ! for the numbering of the other landmasses, see generated map(i,j) character(16) :: variable common /emode_c/ variable logical converged, imask common /emode_l/ converged, imask (-mnisle:mnisle) integer mxscan, mscan, map, nippts, iofs, iperm integer jperm, nisle, imain common /emode_i/ mxscan, mscan common /emode_i/ map(imt,jmt) common /emode_i/ nippts(mnisle), iofs(mnisle), iperm(maxipp) common /emode_i/ jperm(maxipp), nisle, imain real tolrsf, tolrsp, tolrfs, sor, esterr, ptd, res, hr real h, zu, psi, ztd, alph, gam, theta, apgr, uhat, pguess real ps, divf, ubar, ubarm1 common /emode_r/ tolrsf, tolrsp, tolrfs, sor, esterr common /emode_r/ ptd(imt,jmt), res(imt,jmt), hr(imt,jmt) common /emode_r/ h(imt,jmt), zu(imt,jmt,2) common /emode_r/ psi(imt,jmt,2), ztd(imt,jmt) !====================== include file "grdvar.h" ======================== ! variables which are functions of the grid defined by "coord.h" ! dxt = longitudinal width of "t" grid box at the ! equator (in cm) ! dxtr = reciprocal of "dxt" ! dxt2r = reciprocal of "2*dxt" ! dxt4r = reciprocal of "4*dxt" ! dxu = longitudinal width of "u,v" grid box at the ! equator (in cm) ! dxur = reciprocal of "dxu" ! dxu2r = reciprocal of "2*dxu" ! dxu4r = reciprocal of "4*dxu" ! dxmetr = reciprocal of "(dxt(i)+dxt(i+1))" ! duw = xu(i) - xt(i) ! due = xt(i+1) - xu(i) ! dus = yu(jrow) - yt(jrow) ! dun = yt(jrow+1) - yu(jrow) ! dyt = latitudinal height of "t" grid box (in cm) ! dytr = reciprocal of "dyt" ! dyt2r = reciprocal of "2*dyt" ! dyt4r = reciprocal of "4*dyt" ! dyu = latitudinal height of "u,v" grid box (in cm) ! dyur = reciprocal of "dyu" ! dyu2r = reciprocal of "2*dyu" ! dyu4r = reciprocal of "4*dyu" ! csu = cosine of "u,v" grid point latitude ! csur = reciprocal of "csu" ! cst = cosine of "t" grid point latitude ! cstr = reciprocal of "cst" ! phi = latitude of "u,v" grid point in radians ! phit = latitude of "t" grid point in radians ! sine = sine of "u,v" grid point latitude ! tng = tan of "u,v" grid point latitude ! c2dzt(k)= "2*dzt" ! dztr(k) = reciprocal of dzt ("t" cell vertical resolution) ! dzt2r(k)= reciprocal of "2*dzt" ! dzwr(k) = reciprocal of dzw ("w" cell vertical resolution) ! dzw2r(k)= reciprocal of "2*dzw" ! dztur(k)= upper diffusion grid factor = 1.0/(dzw(k-1)*dzt(k)) ! dztlr(k)= lower diffusion grid factor = 1.0/(dzw(k)*dzt(k)) ! tlat = tracer grid geographic latitude (degrees) ! ulat = velocity grid geographic latitude (degrees) ! tlon = tracer grid geographic longitude (degrees) ! ulon = velocity grid geographic longitude (degrees) ! tgarea = tracer grid geographic area (cm2) ! ugarea = velocity grid geographic area (cm2) real dxt, dxtr, dxt2r, dxu, dxur, dxu2r, dxu4r, dxt4r real dyt, dytr, dyt2r, dyu, dyur, dyu2r, dyu4r, dyt4r real csu, csur, cst, cstr, cstdytr, cstdyt2r real csudyur, csudyu2r, cst_dytr, csu_dyur real phi, phit, sine, tng, c2dzt, dztr, dzt2r real dzwr, dzw2r, dxmetr, duw, due, dun, dus, dztur, dztlr real quick_x, curv_xp, curv_xn, quick_y, curv_yp, curv_yn real quick_z, curv_zp, curv_zn real tlat, tlon, ulat, ulon, tgarea, ugarea common /grdvar/ dxt(imt), dxtr(imt), dxt2r(imt), dxu(imt) common /grdvar/ dxur(imt), dxu2r(imt), dxu4r(imt), dxt4r(imt) common /grdvar/ dyt(jmt), dytr(jmt), dyt2r(jmt), dyu(jmt) common /grdvar/ dyur(jmt), dyu2r(jmt), dyu4r(jmt), dyt4r(jmt) common /grdvar/ csu(jmt), csur(jmt), cst(jmt), cstr(jmt) common /grdvar/ cstdytr(jmt), cstdyt2r(jmt) common /grdvar/ csudyur(jmt), csudyu2r(jmt) common /grdvar/ cst_dytr(jmt), csu_dyur(jmt) common /grdvar/ phi(jmt), phit(jmt), sine(jmt), tng(jmt) common /grdvar/ c2dzt(km), dztr(km), dzt2r(km) common /grdvar/ dzwr(0:km), dzw2r(0:km) common /grdvar/ dxmetr(imt), duw(imt), due(imt) common /grdvar/ dun(jmt), dus(jmt) common /grdvar/ dztur(km), dztlr(km) common /grdvar/ tlat(imt,jmt), tlon(imt,jmt), tgarea(imt,jmt) common /grdvar/ ulat(imt,jmt), ulon(imt,jmt), ugarea(imt,jmt) ! tcella = "t" cell surface area cm**2 (entire ocean) ! ucella = "u" cell surface area cm**2 (entire ocean) ! tcellv = "t" cell volume cm**3 (entire ocean) ! ucellv = "u" cell volume cm**3 (entire ocean) real tcellv, ucellv, tcella, ucella common /grdvar/ tcellv, ucellv, tcella(km), ucella(km) !======================= include file "hmixc.h" ======================== ! horizontal mixing coefficients ! visc_cnu = viscosity coeff for northern face of "u" cell ! visc_ceu = viscosity coeff for eastern face of "u" cell ! diff_cnt = diffusion coeff for northern face of "T" cell ! diff_cet = diffusion coeff for eastern face of "T" cell ! am = constant lateral viscosity coeff for momentum ! ah = constant lateral diffusion coeff for tracers ! am3 = viscosity coeff for metric term on "u" cell ! am4 = another viscosity coeff for metric term on "u" cell ! ambi = constant lateral biharmonic viscosity coeff for momentum ! ahbi = constant lateral biharmonic diffusion coeff for tracers !======================================================================= real am, ambi, am3, am4, ah, ahbi, visc_ceu, visc_cnu, amc_north real amc_south, Ahh(km), diff_cnt, diff_cet, ahc_north, ahc_south real strain, am_lambda, am_phi, smag_metric, diff_c_back real hl_depth, hl_back, hl_max, hl_u, hl_n, hl_e, hl_b real droz, rich_inv common /diffus_r/ am, ambi, am3(jmt), am4(jmt,2) common /diffus_r/ ah, ahbi common /diffus_r/ visc_ceu(imt,km,jmt) common /diffus_r/ visc_cnu(imt,km,jmt) common /diffus_r/ amc_north(imt,km,jmt) common /diffus_r/ amc_south(imt,km,jmt) common /diffus_r/ diff_cnt, diff_cet common /diffus_r/ ahc_north(jmt), ahc_south(jmt) !====================== include file "levind.h" ======================== ! vertical level indicators which define model geometry & bottom ! topography: ! kmt = number of vertical boxes over "t" points ! kmu = number of vertical boxes over "u,v" points integer kmt, kmu real sg_bathy common /levind/ kmt(imt,jmt), kmu(imt,jmt) common /levind/ sg_bathy(imt,jmt,km) !======================= include file "mw.h" =========================== ! M E M O R Y W I N D O W ! Data layout on disk ! On disk, a latitude consists of a row of prognostic quantities. ! Latitudes are ordered from south to north on two logical disk ! units: one for time level "tau" and one for time level "tau-1". ! Newly computed "tau+1" quantities are written to the "tau-1" disk ! which serves double duty as a "tau+1" disk. ! Data layout in memory ! A memory window "MW" is constructed to hold a group of "jmw" ! adjacent latitude rows of prognostic quantities from disks "tau-1" ! and "tau". Parameter "jmw" controls the size of the MW and can be ! set anywhere from a minimum of three latitude rows to all "jmt" ! latitudes. In addition, the MW holds diagnostic quantities ! (density, hydrostatic pressure gradients, and advective ! velocities) along with work arrays for constructing intermediate ! fluxes used in solving the tracer and momentum equations. A ! latitude row in the MW is referred to by index "j" and corresponds ! to the physical latitude row "jrow" on disk. ! Data flow between disk and memory ! The MW is loaded with prognostic data from the first "jww" ! latitude rows on disks "tau" and "tau-1". As the tracer and ! momentum equations are solved for rows j=2 through jmw-1 in the ! MW, the solutions are written to disk "tau+1". When finished, the ! MW is moved northward by moving quantities from the northernmost ! two rows into the southernmost two rows within the MW. The ! remaining MW rows are then loaded with more latitude rows from ! disk. The process continues until all latitude rows 2 through ! jmt-1 on disk "tau+1" have been updated. !======================================================================= ! taum1 = tau-1 time level for variables in MW ! tau = tau time level for variables in MW ! taup1 = tau+1 time level for variables in MW integer taum1, tau, taup1 common /mw_i/ taum1, tau, taup1 !----------------------------------------------------------------------- ! MW arrays for prognostic equations: !----------------------------------------------------------------------- ! u(i,k,j,n,tau) = total velocity where: ! i = index for longitude ! k = index for depth ! j = index for latitude row within MW ! n = component (1 = zonal, 2 = meridional) ! tau = time level (tau-1, tau, tau+1) ! (only internal modes are on disk and at tau+1 in the MW) ! t(i,k,j,n,tau) = tracer where: ! i = index for longitude ! k = index for depth ! j = index for latitude row within MW ! n = component (1 = temperature, 2 = salinity) ! if nt > 2 then other tracers are allowed. ! tau = time level (tau-1, tau, tau+1) ! note: temperature is potential temperature in degrees Celsius and ! salinity is in "model units", the deviation from 0.035 grams ! of salt/cm**3 of water, or, assuming a water density of ! 1 gram/cm**3, the deviation from 0.035 g of salt/g of water. ! one can convert model units to the more common units of parts ! per thousand (ppt) by adding 0.035 grams/cm**3 to the model ! units and then multiplying by 1000. real u,t common /mw_r/ u(imt,km,jmw,2,-1:1), t(imt,km,jmw,nt,-1:1) ! indicies for ocean tracer array ! itemp = index for temperature ! isalt = index for salinity ! idic = index for dissolved inorganic carbon ! ic14 = index for carbon 14 ! icfc11 = index for cfc11 ! icfc12 = index for cfc12 ! ialk = index for alkalinity ! ipo4 = index for phosphate ! iphyt = index for phytoplankton ! izoop = index for zooplankton ! idetr = index for detritus ! idfe = index for iron ! idetrfe= index for particulate iron ! io2 = index for oxygen ! ino3 = index for nitrate ! idiaz = index for diazotrophs ! mapt = map for ocean tracer names ! { wyao 20170627 ! wyao 20170627} character(10) :: mapt common /mw_c/ mapt(nt) integer itemp, isalt, idic, ic14, icfc11, icfc12, ialk, ipo4 integer iphyt, izoop, idetr, io2, ino3, idiaz common /mw_i/ itemp, isalt, idic, ic14, icfc11, icfc12, ialk common /mw_i/ ipo4, iphyt, izoop, idetr, io2, ino3, idiaz integer idfe, idetrfe common /mw_i/ idfe, idetrfe ! { wyao 20170627 ! wyao 20170627} ! indicies for ocean tracer source array ! istemp = index for temperature ! issalt = index for salinity ! isdic = index for carbon ! isc14 = index for carbon 14 ! isalk = index for alkalinity ! ispo4 = index for phosphate ! isphyt = index for phytoplankton ! iszoop = index for zooplankton ! isdetr = index for detritus ! isdfe = index for iron ! isdetrfe= index for particulate iron ! { wyao 20170627 ! wyao 20170627} ! itrc = index of tracer sources for all tracers (0 = no source) ! mapst = map for ocean tracer source names character(10) :: mapst common /mw_c/ mapst(nt) integer itrc common /mw_i/ itrc(nt) integer istemp, issalt, isdic, isc14, isalk, ispo4, isphyt integer iszoop, isdetr, iso2, isno3, isdiaz common /mw_i/ istemp, issalt, isdic, isc14,isalk, ispo4, isphyt common /mw_i/ iszoop, isdetr, iso2, isno3, isdiaz integer isdfe, isdetrfe common /mw_i/ isdfe, isdetrfe ! { wyao 20170627 ! wyao 20170627} !----------------------------------------------------------------------- ! MW arrays for diagnostic equations and workspace: !----------------------------------------------------------------------- ! diagnostic advective velocities are in units of cm/sec ! adv_vet = advective velocity on the eastern face of a "T" cell ! adv_vnt = advective velocity on the northern face of a "T" cell ! adv_veu = advective velocity on the eastern face of a "u" cell ! adv_vnu = advective velocity on the northern face of a "u" cell ! adv_vbt = advective velocity on the bottom face of a "T" cell ! adv_vbu = advective velocity on the bottom face of a "u" cell ! rho = density at centre of a "T" cell in units of gm/cm**3 ! note: there is an arbitrary constant which is only a ! function of depth in "rho". It is related to ! subtracting a reference level density for purposes of ! accuracy. ! grad_p = hydrostatic pressure gradient for "u" cell. There are ! two components: (1,2) is for (dp/dx, dp/dy) real adv_vet, adv_vnt, adv_veu, adv_vnu, adv_vbt, adv_vbu, rho real rhotaum1, rhotaup1, rhotilde, grad_p common /mw_r/ adv_vet(imt,km,jsmw:jmw), adv_vnt(imt,km,1:jmw) common /mw_r/ adv_veu(imt,km,jsmw:jemw) common /mw_r/ adv_vnu(imt,km,1:jemw) common /mw_r/ adv_vbt(imt,0:km,jsmw:jmw) common /mw_r/ adv_vbu(imt,0:km,jsmw:jemw) common /mw_r/ rho(imt,km,jsmw:jmw) common /mw_r/ grad_p(imt,km,jsmw:jemw,2) ! tmask = tracer cell land/sea mask = (0.0, 1.0) on (land, sea) ! umask = velocity cell land/sea mask = (0.0, 1.0) on (land, sea) real tmask, umask common /mw_r/ tmask(imt,km,1:jmw), umask(imt,km,1:jmw) ! these workspace arrays are recalculated for each component of the ! equations so do not have to be moved as the MW moves northward. ! adv_fe = advective flux across the eastern face of a cell ! adv_fn = advective flux across the northern face of a cell ! (removed in most cases and put directly into the ! statement functions for speed optimization.) ! adv_fb = advective flux across the bottom face of a cell ! diff_fe = diffusive flux across the eastern face of a cell ! diff_fn = diffusive flux across the northern face of a cell ! diff_fb = diffusive flux across the bottom face of a cell ! source = source term real adv_fe common /mw_r/ adv_fe(imt,km,jsmw:jemw) real adv_fn common /mw_r/ adv_fn(imt,km,1:jemw) real adv_fb, diff_fe common /mw_r/ adv_fb(imt,0:km,jsmw:jemw) common /mw_r/ diff_fe(imt,km,jsmw:jemw) real diff_fn common /mw_r/ diff_fn(imt,km,1:jemw) real diff_fb common /mw_r/ diff_fb(imt,0:km,jsmw:jemw) real diff_fbiso common /mw_r/ diff_fbiso(imt,0:km,jsmw:jemw) real source common /mw_r/ source(imt,km,jsmw:jemw) real zzi common /mw_r/ zzi(imt,km,jsmw:jemw) ! these grid factors are for optimizations and are recalculated as ! the MW moves northward so they do not have to be moved. real cstdxtr, cstdxur, cstdxt2r, ah_cstdxur, csudxur, csudxtr real csudxu2r, am_csudxtr common /mw_r/ cstdxtr(imt,jsmw:jmw), cstdxur(imt,jsmw:jmw) common /mw_r/ cstdxt2r(imt,jsmw:jmw), ah_cstdxur(imt,jsmw:jmw) common /mw_r/ csudxur(imt,jsmw:jmw) common /mw_r/ csudxtr(imt,jsmw:jmw) common /mw_r/ csudxu2r(imt,jsmw:jmw) common /mw_r/ am_csudxtr(imt,km,jsmw:jmw) ! these variables are either constant or globally dimensioned by ! "jmt", so they do not need to be moved as the MW moves northward ! advmet = coeff for metric advection. ! cori = Coriolis parameter for velocity component "n" real advmet, cori, cf common /advec_r/ advmet(jmt,2) common /coriol_r/ cori(imt,jmt,2) common /extwrk_r/ cf(imt,jmt,-1:1,-1:1) ! smf = surface momentum flux ! 1 => zonal wind stress (dynes/cm**2) ! 2 => meridional wind stress (dynes/cm**2) ! bmf = bottom momentum flux ! 1 => zonal bottom drag (dynes/cm**2) ! 2 => meridional bottom drag (dynes/cm**2) ! stf = surface tracer flux ! 1 => surface heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec) ! 2 => surface salt flux (grams of salt/cm**2/sec) ! btf = bottom tracer flux (for consistency but normally zero!) ! 1 => bottom heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec) ! 2 => bottom salt flux (grams of salt/cm**2/sec) real smf, bmf, stf, btf common /mw_r/ smf(imt,1:jmw,2), bmf(imt,1:jmw,2) common /mw_r/ stf(imt,1:jmw,nt), btf(imt,1:jmw,nt) ! "antidiffusive" flux as in Zalesak, 1989( see FCTstm.h) ! same for R+, R- ! anti_fe = antidiffusive flux across the eastern face of a T cell ! anti_fn = antidiffusive flux across the northern face of a T cell ! anti_fb = antidiffusive flux across the bottom face of a T cell ! R_plusY = ratio of maximal feasible to maximal possible change ! of tracer T in subroutine tracer.F, N-S dimension delimiter ! R_minusY = ratio of minimal feasible to minimal possible change ! of tracer T in subroutine tracer.F, N-S dimension delimiter real anti_fe, anti_fn, anti_fb, R_plusY, R_minusY common /mw_r/ anti_fe(imt,km,jsmw:jmw,nt) common /mw_r/ anti_fn(imt,km,1:jmw-1+jmw/jmt,nt) common /mw_r/ anti_fb(imt,0:km,jsmw:jmw,nt) common /mw_r/ R_plusY(imt,km,1:jmw-1+jmw/jmt,nt) common /mw_r/ R_minusY(imt,km,1:jmw-1+jmw/jmt,nt) !====================== include file "scalar.h" ======================== ! various scalar quantities: ! dtts = time step for density & tracers (in seconds) ! dtuv = time step for baroclinic velocity (in seconds) ! dtsf = time step for barotropic velocity (in seconds) ! c2dtts = 2*dtts ! c2dtuv = 2*dtuv ! c2dtsf = 2*dtsf ! acor = (>0, 0) = (implicit, explicit) treatment of Coriolis ! term for internal and external modes. ! rho0 = mean density for Bousinessq approximation ! rho0r = 1/rho0 ! omega = earth`s rotation rate (radians/sec) ! radius = earth`s radius (cm) ! grav = earth`s gravitational acceleration (cm/sec**2) ! cdbot = bottom drag coefficient ! ncon = number of passes through convective code in tracer ! gcor = time centring for Coriolis term ! taux0 = constant zonal windstress (dynes/cm**2) for idealized ! equatorial studies ! tauy0 = constant meridional windstress (dynes/cm**2) for ! idealized equatorial studies integer ncon common /scalar_i/ ncon real dtts, dtuv, dtsf, c2dtts, c2dtuv, c2dtsf, acor, rho0 real rho0r, omega, radius, grav, cdbot, gcor, taux0, tauy0 common /scalar_r/ dtts, dtuv, dtsf, c2dtts, c2dtuv, c2dtsf, acor common /scalar_r/ rho0, rho0r, omega, radius, grav, cdbot, gcor common /scalar_r/ taux0, tauy0 ! various non dimensional quantities: ! radian = degrees per radian ! pi = something good to eat real radian, pi common /ndcon/ radian, pi !====================== include file "switch.h" ======================== ! all time dependent decisions are made by time manager "tmngr.F" ! and communicated elsewhere to the model via logical switches. ! inputs: (defaulted in "blkdta.F", optionally reset via namelist) ! runlen = integration period (see rununits). note "runlen" should ! be an integral number of density time steps. if not, ! then "runlen" is automatically adjusted to insure this. ! fractional days are supported but not fractional months ! or years. ! rununits= units of "runlen". may be "days", "months", or "years". ! tmngr will convert "runlen" which is in "rununits" ! to "rundays" in units of days. ! segtim = the integration time "runlen" is broken into a number of ! segments each of length "segtim" days. updated surface ! boundary conditions are applied to MOM every "segtim" ! days. this is useful when coupling to atmospheric models ! in which case both models exchange surface boundary ! conditions every "segtim" days where "segtim" ! is 1/(coupling frequency). without an atmospheric model, ! when getting surface boundary conditions from data, ! "segtim" is set to the time step (in days) by mom.F. in ! either case, "runlen" (in days) should be an integral ! number of "segtim". ! nmix = number of time steps between mixing timesteps. used ! to damp timestep splitting due to centred leapfrog. ! init = (true,false) indicates that this run is a ! (start from initial conditions, restart) ! initbgc = (true,false) indicates that this run is a ! (start from initial BGC conditions, restart) ! restrt = (true,false) = (do,don`t) write a restart at the end ! of the run ! eb = (true,false) configures for the use of a ! (euler backward,forward) type mixing timestep ! init_time = (true,false) sets restarts to initial time ! init_time_in = (true,false) sets input restart to initial time ! init_time_out = (true,false) sets output restart to initial time !----------------------------------------------------------------------- ! inputs to tmngr.F: diagnostic intervals !----------------------------------------------------------------------- ! note: switches are used to control the interval between doing ! diagnostics. units for all switches are in days. ! setting a switch < 0.0 disables whatever the switch is ! controlling. setting it = 0.0 causes the diagnostic to be ! done every time step, and setting it > 0.0 causes the ! diagnostic to be done repeatedly on the specified interval. ! cmixint = number of days between writing estimated mixing coeffs ! on faces of T cells and U cells ! crossint = number of days between writing diapycnal and isopycnal ! components of flow ! fctint = number of days between writing difference between ! FCT and leapfrog advection ! densityint = number of days between writing density ! exconvint = number of days between writing temperature rate of ! change due to explicit convection ! glenint = number of days between global energetics integrals. ! trmbint = number of days between momentum and tracer term ! balances (global and regional). ! itrmb = (true,false) = (do,don`t) write regional mask info for ! the term balance diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! gyreint = number of days between calculation of tracer northward ! transport. ! igyre = (true,false) = (do,don`t) write regional mask info for ! the gyre diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! vmsfint = number of days between calculation of vertical and ! meridional stream function. ! tyzint = number of days between calculation of zonally averaged ! tracer components. ! prxzint = number of days between printouts of x-z data. ! extint = number of days between printouts of external mode. ! dspint = number of days between surface pressure calculation. ! Note: only when "diagnostic_surface_height" is enabled. ! dspper = averaging period for "diagnostic_surface_height" ! tavgint = number of days between regional tracer averages (under ! horizontal regions). ! itavg = (true,false) = (do,don`t) write regional mask info for ! the tracer average diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! tmbint = number of days over which tracer equation in averaged ! in depth and longitude to determine the meridional ! balance among storage, divergence, dissipation and ! forcing. ! tmbper = averaging period for "meridional_tracer_balance" ! itmb = (true,false) = (do,don`t) write "msktmb" for tracer ! the meridional balance diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! tsiint = number of days between printing of time step integrals. ! tsiper = averaging period for "time_step_monitor" ! stabint = number of days between sampling for various stability ! criteria. ! timavgint= interval (days) for writing time mean data from ! the "averaging" grid (only when "time_averages" is ! enabled). if "timavgint" is not an integral number of ! density time steps,"timavgint" is automatically adjusted ! to insure this. if the number of days to integrate is ! not an integral number of "timavgint" then the last ! averaging period will be less than "timavgint" days.this ! may lead to one more averaging period than expected. ! see "iounit.h" for more details. ! timavgper= averaging period for "time_averages" ! xbtint = averaging period (days) for writing XBT data (only when ! "xbts" is enabled). if "xbtint" is not an integral ! number of density time steps, "xbtint" is automatically ! adjusted to insure this. if the number of days to ! integrate is not an integral number of "xbtint" then the ! last averaging period will be less than "xbtint" days. ! this may lead to one more averaging period than ! expected. see "iounit.h" for more details. ! xbtper = averaging period for "xbts" ! zmbcint = number of days between calculation of zonal mean ! surface boundary conditions (and related quantities) ! restint = number of days between saving restarts ! tbtint = averaging period (days) for writing term balances ! tbtper = averaging period for tracer term balances !----------------------------------------------------------------------- ! outputs from tmngr.F: logical switches !----------------------------------------------------------------------- ! rundays = integration time in days (from "runlen") ! the following are logical counterparts to the above switches are ! set within "tmngr" every time step. logical switches control all ! decisions about when to do things in MOM. ! cmixts = (false,true) = (don`t, do) do write estimated mixing ! coefficients on this time step. ! based on "cmixint". ! crossts = (false,true) = (don`t, do) write diapycnal and ! isopycnal components of flow on this time step. ! based on "crossint". ! fctts = (false,true) = (don`t, do) write difference between ! FCT and leapfrog advection on this time step. ! based on "fctint". ! densityts = (false,true) = (don`t, do) write density on this time ! step. based on "densityint". ! exconvts = (false,true) = (don`t, do) do write temperature change ! due to explicit convection on this time step. ! based on "exconvint". ! glents = (false,true) = (don`t, do) do calculation of global ! energy integrals on this time step. based on "glenint". ! trmbts = (false,true) = (don`t, do) do calculation of momentum & ! tracer term balance on this timestep. based on "trmbint" ! gyrets = (false,true) = (don`t, do) do calculation of tracer ! northward transport on this timestep. based on "gyreint" ! vmsfts = (false,true) = (don`t, do) do calculation of vertical ! and meridional stream function on this time step. ! based on "vmsfint" ! tyzts = (false,true) = (don`t, do) do calculation of zonally ! averaged tracer components on this time step. ! based on "tyzint" ! prxzts = (false,true) = (don`t, do) do printouts of x-z data ! on this time step. based on "prxzint" ! extts = (false,true) = (don`t, do) do printout of external mode ! on this time step. based on "extint" ! dspts = (false,true) = (don`t, do) do calculation of diagnostic ! surface pressure on this time step. based on "dspint" ! stabts = (false,true) = (don`t, do) test for stability on this ! time step. based on "stabint" ! tavgts = (false,true) = (don`t do) do tracer averages on this ! time step. based on "tavgint" ! tmbts = (false,true) = (don`t, do) write out tracer meridional . ! balance on this time step. based on "tmbint" ! tsits = (false,true) = (don`t, do) print time step integrals ! on this time step. based on "tsiint" ! zmbcts = (false,true) = (don`t, do) print zonal mean boundary ! conditions on this time step. based on "zmbcint" ! timavgts = (false,true) = (don`t, do) write time mean data ! on this time step. based on "timavgint" ! xbtts = (false,true) = (don`t, do) write averaged XBT data on ! this time step based on "xbtint" ! restts = (false,true) = (don`t, do) save a restart on this time ! step based on "restint" ! tbtts = (false,true) = (don`t, do) write averaged tracer term ! balance data on this time step based on "tbint" ! leapfrog= (false,true) on a (mixing, normal leapfrog) time step ! based on "nmix" ! euler1 = true on the 1st pass of an euler backward time step ! otherwise false. (applies when "eb" = true) ! euler2 = true on the 2nd pass of an euler backward time step ! otherwise false. (applies when "eb" = true) ! forward = true on a forward time step. otherwise false ! (applies when "eb" = false) ! the following logical switches are based on the model time step. ! first = (true,false) = when it`s (the first, not the first) ! time step of a run ! eots = end of a time step. always true except for first ! pass of an euler backward time step ! eorun = last time step of a run. always false except during the ! last time step of the run. ! eoday = true when within 1/2 time step of the end of a day ! else ... false ! eoweek = true when within 1/2 time step of the end of a 7 day ! week (referenced to the start of a year) else ...false ! eo2wks = true when within 1/2 time step of the end of two weeks ! (referenced to the start of a year) else ... false ! midmon = true when within 1/2 time step of the middle of a month ! else ... false ! eomon = true when within 1/2 time step of the end of a month ! else ... false ! eoyear = true when within 1/2 time step of the end of a year ! else ... false ! osegs = true on the 1st time step of an ocean segment in mom.F ! otherwise false. ! osege = true on the last time step of an ocean segment in mom.F ! otherwise false. character(8) :: rununits common /switc_c/ rununits integer nmix, ieoday,ieoweek,ieo2wks integer ieomon,imidmon,ieoyear,ieorun common /switc_i/ nmix, ieoday,ieoweek,ieo2wks common /switc_i/ ieomon,imidmon,ieoyear,ieorun logical eb, leapfrog, euler1, euler2, forward, eots logical init, first, restrt, itavg, itmb, itrmb, igyre logical eoday, eoweek, eo2wks, eomon, midmon, eoyear, eorun common /switc_l/ eb, leapfrog, euler1, euler2, forward, eots common /switc_l/ init, first, restrt common /switc_l/ itavg, itmb, itrmb, igyre common /switc_l/ eoday, eoweek, eo2wks common /switc_l/ eomon, midmon, eoyear, eorun logical init_time, init_time_in, init_time_out, initbgc common /switc_l/ init_time, init_time_in, init_time_out, initbgc real runlen, rundays common /switc_r/ runlen, rundays !----------------------------------------------------------------------- ! S W I T C H E S B A S E D O N A N I N T E R V A L ! each interval switch needs three variables in common. The ! following naming convention is used. ! 1) an interval (real) for diagnostic output (e.g,. glenint) ! 2) a switch (logical) for the interval (e.g., glents ) ! the third is an internal variable needed by the time manager ! to support calculation of the logical switch ! 3) an index (integer) (e.g., iglenint) ! the user must specify the interval [e.g., glenint] for diagnostic ! output in units of days. tmngr sets the corresponding logical ! switch [e.g., glents] every time step. It is set to true when ! within half a time step of the requested interval, otherwise it is ! false. All decisions relating to the interval [e.g., glenint] ! are based on the logical switch [e.g., glents]. ! internal time structures ! The switch index [e.g., iglenint] is used to subsrcipt into ! internal arrays maintained by tmngr.F. The switch index is ! allocated on the first call to function "alarm". ! The array entry [e.g., iinterval(iglenint)] is a time index to the ! internal representation of the interval [e.g., glenint]. ! The array entry [e.g., ialarm(iglenint)] is a time index to the ! next time the alarm will be true. !----------------------------------------------------------------------- integer itavgint, iglenint, itrmbint, iprxzint common /switc_i/ itavgint, iglenint, itrmbint, iprxzint logical tavgts, glents, trmbts, prxzts common /switc_l/ tavgts, glents, trmbts, prxzts real tavgint, glenint, trmbint, prxzint common /switc_r/ tavgint, glenint, trmbint, prxzint integer iextint, iexconvint, icmixint common /switc_i/ iextint, iexconvint, icmixint logical extts, exconvts, cmixts common /switc_l/ extts, exconvts, cmixts real extint, exconvint, cmixint common /switc_r/ extint, exconvint, cmixint integer ivmsfint, igyreint, ityzint, ifctint common /switc_i/ ivmsfint, igyreint, ityzint, ifctint logical vmsfts, gyrets, tyzts, fctts common /switc_l/ vmsfts, gyrets, tyzts, fctts real vmsfint, gyreint, tyzint, fctint common /switc_r/ vmsfint, gyreint, tyzint, fctint integer istabint, izmbcint, icrossint, idensityint common /switc_i/ istabint, izmbcint, icrossint, idensityint logical stabts, zmbcts, crossts, densityts common /switc_l/ stabts, zmbcts, crossts, densityts real stabint, zmbcint, crossint, densityint common /switc_r/ stabint, zmbcint, crossint, densityint integer iosegs, iosege common /switc_i/ iosegs, iosege logical osegs, osege common /switc_l/ osegs, osege real segtim common /switc_r/ segtim integer irestint common /switc_i/ irestint logical restts common /switcl/ restts real restint common /switcr/ restint !----------------------------------------------------------------------- ! S W I T C H E S B A S E D O N A N I N T E R V A L ! A N D A V E R A G I N G P E R I O D ! each averaging period switch needs five variables in common. The ! following naming convention is used. ! 1) an interval (real) for diagnostic output (e.g. xbtint ) ! 2) a switch (logical) for the interval (e.g. xbtts ) ! 3) an averaging period (real) (e.g. xbtper ) ! 4) a switch (logical) for accumulating (e.g. xbtperts) ! the third is an internal variable needed by the time manager ! to support calculation of the logical switches ! 5) an index (integer) (e.g. ixbtint ) ! The user must specify the interval [e.g., xbtint] for diagnostic ! output in units of days and the averaging period [e.g., xbtper] ! in units of days. The averaging period may be less than or equal ! to the interval. For example, if the interval is 30.0 days and the ! averaging period is 5.0 days, results will be averaged over all ! time steps within days 26, 27, 28, 29, and 30. An averaging period ! of 0.0 days averages over the last time step of the interval (as ! does xbtper = dt), and an averaging period less than zero turns ! the switches off for all time steps. ! The logical switch for writing output at the specified interval ! [e.g., xbtts] is set to true on the last time step of the ! averaging period. The logical switch for accumulating results ! [e.g., xbtperts] is true for all time steps within the averaging ! period, otherwise it is false. ! internal time structures ! The index [e.g., ixbtint] is allocated on the first call to ! function "avg_alarm". The array element iperiod(ixbtint) is an ! index to the time structure for the internal representation of ! "xbtper", and ilastsw(ixbtint) is the index of the switch that ! flags the last time step of the accumulation period. ! Depending on use, ilastsw(ixbtint) may either be the index ! of another "named" switch or the index of a new switch ! allocated on the first time step. ! In the latter case, iinterval(ilastsw(ixbtint)) is the index of ! the time structure where "xbtint" is stored in internal form, ! and ialarm(ilastsw(ixbtint)) is the index of the time when an ! accumulation period will next end. ! The variable nextts(ixbtint) is true whenever the next ! time step will begin the accumulation period. !----------------------------------------------------------------------- integer ixbtint, idspint, itmbint, itimavgint common /switc_i/ ixbtint, idspint, itmbint, itimavgint logical xbtts, dspts, tmbts, timavgts logical xbtperts, dspperts, tmbperts, timavgperts common /switc_l/ xbtts, dspts, tmbts, timavgts common /switc_l/ xbtperts, dspperts, tmbperts, timavgperts real xbtint, dspint, tmbint, timavgint real xbtper, dspper, tmbper, timavgper common /switc_r/ xbtint, dspint, tmbint, timavgint common /switc_r/ xbtper, dspper, tmbper, timavgper integer itbtint, itsiint common /switc_i/ itbtint, itsiint logical tbtts, tsits logical tbtperts, tsiperts common /switc_l/ tbtts, tsits common /switc_l/ tbtperts, tsiperts real tbtint, tsiint real tbtper, tsiper common /switc_r/ tbtint, tsiint common /switc_r/ tbtper, tsiper !----------------------------------------------------------------------- ! S W I T C H E S B A S E D O N ! C A L E N D A R O R P R E V I O U S S W I T C H ! A N D A V E R A G I N G P E R I O D ! the following logical switches are based on any calendar or ! interval switch and an averaging period (in days). The averaging ! period must be less than or equal to the interval. The last ! time step of the averaging period is at the end of the interval. ! If the averaging period is set to zero, the averaging period ! consists only of the last time period of the interval. If ! the averaging period is less than zero, these switches are always ! false. ! each averaging period switch needs four variables in common. For ! example, if the averaging period is before the end of each month ! then the calendar switch (eomon), and index (ieomon) are presumed ! to exist in common and need not be added. ! Additionally, four items are needed. ! 1) an averaging period (real) (e.g. testper ) ! 2) a switch (logical) for accumulating results (e.g. testperts) ! 3) a switch (logical) for the end of interval (e.g. testts ) ! the fourth is an internal variable needed by the time manager ! to support calculation of the logical switch ! 4) an index (integer) (e.g. itestper ) ! Suppose it is required to produce averages over all time steps ! during the last 5 days of each month. Then "testper" = 5.0 and ! the following will calculate the accumulating switch. ! testts = avg_alarm(itestper, ihalfstep, 0, testper, iref, ieomon) ! testperts = on(itestper) ! Note the use of "ieomon" to key off the months. The switch ! "testts" will be true whenever "eomon" is true. ! Also note that when an averaging switch is keyed off another ! switch, the switch interval argument is not used, but is ! retained for consistency with the form of other averaging ! switches. !----------------------------------------------------------------------- integer maxsw parameter (maxsw=100) integer itestper, nsw, ialarm, iinterval, iperiod, ilastsw common /switc_i/ itestper, nsw, ialarm(maxsw), iinterval(maxsw) common /switc_i/ iperiod(maxsw), ilastsw(maxsw) logical testts, testperts, cplts, tsicpl, tavgcpl, on, lastts logical nextts common /switc_l/ testts, testperts, cplts, tsicpl, tavgcpl common /switc_l/ on(maxsw), lastts(maxsw), nextts(maxsw) real testint, testper common /switc_r/ testint, testper ! source file: /Users/dkeller/Desktop/UVic_ESCM/2.9/source/mom/timeavgs.h !====================== include file "timeavgs.h" ====================== ! imtav = # of longitudes used for the time averages grid ! jmtav = # of latitudes used for the time averages grid ! kmav = # of levels used for the time averages grid integer imtav, jmtav, kmav parameter (imtav=imt, jmtav=jmt-2, kmav=km) real ta_vetiso, ta_vntiso, ta_vbtiso common /ta_gm_r/ ta_vetiso(imt,km,jmt), ta_vntiso(imt,km,jmt) common /ta_gm_r/ ta_vbtiso(imt,km,jmt) real ta_dc14 common /ta_dc14/ ta_dc14(imt,km,jmt) real ta_rnpp, ta_rgraz, ta_rmorp, ta_rmorpt, ta_rmorz, ta_rexcr real ta_rremi, ta_rexpo, ta_rnpp_D, ta_rgraz_D, ta_rmorp_D real ta_rnfix, ta_rdeni, ta_rexpocal, ta_rprocal, ta_rgraz_Z real ta_rgraz_Det, ta_ravej, ta_ravej_D, ta_rgmax, ta_rno3P real ta_rpo4P, ta_rpo4_D real ta_rfeorgads, ta_rfecol real ta_rchl, ta_rchl_D real ta_rdeffe, ta_rremife, ta_rexpofe, ta_rfeprime real ta_rfesed, ta_rbfe ! wyao 20170427{ real ta_rfehydro ! wyao 20170427} ! wyao 20180627 { ! wyao 20180627} ! wyao 20170628{ ! wyao 20170628} common /ta_npzd_r/ ta_rnpp(imt,kpzd,jmt) common /ta_npzd_r/ ta_rgraz(imt,kpzd,jmt) common /ta_npzd_r/ ta_rmorp(imt,kpzd,jmt) common /ta_npzd_r/ ta_rmorpt(imt,kpzd,jmt) common /ta_npzd_r/ ta_rmorz(imt,kpzd,jmt) common /ta_npzd_r/ ta_rexcr(imt,kpzd,jmt) common /ta_npzd_r/ ta_rremi(imt,km,jmt) common /ta_npzd_r/ ta_rexpo(imt,km,jmt) common /ta_npzd_r/ ta_rexpocal(imt,km,jmt) common /ta_npzd_r/ ta_rprocal(imt,jmt) common /ta_npzd_r/ ta_rnpp_D(imt,kpzd,jmt) common /ta_npzd_r/ ta_rgraz_D(imt,kpzd,jmt) common /ta_npzd_r/ ta_rmorp_D(imt,kpzd,jmt) common /ta_npzd_r/ ta_rnfix(imt,kpzd,jmt) common /ta_npzd_r/ ta_rdeni(imt,km,jmt) common /ta_npzd_r/ ta_rgraz_Det(imt,kpzd,jmt) common /ta_npzd_r/ ta_rgraz_Z(imt,kpzd,jmt) common /ta_npzd_r/ ta_ravej(imt,kpzd,jmt) common /ta_npzd_r/ ta_ravej_D(imt,kpzd,jmt) common /ta_npzd_r/ ta_rgmax(imt,kpzd,jmt) common /ta_npzd_r/ ta_rno3P(imt,kpzd,jmt) common /ta_npzd_r/ ta_rpo4P(imt,kpzd,jmt) common /ta_npzd_r/ ta_rpo4_D(imt,kpzd,jmt) ! wyao 20170628{ ! wyao 20170628} common /ta_npzd_r/ ta_rfeorgads(imt,kpzd,jmt) common /ta_npzd_r/ ta_rfecol(imt,kpzd,jmt) common /ta_npzd_r/ ta_rdeffe(imt,kpzd,jmt) common /ta_npzd_r/ ta_rremife(imt,kpzd,jmt) common /ta_npzd_r/ ta_rexpofe(imt,kpzd,jmt) common /ta_npzd_r/ ta_rfeprime(imt,kpzd,jmt) common /ta_npzd_r/ ta_rfesed(imt,kpzd,jmt) common /ta_npzd_r/ ta_rbfe(imt,kpzd,jmt) common /ta_npzd_r/ ta_rchl(imt,kpzd,jmt) common /ta_npzd_r/ ta_rchl_D(imt,kpzd,jmt) ! wyao 20170427{ common /ta_npzd_r/ ta_rfehydro(imt,kpzd,jmt) ! wyao 20170427} integer nta_conv common /ta_conv_i/ nta_conv real ta_totalk, ta_vdepth, ta_pe common /ta_conv_r/ ta_totalk(imt,jmt), ta_vdepth(imt,jmt) common /ta_conv_r/ ta_pe(imt,jmt) integer nta_sscar common /ta_car_i/ nta_sscar real ta_sspH, ta_ssCO3, ta_ssOc, ta_ssOa, ta_sspCO2 common /ta_car_r/ ta_sspH(imt,jmt), ta_ssCO3(imt,jmt) common /ta_car_r/ ta_ssOc(imt,jmt), ta_ssOa(imt,jmt) common /ta_car_r/ ta_sspCO2(imt,jmt) !======================= include file "tmngr.h" ======================== ! time manager variables !----------------------------------------------------------------------- ! time manager inputs: !----------------------------------------------------------------------- ! how to choose a reference time: ! refrun = (true,false) to base calculation for diagnostic switches ! on (the start of each job, other reference time) ! example: ! suppose each job submission integrates ! for one month but the number of days per month changes. ! setting "refrun" = true and setting ! "timavgint" = (days in month)/3 will give 3 averaging ! periods per month of approximately 10 days each. the ! only restriction is that "timavgint"is an integral number ! of time steps (if not then "timavgint" is reset to insure ! this condition. other diagnostic switches do not have ! this restriction). ! refinit = (true, false) for basing calculation of logical switches ! on (initial conditions, other reference time) ! example: if term balances are desired every 20 days ! (trmbint=20.0) and refinit = true, then they ! will be done every 20 days starting from initial ! condition time. ! refuser = (true, false) to base calculations of logical switches ! on (user-chosen reference time, other reference time) ! if refuser = true, the user must also supply values for ! ryear, rmonth, rday, rhour, rmin, rsec (integer) ! example: if term balances are desired every 20 days ! (trmbint=20.0) and refuser = true, then they will be done ! every 20 days counting from reference time, ignoring the ! initial condition time. for comparing diagnostics from ! various experiments with different initial condition ! times, refuser = true will be more appropriate. setting ! refuser = true and choosing the reference time to be ! the initial condition time is the same as refinit = true. ! summary of how to choose the time for referencing calculations ! of logical switches ! refrun = T ==> referenced to the start of each run ! refinit = T ==> referenced to initial condition time given by: ! year0, month0, day0, hour0, min0, sec0 ! refuser = T ==> referenced to user specified reference time so ! must set: ryear, rmonth, rday, rhour, rmin, rsec !----------------------------------------------------------------------- ! time variable arrays ! arrays "iday" and "msday" contain the primary internal ! representation of all times within the time manager. they are ! referenced by using a subscript to indicate which time. ! iday = integer days (since Dec 31, 1899 when specifying a date) ! msday = non-negative integer milliseconds after midnight ! it is desirable to have time information expanded to include the ! following secondary time fields: ! year = ! month = ! day = ! hour = ! minute = ! second = ! tstamp = 32 character date and time stamp m/d/y h:m:s ! dayofyear = integer day of the year (1..yrlen) ! dayofweek = 1=sun - 7=sat ! daysinmon = days in the month ! daysinyear = days in the year ! those times for which primary and secondary information is ! maintained by the time manager are called "full times". those for ! which only primary information is kept are called "short times" ! indices to "full times" (including year, month ,day, etc). ! itime = simulation time corresponding to "itt" ! initial = time of the initial conditions ! irunstart = time of the start of the run ! iuser = user defined reference time ! iref = one of the three above selected by logicals ! (refinit, refrun, refuser) ! indices to "short times". ("iday", "msday" only) ! isunday = time of a sunday for week and two week switches ! ihalfstep = dt/2 beyond itime ! imodeltime = time since initial conditions ! iruntime = time since run start ! iusertime = time since user specified reference time ! idt = integer days and milliseconds of dt ! idtd2 = integer days and milliseconds of dt/2 ! ireftime = time used locally in alarm function ! for any time index (short or full) the internal representation ! may be converted to either real days or real seconds using ! the functions: ! realdays(index) ! realsecs(index) ! dayoyr = relative day number referenced to the beginning ! of the current year. (real) ! relyr = number of years (and fractional years) of model ! integration (for time tau+1 {itt}) relative to ! initial condition ! prelyr = relyr for previous time step ! stamp = 32 character date and time for current model timestep ! pstamp = 32 character date and time for previous model timestep ! calendar = calendar name ! itt = current time step counter (from initial cond.) ! itt0 = time step at start of current run ! variables used for initialization ! irstdy = integer number of days at start of run ! msrsdy = fractional day in millisec at start of run ! year0 = year of initial conditions ! month0 = month of initial conditions ! day0 = day of initial conditions ! hour0 = hour of initial conditions ! min0 = minute of initial conditions ! sec0 = second of initial conditions ! ryear = year of user specified reference time ! rmonth = month of user specified reference time ! rday = day of user specified reference time ! rhour = hour of user specified reference time ! rmin = minute of user specified reference time ! rsec = second of user specified reference time !----------------------------------------------------------------------- integer ntimes, nfulltimes parameter (ntimes = 100, nfulltimes = 30) character(32) :: tstamp, stamp, pstamp character(120) :: calendar common /tmngrc/ tstamp(nfulltimes), stamp, pstamp, calendar integer nextfulltime, nexttime integer initial, iref, irunstart, itime, iuser integer iruntime, imodeltime, ireftime, iusertime integer ihalfstep, isunday integer itemptime,itemptime2,itmptime,itmptime2,itmptime3 integer idt, idtd2 integer iday, msday, year, month, day, hour, minute, second integer dayofyear, dayofweek, daysinmon, daysinyear integer itt0, itt, irstdy, msrsdy integer year0, month0, day0, hour0, min0, sec0 integer ryear, rmonth, rday, rhour, rmin, rsec common /tmngri_i/ nextfulltime, nexttime common /tmngri_i/ initial, iref, irunstart, itime, iuser common /tmngri_i/ iruntime, imodeltime, ireftime, iusertime common /tmngri_i/ ihalfstep, isunday common /tmngri_i/ itemptime, itemptime2, itmptime, itmptime2 common /tmngri_i/ itmptime3 common /tmngri_i/ idt, idtd2 common /tmngri_i/ iday(ntimes), msday(ntimes) common /tmngri_i/ year(nfulltimes), month(nfulltimes) common /tmngri_i/ day(nfulltimes), hour(nfulltimes) common /tmngri_i/ minute(nfulltimes), second(nfulltimes) common /tmngri_i/ dayofyear(nfulltimes), dayofweek(nfulltimes) common /tmngri_i/ daysinmon(nfulltimes), daysinyear(nfulltimes) common /tmngri_i/ itt0, itt, irstdy, msrsdy common /tmngri_i/ year0, month0, day0, hour0, min0, sec0 common /tmngri_i/ ryear, rmonth, rday, rhour, rmin, rsec logical refrun, refinit, refuser common /tmngr_l/ refrun, refinit, refuser real dayoyr, relyr, prelyr common /tmngr_r/ dayoyr, relyr, prelyr !----------------------------------------------------------------------- ! acceleration parameters for long term forcing !----------------------------------------------------------------------- ! an acceleration factor can be used to accelerate long term ! forcing. forcing and output at time intervals of less than one ! year are not affected (eg. seasonal cycle of insolation). ! averaging intervals are also not changed. only longer term ! forcings are accelerated (eg. CO2 or orbital changes). ! runlen, runstep and most output intervals will will be divided by ! this factor (accel). output intervals are divided by accel if the ! interval divided by accel is larger or equal to a year. both ! runlen and runstep divided accel must still be evenly divisible ! by the coupling time (segtim). output intervals divided by accel ! should still be evenly divisible by the coupling time and the ! averaging period. ! times in various output files may appear to be inconsistent. ! for example the year may switch part way through a seasonal cycle. ! be cautious when looking at output averaged over less than a year. ! it is probably best to use an integer value for accel that will ! leave the runlen and runstep as integer multiples of a year. the ! initial year can be set to offset any existing time in a restart. ! choosing an inappropriate acceleration factor may mess up output ! intervals and averaging periods. think before you use. ! accel = factor for accelerating forcing transient data ! accel_yr0 = relyr to start acceleration real accel, accel_yr0 common /tmngr_r/ accel, accel_yr0 !====================== include file "vmixc.h" ========================= ! vertical mixing coefficients and related variables ! kappa_h = constant vertical diffusion coefficient (cm**2/sec) ! kappa_m = constant vertical viscosity coefficient (cm**2/sec) ! visc_cbu = viscosity coeff at bottom of "u" cell ! diff_cbt = diffusion coeff at bottom of "T" cell ! visc_cbu_limit = largest allowable "visc_cbu" ! diff_cbt_limit = largest allowable "diff_cbt" ! aidif = coefficient for implicit time differencing for ! vertical diffusion. aidif=1 gives the fully implicit ! case. aidif=0 gives the fully explicit case ! note: not used unless "implicitvmix" or "isopycmix" ! is enabled !======================================================================= real visc_cbu_limit, diff_cbt_limit, aidif, kappa_h, kappa_m real wndmix, fricmx, diff_cbt_back, visc_cbu_back, rhom1z real uzsq, diff_cbt, visc_cbt, visc_cbu, Ahv common /vmix_r/ visc_cbu_limit, diff_cbt_limit, aidif common /vmix_r/ kappa_h, kappa_m common /vmix_r/ visc_cbu(imt,km,jsmw:jemw) common /vmix_r/ diff_cbt(imt,km,jsmw:jemw) !====================== include file "diaga.h" ========================= ! variables used for computing diagnostics: ! totalk = total number of levels involved in convection ! vdepth = ventilation depth (cm) ! pe = potential energy lost due to explicit convection (g/s2) real totalk, vdepth, pe common /cdiaga_r/ totalk(imt,jsmw:jemw), vdepth(imt,jsmw:jemw) common /cdiaga_r/ pe(imt,jsmw:jemw) ! sspH = sea surface pH ! ssCO3 = sea surface CO3 concentration ! ssOc = sea surface Omega_c ! ssOa = sea surface Omega_a ! sspCO2 = sea surface pCO2 real sspH, ssCO3, ssOc, ssOa, sspCO2 common /cdiaga_r/ sspH(imt,jmt), ssCO3(imt,jmt), ssOc(imt,jmt) common /cdiaga_r/ ssOa(imt,jmt), sspCO2(imt,jmt) !======================== include file "ice.h" ========================= ! arrays for the ice model ! hice = ice thickness (tau-1, tau, tau+1) ! aice = ice area (tau-1, tau, tau+1) ! tice = ice surface temperature ! hsno = effective thickness of snow ! psno = precipitation as snow (+ accumulation, - melt) ! frzpt = freezing temperature of sea water ! cpts ice model ("ice_cpts") ! A = area ! heff = effective thickness of ice ! Ts = ice surface temperature ! groice = thermo ice thickness change from diagnostic ! E = enthalpy of each ice layer ! strength = pressure (from Rothrock), must have ncat > 3 ! hseff = effective thickness of snow ! ice dynamics ("ice_evp") ! uice = u ice velocity ! vice = v ice velocity ! pice = pressure due to internal stress ! xint = x component of ice interaction ! yint = y component of ice interaction ! del = delta for ice dynamics models ! eI = divergence of ice velocity ! nseg = number of domains needed for ice calculations ! jsi = start of limited domain for ice calculations ! jei = end of limited domain for ice calculations ! brine convection ("convect_brine") ! cbf = flux from rejected brine over ice categories ! cba = area of brine rejection for each categories ! cba0 = area of the cell without brine rejection ! land ice ("landice_data") ! hicel = land ice thickness ! aicel = land ice area integer nseg, jsi, jei real hice, aice, tice, hsno, psno, frzpt real A, heff, Ts, groice, E, hseff, strength real uice, vice, pice, xint, yint, del, eI real cbf, cba, cba0, hicel, aicel common /ice_r/ hice(imt,jmt,2), aice(imt,jmt,2), tice(imt,jmt) common /ice_r/ hsno(imt,jmt,2), psno(imt,jmt), frzpt(imt,jmt) integer ncat parameter (ncat=1) common /ice_r/ uice(imt,jmt), vice(imt,jmt), pice(imt,jmt) common /ice_r/ xint(imt,jmt), yint(imt,jmt) common /ice_r/ del(imt,jmt), eI(imt,jmt) common /ice_i/ nseg, jsi(jmt), jei(jmt) !======================== include file "atm.h" ========================= ! arrays for the energy-moisture balance model ! note: units for heat flux are: g s-3 (or mW m-2) ! units for fresh water flux are in g cm-2 s-1 ! downward is into the top surface (ocean, ice or land) ! upward is into the bottom of the atmosphere ! outward is out of the top of the atmosphere ! inward is into the top of the atmosphere ! at = tracers (previous and most recent) ! elev = elevations (cm) ! flux = flux accumulator (1=heat, 2=water, 3=carbon, ...) ! rh = relative humidity (%/100) ! tmsk = tracer grid ocean mask (0 = land, 1 = ocean) ! umsk = velocity grid ocean mask (0 = land, 1 = ocean) ! dn = northward atmospheric tracer diffusivity ! de = eastward atmospheric tracer diffusivity ! fcor = Coriolis factor ! heat fluxes ! solins = incoming short wave radiation (g s-3) ! dnswr = downward surface shortwave flux absorbed (g s-3) ! uplwr = upward surface longwave flux (g s-3) ! upsens = upward surface sensible heat flux (g s-3) ! upltnt = upward surface latent heat flux (g s-3) ! outlwr = outgoing atmosphere longwave flux (g s-3) ! fresh water fluxes ! precip = precipitation (g cm-2 s-1) ! evap = evaporation (g cm-2 s-1) ! disch = discharge from continents (g cm-2 s-1) ! vflux = normalized virtual flux ! land model ! soilm = soil moisture, as depth in bucket (g cm -2) ! runoff = water runoff from continents (g cm-2 s-1) ! surf = land surface temperature (C) ! annual average temperature ("embm_awind" and "embm_adiff") ! rtbar = running average atmospheric temperature (C) ! atbar = average temperature accumulator (C) ! tbar = climatological average atmospheric temperature (C) ! wind feedback ("embm_awind") ! awx = anomalous x component of wind (cm s-1) ! awy = anomalous y component of wind (cm s-1) ! apress = anomalous pressure (g cm-2 s-2) ! sulphates ("sulphate_data_transient") ! sulph = sulphate forcing ! distributed CO2 emissions ("co2emit_data" and "carbon_co2_2d") ! co2dist = distributed co2 forcing ! sea level ("sealev" or "sealev_data") ! elev_sealev = elevation anomaly due to sea level changes ! flux adjustment ("save_flxadj") ! flxadj = flux adjustment ! co2 emissions from tracking co2 ("co2emit_track_co2") ! track_co2 = array of global average co2 values ! co2 emissions from tracking sat("co2emit_track_sat" or "embm_vcs") ! track_sat = array of global average sat values ! indicies for atmospheric tracer array ! isat = index for surface air temperature ! ishum = index for surface specific humidity ! ico2 = index for atmospheric co2 ! mapat = map for atmospheric tracer names ! mtrack = maximum size for tracking arrays character(10) :: mapat(nat) common /embm_c/ mapat integer isat, ishum, ico2, mtrack parameter (mtrack=3650) common /embm_i/ isat, ishum, ico2 real at, elev, flux, rh, tmsk, umsk, dn, de, fcor, solins, dnswr real uplwr, upsens, upltnt, outlwr, precip, evap, disch, vflux real soilm, runoff, surf, rtbar, atbar, tbar, awx, awy, apress real sulph, co2dist, elev_sealev, flxadj, track_co2, track_sat common /embm_r/ at(imt,jmt,2,nat), elev(imt,jmt) common /embm_r/ flux(imt,jmt,nat+2), rh(imt,jmt) common /embm_r/ tmsk(imt,jmt), umsk(imt,jmt), dn(imt,jmt,nat) common /embm_r/ de(imt,jmt,nat), fcor(imt,jmt), solins(imt,jmt) common /embm_r/ dnswr(imt,jmt), uplwr(imt,jmt), upsens(imt,jmt) common /embm_r/ upltnt(imt,jmt), outlwr(imt,jmt), precip(imt,jmt) common /embm_r/ evap(imt,jmt), disch(imt,jmt), vflux(imt,jmt) common /embm_r/ soilm(imt,jmt,2), runoff(imt,jmt) common /embm_r/ surf(imt,jmt) common /embm_r/ rtbar(imt,jmt), atbar(imt,jmt), tbar(imt,jmt) !====================== include file "npzd.h" ========================= ! variables for npzd model ! ntnpzd = number of npzd tracers ! nbio = number of npzd timesteps per ocean timestep ! trcmin = minimum tracer for flux calculations ! alpha = initial slope P-I curve [(W/m^2)^(-1)/day] ! kw = light attenuation due to water [1/m] ! kc = light attenuation by phytoplankton [1/(m*mmol m-3)] ! ki = light attenuation through sea ice & snow ! abio = maximum growth rate parameter [1/day] ! bbio = b ! cbio = [1/deg_C] ! k1n = half saturation constant for N uptake [mmol m-3] ! gamma1 = assimilation efficiency (zpk) ! gbio = maximum grazing rate at 0 deg C [day-1] ! nuz = quadratic mortality (zpk) ! nud0 = remineralization rate [day-1] !{ wyao 20170707 ! wyao 20170707} ! LFe = Iron limitation ! wd = sinking speed of detritus [m day-1] ! ztt = depth to top of grid cell [cm] ! rkwz = reciprical of light attenuation times grid depth ! par = fraction of photosythetically active radiation ! dtnpzd = time step of biology ! capr = carbonate to carbon production ratio ! dcaco3 = remineralisation depth of calcite [cm] ! rcak = array used in calculating calcite remineralization ! rcab = array used in calculating bottom calcite remineralization !{ wyao 20170707 ! nup = quick recycling rate(Phytoplankton)[1/day] ! nupt0 = quick recycling rate(Temprtr depnd.Phytoplankton)[1/day] ! wyao 20170707} ! wd0 = sinking speed of detritus at surface [m/day] ! k1p = half saturation constant for P uptake ! jdiar = factor reducing the growth rate of diazotrophs ! redctn = C/N Redfield ratio (includes mol to mmol conversion) ! redctp = C/P Redfield ratio (includes mol to mmol conversion) ! redptn = P/N Redfield ratio ! redntp = N/P Redfield ratio ! redotn = O/N Redfield ratio (includes mol to mmol conversion) ! redotp = O/P Redfield ratio (includes mol to mmol conversion) ! rnbio = reciprical of nbio ! rdtts = reciprical of dtts [s-1] ! dtbio = npzd time step [s] ! rnpp = rate of net primary production [nmol cm-3 s-1] ! rgraz = rate of grazing [nmol cm-3 s-1] ! rmorp = rate of mortality of phytoplankton [nmol cm-3 s-1] !{ wyao 20170707 ! rmorpt = rate of quick recy of phytoplankton [nmol cm-3 s-1] ! wyao 20170707} ! rmorz = rate of mortality of zooplankton [nmol cm-3 s-1] ! rremi = rate of remineralization [nmol cm-3 s-1] ! rexcr = rate of excretion [nmol cm-3 s-1] ! rexpo = rate of export through the bottom [nmol cm-3 s-1] ! rnpp_D = npp for diazotraphs [nmol cm-3 s-1] ! rgraz_D = rgraz for diazotraphs [nmol cm-3 s-1] !{ wyao 20170707 ! rmorp_D = rate of mortality for diazotraphs [nmol cm-3 s-1] ! wyao 20170707} ! rnfix = rate of nitrogen fixation [nmol cm-3 s-1] ! rdeni = rate of denitrification [nmol cm-3 s-1] ! kzoo = half saturation constant for Z grazing ! zprefP = Z preference for grazing on P ! zprefD = Z preference for grazing on Diaz ! zprefZ = Z preference for grazing on other Z ! zprefDet = Z preference for grazing on Detritus ! rgraz_Det = rate of grazing on Detritus [nmol cm-3 s-1] ! rgraz_Z = rate of grazing on other Zooplankton [nmol cm-3 s-1] ! geZ = growth efficiency of zooplankton !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! New diagnostic output ! ravej = light-dependant growth rate of P ! ravej_D = light-dependant growth rate of Diaz ! rgmax = temp-dependant growth rate of zoo ! rno3P = nitrate-dependant growth rate of P ! rpo4P = phosphate-dependant growth rate of P ! rpo4_D = phosphate-dependant growth rate of D ! ! fe_dissolved = dissolved iron concentration ! kfe = Fe limitation half saturation parameter ! kfe_D = Fe limitation half sat. param. for diaz. ! kmfe = number of depth levels for the iron mask ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Dynamic iron cycle ! kfeleq = Fe-ligand stability constant [m^3/(mmol ligand)] ! lig = Ligand concentration [mmol/m^3] ! thetamaxhi = Maximum Chl:C ratio, abundant iron [gChl/(gC)] ! thetamaxlo = Maximum Chl:C ratio, extreme iron limitation [gChl/(gC)] ! alphamax = Maximum initial slope in PI-curve [day^-1 (W/m^2)^-1 (gChl/(gC))^-1] ! alphamin = Minimum intital slope in PI-curve [day^-1 (W/m^2)^-1 (gChl/(gC))^-1] ! mc = Molar mass of carbon [g/mol] ! fetopsed = Fe:P ratio for sedimentary iron source [molFe/molP] ! o2min = Minimum O2 concentration for aerobic respiration [mmolO_2/m^3] ! kfeorg = Organic-matter dependent scavenging rate [(m^3/(gC s))^0.58] ! kfeorg_ps = Organic-matter dependent scavenging power lay scale parameter default 0.58 ! kfecol = Colloidal Fe production rate (inorganic scavenging) [s^-1] ! fehydro = hydrothermal flux of iron [mol/m^3*s^-1] !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! DOP and DON guessing by Ben ! option O_npzd_dop & O_npzd_don ! dfr = refractory/semi-labile DOM fraction from mortality. ! dfrt = refractory/semi-labile DOM fraction from quick recy. ! redotc = o/c redfield ratio ! redntc = n/c redfield ratio ! redotc = b/c redfield ratio ! redptc = p/c redfield ratio ! hdop = DOP growth rate handicap ! rnpp_dop = rate of net primary production from DOP [nmol cm-3 s-1] ! rnpp_D_dop = rate of net primary production for diaz from DOP [nmol cm-3 s-1] ! diazntp = n/p diazo ratio ! diazptn = p/N diazo ratio ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Sinking speed parameters ! option O_npzd_caped_sinkspeed ! wdi the dependency on depth ! wdit the handicaped depth, fixed speed after certain depth. ! ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! non-redfield iron (diaz and phyt have different uptake fe to N ratio) ! option O_npzd_iron_diazratio ! ! rfeton_D Uptake ratio of iron to nitrogen for Diaz integer ntnpzd, nbio, kmfe parameter (ntnpzd = 4 $ + 2 $ + 2 $ ) common /npzd_i/ nbio(km) real trcmin parameter (trcmin=5e-12) real alpha, kw, kc, ki, abio, bbio, cbio, k1n, nup, gamma1, gbio real epsbio, nuz, gamma2, nud0, LFe, wd, ztt, rkwz, par, dtnpzd real capr, dcaco3, rcak, rcab, nupt0, wd0, k1p, jdiar, redctn real redctp, redptn, redntp, redotn, redotp, rnbio, rdtts, dtbio real rnpp, rgraz, rmorp, rmorpt, rmorz, rremi, rexcr, rexpo !{ wyao 20170707 real wdi, wdit ! wyao 20170707} real rnpp_D, rgraz_D, rmorp_D, rnfix, rdeni, kzoo, zprefP real zprefD, zprefZ, zprefDet, rgraz_Det, rgraz_Z, geZ, kfe real ravej, ravej_D, rgmax, rno3P, rpo4P, rpo4_D, kfe_D real kfemax, kfemin, pmax real kfeleq, lig, thetamaxhi, alphamax, alphamin real thetamaxlo, mc, fetopsed, o2min, kfeorg, rfeton real kfecol real kfeorg_ps real rfeorgads, rdeffe, rremife, rexpofe, rfeprime real rfesed, rbfe, rfecol real fe_hydr real rchl, rchl_D real fehydro, rfehydro common /npzd_r/ alpha, kw, kc, ki, abio, bbio, cbio, k1n, nup common /npzd_r/ gamma1, gbio, epsbio, nuz, nud0, LFe common /npzd_r/ wd(km), ztt(km), rkwz(km), par, dtnpzd, capr common /npzd_r/ dcaco3, rcak(km), rcab(km), nupt0, wd0, k1p common /npzd_r/ jdiar, redctn, redctp, redptn, redntp, redotn common /npzd_r/ redotp, rnbio(km), rdtts(km), dtbio(km), geZ common /npzd_r/ kzoo, zprefP, zprefD, zprefZ, zprefDet, kfe, kfe_D common /npzd_r/ kfeleq, alphamax, alphamin common /npzd_r/ thetamaxhi, thetamaxlo, lig, fetopsed, o2min common /npzd_r/ mc, kfeorg, rfeton, kfecol, kfeorg_ps common /npzd_r/ kfemax, kfemin, pmax, fe_hydr(imt,jmt,km) !{ wyao 20170711 ! wyao 20170711} !{ wyao 20170707 common /npzd_r/ wdi ! wyao 20170707} common /npzd_r/ rnpp(imt,kpzd,jsmw:jemw) common /npzd_r/ rgraz(imt,kpzd,jsmw:jemw) common /npzd_r/ rmorp(imt,kpzd,jsmw:jemw) common /npzd_r/ rmorpt(imt,kpzd,jsmw:jemw) common /npzd_r/ rmorz(imt,kpzd,jsmw:jemw) common /npzd_r/ rexcr(imt,kpzd,jsmw:jemw) common /npzd_r/ rremi(imt,km,jsmw:jemw) common /npzd_r/ rexpo(imt,km,jsmw:jemw) common /npzd_r/ rgraz_Det(imt,kpzd,jsmw:jemw) common /npzd_r/ rgraz_Z(imt,kpzd,jsmw:jemw) common /npzd_r/ rnpp_D(imt,kpzd,jsmw:jemw) common /npzd_r/ rgraz_D(imt,kpzd,jsmw:jemw) common /npzd_r/ rmorp_D(imt,kpzd,jsmw:jemw) common /npzd_r/ rnfix(imt,kpzd,jsmw:jemw) common /npzd_r/ rdeni(imt,km,jsmw:jemw) common /npzd_r/ rremife(imt,km,jsmw:jemw) common /npzd_r/ rexpofe(imt,km,jsmw:jemw) ! wyao 20170427{ common /npzd_r/ fehydro(imt,jmt,km) common /npzd_r/ rfehydro(imt,km,jsmw:jemw) ! wyao 20170427} common /npzd_r/ rfeorgads(imt,km,jsmw:jemw) common /npzd_r/ rdeffe(imt,km,jsmw:jemw) common /npzd_r/ rfeprime(imt,km,jsmw:jemw) common /npzd_r/ rfesed(imt,km,jsmw:jemw) common /npzd_r/ rbfe(imt,km,jsmw:jemw) common /npzd_r/ rfecol(imt,km,jsmw:jemw) common /npzd_r/ rchl(imt,km,jsmw:jemw) common /npzd_r/ rchl_D(imt,km,jsmw:jemw) parameter (istrt=1, iend=imt) real twodt(km) real snpzd(ntnpzd), tnpzd(ntnpzd) real src(imt,km,jsmw:jemw,nsrc) common/uvokcomm/src !======================== include file "isopyc.h" ====================== ! isopycnal diffusion variables: ! ahisop = isopycnal tracer mixing coefficient (cm**2/sec) ! drodx = d(rho)/dx local to east face of T cell ! drody = d(rho)/dy local to north face of T cell ! drodz = d(rho)/dz local to bottom face of T cell ! Ai_e = diffusion coefficient on eastern face of T cell ! Ai_n = diffusion coefficient on northern face of T cell ! Ai_bx = diffusion coefficient on bottom face of T cell ! Ai_by = diffusion coefficient on bottom face of T cell ! fisop = structure function for isopycnal diffusion coefficient. ! slmxr = reciprocal of maximum allowable slope of isopycnals for ! small angle approximation real alphai, betai common /cisop_r/ alphai(imt,km,jmw), betai(imt,km,jmw) real ddxt, ddyt, ddzt, Ai_ez, Ai_nz, Ai_bx, Ai_by, K11, K22, K33 real ahisop, fisop, slmxr, del_dm, s_dmr common /cisop_r/ ddxt(imt,km,jsmw:jemw,2) common /cisop_r/ ddyt(imt,km,1:jemw,2) common /cisop_r/ ddzt(imt,0:km,jmw,2) common /cisop_r/ Ai_ez(imt,km,jsmw:jemw,0:1,0:1) common /cisop_r/ Ai_nz(imt,km,1:jemw,0:1,0:1) common /cisop_r/ Ai_bx(imt,km,jsmw:jemw,0:1,0:1) common /cisop_r/ Ai_by(imt,km,jsmw:jemw,0:1,0:1) common /cisop_r/ K11(imt,km,jsmw:jemw) common /cisop_r/ K22(imt,km,1:jemw) common /cisop_r/ K33(imt,km,jsmw:jemw) common /cisop_r/ ahisop, fisop(imt,jmt,km), slmxr real delta_iso, s_minus, s_plus common /cisop_r/ delta_iso, s_minus, s_plus ! adv_vetiso = zonal isopycnal mixing velocity computed at the ! center of the eastern face of the "t" cells ! adv_vntiso = meridional isopycnal mixing velocity computed at ! the center of the northern face of the "t" cells ! (Note: this includes the cosine as in "adv_vnt") ! adv_vbtiso = vertical isopycnal mixing velocity computed at the ! center of the top face of the "t" cells ! adv_fbiso = "adv_vbtiso" * (tracer) evaluated at the center of ! the bottom face of the "t" cells ! athkdf = isopycnal thickness diffusivity (cm**2/sec) real athkdf, adv_vetiso, adv_vntiso, adv_vbtiso, adv_fbiso common /cisop_r/ athkdf common /cisop_r/ adv_vetiso(imt,km,jsmw:jemw) common /cisop_r/ adv_vntiso(imt,km,1:jemw) common /cisop_r/ adv_vbtiso(imt,0:km,jsmw:jemw) common /cisop_r/ adv_fbiso(imt,0:km,jsmw:jemw) real drodxe, drodze, drodyn, drodzn, drodxb, drodyb, drodzb real drodye, drodxn ! statement functions drodxe(i,k,j,ip) = alphai(i+ip,k,j)*ddxt(i,k,j,1) + & betai(i+ip,k,j)*ddxt(i,k,j,2) drodze(i,k,j,ip,kr) = alphai(i+ip,k,j)*ddzt(i+ip,k-1+kr,j,1) + & betai(i+ip,k,j)*ddzt(i+ip,k-1+kr,j,2) drodyn(i,k,j,jq) = alphai(i,k,j+jq)*ddyt(i,k,j,1) + & betai(i,k,j+jq)*ddyt(i,k,j,2) drodzn(i,k,j,jq,kr) = alphai(i,k,j+jq)*ddzt(i,k-1+kr,j+jq,1) + & betai(i,k,j+jq)*ddzt(i,k-1+kr,j+jq,2) drodxb(i,k,j,ip,kr) = alphai(i,k+kr,j)*ddxt(i-1+ip,k+kr,j,1) + & betai(i,k+kr,j)*ddxt(i-1+ip,k+kr,j,2) drodyb(i,k,j,jq,kr) = alphai(i,k+kr,j)*ddyt(i,k+kr,j-1+jq,1) + & betai(i,k+kr,j)*ddyt(i,k+kr,j-1+jq,2) drodzb(i,k,j,kr) = alphai(i,k+kr,j)*ddzt(i,k,j,1) + & betai(i,k+kr,j)*ddzt(i,k,j,2) !====================== include file "fdift.h" ========================= ! finite difference numerics for tracers !======================================================================= T_i(i,k,j,n,ip) = t(i+ip,k,j,n,taum1) T_j(i,k,j,n,jp) = t(i,k,j+jp,n,taum1) dz_t2r(i,k,j) = dzt2r(k) dz_tr(i,k,j) = dztr(k) dz_wtr(i,k,j) = dzwr(k) dx_t2r(i,k,j) = cstdxt2r(i,j) dx_tr(i,k,j) = cstdxtr(i,j) dy_t2r(i,k,j) = cstdyt2r(jrow) dy_tr(i,k,j) = cstdytr(jrow) !----------------------------------------------------------------------- ! advective terms !----------------------------------------------------------------------- ADV_Tx(i,k,j) = (adv_fe(i,k,j) - adv_fe(i-1,k,j))*cstdxt2r(i,j) ADV_Ty(i,k,j,jrow,n) = (adv_fn(i,k,j) - adv_fn(i,k,j-1)) & *cstdyt2r(jrow) ADV_Tz(i,k,j) = (adv_fb(i,k-1,j) - adv_fb(i,k,j))*dzt2r(k) ! gent_mcwilliams isopycnal advective terms simulating the effect ! of eddies on the isopycnals ADV_Txiso(i,k,j,n) = cstdxt2r(i,j)*(adv_vetiso(i,k,j) & *(t(i+1,k,j,n,taum1) + t(i,k,j,n,taum1)) - adv_vetiso(i-1,k,j) & *(t(i,k,j,n,taum1) + t(i-1,k,j,n,taum1))) ADV_Tyiso(i,k,j,jrow,n) = cstdyt2r(jrow)*(adv_vntiso(i,k,j) & *(t(i,k,j+1,n,taum1) + t(i,k,j,n,taum1)) - adv_vntiso(i,k,j-1) & *(t(i,k,j,n,taum1) + t(i,k,j-1,n,taum1))) ADV_Tziso(i,k,j) = dzt2r(k)*(adv_fbiso(i,k-1,j)-adv_fbiso(i,k,j)) !----------------------------------------------------------------------- ! diffusive terms !----------------------------------------------------------------------- ! zonal component DIFF_Tx(i,k,j) = (diff_fe(i, k,j)*tmask(i+1,k,j) & - diff_fe(i-1,k,j)*tmask(i-1,k,j))*cstdxtr(i,j) ! meridional component DIFF_Ty(i,k,j,jrow,n) = (diff_fn(i,k,j )*tmask(i,k,j+1) & - diff_fn(i,k,j-1)*tmask(i,k,j-1))*cstdytr(jrow) ! vertical component DIFF_Tz(i,k,j) = (diff_fb(i,k-1,j) - diff_fb(i,k,j))*dztr(k) & *(c1-aidif) & + (diff_fbiso(i,k-1,j) - diff_fbiso(i,k,j))*dztr(k) !----------------------------------------------------------------------- ! calculation of biological interactions !----------------------------------------------------------------------- declin = sin((mod(relyr,1.) - 0.22)*2.*pi)*0.4 ! declination do k=1,km twodt(k) = c2dtts*dtxcel(k) nbio(k) = twodt(k)/dtnpzd dtbio(k) = twodt(k)/nbio(k) rdtts(k) = 1./twodt(k) rnbio(k) = 1./nbio(k) enddo tap = 2.*par ! alpha is iron dependent and calculated in npzd_src do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .gt. 0) then ai = aice(i,jrow,2) hi = hice(i,jrow,2) hs = hsno(i,jrow,2) ! calculate day fraction and incoming solar ! angle of incidence = lat - declin, refraction index = 1.33 rctheta = max(-1.5, min(1.5, tlat(i,jrow)/radian - declin)) rctheta = kw/sqrt(1. - (1. - cos(rctheta)**2.)/1.33**2.) dayfrac = min( 1., -tan(tlat(i,jrow)/radian)*tan(declin)) dayfrac = max(1e-12, acos(max(-1., dayfrac))/pi) swr = dnswr(i,jrow)*1e-3*(1. + ai*(exp(-ki*(hi + hs)) - 1.)) expo = 0.0 impo = 0.0 expofe = 0.0 impofe = 0.0 phin = 0.0 ! integrated phytoplankton prca = 0.0 ! integrated production of calcite kmx = min(kmt(i,jrow), kpzd) do k=1,kmx !----------------------------------------------------------------------- ! Initialize tracers. ! In the old code tracers were limited to positive values here ! this should not be done here! Only in npzd_src.F !----------------------------------------------------------------------- tnpzd(1) = t(i,k,j,ipo4,taum1) tnpzd(2) = t(i,k,j,iphyt,taum1) tnpzd(3) = t(i,k,j,izoop,taum1) tnpzd(4) = t(i,k,j,idetr,taum1) tnpzd(5) = t(i,k,j,ino3,taum1) tnpzd(6) = t(i,k,j,idiaz,taum1) tnpzd(7) = t(i,k,j,idfe,taum1) tnpzd(8) = t(i,k,j,idetrfe,taum1) !{wyao 20171018 ! if ((k.eq.1).and.(i.eq.50).and.(j.eq.50)) then ! print*, 'ironvalue20170428.', tnpzd(7) ! endif ! wyao 20171018} !{ wyao 20170713 ! wyao 20170713} swr = swr*exp(-kc*phin) ! phin has been removed from the following code to fix a bug phin = (max(tnpzd(6),trcmin) + max(tnpzd(2),trcmin)) & *dzt(k) gl = tap*swr*exp(ztt(k)*rctheta) impo = expo*dztr(k) impofe = expofe*dztr(k) bct = bbio**(cbio*t(i,k,j,itemp,taum1)) if (t(i,k,j,itemp,taum1).gt.20) then bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1)) & *bbio**(cbio*20) else bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1)) & *bct endif ! Decrease remineralisation rate in oxygen minimum zone and when there is no ! nitrate. if (t(i,k,j,ino3,taum1).le.0.) then nud = 0 else nud = nud0*(0.65+0.35*tanh(t(i,k,j,io2,taum1) & *1000.-6.)) endif !{ wyao 20170713 ! wyao 20170713} !----------------------------------------------------------------------- ! call the npzd model !----------------------------------------------------------------------- call npzd_src (tnpzd, nbio(k), dtbio(k), gl, bct, impo &, dzt(k), dayfrac, wd(k), rkwz(k), nud &, snpzd, expo, graz, morp, morz, graz_Det &, graz_Z &, npp, morpt, remi, excr &, npp_D, graz_D, morp_D, nfix &, bctz &, expofe, impofe, remife &, feorgads, thetamax, deffe, feprime &, fecol &, t(i,k,j,io2,taum1) &, chl, thetachl &, chl_D & ) ! These are source/sink terms snpzd(1:4) = snpzd(1:4)*rdtts(k) snpzd(5:6) = snpzd(5:6)*rdtts(k) snpzd(7:8) = snpzd(7:8)*rdtts(k) expofe=expofe*rnbio(k) !{ wyao 20170713 ! wyao 20170713} expo = expo*rnbio(k) rexpo(i,k,j) = expo rgraz(i,k,j) = graz*rnbio(k) rgraz_Det(i,k,j) = graz_Det*rnbio(k) rgraz_Z(i,k,j) = graz_Z*rnbio(k) rmorp(i,k,j) = morp*rnbio(k) rmorz(i,k,j) = morz*rnbio(k) rnpp(i,k,j) = npp*rnbio(k) !{ wyao 20170713 rmorpt(i,k,j) = morpt*rnbio(k) rremi(i,k,j) = remi*rnbio(k) rexcr(i,k,j) = excr*rnbio(k) rnpp_D(i,k,j) = npp_D*rnbio(k) ! wyao 20170713} rgraz_D(i,k,j) = graz_D*rnbio(k) rmorp_D(i,k,j) = morp_D*rnbio(k) rnfix(i,k,j) = nfix*rnbio(k) rexpofe(i,k,j) = expofe rremife(i,k,j) = remife*rnbio(k) rfeorgads(i,k,j) = feorgads*rnbio(k) rdeffe(i,k,j) = deffe*rnbio(k) rfeprime(i,k,j) = feprime*rnbio(k) rfecol(i,k,j) = fecol*rnbio(k) rchl(i,k,j) = chl*rnbio(k) rchl_D(i,k,j) = chl_D*rnbio(k) !----------------------------------------------------------------------- ! calculate detritus at the bottom and remineralize !----------------------------------------------------------------------- sgb = sg_bathy(i,j,k) if (sgb .gt. 0) then rremi(i,k,j) = rremi(i,k,j) + expo*sgb snpzd(1) = snpzd(1) + redptn*expo*sgb snpzd(5) = snpzd(5) + expo*sgb fesed=fetopsed*bct*redptn*expo*sgb bfe=fesed snpzd(7) = snpzd(7) + fesed if (t(i,k,j,io2,taum1)*1e3 .lt. o2min) then snpzd(7) = snpzd(7) + expofe*sgb bfe=bfe+expofe*sgb expofe = expofe - sgb * expofe endif rremife(i,k,j) = rremife(i,k,j) + bfe rfesed=fesed rbfe=bfe expo = expo - sgb * expo endif !----------------------------------------------------------------------- ! set source/sink terms !----------------------------------------------------------------------- src(i,k,j,ispo4) = snpzd(1) src(i,k,j,isphyt) = snpzd(2) src(i,k,j,iszoop) = snpzd(3) src(i,k,j,isdetr) = snpzd(4) src(i,k,j,isno3) = snpzd(5) src(i,k,j,isdiaz) = snpzd(6) if ( k .eq. 1 ) then snpzd(7) = snpzd(7) + ! comes in as mol/m2/s gets converted to mmol/m3/s & max(0.0,sbc(i,j,idfeadep)*1000./(dzt(1)/100.)) !dzt is in cm endif ! wyao 20170427{ ! add hydrothermal flux of iron snpzd(7) = snpzd(7) + & fehydro(i,j,k) !dzt is in cm rfehydro(i,k,j)=fehydro(i,j,k) ! print*, 'hydrotest20170428.', fehydro(i,j,k) ! wyao 20170427} src(i,k,j,isdfe) = snpzd(7) src(i,k,j,isdetrfe) = snpzd(8) !{ wyao 20170713 ! wyao 20170713} ! production of calcite dprca = (morp+morz+(graz+graz_Z)*(1.-gamma1)) & *capr*redctn*rnbio(k) prca = prca + dprca*dzt(k) ! These are sources and sinks of DIC (i.e. remin - pp) ! all are based on po4 uptake and remineralization ! dprca is a correction term src(i,k,j,isdic) = (snpzd(1)*redctp - dprca) src(i,k,j,isalk) = (-snpzd(1)*redntp*1.e-3 - 2.*dprca) !----------------------------------------------------------------------- ! accumulate time averages !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then ta_rnpp(i,k,jrow) = ta_rnpp(i,k,jrow) + rnpp(i,k,j) !{ wyao 20170713 ! wyao 20170713} ta_rgraz(i,k,jrow) = ta_rgraz(i,k,jrow) + rgraz(i,k,j) ta_rgraz_Z(i,k,jrow) = ta_rgraz_Z(i,k,jrow) & + rgraz_Z(i,k,j) ta_rgraz_Det(i,k,jrow) = ta_rgraz_Det(i,k,jrow) & + rgraz_Det(i,k,j) ta_rmorp(i,k,jrow) = ta_rmorp(i,k,jrow) + rmorp(i,k,j) ta_rmorpt(i,k,jrow)= ta_rmorpt(i,k,jrow) + rmorpt(i,k,j) ta_rmorz(i,k,jrow) = ta_rmorz(i,k,jrow) + rmorz(i,k,j) ta_rexcr(i,k,jrow) = ta_rexcr(i,k,jrow) + rexcr(i,k,j) ta_rnpp_D(i,k,jrow) = ta_rnpp_D(i,k,jrow) & + rnpp_D(i,k,j) !{ wyao 20170713 ! wyao 20170713} ta_rgraz_D(i,k,jrow) = ta_rgraz_D(i,k,jrow) & + rgraz_D(i,k,j) ta_rmorp_D(i,k,jrow) = ta_rmorp_D(i,k,jrow) & + rmorp_D(i,k,j) ta_rnfix(i,k,jrow) = ta_rnfix(i,k,jrow) + rnfix(i,k,j) ! wyao 20160719{ ta_rfehydro(i,k,jrow)=ta_rfehydro(i,k,jrow) & +rfehydro(i,k,j) ! wyao 20160719} ta_rfeorgads(i,k,jrow) = ta_rfeorgads(i,k,jrow) & + rfeorgads(i,k,j) ta_rdeffe(i,k,jrow) = ta_rdeffe(i,k,jrow) & + rdeffe(i,k,j) ta_rfeprime(i,k,jrow) = ta_rfeprime(i,k,jrow) & + rfeprime(i,k,j) ta_rfesed(i,k,jrow) = ta_rfesed(i,k,jrow) & + rfesed(i,k,j) ta_rbfe(i,k,jrow) = ta_rbfe(i,k,jrow) & + rbfe(i,k,j) ta_rfecol(i,k,jrow) = ta_rfecol(i,k,jrow) & + rfecol(i,k,j) ta_rchl(i,k,jrow) = ta_rchl(i,k,jrow) & + rchl(i,k,j) ta_rchl_D(i,k,jrow) = ta_rchl_D(i,k,jrow) & + rchl_D(i,k,j) endif ! calculate total export to get total import for next layer expo = expo*dzt(k) expofe = expofe*dzt(k) enddo kmx = kmt(i,jrow) do k=1,kmx ! limit oxygen consumption below concentrations of ! 5umol/kg as recommended in OCMIP fo2 = 0.5*tanh(t(i,k,j,io2,taum1)*1000. - 5.) ! sink of oxygen ! O2 is needed to generate the equivalent of NO3 from N2 during N2 fixation ! 0.5 H2O + 0.5 N2+1.25O2 -> HNO3 ! note that so2 is -dO2/dt so2 = src(i,k,j,ispo4)*redotp & + rnfix(i,k,j)*1.25e-3 src(i,k,j,iso2) = -so2*(0.5 + fo2) ! add denitrification as source term for NO3 no3flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-trcmin) ! 800 = 0.8*1000 = (elec/mol O2)/(elec/mol NO3)*(mmol/mol) deni = max(0., 800.*no3flag*so2*(0.5 - fo2)) src(i,k,j,isno3) = src(i,k,j,isno3) - deni ! Correct the ALK stoichiometry to account for N2 fixation ! New stoichiometric model parameters formulated as in Paulmier et al. 2009 BG Paulmier_a = 1.e3*redctn Paulmier_z = 4*1.e3*redotp - 4.*Paulmier_a - 8.*redntp Paulmier_R0 = Paulmier_a + 0.25*Paulmier_z src(i,k,j,isalk) = src(i,k,j,isalk) & + no3flag*src(i,k,j,ispo4)*(0.5-fo2) & *(4./5.*Paulmier_R0 + (3./5. +1.)*redntp) * 1.e-3 ! Now account for N2 fixation (ALK production is tied to PO4 change and ! thus in the case of N2 fixation was not correct without this fix). src(i,k,j,isalk) = src(i,k,j,isalk) & - rnfix(i,k,j)*1.e-3 rdeni(i,k,jrow) = deni enddo !----------------------------------------------------------------------- ! remineralize calcite !----------------------------------------------------------------------- kmx = kmt(i,jrow) do k=1,kmx-1 src(i,k,j,isdic) = src(i,k,j,isdic) + prca*rcak(k) src(i,k,j,isalk) = src(i,k,j,isalk) + 2.*prca*rcak(k) enddo src(i,kmx,j,isdic) = src(i,kmx,j,isdic) + prca*rcab(kmx) src(i,kmx,j,isalk) = src(i,kmx,j,isalk) + 2.*prca*rcab(kmx) !----------------------------------------------------------------------- ! accumulate time averages for full depth variables !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then kmx = kmt(i,jrow) expo = prca ta_rprocal(i,jrow) = ta_rprocal(i,jrow) + prca do k=1,kmx expo = expo*dztr(k) ta_rremi(i,k,jrow) = ta_rremi(i,k,jrow) + rremi(i,k,j) ta_rexpo(i,k,jrow) = ta_rexpo(i,k,jrow) + rexpo(i,k,j) expo = expo - prca*rcak(k) ta_rexpocal(i,k,jrow) = ta_rexpocal(i,k,jrow) + expo expo = expo*dzt(k) ta_rdeni(i,k,jrow) = ta_rdeni(i,k,jrow) + rdeni(i,k,j) ta_rremife(i,k,jrow) = ta_rremife(i,k,jrow) & + rremife(i,k,j) ta_rexpofe(i,k,jrow) = ta_rexpofe(i,k,jrow) & + rexpofe(i,k,j) enddo endif endif enddo enddo return end