! source file: /Users/nmengis/Documents/UVic_ESCM/2.10.backup/updates/02_CE_permafrost_merge_CMIP6forcing/source/common/UVic_ESCM.F program UVic_ESCM !======================================================================= ! UNIVERSITY OF VICTORIA EARTH SYSTEM CLIMATE MODEL (UVic ESCM) ! The UVic ESCM is a climate model developed by researchers in the ! Climate Modelling Group, within the School of Earth and Ocean ! Sciences, located at the University of Victoria, Victoria, ! British Columbia, Canada. ! Many people have contributed to the development of this code. ! This is a collective effort and although individual contributions ! are appreciated, individual authorship is not indicated. ! Some of the people who have contributed to the development of the ! model are: ! D. Archer, C. Avis, A. Berger, R. Betts, A. Biastoch, C. Bitz, ! J. Brauch, C. Brennan, B. Bryan, J. Burton, S. Carto, ! M. Cottet-Puinel, M. Cox, P. Cox, G. Danabasoglu, K. Dixon, ! J. Dukowicz, M. Eby, R. Essery, T. Ewen, A. Fanning, J. Fyke, ! R. Gerdes, A. Gnanadesikan, C. Goldberg, D. Goldberg, ! D. Gregory, J. Gregory, S. Griffies, R. Hanson, R. Hetherington, ! H. Hickey, M. Holland, G. Holloway, T. Huck, T. Hughes, ! E. Hunke, W. Hurlin, W. Ingram, R. Key, E. Kluzek, C. Koeberle, ! K. Kvale, C. Lawson, J. Lewis, M. Loutre, R. Malone, ! D. Matthews, K. Meissner, A. Montenegro, A. Mouchet, T. Murdock, ! P. Myers, A. Oschlies, R. Pacanowski, M. Pahlow, P. Poussart ! S. Rahmstorf, R. Redler,, D. Robitaille, A. Rosati, M. Roth ! C. Sabine, O. Saenko, A. Schmittner, B. Semtner, W. Sijp ! T. Silva, H. Simmons, A. Skvortsov, R. Smith, P. Spence ! D. Stone, R. Tonkonojenkov, C. Tricot, S. Valcke, A. Weaver ! E. Wiebe, J. Willebrandt, M. Yoshimori, K. Zickfeld ! Please direct problems or questions to the code contact ! person at: http://climate.uvic.ca/model ! Requirements: ! Standard fortran 90 is used ! Disclaimer: ! The UVic Earth System Climate Model is a climate modelling ! research tool developed at the University of Victoria. Others may ! use it freely but we assume no responsibility for problems or ! incorrect use. It is left to the user to ensure that a particular ! configuration is working correctly. !======================================================================= implicit none integer i, j, numots, numats, numseg, n, loop real relday, hst, relyr0 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "coord.h" include "calendar.h" include "csbc.h" include "iounit.h" include "emode.h" include "levind.h" include "scalar.h" include "switch.h" include "tmngr.h" include "cembm.h" include "atm.h" include "mw.h" print*, '== UNIVERSITY OF VICTORIA EARTH SYSTEM CLIMATE MODEL ==' print*, ' ' !----------------------------------------------------------------------- ! initialize i/o units !----------------------------------------------------------------------- call ioinit !----------------------------------------------------------------------- ! setup file renaming !----------------------------------------------------------------------- call file_names !----------------------------------------------------------------------- ! Initialize S.B.C. indices !----------------------------------------------------------------------- call sbc_init !----------------------------------------------------------------------- ! Initialize tracers !----------------------------------------------------------------------- call tracer_init !----------------------------------------------------------------------- ! read namelist variables !----------------------------------------------------------------------- call read_namelist !----------------------------------------------------------------------- ! read grid !----------------------------------------------------------------------- call grids !----------------------------------------------------------------------- ! read topography !----------------------------------------------------------------------- call topog (kmt, kmu, map, xt, yt, zt, xu, yu, zw, imt, jmt, km) call isleperim (kmt, map, iperm, jperm, iofs, nippts, nisle, imt &, jmt, km, mnisle, maxipp, xu, yu, zw) !----------------------------------------------------------------------- ! common setup !----------------------------------------------------------------------- call setcom (1, imt, 1, jmt) !----------------------------------------------------------------------- ! ocean setup !----------------------------------------------------------------------- call setmom (1, imt, 1, jmt) !----------------------------------------------------------------------- ! sediment setup (requires mom) !----------------------------------------------------------------------- call setsed (1, imt, 1, jmt) !----------------------------------------------------------------------- ! atmosphere setup !----------------------------------------------------------------------- call setembm (1, imt, 1, jmt) !----------------------------------------------------------------------- ! land setup (requires embm) !----------------------------------------------------------------------- call setmtlm (1, imt, 1, jmt) !----------------------------------------------------------------------- ! initialize data reading routine !----------------------------------------------------------------------- call setdata (1, imt, 1, jmt) !----------------------------------------------------------------------- ! compute the number of ocean time steps "numots" for this run and ! the number of ocean time steps per ocean segment "ntspos". ! compute the number of atmos time steps "numats" for this run and ! the number of atmos time steps per atmos segment "ntspas". ! divide the integration time "days" into "numseg" segments. ! each will be length "segtim" days. Surface boundary conditions ! are supplied every "segtim" days. !----------------------------------------------------------------------- numots = nint(rundays/(dtocn*secday)) ntspos = nint(segtim/(dtocn*secday)) numats = nint(rundays/(dtatm*secday)) ntspas = nint(segtim/(dtatm*secday)) numseg = numots/ntspos if (segtim .gt. 1.) then ntspls = nint(c1/(dtlnd*secday)) else ntspls = nint(segtim/(dtlnd*secday)) endif !----------------------------------------------------------------------- ! check for consistency in the S.B.C. setup !----------------------------------------------------------------------- call chkcpl !----------------------------------------------------------------------- ! get global sums at the start of the run !----------------------------------------------------------------------- call globalsum (1) hst = 0.5*segtim relyr0 = 0. if (refrun) relyr0 = relyr !----------------------------------------------------------------------- ! S T A R T S E G M E N T L O O P !----------------------------------------------------------------------- do n=1,numseg !----------------------------------------------------------------------- ! set output flags !----------------------------------------------------------------------- relday = (relyr-relyr0)*yrlen + hst if (mod(relday+segtim,tsiint) .lt. segtim) then tsits = .true. else tsits = .false. endif if (mod(relday+tsiper,tsiint) .lt. tsiper) then tsiperts = .true. else tsiperts = .false. endif if (mod(relday+segtim,timavgint) .lt. segtim) then timavgts = .true. else timavgts = .false. endif if (mod(relday+timavgper,timavgint) .lt. timavgper) then timavgperts = .true. else timavgperts = .false. endif if (mod(relday+segtim,restint) .lt. segtim) then restts = .true. else restts = .false. endif if (n .eq. numseg) then eorun = .true. if (timavgperts) timavgts = .true. if (tsiperts) tsits = .true. endif !----------------------------------------------------------------------- ! get the atmospheric S.B.C. !----------------------------------------------------------------------- call gasbc (1, imt, 1, jmt) !----------------------------------------------------------------------- ! call the atmospheric model once for each time step until one ! segment of "segtim" days is complete. hold atmos S.B.C. fixed ! during each segment and predict average S.B.C. for ocean !----------------------------------------------------------------------- do loop=1,ntspas call embm (1, imt, 1, jmt) enddo !----------------------------------------------------------------------- ! get land S.B.C.s !----------------------------------------------------------------------- call glsbc (1, imt, 1, jmt) !---------------------------------------------------------------------- ! call the land-surface and vegetation model once for each time ! step until one segment of "segtim" days is complete. !----------------------------------------------------------------------- do loop=1,ntspls call mtlm (1, imt, 1, jmt) enddo !----------------------------------------------------------------------- ! get ocean S.B.C.s !----------------------------------------------------------------------- call gosbc (1, imt, 1, jmt) !----------------------------------------------------------------------- ! call the ocean model !----------------------------------------------------------------------- eoseg = .false. do loop=1,ntspos if (loop .eq. ntspos) eoseg = .true. call mom enddo !----------------------------------------------------------------------- ! call the sediment model if needed !----------------------------------------------------------------------- if (mod(relyr,dtsedyr) .lt. 1.e-6) call sed (1, imt, 1, jmt) !----------------------------------------------------------------------- ! write output !----------------------------------------------------------------------- call embmout (1, imt, 1, jmt) call mtlmout (1, imt, 1, jmt) call sedout (1, imt, 1, jmt) !----------------------------------------------------------------------- ! write change in global sums for heat and fresh water !----------------------------------------------------------------------- if (tsits) call globalsum (2) !----------------------------------------------------------------------- ! close any open netcdf files to flush buffers !----------------------------------------------------------------------- call closeall enddo !----------------------------------------------------------------------- ! E N D S E G M E N T L O O P !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! get global sums at the end of the run !----------------------------------------------------------------------- call globalsum (3) print*, ' ==> UVIC_ESCM integration is complete.' call closeall call release_all stop end subroutine chkcpl implicit none integer jrow, j, i, k logical errorc real critv, tmp, t1, r1, r2, r3, r4, r5, r6 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "csbc.h" include "cembm.h" include "calendar.h" include "switch.h" include "tmngr.h" !----------------------------------------------------------------------- ! adjust output intervals for forcing acceleration !----------------------------------------------------------------------- if (accel .le. 1) then accel = 1. else print*, ' ' print*, '==> Warning: any forcing that changes on time scales' & , ' longer than 1 year has been accelerated' print*, ' acceleration factor: ', accel print*, ' starts at relyr: ', accel_yr0 print*, ' current year: ', year0 + relyr tmp = year0 + accel_yr0 + (relyr - accel_yr0)*accel print*, ' current accelerated year: ', tmp print*, ' ' print*, '==> Warning: adjusting output intervals for' &, ' long term forcing acceleration' print*, ' ' if (tsiint .ge. yrlen*accel) tsiint = tsiint/accel if (tavgint .ge. yrlen*accel) tavgint = tavgint/accel if (tmbint .ge. yrlen*accel) tmbint = tmbint/accel if (stabint .ge. yrlen*accel) stabint = stabint/accel if (zmbcint .ge. yrlen*accel) zmbcint = zmbcint/accel if (glenint .ge. yrlen*accel) glenint = glenint/accel if (trmbint .ge. yrlen*accel) trmbint = trmbint/accel if (vmsfint .ge. yrlen*accel) vmsfint = vmsfint/accel if (gyreint .ge. yrlen*accel) gyreint = gyreint/accel if (extint .ge. yrlen*accel) extint = extint/accel if (prxzint .ge. yrlen*accel) prxzint = prxzint/accel if (exconvint .ge. yrlen*accel) exconvint = exconvint/accel if (dspint .ge. yrlen*accel) dspint = dspint/accel if (timavgint .ge. yrlen*accel) timavgint = timavgint/accel if (cmixint .ge. yrlen*accel) cmixint = cmixint/accel if (xbtint .ge. yrlen*accel) xbtint = xbtint/accel if (crossint .ge. yrlen*accel) crossint = crossint/accel if (densityint .ge. yrlen*accel) densityint = densityint/accel if (fctint .ge. yrlen*accel) fctint = fctint/accel if (tyzint .ge. yrlen*accel) tyzint = tyzint/accel if (restint .ge. yrlen*accel) restint = restint/accel if (tbtint .ge. yrlen*accel) tbtint = tbtint/accel endif if (timavgint .lt. timavgper) then write (stdout,'(/,(1x,a))') & '==> Warning: averaging period "timavgper" exceeds interval ' &, ' "timavgint". Setting timavgper = timavgint ' timavgper = timavgint endif write (stdout,'(/,1x,a,f10.2,a,/,1x,a,f10.2,a)') & '==> Time averages written every ', timavgint, ' days, ' &, ' with an averaging period of ', timavgper, ' days. ' if (tsiint .lt. tsiper) then write (stdout,'(/,(1x,a))') & '==> Warning: averaging period "tsiper" exceeds interval ' &, ' "tsiint". Setting tsiper = tsiint ' tsiper = tsiint endif write (stdout,'(/,1x,a,f10.2,a,/,1x,a,f10.2,a)') & '==> Time step integrals written every ', tsiint, ' days, ' &, ' with an averaging period of ', tsiper, ' days. ' if (mod(namix,2) .ne. mod(nint(segtim*daylen/dtatm),2)) then write(*,*) '==> Error: time steps between mixing and coupling' write(*,*) ' in the atmosphere must both be even to' write(*,*) ' use even_fluxes.' stop '=>UVic_ESCM' endif if (mod(nmix,2) .ne. mod(nint(segtim*daylen/dtocn),2)) then write(*,*) '==> Error: time steps between mixing and coupling' write(*,*) ' in the ocean must both be even to use' write(*,*) ' even_fluxes.' stop '=>UVic_ESCM' endif if (mod(namix,2) .ne. mod(nats,2)) then write(*,*) '==> Warning: restart was not saved with even flux' write(*,*) ' averaging. Starting with a mixing time' write(*,*) ' step in atmosphere.' nats = namix endif if (mod(nmix,2) .ne. mod(nots,2)) then write(*,*) '==> Warning: restart was not saved with even flux' write(*,*) ' averaging. Starting with a mixing time' write(*,*) ' step in ocean.' nots = nmix endif !----------------------------------------------------------------------- ! do consistency checks before allowing model to continue !----------------------------------------------------------------------- errorc = .false. write (stdout,*) ' ' write (stdout,*) ' (checking setup)' if (dtatm .eq. c0) then write (stdout,9000) & '==> Error: the atmospheric time step must be set in "setatm"' errorc = .true. dtatm = 1.e-6 endif ! critv = 1.e-6 critv = 1.e-4 if (segtim .ne. c0) then r1 = rundays/segtim else r1 = 0.5 endif if (segtim .eq. c0) then write (stdout,9000) & '==> Error: coupling period "segtim" must be specified' errorc = .true. endif if (abs(r1-nint(r1)) .gt. critv) then write (stdout,9000) & '==> Error: there must be an integral number of segments ' &,'"segtim" within "rundays" (the length of the run)' errorc = .true. endif r2 = segtim/(dtocn*secday) if (abs(r2-nint(r2)) .gt. critv) then write (stdout,9000) & '==> Error: there must be an integral number of density time ' &,'steps "dtocn" within "segtim" (the segment time)' errorc = .true. endif r3 = segtim/(dtatm*secday) if (abs(r3-nint(r3)) .gt. critv) then write (stdout,9000) & '==> Error: there must be an integral number of atmos time ' &,'steps "dtatm" within "segtim" (the segment time)' errorc = .true. endif r5 = segtim/(dtlnd*secday) if (abs(r5-nint(r5)) .gt. critv) then write (stdout,9000) & '==> Error: there must be an integral number of land time ' &,'steps "dtlnd" within "segtim" (the segment time)' errorc = .true. endif r6 = (dtsed*secday)/segtim if (abs(r6-nint(r6)) .gt. critv) then write (stdout,9000) & '==> Error: there must be an integral number of "segtim" steps ' &,'within "dtsed"' errorc = .true. endif tmp = 0 & + 1 if (tmp .gt. 1) then write (stdout,9000) & '==> Error: use only one of O_co2ccn_data, O_co2emit_data, ' &,'or O_co2ccn_user' errorc = .true. endif write (stdout,*) ' ' write (stdout,*) ' (End of checks) ' write (stdout,*) ' ' if (errorc) stop '=> ERRORS found in chkcpl' 9000 format (/,(1x,a)) return end 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 getst (jrow, ocnout, ntabc) !----------------------------------------------------------------------- ! read surface tracers from disk row "jrow" !----------------------------------------------------------------------- implicit none integer i, jrow, ntabc include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "iounit.h" include "mw.h" include "tmngr.h" real ocnout(imt,jmt) call getrow (latdisk(taup1disk), nslab, jrow &, u(1,1,jmw,1,taup1), t(1,1,jmw,1,taup1)) do i=1,imt ocnout(i,jrow) = t(i,1,jmw,ntabc,taup1) enddo return end subroutine sbc_init !----------------------------------------------------------------------- ! Initialize S.B.C. indices !----------------------------------------------------------------------- implicit none integer m include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "csbc.h" 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 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 (iwa, m, mapsbc(m), 'wa', 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 (issc14, m, mapsbc(m), 'ssc14', m) call set (ic14flx, m, mapsbc(m), 'c14flx', 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) call set (issno3, m, mapsbc(m), 'ssno3', m) call set (ino3flx, m, mapsbc(m), 'no3flx', m) call set (iat, m, mapsbc(m), 'at', m) call set (irh, m, mapsbc(m), 'rh', m) call set (ipr, m, mapsbc(m), 'pr', m) call set (ips, m, mapsbc(m), 'ps', m) call set (iaws, m, mapsbc(m), 'aws', m) call set (iswr, m, mapsbc(m), 'swr', m) call set (ilwr, m, mapsbc(m), 'lwr', m) call set (isens, m, mapsbc(m), 'sens', m) call set (ievap, m, mapsbc(m), 'evap', m) call set (idtr, m, mapsbc(m), 'dtr', m) call set (isr, m, mapsbc(m), 'sr', m) call set (inpp, m, mapsbc(m), 'npp', m) call set (iburn, m, mapsbc(m), 'burn', m) call set (ircal, m, mapsbc(m), 'rcal', m) call set (irorg, m, mapsbc(m), 'rorg', m) call set (ibtemp, m, mapsbc(m), 'btemp', m) call set (ibsalt, m, mapsbc(m), 'bsalt', m) call set (ibo2, m, mapsbc(m), 'bo2', m) call set (ibalk, m, mapsbc(m), 'balk', m) call set (ibdic, m, mapsbc(m), 'bdic', m) call set (ibdicfx, m, mapsbc(m), 'bdicfx', m) call set (ibalkfx, m, mapsbc(m), 'balkfx', m) 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 "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "atm.h" include "mw.h" !----------------------------------------------------------------------- ! Initialize ocean tracer names !----------------------------------------------------------------------- mapt(:) = " " itemp = 0 isalt = 0 idic = 0 ic14 = 0 icfc11 = 0 icfc12 = 0 io2 = 0 ialk = 0 ipo4 = 0 ino3 = 0 iphyt = 0 izoop = 0 idetr = 0 idiaz = 0 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 (ic14, m, mapt(m), 'c14', m) call set (io2, m, mapt(m), 'o2', m) call set (ialk, m, mapt(m), 'alk', m) call set (ipo4, m, mapt(m), 'po4', m) 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) 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 (isc14, m, mapst(m), 'sc14', m) itrc(ic14) = m-1 call set (isdic, m, mapst(m), 'sdic', m) itrc(idic) = m-1 call set (iso2, m, mapst(m), 'so2', m) itrc(io2) = m-1 call set (isalk, m, mapst(m), 'salk', m) itrc(ialk) = m-1 call set (ispo4, m, mapst(m), 'spo4', m) itrc(ipo4) = m-1 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 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 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "calendar.h" include "cembm.h" include "coord.h" include "cnep.h" include "csbc.h" include "cprnts.h" include "diag.h" include "emode.h" include "fwa.h" include "hmixc.h" include "insolation.h" include "iounit.h" include "isopyc.h" include "mtlm.h" include "npzd.h" include "scalar.h" include "sed.h" include "stab.h" include "switch.h" include "tmngr.h" include "veg.h" include "vmixc.h" include "transpiration_sens.h" include "CO2_fertilization_scaling.h" !----------------------------------------------------------------------- ! read all model namelist variables !----------------------------------------------------------------------- namelist /contrl/ init, runlen, rununits, restrt, initpt &, num_processors, runstep 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 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 &, aero_scal namelist /carbon/ co2ccn, co2emit, co2for, co2_yr, c14_yr, dc14ccn &, c14prod, co2pi 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, sens_fact, co2_fert_scal namelist /pfrost/ HCON_D_P, HCON_SNOW, KAPSp, KOVCON, KOVSAT &, PPtm, iAF ! 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 ! 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 kfe = 0.1 ! Half saturation constant for Fe limitation kfe_D = 0.1 ! Half saturation constant for Diaz Fe limitation ! 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 aero_scal = 1. ! set defaults namelist carbon co2ccn = 280. co2emit = 0. co2for = 5.35e03 co2_yr = 2.e20 c14_yr = 2.e20 dc14ccn = 0. c14prod = 0. co2pi = 278.05158 ! 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. sens_fact = 1 co2_fert_scal = 1 ! set defaults namelist pfrost ! 1 defines HCON_D_P and HCON_SNOW here HCON_D_P = 0.05 HCON_SNOW = 0.25 KAPSp = 1.70E-8 KOVCON = 0.001 KOVSAT = 0.2587 PPtm = 8.8E-11 iAF = 0.154 !----------------------------------------------------------------------- ! 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) call getunit (ioun, 'control.in', 'f s r') read (ioun, pfrost, end=123) 123 continue write (stdout, pfrost) 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 (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