subroutine set (index, num, name, text, inc) !----------------------------------------------------------------------- ! increment counter, set index and text !----------------------------------------------------------------------- character(*) :: name, text name = text index = num inc = index + 1 return end subroutine sbc_init !----------------------------------------------------------------------- ! Initialize S.B.C. indices !----------------------------------------------------------------------- implicit none integer m !======================= 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 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +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 "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 & ) 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 sbc(:,:,:) = 0.0 mapsbc(:) = " " itaux = 0 itauy = 0 iws = 0 iaca = 0 isca = 0 ihflx = 0 isflx = 0 isst = 0 isss = 0 iro = 0 iwa = 0 iwxq = 0 iwyq = 0 iwxt = 0 iwyt = 0 iwxc = 0 iwyc = 0 ipsw = 0 isu = 0 isv = 0 igu = 0 igv = 0 issdic = 0 idicflx = 0 issalk = 0 ialkflx = 0 isso2 = 0 io2flx = 0 isspo4 = 0 !{ wyao 20170714 ! wyao 20170714} ipo4flx = 0 issphyt = 0 iphytflx = 0 isszoop = 0 izoopflx = 0 issdetr = 0 idetrflx = 0 issno3 = 0 ino3flx = 0 issdiaz = 0 idiazflx = 0 issc14 = 0 ic14flx = 0 isscfc11 = 0 icfc11flx = 0 isscfc12 = 0 icfc12flx = 0 iat = 0 irh = 0 ipr = 0 ips = 0 iaws = 0 iswr = 0 ilwr = 0 isens = 0 ievap = 0 idtr = 0 isr = 0 inpp = 0 iburn = 0 ibtemp = 0 ibsalt = 0 ibdic = 0 ibdicfx = 0 ibalk = 0 ibalkfx = 0 ibo2 = 0 ircal = 0 irorg = 0 ibtemp = 0 ibsalt = 0 ibo2 = 0 ibalk = 0 ibdic = 0 ibdicfx = 0 ibalkfx = 0 m = 1 call set (itaux, m, mapsbc(m), 'taux', m) call set (itauy, m, mapsbc(m), 'tauy', m) call set (iws, m, mapsbc(m), 'ws', m) call set (iaca, m, mapsbc(m), 'a_calb', m) call set (isca, m, mapsbc(m), 's_calb', m) call set (ihflx, m, mapsbc(m), 'hflx', m) call set (isflx, m, mapsbc(m), 'sflx', m) call set (isst, m, mapsbc(m), 'sst', m) call set (isss, m, mapsbc(m), 'sss', m) call set (iro, m, mapsbc(m), 'ro', m) call set (iwxq, m, mapsbc(m), 'wx_q', m) call set (iwyq, m, mapsbc(m), 'wy_q', m) call set (iwxt, m, mapsbc(m), 'wx_t', m) call set (iwyt, m, mapsbc(m), 'wy_t', m) call set (isu, m, mapsbc(m), 'su', m) call set (isv, m, mapsbc(m), 'sv', m) call set (igu, m, mapsbc(m), 'gu', m) call set (igv, m, mapsbc(m), 'gv', m) call set (issdic, m, mapsbc(m), 'ssdic', m) call set (idicflx, m, mapsbc(m), 'dicflx', m) call set (issalk, m, mapsbc(m), 'ssalk', m) call set (ialkflx, m, mapsbc(m), 'alkflx', m) call set (isso2, m, mapsbc(m), 'sso2', m) call set (io2flx, m, mapsbc(m), 'o2flx', m) call set (isspo4, m, mapsbc(m), 'sspo4', m) call set (ipo4flx, m, mapsbc(m), 'po4flx', m) !{ wyao 20170714 ! wyao 20170714} call set (issno3, m, mapsbc(m), 'ssno3', m) call set (ino3flx, m, mapsbc(m), 'no3flx', m) !{ wyao 20170714 ! wyao 20170714} if ( m-1 .gt. numsbc) then print*, '==> Error: increase numsbc in csbc.h to ', m-1 stop '=>UVic_ESCM' endif return end subroutine tracer_init !----------------------------------------------------------------------- ! Initialize ocean tracer names !----------------------------------------------------------------------- implicit none integer m !======================= 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 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +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 "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 "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 ! 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 ! { 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 ! { 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) !----------------------------------------------------------------------- ! Initialize ocean tracer names !----------------------------------------------------------------------- mapt(:) = " " itemp = 0 isalt = 0 idic = 0 ic14 = 0 icfc11 = 0 icfc12 = 0 ialk = 0 io2 = 0 ipo4 = 0 ino3 = 0 iphyt = 0 izoop = 0 idetr = 0 idiaz = 0 !{ wyao 20170714 ! wyao 20170714} m = 1 call set (itemp, m, mapt(m), 'temp', m) call set (isalt, m, mapt(m), 'salt', m) call set (idic, m, mapt(m), 'dic', m) call set (ialk, m, mapt(m), 'alk', m) call set (io2, m, mapt(m), 'o2', m) call set (ipo4, m, mapt(m), 'po4', m) !{ wyao 20170714 ! wyao 20170714} call set (iphyt, m, mapt(m), 'phyt', m) call set (izoop, m, mapt(m), 'zoop', m) call set (idetr, m, mapt(m), 'detr', m) call set (ino3, m, mapt(m), 'no3', m) !{ wyao 20170714 ! wyao 20170714} call set (idiaz, m, mapt(m), 'diaz', m) if ( m-1 .gt. nt) then print*, '==> Error: increase nt for tracers in size.h' stop '=>UVic_ESCM' endif !----------------------------------------------------------------------- ! Initialize ocean tracer sourcr names, must have equivalent tracer !----------------------------------------------------------------------- mapst(:) = " " itrc(:) = 0 m = 1 call set (isdic, m, mapst(m), 'sdic', m) itrc(idic) = m-1 call set (isalk, m, mapst(m), 'salk', m) itrc(ialk) = m-1 call set (iso2, m, mapst(m), 'so2', m) itrc(io2) = m-1 call set (ispo4, m, mapst(m), 'spo4', m) itrc(ipo4) = m-1 !{ wyao 20170714 ! wyao 20170714} call set (isphyt, m, mapst(m), 'sphyt', m) itrc(iphyt) = m-1 call set (iszoop, m, mapst(m), 'szoop', m) itrc(izoop) = m-1 call set (isdetr, m, mapst(m), 'sdetr', m) itrc(idetr) = m-1 call set (isno3, m, mapst(m), 'sno3', m) itrc(ino3) = m-1 !{ wyao 20170714 ! wyao 20170714} call set (isdiaz, m, mapst(m), 'sdiaz', m) itrc(idiaz) = m-1 if ( m-1 .gt. nt) then print*, '==> Error: increase nsrc for tracer sources in size.h' stop '=>UVic_ESCM' endif !----------------------------------------------------------------------- ! Initialize atmosphere tracer names !----------------------------------------------------------------------- mapat(:) = " " isat = 0 ishum = 0 ico2 = 0 m = 1 call set (isat, m, mapat(m), 'sat', m) call set (ishum, m, mapat(m), 'shum', m) if ( m-1 .gt. nat) then print*, '==> Error: increase nat in size.h' stop '=>UVic_ESCM' endif return end subroutine read_namelist !----------------------------------------------------------------------- ! read all model namelist variables !----------------------------------------------------------------------- implicit none character (120) :: fname, new_file_name, logfile integer iotraj, i, j, k, ip, kr, jq, n, ioun, num_processors logical exists, annlevobc, annlev, initpt real dtatms, snapint, snapls, snaple, snapde, ahbkg, runstep real slmx, s_dm, afkph, dfkph, sfkph, zfkph, ahs, ahb, trajint real crops_yr integer isot1, ieot1, isot2, ieot2, jsot, jeot, ksot, keot, mrot !======================= 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 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +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 "calendar.h" ====================== ! calendar specification arrays !----------------------------------------------------------------------- ! eqyear = true to select a calendar in which each year ! has the same number of days (i.e., no leap years) ! false selects a Julian calendar ! eqmon = true to force all months to have the same number of days ! false => the usual 31, 28, 31, 30, ..., days per month. ! only used when eqyear = true ! dayname = character names of days ! monname = character names of months ! monlen = the length of each month (in days) when eqmon is true ! yrlen = the length of a typical (non-leap) year in days ! daypm = array of month lengths in days (non-leap) ! msum = array of cumulative days preceding each month ! (again, non-leap) ! daylen = the length of a day in seconds !----------------------------------------------------------------------- character(10) :: dayname character(12) :: monname common /calen_c/ dayname(7), monname(12) logical eqyear, eqmon common /calen_l/ eqyear, eqmon integer daypm, msum, yrlen, monlen common /calen_i/ daypm(12), msum(12), yrlen, monlen real daylen common /calen_r/ daylen !======================= include file "cembm.h" ======================== ! parameters for use in the energy balance model (also see atm.h) ! namix = time steps between mixing (set in atmos.in) ! lf = time step flag (1=>leapfrog, 2=>forward) ! niats = number of ice advection sub time steps ! nivts = time steps between recalculating ice velocities ! nivc = time step counter for nivts ! ns = number of subcycles for advection and diffusion ! pyear = default paleo calendar year (-/+ = BC/AD) ! dts = time step (2*dtatm=>leapfrog, dtatm=>forward) ! co2ccn = atmospheric CO2 concentration (ppmv) ! co2emit = atmospheric CO2 emissions flux (g cm-2 s-1) ! co2emit_fuel = emissions flux from fossil fuels (g cm-2 s-1) ! co2emit_land = emissions flux from land use change (g cm-2 s-1) ! anthro = radiative forcing by atmospheric CO2 ! co2for = atmospheric CO2 forcing term (in units of heat flux) ! c14ccn = atmospheric C14 concentration (ppmv) ! dc14ccn = atmospheric dC14 concentration (permil) ! dc14ccnn = northern hemisphere atmospheric dC14 concentration ! dc14ccne = equatorial atmospheric dC14 concentration ! dc14ccns = southern hemisphere atmospheric dC14 concentration ! c14prod = atmospheric C14 production (umol cm-2 s-1) ! cfc11ccnn = northern hemisphere atmospheric CFC11 concentration ! cfc11ccns = southern hemisphere atmospheric CFC11 concentration ! cfc12ccnn = northern hemisphere atmospheric CFC12 concentration ! cfc12ccns = southern hemisphere atmospheric CFC12 concentration ! scatter = proportion of solar scattered by the atmosphere ! solarconst = solar constant (g s-3) ! cssh = constant used in calculation of ssh (g g-1) ! cdatm = drag coefficient (dimensionless) ! cpatm = atmospheric heat capacity (cm**2 s-2 K-1) ! sht = scale height for temperature ! shq = scale height for specific humidity ! shc = scale height for carbon ! rhoatm = density of air at sea surface (g cm-3) ! esatm = atmosphere emissivity times Stefan's constant ! rhoocn = representative sea surface density ! esocn = ocean emissivity times Stefan's constant ! vlocn = latent heat of vapourization of water ! cdice = drag coefficient (dimensionless) ! dampice = time scale for freezing first layer under ice (days) ! rhoice = ice density (g cm-3) ! rhosno = snow density (g cm-3) ! esice = ice emissivity times Stefan's constant ! slice = latent heat of sublimation of ice ! flice = latent heat of fusion of ice (cm2 s-2) ! condice = ice conductivity (g*cm s-3 K-1) ! tsno = air temperature for accumulating snow ! hsno_max = maximum snow depth ! totaltime = total time for long term averages ! rlapse = lapse rate (K cm-1) ! soilmax = soil water field capacity (cm) ! eslnd = land emissivity time Stefan's constant ! pass = atmospheric transmission coefficient (%/100) ! ice_calb = ice coalbedo (%/100) ! sno_calb = snow coalbedo (%/100) ! pcfactor = precip - cloud correlation factor ( %/100) ! rf1 = factor used in calculating lapse rate reduction ! rf2 = factor used in calculating lapse rate reduction ! co2_yr = co2 forcing year (-/+ = BC/AD) ! agric_yr = agricutural forcing year (-/+ = BC/AD) ! landice_yr = ice sheet forcing year (-/+ = BC/AD) ! solar_yr = solar forcing year (-/+ = BC/AD) ! orbit_yr = orbital forcing year (-/+ = BC/AD) ! volcano_yr = volcanic forcing year (-/+ = BC/AD) ! sulph_yr = sulphate forcing year (-/+ = BC/AD) ! aggfor_yr = additional greenhouse gas forcing year (-/+ = BC/AD) ! cfcs_yr = cfc forcing year (-/+ = BC/AD) ! c14_yr = c14 forcing year (-/+ = BC/AD) ! ice_yr = redundant (same as landice_yr, kept for compatibility) ! dalt_v = dalton number over vegetation (no vegetation model) ! dalt_o = dalton number over ocean ! dalt_i = dalton number over ice ! rhmax = maximum relative humidity ! volcfor = anomalous volcanic forcing (g/s**3) ! aggfor = anomalous additional greenhouse gas forcing (g s-3) ! aggfor_os = additional greenhouse gas forcing offset (g s-3) ! atmsa = atmospheric surface area (cm 2) ! ocnsa = exposed ocean surface area (cm 2) ! sealev = anomalous sea level change (cm) ! dsealev = anomalous sea level change in segtim (cm) ! sealev_yr = sea level forcing year (-/+ = BC/AD) ! itrack_co2 = index for averaging co2 array ! ntrack_co2 = number of points for average co2 ! itrack_sat = index for averaging sat array ! ntrack_sat = number of points for average sat ! vcsref = climate sensitivity reference temperature (C) ! vcsfac = climate sensitivity factor (mW/m2/C) ! vcsyri = year to start changing the climate sensitivity (year) ! gtoppm = conversion from g cm-2 to ppmv ! carbemit = accumulated co2 emissions (Pg) ! adiff = anomolous diffusion factor (percent/100/C) ! dtbar = global average sat anomoly (C) ! cropf = scaler for crop area ! pastf = scaler for pasture area integer namix, lf, niats, nivts, nivc, ns, itrack_co2 integer ntrack_co2, itrack_sat, ntrack_sat common /cembm_i/ namix, lf, niats, nivts, nivc, ns, itrack_co2 common /cembm_i/ ntrack_co2, itrack_sat, ntrack_sat real pyear, dts, co2ccn, co2emit, co2emit_fuel, co2emit_land real anthro, co2for, c14ccn, dc14ccn, dc14ccnn, dc14ccne, dc14ccns real c14prod, cfc11ccnn, cfc11ccns, cfc12ccnn, cfc12ccns, scatter real solarconst, cssh, cdatm, cpatm, sht, shq, shc, rhoatm, esatm real rhoocn, esocn, vlocn, cdice, dampice, rhoice, rhosno, esice real slice, flice, condice, tsno, hsno_max, totaltime, rlapse real soilmax, eslnd, pass, ice_calb, sno_calb, pcfactor, rf1, rf2 real co2_yr, agric_yr, landice_yr, solar_yr, orbit_yr, volcano_yr real sulph_yr, aggfor_yr, cfcs_yr, c14_yr, ice_yr, dalt_v, dalt_o real dalt_i, rhmax, volcfor, aggfor, aggfor_os, atmsa, ocnsa real sealev, dsealev, sealev_yr, vcsref, vcsfac, vcsyri, gtoppm real carbemit, adiff, dtbar, cropf, pastf, anthro_srm common /cembm_r/ pyear, dts, co2ccn, co2emit, co2emit_fuel common /cembm_r/ co2emit_land, anthro, co2for, c14ccn, dc14ccn common /cembm_r/ dc14ccnn, dc14ccne, dc14ccns, c14prod, cfc11ccnn common /cembm_r/ cfc11ccns, cfc12ccnn, cfc12ccns, cssh, cdatm common /cembm_r/ cpatm, sht, shq, shc, rhoatm, esatm, rhoocn common /cembm_r/ scatter, solarconst, esocn, dampice, rhoice common /cembm_r/ rhosno, esice, slice, flice, condice, vlocn common /cembm_r/ cdice, tsno, hsno_max, totaltime, rlapse, soilmax common /cembm_r/ eslnd, pass, ice_calb, sno_calb, pcfactor, rf1 common /cembm_r/ rf2, co2_yr, agric_yr, landice_yr, solar_yr common /cembm_r/ orbit_yr, volcano_yr, sulph_yr, aggfor_yr common /cembm_r/ cfcs_yr, c14_yr, ice_yr, dalt_v, dalt_o, dalt_i common /cembm_r/ rhmax, volcfor, aggfor, aggfor_os, atmsa, ocnsa common /cembm_r/ sealev, dsealev, sealev_yr, vcsref, vcsfac common /cembm_r/ vcsyri, gtoppm, carbemit, adiff, dtbar, cropf common /cembm_r/ pastf, anthro_srm ! ntatsa = time step counter for time averaging ! ntatia = number of time averaged time step integrals ! tai_sat = average integrated elevated surface air temperature ! tai_shum = average integrated surface specific humidity ! tai_precip = average integrated precipitation ! tai_evap = average integrated evaporation ! tai_ohice = total integrated sea ice volume ! tai_oaice = total integrated sea ice area ! tai_hsno = total integrated snow volume ! tai_lhice = total integrated land ice volume ! tai_laice = total integrated land ice area ! tai_co2ccn = average integrated CO2 concentration ! tai_co2emit = average integrated CO2 emissions ! tai_dc14ccn = average integrated dC14 concentration ! tai_cfc11ccn = average integrated CFC11 concentration ! tai_cfc12ccn = average integrated CFC12 concentration ! tai_maxit = average maximum iterations for atmospheric solver ! tai_nsat = average northern hemisphere surface air temperature ! tai_ssat = average southern hemisphere surface air temperature ! tai_nshum = average northern hemisphere surface specific humidity ! tai_sshum = average southern hemisphere surface specific humidity ! tai_nprecip = average northern hemisphere precipitation ! tai_sprecip = average southern hemisphere precipitation ! tai_nevap = average northern hemisphere evaporation ! tai_sevap = average southern hemisphere evaporation ! tai_nohice = total northern hemisphere sea ice volume ! tai_sohice = total southern hemisphere sea ice volume ! tai_noaice = total northern hemisphere sea ice area ! tai_soaice = total southern hemisphere sea ice area ! tai_nhsno = total northern hemisphere snow volume ! tai_shsno = total southern hemisphere snow volume ! tai_nlhice = total northern hemisphere land ice volume ! tai_slhice = total southern hemisphere land ice volume ! tai_nlaice = total northern hemisphere land ice area ! tai_slaice = total southern hemisphere land ice area ! tai_lsat = average sat over land ! tai_osat = average sat over ocean ! tai_lprecip = average precip over land ! tai_oprecip = average precip over ocean ! tai_levap = average evap over land ! tai_oevap = average evap over ocean ! tai_solins = average incoming solar ! tai_upsens = average surface sensible heat ! tai_uplwr = average surface upward longwave radiation ! tai_outlwr = average outgoing longwave radiation ! tai_dnswr = average downward (absorbed) shortwave ! tai_absswr = net absorbed shortwave radiation ! tai_netrad = net radiation at the top of the atmosphere ! tai_palb = average planetary albedo ! tai_aalb = average atmospheric albedo ! tai_salb = average surface albedo ! tai_lsalb = average surface albedo over land ! tai_osalb = average surface albedo over ocean ! tai_sst = average sea surface temperature ! tai_sss = average sea surface salinity ! tai_ssdic = average sea surface dic ! tai_ssc14 = average sea surface c14 ! tai_ssalk = average sea surface alkalinity ! tai_sso2 = average sea surface o2 ! tai_sspo4 = average sea surface po4 ! tai_ssno3 = average sea surface no3 !{ wyao 20170703 ! wyao 20170703} ! tai_ssdfe = average sea surface dfe ! tai_sscfc11 = average sea surface cfc11 ! tai_sscfc12 = average sea surface cfc12 ! tai_sulph = average tropospheric sulphate forcing ! tai_volc = average volcanic forcing ! tai_agg = average additional greenhouse gas forcing ! tai_catm = average total carbon in atmosphere ! tai_carbemit = average taccumulated co2 emissions integer ntatsa, ntatia common /cembm_i/ ntatsa, ntatia real tai_sat, tai_shum, tai_precip, tai_evap, tai_ohice, tai_oaice real tai_hsno, tai_lhice, tai_laice, tai_co2ccn, tai_co2emit real tai_dc14ccn, tai_cfc11ccn, tai_cfc12ccn, tai_maxit, tai_nsat real tai_ssat, tai_nshum, tai_sshum, tai_nprecip, tai_sprecip real tai_nevap, tai_sevap, tai_nohice, tai_sohice, tai_noaice real tai_soaice, tai_nhsno, tai_shsno, tai_nlhice, tai_slhice real tai_nlaice, tai_slaice, tai_lsat, tai_osat, tai_lprecip real tai_oprecip, tai_levap, tai_oevap, tai_solins, tai_upsens real tai_uplwr, tai_outlwr, tai_dnswr, tai_absswr, tai_netrad real tai_palb, tai_aalb, tai_salb, tai_lsalb, tai_osalb, tai_sst real tai_sss, tai_ssdic, tai_ssc14, tai_ssalk, tai_sso2, tai_sspo4 real tai_ssno3, tai_sscfc11, tai_sscfc12, tai_sulph, tai_volc real tai_agg, tai_catm, tai_ssdfe, tai_carbemit !{ wyao 20170703 ! wyao 20170703} common /cembm_r/ tai_sat, tai_shum, tai_precip, tai_evap common /cembm_r/ tai_ohice, tai_oaice, tai_hsno, tai_lhice common /cembm_r/ tai_laice, tai_co2ccn, tai_co2emit, tai_dc14ccn common /cembm_r/ tai_cfc11ccn, tai_cfc12ccn, tai_maxit, tai_nsat common /cembm_r/ tai_ssat, tai_nshum, tai_sshum, tai_nprecip common /cembm_r/ tai_sprecip, tai_nevap, tai_sevap, tai_nohice common /cembm_r/ tai_sohice, tai_noaice, tai_soaice, tai_nhsno common /cembm_r/ tai_shsno, tai_nlhice, tai_slhice, tai_nlaice common /cembm_r/ tai_slaice, tai_lsat, tai_osat, tai_lprecip common /cembm_r/ tai_oprecip, tai_levap, tai_oevap, tai_solins common /cembm_r/ tai_upsens, tai_uplwr, tai_outlwr, tai_dnswr common /cembm_r/ tai_absswr, tai_netrad, tai_palb, tai_aalb common /cembm_r/ tai_salb, tai_lsalb, tai_osalb, tai_sst common /cembm_r/ tai_sss, tai_ssdic, tai_ssc14, tai_ssalk common /cembm_r/ tai_sso2, tai_sspo4, tai_ssno3, tai_sscfc11 common /cembm_r/ tai_sscfc12, tai_sulph, tai_volc, tai_agg common /cembm_r/ tai_catm, tai_ssdfe, tai_carbemit !{ wyao 20170703 ! wyao 20170703} !======================= 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 "cnep.h" =========================== ! the option "neptune" provides a subgridscale parameterization ! for the interaction of eddies and topography ! reference: ! Holloway, G., 1992: representing topographic stress for large ! scale ocean models, J. Phys. Oceanogr., 22, 1033-1046 ! neptune is calculated as an equilibrium streamfunction given by ! pnep=-f*snep*snep*hnep and is applied through eddy viscosity ! hnep = model streamfunction depth ! snep = spnep + (senep-spnep)*(0.5 + 0.5*cos(2.0*latitude)) ! the neptune length scale snep has a value of senep at the ! equator and smoothly changes to spnep at the poles ! variables used in applying neptune ! spnep = neptune length scale at the pole ! senep = neptune length scale at the equator ! unep = neptune velocity ! pnep = neptune streamfunction real spnep, senep, pnep, unep common /cnep_r/ spnep, senep, pnep(imt,jmt), unep(imt,jmt,2) !======================= 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 & ) 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 "cprnts.h" ======================= ! variables used for controlling matrix printouts during diagnostic ! timesteps: ! nlatpr = maximum number of latitudes for matrix printouts ! on diagnostic time steps ! prlat = latitudes (deg) at which (x,z) printouts are desired ! start & end coordinates for matrix printouts of (x,z) sections ! prslon = starting longitudes (deg) ! prelon = ending longitudes (deg) ! prsdpt = starting depths (cm) ! predpt = ending depths (cm) ! start & end coordinates for matrix printouts of (x,y) sections ! slonxy = starting longitude (deg) ! elonxy = ending longitude (deg) ! slatxy = starting latitude (deg) ! elatxy = ending latitude (deg) ! matrix printouts of (y,z) sections will use above coordinates real prlat, prslon, prelon, prsdpt, predpt, slatxy, elatxy real slonxy, elonxy common /cprnts/ prlat(nlatpr), prslon(nlatpr), prelon(nlatpr) common /cprnts/ prsdpt(nlatpr), predpt(nlatpr) common /cprnts/ slatxy, elatxy, slonxy, elonxy !====================== include file "diag.h" ========================== ! variables used for computing diagnostics: ! engint = volume averaged internal mode energy integral ! components ! engext = volume averaged external mode energy integral ! components ! buoy = volume averaged buoyancy ! tcerr = maximum "t" cell continuity error ! ucerr = maximum "u" cell continuity error ! itcerr = "i" index corresponding to "tcerr" ! jtcerr = "jrow" index corresponding to "tcerr" ! ktcerr = "k" index corresponding to "tcerr" ! iucerr = "i" index corresponding to "ucerr" ! jucerr = "jrow" index corresponding to "ucerr" ! kucerr = "k" index corresponding to "ucerr" ! wtbot = maximum "adv_vbt" error at ocean bottom ! iwtbot = "i" index corresponding to "wtbot" ! jwtbot = "jrow" index corresponding to "wtbot" ! kwtbot = "k" index corresponding to "wtbot" ! wubot = maximum "adv_vbu" at ocean bottom ! iwubot = "i" index corresponding to "wubot" ! jwubot = "jrow" index corresponding to "wubot" ! kwubot = "k" index corresponding to "wubot" ! wtlev = zonally integrated adv_vbt for each level ! wulev = zonally integrated adv_vbu for each level integer itcerr, jtcerr, ktcerr, iucerr, jucerr, kucerr integer iwtbot, jwtbot, kwtbot, iwubot, jwubot, kwubot common /cdiag_i/ itcerr(jmt), jtcerr(jmt), ktcerr(jmt) common /cdiag_i/ iucerr(jmt), jucerr(jmt), kucerr(jmt) common /cdiag_i/ iwtbot(jmt), jwtbot(jmt), kwtbot(jmt) common /cdiag_i/ iwubot(jmt), jwubot(jmt), kwubot(jmt) real buoy, engint, engext, tcerr, ucerr, wtbot, wubot, wtlev real wulev common /cdiag_r/ buoy(0:km,jmt), engint(0:km,8,jmt), engext(8,jmt) common /cdiag_r/ tcerr(jmt), ucerr(jmt) common /cdiag_r/ wtbot(jmt), wubot(jmt) common /cdiag_r/ wtlev(km,0:jmt), wulev(km,0:jmt) ! vmsf = vertical_meridional stream function real vmsf common /cdiag_r/ vmsf(jmt,km) ! term balances are instantaneous breakdowns of all terms in the ! momentum & tracer equations. They are averaged over ocean volumes ! defined by horizontal and vertical regional masks: ! termbt = term balance components for time rate of change of ! tracers within a volume. the total time rate of change ! is broken down into components as follows: ! the form is d( )/dt = terms (2) ... (10) where each ! term has the units of "tracer units/sec" using ! schematic terms for illustration. ! (1) = total time rate of change for the tracer ! (2) = change due to zonal nonlinear term: (UT)x ! (3) = change due to meridional nonlinear term: (VT)y ! (4) = change due to vertical nonlinear term: (WT)z ! (5) = change due to zonal diffusion: Ah*Txx ! (6) = change due to meridional diffusion: Ah*Tyy ! (7) = change due to vertical diffusion: kappa_h*Tzz ! (8) = change due to source term ! (9) = change due to explicit convection ! (10) = change due to filtering ! the nonlinear terms can be broken into two parts: advection and a ! continuity part: The physically meaningful part is advection. ! eg: Zonal advection of tracer "A" is -U(A)x = A(U)x - (UA)x ! (11) = zonal advection U(Ax) ! (12) = meridional advection V(Ay) ! (13) = vertical advection W(Az) ! (14) = change of tracer variance (tracer**2 units) ! (15) = average tracer within volume (tracer units) ! terr = error term = (1) - sum (2) ... (10) ! asst = average sea surface tracer for regional surface areas ! stflx = average surface tracer flux for regional surface areas ! tracer (#1,#2) units = (cal/cm**2/sec, gm/cm**2/sec) ! termbm = term balance components for time rate of change of ! momentum within a volume. the total time rate of change ! is broken down into components as follows: ! the form is d( )/dt = terms (2) ... (13) where each ! term has the units of "cm/sec**2" and "Q" is the ! momentum component {zonal or meridional} using ! schematic terms for illustration. ! (1) = total time rate of change for the momentum ! (2) = change due to the pressure gradient: grad_p ! without the surface pressure gradients ! (i.e., for computing the internal modes) ! (3) = change due to zonal nonlinear term: (UQ)x ! (4) = change due to meridional nonlinear term: (VQ)y ! (5) = change due to vertical nonlinear term: (wQ)z ! (6) = change due to zonal viscosity: Am*Qxx ! (7) = change due to meridional viscosity: Am*Qyy ! (8) = change due to vertical viscosity: kappa_m*Qzz ! (9) = change due to metric diffusion terms ! (10) = change due to coriolis terms: fQ ! (11) = change due to source terms ! (12) = change due to surface pressure gradient ! this is obtained after solving the external mode ! in the stream function technique. It is solved ! directly from the elliptic equation for the ! prognostic surface pressure technique ! (13) = change due to metric advection ! the nonlinear terms can be broken into two parts: advection and a ! continuity part: The physically meaningful part is advection. ! eg: Zonal advection of vel component "Q" is -U(Q)x = Q(U)x - (UQ)x ! (14) = zonal advection U(Qx) ! (15) = meridional advection V(Qy) ! (16) = vertical advection W(Qz) ! (17) = average velocity component ! smflx = average surface momentum flux for regional surf areas ! in dynes/cm**2 ! avgw = average vertical velocity (cm/sec) ! ustf = names & units for surface tracer fluxes integer ntterms, nuterms parameter (ntterms=15, nuterms=17) character(15) :: ustf(nt,2) common /termb_c/ ustf real asst, avgw, termbt, termbm, smflx, stflx, terr common /termb_r/ asst(nt,0:nhreg), avgw(numreg) common /termb_r/ termbt(0:km,ntterms,nt,0:numreg) common /termb_r/ termbm(0:km,nuterms,2,numreg) common /termb_r/ smflx(2,0:nhreg), stflx(nt,0:nhreg), terr(nt) !======================= 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 "fwa.h" ========================= ! parameters for adding anomalous fresh water pulses ! isfwa1 = starting i index of section 1 for freshwater ! iefwa1 = ending i index of section 1 for freshwater ! isfwa2 = starting i index of section 2 for freshwater ! iefwa2 = ending i index of section 2 for freshwater ! jsfwa = starting j index for freshwater ! jefwa = ending j index for freshwater ! mrfwa = regional mask region for freshwater ! fwaflxi = initial fresh water flux (Sv) ! fwayri = year to start fresh water flux ! fwayrf = year to end fresh water flux ! fwarate = rate of increase in flux (Sv year-1) ! areafwa = area of fresh water anomaly ! areafwc = area of fresh water compensation ! compensate = flag to compensate for flux everywhere else ! fwawt = weighting for fresh water anomaly area ! fwcwt = weighting for fresh water compensation area integer isfwa1, iefwa1, isfwa2, iefwa2, jsfwa, jefwa, mrfwa common /fwa_i/ isfwa1, iefwa1, isfwa2, iefwa2, jsfwa, jefwa common /fwa_i/ mrfwa logical compensate common /fwa_l/ compensate real fwaflxi, fwayri, fwayrf, fwarate, areafwa, areafwc real fwawt, fwcwt common /fwa_r/ fwaflxi, fwayri, fwayrf, fwarate, areafwa, areafwc common /fwa_r/ fwawt(imt,jmt), fwcwt(imt,jmt) !======================= 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 "insolation.h" ===================== ! parameters used in calculating solar insolation ! eccen = orbital eccentricity ! obliq = obliquity (degrees) ! mvelp = moving vernal equinox longitude of perihelion (degrees) ! lambm0 = Mean long of perihelion at vernal equinox (radians) ! sindec = sine of solar declination angle in rad ! eccf = Earth-sun distance factor (ie. (1/r)**2) real eccen, obliq , mvelp, lambm0, sindec, eccf common /insolation_r/ eccen, obliq , mvelp, lambm0, sindec, eccf !====================== include file "iounit.h" ======================== ! i/o units and related variables ! taum1disk = disk pointer for tau-1 latitude rows ! taudisk = disk pointer for tau latitude rows ! taup1disk = disk pointer for tau+1 latitude rows ! kflds = disk unit used for two dimensional fields ! latdisk = disk units for latitude rows (alternately pointed to ! by taum1disk, taudisk, and taup1disk) ! wide_open_mw = logical to indicate that the MW is fully opened. ! if .true. then jmw = jmt and there are no latitude rows ! on disk. Instead, they are all in the MW. ! if .false. then jmw < jmt and all latitude rows are on ! disk so they must be transferred between the MW and disk. ! iodoc = unit for documentation ! iostab = unit for stability testing ! iotim = unit for time means ! iotim1 = scratch disk (SSD) unit for accumulating time means ! for the following, a control # < 0 implies that unformatted data ! will be written to a unit selected by the i/o manager "iomngr.F" ! and given a hardwired name (grep getunit *.F to see names) ! and formatted data (to stdout) will be written. if a # > 0 and ! # <> stdout, only unformatted data will be written. ! iotavg = control # for tracer averages ! iotmb = control # for writing tracer meridional budget. ! iotrmb = control # for term balances for tracer and momentum ! ioglen = control # for writing global energetics integrals ! iovmsf = control # for writing meridional stream function ! iogyre = control # for writing gyre transport. ! ioprxz = control # for writing x-z sections from latitudes ! ioext = control # for writing external mode (stream function) ! iodsp = control # for writing diagnostic surface pressure ! iotsi = control # for writing time step integrals ! ioxbt = control # for writing time averaged xbt data ! iozmbc = control # for writing zonal mean surf boundary conditions ! iotext = 120 character text string for describing the details ! of the next unformatted data record. ! expnam = 60 character text string for the experiment name ! runstamp = 120 character text string for the run stamp ! when writing unformatted data records in MOM, each data record is ! proceeded by a header record which was written as: ! write(unit) stamp, iotext, expnam ! where stamp is a 32 character specification of the model date & ! time corresponding to the time step when the data was written and ! iotext is a 120 character description of what is in the ! data record and how it is to be read. expnam is a 60 character ! experiment name which shows which experiment wrote the data. ! this makes it easy to decipher any unformatted output from the ! model by using a program similar to the following: ! program decifr !----------------------------------------------------------------------- ! decipher an unformatted file from MOM by showing the header ! records. the file needs to copied to file "fort.21" !----------------------------------------------------------------------- ! character*32 stamp ! character*120 iotext ! character*60 expnam ! iounit = 21 ! rewind iounit ! do n=1,100000 ! read the header record ! read (iounit, end=110) stamp, iotext, expnam ! write (*,'(1x,a32,1x,a80)') stamp, iotext ! skip the data record ! read (iounit) ! enddo !110 continue ! write (*,*) " => end of file on fort.",iounit ! stop ! end ! note: all unformatted diagnostic MOM data is handled this way. ! to insure that data is read properly, verify that arrays are ! dimensioned correctly by comparing the listed variables against ! those in the *.h files. (grep -i -n "variable" *.h) Also, most ! data from MOM is written IEEE 32bit so it is read directly by ! most workstations. However, when trying to read these IEEE files ! on the CRAY, they must be assigned IEEE before being read. ! Some diagnostic data is averaged over time before being written. ! In these cases, the time "stamp" refers to the last time step ! at the end of the averaging period. An averaging interval is ! also written as part of the data. Averaging periods = zero ! indicate instantaneous data. character iotext*120, expnam*60, runstamp*120 common /iounit_c/ iotext, expnam, runstamp integer taum1disk, taudisk, taup1disk, latdisk, kflds integer iodoc, iostab, iotavg, iotmb, iotrmb, iotim, iotim1 integer ioglen, iovmsf, iogyre, ioprxz, ioext, iodsp integer iotsi, iozmbc, ioxbt common /iounit_i/ taum1disk, taudisk, taup1disk, latdisk(2), kflds common /iounit_i/ iodoc, iostab, iotavg, iotmb, iotrmb, iotim common /iounit_i/ iotim1, ioglen, iovmsf, iogyre, ioprxz, ioext common /iounit_i/ iodsp, iotsi, iozmbc, ioxbt logical wide_open_mw common /iounit_l/ wide_open_mw !======================== 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 "mtlm.h" ======================== !********************************************************************** ! this file is based on code that may have had the following copyright: ! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE. ! Permission has been granted by the authors to the public to copy ! and use this software without charge, provided that this Notice and ! any statement of authorship are reproduced on all copies. Neither the ! Crown nor the U.K. Meteorological Office makes any warranty, express ! or implied, or assumes any liability or responsibility for the use of ! this software. !********************************************************************** !----------------------------------------------------------------------- ! Atmospheric CO2 variables !----------------------------------------------------------------------- ! CO2 = CO2 mass mixing ratio (kg/kg). ! EPCO2 = Ratio of molecular weights of CO2 real EPCO2, CO2 parameter (EPCO2 = 1.5194) common /land_r/ CO2(POINTS) !----------------------------------------------------------------------- ! Driving variables !----------------------------------------------------------------------- ! LYING_SNOW = Snow mass (kg/m2). ! LW = Surface longwave radiation (W/m**2). ! SW = Surface shortwave radiation (W/m**2). ! SWN = Net shortwave radiation (W/m**2). ! PSTAR = Surface pressure (Pa). ! Q = Specific humidity (kg/kg). ! SAT_D = Diurnal surface air temperature (K). ! TS1 = Sub-surface temperature (K). ! TIMEDAY = Time of day (s). ! W = Water incident at the soil surface (mm/day). ! WIND = Wind speed (m/s). ! LW_OUT = total net longwave radiation (W/m2) ! LHC = Latent heat of condensation (J/kg). ! LHF = Latent heat of fusion (J/kg). ! SIGMA = Stefan-Boltzman constant (W/m2/K4). real LYING_SNOW, LW, SW, SWN, PSTAR, Q, SAT_D, TS1, TIMEDAY real W, WIND, LW_OUT, LHC, LHF, SIGMA common /land_r/ LYING_SNOW(POINTS), LW(POINTS), SW(POINTS) common /land_r/ SWN(POINTS), PSTAR(POINTS), Q(POINTS) common /land_r/ SAT_D(POINTS), TS1(POINTS), TIMEDAY common /land_r/ WIND(POINTS), LW_OUT(POINTS), LHC, LHF, SIGMA !----------------------------------------------------------------------- ! Variables for forcing with climatological means !----------------------------------------------------------------------- ! DTEMP_DAY = Diurnal temperature range (K). ! LYING_SNOW_C = Snow mass (kg/m2). ! QS = Saturated specific humidity (kg/kg). ! RAIN = Rainfall rate (kg/m2/s). ! SNOW = Snowfall rate (kg/m2/s). ! SNOWMELT = Snow melt (mm/day). ! SW_C = Surface shortwave radiation (W/m**2). ! SUN = Normalized solar radiation. ! SURF_ROFF = Surface runoff (mm/day). ! TIME_MAX = GMT at which maximum temperature occurs (s). ! T_C = Air temperature (K). ! W_C = Water incident at the soil surface (mm/day). ! RH_C = relative humidity from atmospheric model real DTEMP_DAY, QS, RAIN, SNOW, SNOWMELT, SW_C, SUN, SURF_ROFF real TIME_MAX, T_C, W_C, RH_C common /land_r/ DTEMP_DAY(POINTS), QS(POINTS), RAIN(POINTS) common /land_r/ SNOW(POINTS), SNOWMELT(POINTS), SW_C(POINTS) common /land_r/ SUN(POINTS,STEPSM), SURF_ROFF(POINTS) common /land_r/ TIME_MAX(POINTS), T_C(POINTS), W_C(POINTS) common /land_r/ RH_C(POINTS) !----------------------------------------------------------------------- ! Model parameters !----------------------------------------------------------------------- ! TIMESTEP = Timestep for daily calculations (s). real TIMESTEP common /land_r/ TIMESTEP !----------------------------------------------------------------------- ! Soil parameters !----------------------------------------------------------------------- ! ALBSOIL = Soil albedo. ! ALBSNOW = Snow albedo. ! ALBLAND = Surface albedo. ! Z0S = Roughness length for bare soil (m). ! ROOTDEP = Rootdepth (m). ! HCAP_SOIL = Soil heat capacity (W/m3/K). ! HCON_SOIL = Soil heat conductivity (W/m/K). real ALBSOIL, ALBSNOW, ALBLAND, Z0S, ROOTDEP, HCAP_SOIL, HCON_SOIL ! parameter (ROOTDEP=1.0, HCAP_SOIL=3.3E5, HCON_SOIL=0.23) parameter (ROOTDEP=1.0, HCAP_SOIL=3.3E5, HCON_SOIL=0.75) common /land_r/ ALBSOIL(POINTS), ALBSNOW(POINTS) common /land_r/ ALBLAND(POINTS), Z0S(POINTS) !----------------------------------------------------------------------- ! Vegetation parameters !----------------------------------------------------------------------- ! ALBSNC = Cold deep snow albedo. ! ALBSNF = Snow free albedo. ! CATCH = Canopy capacity (kg/m2). ! Z0 = Vegetative roughness length (m). ! VEG_FRAC = Vegetated fraction. ! FRAC_VS = Total fraction of gridbox covered by vegetation and soil. ! FRAC_MIN = Minimum areal fraction for PFTs. ! FRAC_SEED = "Seed" fraction for PFTs. real ALBSNC, ALBSNF, CATCH, Z0, VEG_FRAC, FRAC_VS, FRAC_MIN real FRAC_SEED parameter (FRAC_MIN=1.0E-6, FRAC_SEED=0.01) common /land_r/ ALBSNC(POINTS,NPFT), ALBSNF(POINTS,NPFT) common /land_r/ CATCH(POINTS,NPFT) common /land_r/ Z0(POINTS,NPFT), VEG_FRAC(POINTS) common /land_r/ FRAC_VS(POINTS) !----------------------------------------------------------------------- ! IVM Prognostics !----------------------------------------------------------------------- ! CS = Soil carbon (kg C/m2). ! FRAC = Areal coverage. ! LAI = Leaf area index. real CS, FRAC, LAI common /land_r/ CS(POINTS), FRAC(POINTS,NTYPE) common /land_r/ LAI(POINTS,NPFT) !----------------------------------------------------------------------- ! IVM Fluxes and diagnostics !----------------------------------------------------------------------- ! C_VEG = Vegetation carbon (kg C/m2). ! CV = Gridbox mean vegetation carbon (kg C/m2). ! FTIME = Weighting factor for accumulations. ! HT = Canopy height (m). ! GPP = Gross Primary Productivity (kg C/m2/s). ! G_LEAF = Leaf turnover rate (/yr). ! G_LEAF_DAY = Daily mean leaf turnover rate (/yr). ! G_LEAF_PHEN = Daily leaf turnover rate including phenology (/yr). ! NPP = Net Primary Productivity (kg C/m2/s). ! RESP_S = Soil respiration rate (kg C/m2/s). ! LIT_C_T = Gridbox mean carbon litter (kg C/m2/yr). ! BF = Burn fraction real C_VEG, CV, FTIME, HT, GPP, G_LEAF, G_LEAF_DAY real G_LEAF_PHEN, NPP, RESP_S, LIT_C_T, BF common /land_r/ C_VEG(POINTS,NPFT), CV(POINTS), FTIME common /land_r/ HT(POINTS,NPFT), GPP(POINTS,NPFT) common /land_r/ G_LEAF(POINTS,NPFT), G_LEAF_DAY(POINTS,NPFT) common /land_r/ G_LEAF_PHEN(POINTS,NPFT), NPP(POINTS,NPFT) common /land_r/ RESP_S(POINTS), LIT_C_T(POINTS), BF !----------------------------------------------------------------------- ! Driving variables for TRIFFID !----------------------------------------------------------------------- ! G_LEAF_DR = Accumulated leaf turnover rate (/yr). ! NPP_DR = Accumulated Net Primary Productivity (kg C/m2/yr). ! RESP_S_DR = Accumulated soil respiration rate (kg C/m2/yr). ! RESP_W_DR = Accumulated wood respiration rate (kg C/m2/yr). real G_LEAF_DR, NPP_DR, RESP_S_DR, RESP_W_DR common /land_r/ G_LEAF_DR(POINTS,NPFT), NPP_DR(POINTS,NPFT) common /land_r/ RESP_S_DR(POINTS), RESP_W_DR(POINTS,NPFT) !----------------------------------------------------------------------- ! Hydrology variables !----------------------------------------------------------------------- ! ESUB = Sublimation (kg/m2/s). ! FSMC = Moisture availability factor. ! MAF = Moisture availability factor PFT/soil dependence ! M = Total soil moisture (kg/m2). ! MNEG = Negative soil moisture (kg/m2). real ESUB, FSMC, MAF, M, MNEG common /land_r/ ESUB(POINTS), FSMC(POINTS), MAF(NTYPE), M(POINTS) common /land_r/ MNEG(POINTS) !----------------------------------------------------------------------- ! Temperatures !----------------------------------------------------------------------- ! TSOIL = Temperature of bare soil (K). ! TSTAR = Surface temperature (K). real TSOIL, TSTAR common /land_r/ TSOIL(POINTS), TSTAR(POINTS,NPFT) !----------------------------------------------------------------------- ! Inputs defining locations !----------------------------------------------------------------------- ! LAND_PTS = Number of land points. ! LAND_INDEX = Indices of land points. ! VEG_PTS = Number of land points which include the nth vegetation ! type. ! VEG_INDEX = Indices of land points which include the nth vegetation ! type. integer LAND_PTS, LAND_INDEX, VEG_PTS, VEG_INDEX common /land_i/ LAND_PTS, LAND_INDEX(POINTS), VEG_PTS(NPFT) common /land_i/ VEG_INDEX(POINTS,NPFT) ! LAT = Latitude (degrees) ! LATMIN,LATMAX = Latitudinal limits of the area (degrees). ! GAREA = grid area (m2) ! LONG = Longitude (degrees) ! LONGMIN,LONGMAX = Longitudinal limits of the area (degrees). real LAT, LATMIN, LATMAX, LONG, LONGMIN, LONGMAX, GAREA common /land_r/ LAT(POINTS),LATMIN, LATMAX, GAREA(POINTS) common /land_r/ LONG(POINTS), LONGMIN, LONGMAX !----------------------------------------------------------------------- ! Time parameters !----------------------------------------------------------------------- ! ISTEP = step counter. ! STEP_DAY = Number of steps in a day. ! LAND_COUNTER = number time steps for the model. integer ISTEP, STEP_DAY, LAND_COUNTER common /land_i/ ISTEP, STEP_DAY, LAND_COUNTER ! SEC_DAY = Number of seconds in a day (s). ! DAY_YEAR = Number of days in a year (days). ! SEC_YEAR = Number of seconds in a year (s). real DAY_YEAR, SEC_DAY, SEC_YEAR common /land_r/ DAY_YEAR, SEC_DAY, SEC_YEAR ! INT_VEG = .T. for interactive vegetation ! VEG_EQUIL = .T. if the vegetation equilibrium logical INT_VEG, VEG_EQUIL common /land_l/ INT_VEG, VEG_EQUIL !----------------------------------------------------------------------- ! Variables defining anthropogenic disturbance !----------------------------------------------------------------------- ! FRACA = Areal fraction of agriculture. ! G_ANTH = Anthropogenic disturbance rate (/yr). real FRACA, G_ANTH common /land_r/ FRACA(POINTS), G_ANTH(POINTS) !----------------------------------------------------------------------- ! New variables required for LAND module !----------------------------------------------------------------------- ! ET = Evapotranspiration (kg/m2/s). ! LE = Latent heat flux (W/m2). ! TSTAR_GB = Grid box average surface temperature (K). ! SH = Sensible heat flux (W/m2). real ET, LE, TSTAR_GB, SH common /land_r/ ET(POINTS), LE(POINTS), TSTAR_GB(POINTS) common /land_r/ SH(POINTS) ! DAY_PHEN = IN Number of days between phenology. ! DAY_TRIF = IN Number of days between TRIFFID. integer DAY_PHEN, DAY_TRIF common /land_i/ DAY_PHEN, DAY_TRIF ! L_PHEN = IN .T. if phenology to be updated. ! L_TRIF = IN .T. if vegetation to be updated. logical L_PHEN, L_TRIF common /land_l/ L_PHEN, L_TRIF !====================== 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 $ ) 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 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 real fe_dissolved ! wyao 20180626 real rchl, rchl_D, mc real thetamaxlo, thetamaxhi common /npzd_r/ thetamaxlo, thetamaxhi, mc common /npzd_r/ rchl(imt,km,jsmw:jemw) common /npzd_r/ rchl_D(imt,km,jsmw:jemw) common /fe_dissolved_r/ fe_dissolved(imt,jmt,km,12) common /fe_dissolved_i/ kmfe !{ 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) real rbct, rjmax, rjmax_D, rgd, rgd_D common /npzd_r/ rbct(imt,kpzd,jsmw:jemw) common /npzd_r/ rjmax(imt,kpzd,jsmw:jemw) common /npzd_r/ rjmax_D(imt,kpzd,jsmw:jemw) common /npzd_r/ rgd(imt,kpzd,jsmw:jemw) common /npzd_r/ rgd_D(imt,kpzd,jsmw:jemw) common /npzd_r/ ravej(imt,kpzd,jsmw:jemw) common /npzd_r/ ravej_D(imt,kpzd,jsmw:jemw) common /npzd_r/ rgmax(imt,kpzd,jsmw:jemw) common /npzd_r/ rno3P(imt,kpzd,jsmw:jemw) common /npzd_r/ rpo4P(imt,kpzd,jsmw:jemw) common /npzd_r/ rpo4_D(imt,kpzd,jsmw:jemw) !====================== 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 ! sed include file ! ipmax = maximum number of sediment points ! nzmax = maximum number of mixed layers ! ibmax = maximum number of buried layers integer ipmax, nzmax, ibmax, imap, jmap, map_sed, ipsed, kmax integer kmin parameter (ipmax=imt*jmt, nzmax=8, ibmax=20) ! imap = 1d maping index for 2d map (i index) ! jmap = 1d maping index for 2d map (j index) ! map_sed = 2d map of sediment indices ! ipsed = number of sediment points (must be >= to ipmax) ! kmax = number of mixed layers (must be less than nzmax) ! kmin = minimum kmt level for sediments common /sed_i/ imap(ipmax), jmap(ipmax), map_sed(imt,jmt) common /sed_i/ ipsed, kmax, kmin ! mixed layer fields ! carb = forms of carbon, 1=co2,2=hco3,3=co3 (mol l-1) ! dcpls = diffusion coefficient above, 1=co2,2=hco3,3=co3 ! dcmin = diffusion coefficient below, 1=co2,2=hco3,3=co3 ! pore = porosity ! form = pore**3. ! o2 = bottom water oxygen (mol l-1) ! orggg = organic carbon fraction ! orgml = organic carbon mass (g cm-2) ! calgg = calcite fraction ! calml = calcite mass (g CaCO3 cm-2) ! dopls = diffusion coefficient for oxygen above ! domin = diffusion coefficient for oxygen below ! dbpls = diffusion coefficient for organic carbon above ! dbmin = diffusion coefficient for organic carbon below ! buried layer fields ! buried_mass = buried mass (g CaCO3 cm-2) ! buried_calfrac = buried calcite fraction ! depth_age = age of buried level (years) ! sediment surface or integrated fields ! zrct = maximum respiration depth (cm) ! water_z_p = ocean depth (m) ! k1 = first dissociation constant for carbonic acid ! k2 = second dissociation constant for carbonic acid ! k3 = dissociation constant for boric acid ! csat = carbonate saturation at sediment surface (mol l-1) ! rc = respiration coefficient ! sed_ml_mass = mixed layer mass (g cm-2) ! ttrorg = dissolution rate organic carbon (mol cm-2 dtsed-1) ! ttrcal = dissolution rate of calcite (mol cm-2 dtsed-1) ! c_advect = burial flux [mol cm-2 dtsed-1] ! zsed = mixed layer grid depth (cm) ! delz = mixed layer grid layer thickness (cm) ! rain_org_p = rain rate of organic carbon (mol cm-2 dtsed-1) ! rain_cal_p = rain rate of calcite (mol cm-2 dtsed-1) ! co3_p = carbonate at sediment surface (mol l-1) ! scalars ! dissc = calcite dissolution constant ! dissn = calcite dissolution constant ! weath = weathering flux (kg s-1) ! weathflx = weathering flux (umol s-1) ! sed_year = year for sediment profile ! sedsa = surface area of potential deep sediments (cm2) ! carblith = change in lithosphere carbon (umol) ! set to zero if time is initialized real carb, dcpls, dcmin, pore, form, o2, orggg, orgml, calgg real calml, dopls, domin, dbpls, dbmin, buried_mass real buried_calfrac, depth_age, zrct, water_z_p, k1, k2, k3, csat real rc, sed_ml_mass, ttrorg, ttrcal, c_advect, zsed, delz real rain_org_p, rain_cal_p, co3_p, dissc, dissn, weath, weathflx real dicwflx, sed_year, sedsa, carblith common /sed_r/ carb(nzmax,3,ipmax), dcpls(nzmax,3,ipmax) common /sed_r/ dcmin(nzmax,3,ipmax) common /sed_r/ pore(nzmax,ipmax), form(nzmax,ipmax) common /sed_r/ o2(nzmax,ipmax), orggg(nzmax,ipmax) common /sed_r/ orgml(nzmax,ipmax), calgg(nzmax,ipmax) common /sed_r/ calml(nzmax,ipmax), dopls(nzmax,ipmax) common /sed_r/ domin(nzmax,ipmax), dbpls(nzmax,ipmax) common /sed_r/ dbmin(nzmax,ipmax), buried_mass(ibmax,ipmax) common /sed_r/ buried_calfrac(ibmax,ipmax), depth_age(ibmax) common /sed_r/ zrct(ipmax), water_z_p(ipmax), k1(ipmax), k2(ipmax) common /sed_r/ k3(ipmax), csat(ipmax), rc(ipmax), ttrorg(ipmax) common /sed_r/ ttrcal(ipmax), sed_ml_mass(ipmax), c_advect(ipmax) common /sed_r/ zsed(nzmax), delz(nzmax), rain_org_p(ipmax) common /sed_r/ rain_cal_p(ipmax), co3_p(ipmax) common /sed_r/ dissc, dissn, weath, weathflx, dicwflx, sed_year common /sed_r/ sedsa, carblith ! ntatss = time step counter for time averaging ! ta_ttrcal = time average dissolution rate of calcite ! ta_rain_cal = time average rain rate of calcite ! ta_cal = time average mixed layer calcite fraction ! ta_calmass = time average mixed layer calcite mass ! ta_calmass_bur = time average buried calcite mass ! ta_co3 = time average carbonate at sediment surface ! ta_co3sat = time average carbonate saturation at surface ! ta_rainr = time average rain ratio integer ntatss common /sed_i/ ntatss real ta_ttrcal, ta_rain_cal, ta_cal, ta_calmass, ta_calmass_bur real ta_co3, ta_co3sat, ta_rainr common /sed_r/ ta_ttrcal(imt,jmt), ta_rain_cal(imt,jmt) common /sed_r/ ta_cal(imt,jmt), ta_calmass(imt,jmt) common /sed_r/ ta_calmass_bur(imt,jmt), ta_co3(imt,jmt) common /sed_r/ ta_co3sat(imt,jmt), ta_rainr(imt,jmt) ! ntatis = number of time averaged time step integrals ! tai_ttrcal = time step integral of ttrcal ! tai_rain_cal = time step integral of rain_cal ! tai_cal = time step integral of cal ! tai_weathflx = time step integral of weathflx ! tai_calmass = time step integral of mixed layer calcite mass ! tai_calmass_bur = time step integral of buried calcite mass ! tai_co3 = time step integral of co3 ! tai_co3sat = time step integral of co3sat ! tai_rainr = time step integral of rain ratio ! tai_carblith = time step integral of carblith ! tai_cfo2s = average total flux ocean to sediments ! tai_cfl2o = average total flux land to ocean (weathering) integer ntatis common /sed_i/ ntatis real tai_ttrcal, tai_rain_cal, tai_cal, tai_weathflx, tai_calmass real tai_calmass_bur, tai_co3, tai_co3sat, tai_rainr, tai_carblith real tai_cfo2s, tai_cfl2o common /sed_r/ tai_ttrcal, tai_rain_cal, tai_cal, tai_weathflx common /sed_r/ tai_calmass, tai_calmass_bur, tai_co3, tai_co3sat common /sed_r/ tai_rainr, tai_carblith, tai_cfo2s, tai_cfl2o !====================== include file "stab.h" ========================== ! CFL and other stability criteria information ! cflons = starting longitude (degrees) for stability tests ! cflone = ending longitude (degrees) for stability tests ! cflats = starting latitude (degrees) for stability tests ! cflate = ending latitude (degrees) for stability tests ! cfldps = starting depth (cm) for stability tests ! cfldpe = ending depth (cm) for stability tests ! iscfl = index corresponding to "cflons" ! iecfl = index corresponding to "cflone" ! jscfl = index corresponding to "cflats" ! jecfl = index corresponding to "cflate" ! kscfl = index corresponding to "cfldps" ! kecfl = index corresponding to "cfldpe" ! cflcrt = factor by which the cfl criteria must be exceeded in ! order to show local values (see blkdta.F) ! maxcfl = maximum number of times the "cflcrt" factor can be ! exceeded before stopping. ! numcfl = counter for number of times the "cflcrt" factor was ! exceeded. ! cflum = zonal velocity which comes closest to its cfl criteria ! cflup = percent of cfl criteria reached by "cflum" ! icflu = "i" coordinate of "cflum" ! jcflu = "j" coordinate of "cflum" ! kcflu = "k" coordinate of "cflum" ! cflvm = meridional velocity which comes closest to its cfl ! criteria ! cflvp = percent of cfl criteria reached by "cflvm" ! icflv = "i" coordinate of "cflvm" ! jcflv = "j" coordinate of "cflvm" ! kcflv = "k" coordinate of "cflvm" ! cflwtm = vertical velocity at "t" box bottom ! which comes closest to its cfl criteria ! cflwtp = percent of cfl criteria reached by "cflwtm" ! icflwt = "i" coordinate of "cflwtm" ! jcflwt = "j" coordinate of "cflwtm" ! kcflwt = "k" coordinate of "cflwtm" ! cflwum = vertical velocity at "u,v" box bottom ! which comes closest to its cfl criteria ! cflwup = percent of cfl criteria reached by "cflwum" ! icflwu = "i" coordinate of "cflwum" ! jcflwu = "j" coordinate of "cflwum" ! kcflwu = "k" coordinate of "cflwum" ! reynx = maximum reynolds number in the zonal direction ! ireynx = "i" coordinate of "reynx" ! jreynx = "j" coordinate of "reynx" ! kreynx = "k" coordinate of "reynx" ! reynu = "u" for computing "reynx" ! reynmu = zonal mixing of momentum for computing "reynx" ! reyny = maximum reynolds number in the meridional direction ! ireyny = "i" coordinate of "reyny" ! jreyny = "j" coordinate of "reyny" ! kreyny = "k" coordinate of "reyny" ! reynv = "v" for computing "reyny" ! reynmv = meridional mixing of momentum for computing "reyny" ! reynz = maximum reynolds number in the vertical direction ! ireynz = "i" coordinate of "reynz" ! jreynz = "j" coordinate of "reynz" ! kreynz = "k" coordinate of "reynz" ! reynw = "w" for computing "reynz" ! reynmw = vertical mixing of momentum for computing "reynz" ! peclx = maximum peclet number in the zonal direction ! ipeclx = "i" coordinate of "peclx" ! jpeclx = "j" coordinate of "peclx" ! kpeclx = "k" coordinate of "peclx" ! peclu = "u" for computing "peclx" ! peclmu = zonal mixing of tracer for computing "peclx" ! pecly = maximum peclet number in the meridional direction ! ipecly = "i" coordinate of "pecly" ! jpecly = "j" coordinate of "pecly" ! kpecly = "k" coordinate of "pecly" ! peclv = "v" for computing "pecly" ! peclmv = meridional mixing of tracer for computing "pecly" ! peclz = maximum peclet number in the vertical direction ! ipeclz = "i" coordinate of "peclz" ! jpeclz = "j" coordinate of "peclz" ! kpeclz = "k" coordinate of "peclz" ! peclw = "w" for computing "peclz" ! peclmw = vertical mixing of tracer for computing "peclz" ! tdig = factor by which local tracer extremum must be exceeded ! before showing ficticious creation of tracer integer iscfl, iecfl, jscfl, jecfl, kscfl, kecfl integer icflu, jcflu, kcflu, icflv, jcflv, kcflv integer icflwt, jcflwt, kcflwt, icflwu, jcflwu, kcflwu integer numcfl, maxcfl integer ireynx, jreynx, kreynx, ireyny, jreyny, kreyny integer ireynz, jreynz, kreynz, ipeclx, jpeclx, kpeclx integer ipecly, jpecly, kpecly, ipeclz, jpeclz, kpeclz common /stab_i/ iscfl, iecfl, jscfl, jecfl, kscfl, kecfl common /stab_i/ icflu, jcflu, kcflu, icflv, jcflv, kcflv common /stab_i/ icflwt, jcflwt, kcflwt, icflwu, jcflwu, kcflwu common /stab_i/ numcfl, maxcfl common /stab_i/ ireynx, jreynx, kreynx, ireyny, jreyny, kreyny common /stab_i/ ireynz, jreynz, kreynz, ipeclx, jpeclx, kpeclx common /stab_i/ ipecly, jpecly, kpecly, ipeclz, jpeclz, kpeclz real cflons, cflone, cflats, cflate, cfldps, cfldpe real cflup, cflum, cflvp, cflvm, cflwtp, cflwtm real cflwup, cflwum, cflcrt, tdig real reynx, reynu, reynmu, reyny, reynv, reynmv real reynz, reynw, reynmw, peclx, peclu, peclmu real pecly, peclv, peclmv, peclz, peclw, peclmw common /stab_r/ cflons, cflone, cflats, cflate, cfldps, cfldpe common /stab_r/ cflup, cflum, cflvp, cflvm, cflwtp, cflwtm common /stab_r/ cflwup, cflwum, cflcrt, tdig common /stab_r/ reynx, reynu, reynmu, reyny, reynv, reynmv common /stab_r/ reynz, reynw, reynmw, peclx, peclu, peclmu common /stab_r/ pecly, peclv, peclmv, peclz, peclw, peclmw !====================== 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 !======================= 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 "veg.h" ======================== ! variables for the vegetation data ! Vegetation types ! 1 = tropical forest ! 2 = temperate/boreal forest ! 3 = grass ! 4 = shrub ! 5 = tundra ! 6 = desert ! 7 = ice ! nveg = number of vegetation classes ! iveg = vegetation class ! veg_rl = roughness length ! veg_alb = albedo ! veg_rs = stomatal resistance ! veg_smd = snow masking depth (m) ! veg_dalt = dalton number over land ! idesert = index for desert ! iagric = index for agricultural land ! icrops = this variable is obsolete and its value ignored ! land use ("crop_data" "pasture_data" "agric_data") ! agric = agricultural extent as percentage of grid cell integer nveg, iveg, idesert, iagric, icrops, iice real veg_rl, veg_alb, veg_rs, veg_smd, veg_dalt, agric parameter (nveg=7) common /veg_i/ iveg(imt,jmt), idesert, iagric, icrops, iice common /veg_r/ veg_rl(nveg), veg_alb(nveg), veg_rs(nveg) common /veg_r/ veg_smd(nveg) common /veg_r/ veg_dalt(nveg) !====================== 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) !----------------------------------------------------------------------- ! read all model namelist variables !----------------------------------------------------------------------- namelist /contrl/ init, runlen, rununits, restrt, initpt &, num_processors, runstep, initbgc namelist /tsteps/ dtts, dtuv, dtsf, dtatm, dtatms, namix, segtim &, daylen namelist /riglid/ mxscan, sor, tolrsf, tolrsp, tolrfs namelist /mixing/ am, ah, ahbkg, ambi, ahbi, kappa_m, kappa_h &, cdbot, spnep, senep, aidif, ncon, nmix, eb &, acor, dampts, dampdz, annlev, annlevobc namelist /diagn/ tsiint, tsiper, tavgint, itavg, tmbint, tmbper &, itmb, stabint, zmbcint, glenint, trmbint, itrmb &, vmsfint, gyreint,igyre, extint, prxzint, trajint &, exconvint, dspint, dspper, snapint, snapls &, snaple, snapde, timavgint, timavgper, cmixint &, prlat, prslon, prelon, prsdpt, predpt &, slatxy, elatxy, slonxy, elonxy &, cflons, cflone, cflats, cflate, cfldps, cfldpe &, maxcfl, xbtint, xbtper, crossint, densityint &, fctint, tyzint, restint, tbtint, tbtper &, accel, accel_yr0 namelist /io/ expnam, logfile, runstamp, iotavg, iotmb, iotrmb &, ioglen, iovmsf, iogyre, ioprxz, ioext, iodsp &, iotsi, iozmbc, iotraj, ioxbt, isot1, ieot1 &, isot2, ieot2, jsot, jeot, ksot, keot, mrot namelist /ictime/ year0, month0, day0, hour0, min0, sec0, ryear &, rmonth, rday, rhour, rmin, rsec, refrun, refinit &, refuser, eqyear, eqmon, monlen, init_time &, init_time_in, init_time_out, omega namelist /npzd/ alpha, kw, kc, abio, bbio, cbio, k1n, nup &, gamma1, gbio, nuz, nud0 &, wd0, par, dtnpzd, redctn, ki, redptn, capr &, dcaco3, nupt0, redotn, jdiar, kzoo, geZ &, zprefP, zprefZ, zprefDet, zprefD, kfe, kfe_D &, wdi namelist /fwa/ isfwa1, iefwa1, isfwa2, iefwa2, jsfwa, jefwa &, mrfwa, fwaflxi, fwayri, fwayrf, fwarate &, compensate namelist /blmix/ Ahv, Ahh namelist /hlmix/ hl_depth, hl_back, hl_max namelist /isopyc/ slmx, ahisop, athkdf, del_dm, s_dm namelist /ppmix/ wndmix, fricmx, diff_cbt_back, visc_cbu_back &, visc_cbu_limit, diff_cbt_limit namelist /smagnl/ diff_c_back namelist /sed/ dtsedyr, nsedacc, weath, kmin namelist /embm/ rlapse, rf1, rf2, scatter, rhmax, vcsref, vcsfac &, vcsyri, aggfor_os, adiff, cropf, pastf namelist /carbon/ co2ccn, co2emit, co2for, co2_yr, c14_yr, dc14ccn &, c14prod namelist /paleo/ orbit_yr, pyear, eccen, obliq, mvelp, sealev &, sealev_yr, volcano_yr, sulph_yr, aggfor_yr &, cfcs_yr namelist /ice/ niats, nivts, dampice, tsno, hsno_max, ice_yr &, landice_yr, ice_calb, sno_calb namelist /veg/ veg_alb, veg_rl, veg_rs, veg_smd, agric_yr &, crops_yr, iagric, icrops, idesert, iice namelist /solar/ solarconst, solar_yr namelist /mtlm/ TIMESTEP, INT_VEG, VEG_EQUIL, DAY_TRIF, DAY_PHEN &, BF ! physical constants rho0 = 1.035 rho0r = c1/rho0 grav = 980.6 radius = 6370.0e5 pi = atan(1.0) * 4.0 ! set defaults for namelist contrl init = .true. runlen = 365.0 rununits = 'days' restrt = .true. initpt = .false. num_processors = 1 runstep = -1.0 initbgc = .true. ! set defaults for namelist tsteps dtts = 43200.0 dtuv = 600.0 dtsf = 600.0 dtatm = 43200.0 dtatms = 43200.0 namix = 16 segtim = 1.0 daylen = 86400.0 ! set defaults for namelist riglid mxscan = 300 sor = 1.60 tolrsf = 5.0e8 tolrsp = 1.0e4 tolrfs = 1.0e4 ! set defaults for namelist mixing am = 2.0e9 ah = 2.0e7 ahbkg = 0. ambi = 1.0e23 ahbi = 5.0e22 kappa_m = 10.0 kappa_h = 1.0 cdbot = 1.3e-3 spnep = 3.0e5 senep = 12.0e5 aidif = 1.0 ncon = 1 nmix = 16 eb = .false. acor = 0.0 do n=1,nt dampts(n) = 50.0 dampdz(n) = 26.575e2 enddo annlev = .false. annlevobc = .false. ! set defaults for namelist diagn tsiint = 1.0 tsiper = 1.0 tavgint = -36500.0 itavg = .true. tmbint = -36500.0 tmbper = 365.0 itmb = .true. stabint = -36500.0 zmbcint = -36500.0 glenint = -36500.0 trmbint = -36500.0 itrmb = .true. vmsfint = -36500.0 gyreint = -36500.0 igyre = .true. extint = -36500.0 prxzint = -36500.0 trajint = 0.0 exconvint = -36500.0 dspint = -36500.0 dspper = -365.0 snapint = 0.0 snapls = 0.0 snaple = 0.0 snapde = 0.0 timavgint = 36500.0 timavgper = 365.0 cmixint = -36500.0 do n=1, nlatpr prlat(n) = 100.0 prslon(n) = 0.0 prelon(n) = 0.0 prsdpt(n) = 0.0 predpt(n) = 6000.0e2 if (n. le. 4) then prslon(n) = 180.0 prelon(n) = 250.0 endif enddo prlat(1) = -60.0 prlat(2) = 0.0 prlat(3) = 27.0 prlat(4) = 55.0 slatxy = -90.0 elatxy = 90.0 slonxy = 3.0 elonxy = 357.0 cflons = 0.0 cflone = 360.0 cflats = -90.0 cflate = 90.0 cfldps = 0.0 cfldpe = 6000.0e2 maxcfl = 3 xbtint = -36500.0 xbtper = -365.0 crossint = 365000.0 densityint = -36500.0 fctint = -36500.0 tyzint = -36500.0 restint = 36500.0 tbtint = -36500.0 tbtper = -365.0 accel = 1. accel_yr0 = 0. ! set defaults for namelist io expnam = ' ' logfile = 'machine.log' runstamp = ' ' restrt = .false. iotavg = -1 iotmb = -1 iotrmb = -1 ioglen = -1 iovmsf = -1 iogyre = -1 ioprxz = -1 ioext = -1 iodsp = -1 iotsi = -1 iozmbc = -1 iotraj = -1 ioxbt = -1 isot1 = 2 ieot1 = imtm1 isot2 = 2 ieot2 = 1 jsot = 2 jeot = jmtm1 ksot = 1 keot = km mrot = 0 ! set defaults for namelist ictime year0 = 1 month0 = 1 day0 = 1 hour0 = 0 min0 = 0 sec0 = 0 ryear = 1 rmonth = 1 rday = 1 rhour = 0 rmin = 0 rsec = 0 refrun = .false. refinit = .true. refuser = .false. eqyear = .true. eqmon = .false. monlen = 30 init_time = .false. init_time_in = .false. init_time_out = .false. omega = pi/43082.0 ! set defaults for namelist npzd alpha = 0.1 ! Initial slope P-I curve [(W/m^2)^(-1)/day] kw = 0.04 ! Light attenuation due to water [1/m] kc = 0.047 ! Light atten. by phytoplankton [1/(m*mmol/m^3)] ki = 5.0 ! Light attenuation through sea ice & snow abio = 0.18 ! a; Maximum growth rate parameter [1/day] bbio = 1.066 ! b cbio = 1.0 ! c [1/deg_C] k1n = 0.7 ! Half saturation constant for N uptake [mmol/m^3] nup = 0.025 ! Specific mortality rate (Phytoplankton) [1/day] nupt0 = 0.02 ! Specific mortality rate (Phytoplankton) [1/day] gamma1 = 0.70 ! gama1; Assimilation efficiency (zpk) gbio = 0.18 ! Maximum grazing rate at 0 deg C [1/day] nuz = 0.34 ! Quadratic mortality (zpk) [(mmol/m^3)^(-2)day^(-1)] nud0 = 0.048 ! remineralization rate [1/day] wd0 = 6.0 ! Sinking speed of detritus at surface [m/day] par = 0.43 ! fraction of photosythetically active radiation dtnpzd = dtts ! time step of biology [s] redctn = 7. ! C/N Redfield ratio redptn = 1./16. ! P/N Redfield ratio capr = 0.018 ! carbonate to carbon production ratio dcaco3 = 650000.0 ! remineralisation depth of calcite [cm] redotn = 10.6 ! O2/N Redfield ratio jdiar = 0.5 ! factor reducing the growth rate of diazotrophs kzoo = 0.15 ! half sat. constant for Z grazing [mmol N m^3] geZ = 0.50 ! Zooplankton growth efficiency zprefP = 0.30 ! Zooplankton preference for P zprefZ = 0.30 ! Zooplankton preference for other Z zprefDet = 0.30 ! Zooplankton preference for Detritus zprefD = 0.10 ! Zooplankton preference for diazotrophs !{ wyao 20170714 wdi = 0.06 ! Sinking speed of detritus with depth [m-1] ! wyao 20170714} kmfe = 3 ! number of depth levels used for Fe limitation kfe = 0.1e-3 ! Half saturation constant for Fe limitation [mmol Fe / m-3] kfe_D = 0.1e-3 ! Half saturation constant for Diaz Fe limitation [mmol Fe / m-3] !{wyao 20171018 after unify the iron currency unit in the fe_limit version, the default value is now e-3} ! wyao 20180626 thetamaxhi = 0.04 ! Maximum Chl:C ratio, abundant iron [gChl/(gC)] thetamaxlo = 0.01 ! Maximum Chl:C ratio, extreme iron limitation [gChl/(gC)] mc = 12.011 ! Molar mass of carbon [g/mol] ! set defaults for namelist fwa isfwa1 = 1 iefwa1 = imt isfwa2 = 0 iefwa2 = 0 jsfwa = 1 jefwa = jmt mrfwa = 1 fwaflxi = 0. fwayri = -1.e20 fwayrf = 1.e20 fwarate = 0. compensate = .false. ! set defaults for namelist blmix ! Reference: ! A Water Mass Model of the World Ocean K. Bryan, L.J. Lewis ! JGR, vol 84, No. C5, May 20, 1979 afkph = 0.8 dfkph = 1.05 sfkph = 4.5e-5 zfkph = 2500.0e2 ! Use Bryan & Lewis values for vertical tracer diffusion ! Ahv range of 0.3 to 1.3, crossover at 2500m. ! compute depth dependent vertical diffusion coefficients for ! tracers using the relationship of Bryan and Lewis pi = 4.0 * atan(1.0) ! compute depth dependent horizontal diffusion coefficients for ! tracers using the relationship of Bryan and Lewis ahs = 5.0e+3 ahb = 1.0e+3 ! set defaults for namelist hlmix hl_depth = 500.0e2 hl_back = 1.e4 hl_max = 1.e9 ! set defaults for namelist isopyc slmx = 1.0/100.0 ! maximum isopycnal slope ahisop = 1.e7 ! isopycnal diffusion coefficient athkdf = 1.0e7 ! isopycnal thickness diffusion coefficient del_dm = 4.0/1000.0 ! transition for scaling diffusion coefficient s_dm = 1.0/1000.0 ! half width scaling for diffusion coefficient ! set defaults for namelist ppmix wndmix = 10.0 fricmx = 50.0 diff_cbt_back = 0.1 visc_cbu_back = 1.0 ! in regions of gravitational instability set mixing limits to the ! maximum consistant with the "cfl" criterion. convective adjustment ! will also act on the instability. visc_cbu_limit = fricmx diff_cbt_limit = fricmx ! set defaults for namelist smagnl diff_c_back = c0 ! set defaults for namelist sed dtsedyr = 1. ! time step of sediment model in years nsedacc = 1 ! number of steps for accelerating sediments weath = 2.e20 ! weathering flux in Pg year-1 ! set high to check for constant namelist input kmin = 8 ! minimum ocean level for sediments ! set defaults namelist embm rlapse = 5.e-5 rf1 = 0.3 rf2 = 3.e5 scatter = 0.23 rhmax = 0.85 vcsref = 0. vcsfac = 0. vcsyri = 0. aggfor_os = 0. adiff = 0.03 cropf = 1. pastf = 0.5 ! set defaults namelist carbon co2ccn = 280. co2emit = 0. co2for = 5.35e03 co2_yr = 2.e20 c14_yr = 2.e20 dc14ccn = 0. c14prod = 0. ! set defaults namelist paleo pyear = 2.e20 orbit_yr = 2.e20 ! default user orbit is for 1950 eccen = 1.672393E-02 obliq = 23.44627 mvelp = 102.0391 sealev = 0. sealev_yr = 2.e20 volcano_yr = 2.e20 sulph_yr = 2.e20 aggfor_yr = 2.e20 cfcs_yr = 2.e20 ! set defaults namelist ice niats = 1 nivts = 1 dampice = 5. tsno = 0. hsno_max = 1000. ice_yr = 2.e20 landice_yr = 2.e20 ice_calb = 0.25 sno_calb = 0.2 ! set defaults namelist veg veg_alb(1:7) = (/0.17, 0.17, 0.22, 0.22, 0.22, 0.30, 0.60/) veg_rl(1:7) = (/1.5, 1.0, 0.1, 0.3, 0.01, 0.01, 0.005/) veg_rs(1:7) = (/0.75, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/) veg_smd(1:7) = (/2.0, 3.0, 0.1, 0.1, 0.5, 0.01, 0.01/) agric_yr = 2.e20 crops_yr = 2.e20 iagric = 3 icrops = 3 idesert = 6 iice = 7 ! set defaults namelist solar solarconst = 1.368e6 ! mw m-2 solar_yr = 2.e20 ! set defaults namelist mtlm TIMESTEP = 3600. INT_VEG = .true. VEG_EQUIL = .false. DAY_TRIF = 30 DAY_PHEN = 1 BF = 1. !----------------------------------------------------------------------- ! provide for change in above presets using "namelist" !----------------------------------------------------------------------- call getunit (ioun, 'control.in', 'f s r') read (ioun, contrl, end=101) 101 continue runstep = float(int(runstep)) write (stdout, contrl) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, tsteps, end=102) 102 continue write (stdout, tsteps) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, riglid, end=103) 103 continue write (stdout, riglid) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, mixing, end=104) 104 continue write (stdout, mixing) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, diagn, end=105) 105 continue write (stdout, diagn) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, io, end=106) 106 continue write (stdout, io) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, ictime, end=107) 107 continue write (stdout, ictime) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, npzd, end=108) 108 continue write (stdout, npzd) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, fwa, end=109) 109 continue write (stdout, fwa) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, blmix, end=110) 110 continue write (stdout, blmix) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, hlmix, end=111) write (stdout, hlmix) call relunit (ioun) 111 continue call getunit (ioun, 'control.in', 'f s r') read (ioun, isopyc, end=112) 112 continue write (stdout, isopyc) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, ppmix, end=113) 113 continue write (stdout, ppmix) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, smagnl, end=114) 114 continue write (stdout, smagnl) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, sed, end=115) 115 continue write (stdout, sed) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, embm, end=116) 116 continue write (stdout, embm) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, carbon, end=117) 117 continue write (stdout, carbon) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, paleo, end=118) 118 continue write (stdout, paleo) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, ice, end=119) 119 continue write (stdout, ice) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, veg, end=120) 120 continue write (stdout, veg) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, solar, end=121) 121 continue write (stdout, solar) call relunit (ioun) call getunit (ioun, 'control.in', 'f s r') read (ioun, mtlm, end=122) 122 continue write (stdout, mtlm) call relunit (ioun) ! limit or convert some variables ! ensure pass is between zero and one. pass = min(max((1. - scatter), 0.), 1.) dtocn = dtts nots = nmix isot1 = max(isot1, 1) ieot1 = min(ieot1, imt) isot2 = max(isot1, 1) ieot2 = min(ieot1, imt) jsot = max(jsot, 1) jeot = min(jeot, jmt) ksot = max(ksot, 1) keot = min(keot, km) do k=1,km dtxcel(k) = 1.0 enddo ah = ahbkg slmxr = 0. if (slmx .gt. 0.) slmxr = c1/slmx s_dm = 0. if (s_dm .gt. 0.) s_dmr = c1/s_dm ! sort out forcing years if (pyear .gt. 1.e20) pyear = 1800. if (co2_yr .gt. 1.e20) co2_yr = pyear if (c14_yr .gt. 1.e20) c14_yr = pyear if (orbit_yr .gt. 1.e20) orbit_yr = pyear if (volcano_yr .gt. 1.e20) volcano_yr = pyear if (sulph_yr .gt. 1.e20) sulph_yr = pyear if (aggfor_yr .gt. 1.e20) aggfor_yr = pyear if (cfcs_yr .gt. 1.e20) cfcs_yr = pyear if (sealev_yr .gt. 1.e20) sealev_yr = pyear if (solar_yr .gt. 1.e20) solar_yr = pyear if (agric_yr .gt. 1.e20 .and. crops_yr .lt. 1.e20) & agric_yr = crops_yr if (agric_yr .gt. 1.e20) agric_yr = pyear if (landice_yr .gt. 1.e20 .and. ice_yr .lt. 1.e20) & landice_yr = ice_yr if (landice_yr .gt. 1.e20) landice_yr = pyear ice_yr = landice_yr tsiper = min(tsiper, tsiint) tbtper = min(tbtper, tbtint) if (accel .le. 1) accel = 1. runlen = runlen/accel runstep = runstep/accel if (runstep .gt. 0.0) runlen = min(runstep, runlen) if (init_time) then init_time_in = .true. init_time_out = .true. endif monlen = 30 eqyear = .true. eqmon = .false. calendar = 'noleap' !----------------------------------------------------------------------- ! set runstamp !----------------------------------------------------------------------- if (runstamp .eq. ' ') then ! read runstamp from last line of logfile if it exists fname = new_file_name (logfile) inquire (file=trim(fname), exist=exists) if (exists) then call getunit (ioun, fname, 'f s r') do while (exists) read (ioun, '(a130)', end=130) runstamp enddo 130 continue call relunit (ioun) print*, " " print*, "runstamp: ", trim(runstamp) endif endif return end