subroutine mom c c======================================================================= c c GFDL Modular Ocean Model c c A primitive equation ocean model developed by researchers at c the Geophysical Fluid Dynamics Laboratory /NOAA. c c For details, refer to the manual. c c author: r.c.pacanowski e-mail rcp@gfdl.gov c (mpp items) v. balaji e-mail=> vb@gfdl.gov c======================================================================= c use arrays_2d use datamod use memory_window use asm_x_mod use mom3_mpp_mod c external blkdta c include 'mpif.h' !====================== include file "param.h" ========================= ! !====================== include file "size.h" ========================== ! !----------------------------------------------------------------------- ! USER INPUT: !----------------------------------------------------------------------- ! ! imt = number of grid cells in the longitudinal direction ! (calculated cells are from 2 through imt-1. end cells ! are boundaries) ! ! jmt = number of grid cells (latitude rows) in the latitudinal ! direction (calculated cells are from 2 through jmt-1. ! end cells are boundaries) ! ! km = number of grid cells in the vertical direction ! (calculated cells are from 1 through km. no boundary ! cells) ! ! nt = number of tracers (temperature, salinity, ...) ! ! mnisle = maximum number of islands (unconnected land masses) ! ! maxipp = maximum number of all island perimeter points !----------------------------------------------------------------------- ! ! window size must be appropriate for options with highest order numerics ! specify lomem or himem on 3D diagnostics. default is himem !these are used to determine hardware/OS/compiler (V. Balaji vb@gfdl.gov) !compilers that support Cray pointers !values of kind: double and long are 8-byte, float and int are 4-byte !these might be different on non-SGICRAY, I believe !parallel machine types parameter (imt= 362, jmt= 202, km= 40) parameter (nt=2) parameter (mnisle=50, maxipp=10000) parameter (kbot = km) ! ! !----------------------------------------------------------------------- ! END USER INPUT !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! "jmw" = number of latitude rows within the MW (memory window). ! "jbuf" = number of northermost and southernmost buffer rows. ! The MW is symmetric because it has the same number of ! northermost and southernmost buffer rows "jbuf". ! total number of buffer rows = 2*jbuf !----------------------------------------------------------------------- ! parameter (jbuf=2) parameter (jmw=1 + 2*jbuf) ! !----------------------------------------------------------------------- ! num_rows_to_copy = number of northernmost rows to copy into ! the southernmost rows of the MW when moving ! the MW northward !----------------------------------------------------------------------- ! parameter (num_rows_to_copy = 2*jbuf) ! ! 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" isn`t used for tracer averages ! numreg = total number of regions ( = product of nhreg & nvreg) ! used for term balance calculations ! parameter (lseg=5, nlatpr=10) parameter (nhreg=4, nvreg=6, numreg=nhreg*nvreg) parameter (nvar=nt+2) ! parameter (imtm1=imt-1, imtm2=imt-2) parameter (kmm1=km-1, kmp1=km+1) parameter (jmtm1=jmt-1, jmtm2=jmt-2) parameter (imtkm=imt*km, nwds=imt*jmt, nkflds=2) parameter (nslab=imt*nvar*km, ntmin2=nt+1/nt) ! ! parameter constants ! !====================== 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) ! 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)) ! minimum allowable thickness (cm) of the top cell when using the ! nonlinear free surface. parameter (thickmin=10.0) ! ! standard I/O units ! !====================== 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, stderr = 0) common /stdunit/ stdout c====================== include file "emode.h" ========================= c c c mxscan = max number of allowable scans for poisson solvers c mscan = actual number of scans taken by poisson solvers c tolrsf = tolerance for stream function calculation. c the solution is halted when it is within "tolrsf" c of the "true" solution assuming geometric convergence. c tolrsp = tolerance for surface pressure calculation c the solution is halted when it is within "tolrsp" c of the "true" solution assuming geometric convergence. c tolrfs = tolerance for implicit free surface calculation c the solution is halted when it is within "tolrfs" c of the "true" solution assuming geometric convergence. c esterr = estimated maximum error in elliptic solver assuming c geometric convergence c c c nisle = number of land masses c nippts= number of land mass perimeter points c iperm = "i" coordinate for the land mass perimeter point c jperm = "j" coordinate for the land mass perimeter point c iofs = offset for indexing into the land mass perimeter points c imask = controls whether calculations get done on perimeters c set mask for land mass perimeters on which to perform calculations c imask(-n) = .false. [no equations ever on dry land mass n] c imask(0) = .true. [equations at all mid ocean points] c imask(n) = .true./.false [controls whether there will be c equations on the ocean perimeter of c land mass n] c note: land mass 1 is the northwest-most land mass c for the numbering of the other landmasses, see generated map(i,j) logical imask common /emodei/ nippts(mnisle), iofs(mnisle), iperm(maxipp) common /emodei/ jperm(maxipp), nisle, imain common /emodel/ imask (-mnisle:mnisle) c character*16 variable logical converged common /emodel/ converged common /emodec/ variable c c====================== include file "iounit.h" ======================== c c i/o units and related variables c c taum1disk = disk pointer for tau-1 latitude rows c taudisk = disk pointer for tau latitude rows c taup1disk = disk pointer for tau+1 latitude rows c kflds = disk unit used for two dimensional fields c latdisk = disk units for latitude rows (alternately pointed to c by taum1disk, taudisk, and taup1disk) c c iodoc = unit for documentation c iostab = unit for stability testing c iotim = unit for time means c iotim1 = scratch disk (SSD) unit for accumulating time means c ionew1 = unit for reading sponge layer data c ionew2 = mirror unit of sponge layer data on SSD c c for the following, a control # < 0 implies that unformatted data c will be written to a unit selected by the i/o manager "iomngr.F" c and given a hardwired name (grep getunit *.F to see names) c and formatted data (to stdout) will be written. if a # > 0 and c # <> stdout, only unformatted data will be written. c c iotavg = control # for tracer averages c iotmb = control # for writing tracer meridional budget. c iotrmb = control # for term balances for tracer and momentum c ioglen = control # for writing global energetics integrals c iovmsf = control # for writing meridional stream function c iogyre = control # for writing gyre transport. c ioprxz = control # for writing x-z sections from latitudes c ioext = control # for writing external mode (stream function) c iodsp = control # for writing diagnostic surface pressure c iotsi = control # for writing time step integrals c ioxbt = control # for writing time averaged xbt data c iozmbc = control # for writing zonal mean surf boundary conditions c integer taum1disk, taudisk, taup1disk common /iounit/ taum1disk, taudisk, taup1disk common /iounit/ latdisk(3), kflds common /iounit/ iodoc, iostab, iotavg, iotmb, iotrmb common /iounit/ iotim, iotim1 common /iounit/ ioglen, iovmsf, iogyre, ioprxz, ioext, iodsp common /iounit/ iotsi, iozmbc, ionew1, ionew2, ioxbt c c iotext = 80 character text string for describing the details c of the next unformatted data record. c expnam = 60 character text string for the experiment name c character iotext*80, expnam*60 common /iounitc/ iotext, expnam c c when writing unformatted data records in MOM, each data record is c preceded by a header record which was written as: c write(unit) stamp, iotext, expnam c where stamp is a 32 character specification of the model date & c time corresponding to the time step when the data was written and c iotext is a 80 character description of what is in the c data record and how it is to be read. expnam is a 60 character c experiment name which shows which experiment wrote the data. c this makes it easy to decipher any unformatted output from the c model by using a program similar to the following: c c program decifr c c----------------------------------------------------------------------- c decifer an unformatted file from MOM by showing the header c records. the file needs to copied to file "fort.21" c----------------------------------------------------------------------- c c character*32 stamp c character*80 iotext c character*60 expnam c c iounit = 21 c rewind iounit c do n=1,100000 c c read the header record c c read (iounit, end=110) stamp, iotext, expnam c write (*,'(1x,a32,1x,a80)') stamp, iotext c c skip the data record c c read (iounit) c enddo c110 continue c write (*,*) " => end of file on fort.",iounit c stop c end c c note: all unformatted diagnostic MOM data is handled this way. c to insure that data is read properly, verify that arrays are c dimensioned correctly by comparing the listed variables against c those in the *.h files. (grep -i -n "variable" *.h) Also, most c data from MOM is written IEEE 32bit so it is read directly by c most workstations. However, when trying to read these IEEE files c on the CRAY, they must be assigned IEEE before being read. c Some diagnostic data is averaged over time before being written. c In these cases, the time "stamp" refers to the last time step c at the end of the averaging period. An averaging interval is c also written as part of the data. Averaging periods = zero c indicate instantaneous data. c c======================= include file "mw.h" =========================== c c Refer to the manual for a description of the memory window (MW) c c author: r.c.pacanowski e-mail rcp@gfdl.gov c======================================================================= c c taum1 = tau-1 time level for variables in MW c tau = tau time level for variables in MW c taup1 = tau+1 time level for variables in MW c integer taum1, tau, taup1 common /mwi/ taum1, tau, taup1 c c----------------------------------------------------------------------- c MW arrays for prognostic equations: c----------------------------------------------------------------------- c c u(i,k,j,n,tau) = total velocity where: c i = index for longitude c k = index for depth c j = index for latitude row within MW c n = component (1 = zonal, 2 = meridional) c tau = time level (tau-1, tau, tau+1) c c (only internal modes are on disk and at tau+1 in the MW) c c t(i,k,j,n,tau) = tracer where: c i = index for longitude c k = index for depth c j = index for latitude row within MW c n = component (1 = temperature, 2 = salinity) c if nt > 2 then other tracers are allowed. c tau = time level (tau-1, tau, tau+1) c c note: temperature is potential temperature in degrees Celsius and c salinity is in "model units", the deviation from 0.035 grams c of salt/cm**3 of water, or, assuming a water density of c 1 gram/cm**3, the deviation from 0.035 g of salt/g of water. c one can convert model units to the common units of parts c per thousand (ppt) by adding 0.035 grams/cm**3 to the model c units and then multiplying by 1000. c common /mw/ u(imt,km,jmw,2,-1:1), t(imt,km,jmw,nt,-1:1) c c----------------------------------------------------------------------- c MW arrays for diagnostic equations and workspace: c----------------------------------------------------------------------- c c diagnostic advective velocities are in units of cm/sec c c adv_vet = advective velocity on the eastern face of a T cell c adv_vnt = advective velocity on the northern face of a T cell c adv_veu = advective velocity on the eastern face of a U cell c adv_vnu = advective velocity on the northern face of a U cell c adv_vbt = advective velocity on the bottom face of a T cell c adv_vbu = advective velocity on the bottom face of a U cell c c rho = density at center of a T cell in units of gm/cm**3 c note: there is an arbitrary constant which is only a c function of depth in "rho". It is related to c subtracting a reference level density for purposes of c accuracy. c c pressure = hydrostatic pressure at the grid point within a T cell c grad_p = hydrostatic pressure gradient for U cell. There are c two components: (1,2) is for (dp/dx, dp/dy) c common /mw/ adv_vet(imt,km,2:jmw), adv_vnt(imt,km,1:jmw) common /mw/ adv_veu(imt,km,2:jmw-1) common /mw/ adv_vnu(imt,km,1:jmw-1) common /mw/ adv_vbt(imt,0:km,2:jmw) common /mw/ adv_vbu(imt,0:km,2:jmw-1) common /mw/ rho(imt,km,1:jmw) common /mw/ pressure(imt,km,1:jmw) common /mw/ grad_p(imt,km,1:jmw-1,2) c c tmask = tracer cell land/sea mask = (0.0, 1.0) on (land, sea) c umask = velocity cell land/sea mask = (0.0, 1.0) on (land, sea) c common /mw/ tmask(imt,km,1:jmw), umask(imt,km,1:jmw) c common /mw/ dht(imt,km,jmw+1), dhu(imt,km,jmw) common /mw/ dhwt(imt,0:km,jmw) common /mw/ pc_sink(imt,km,2:jmw-1) common /mw/ ztp(imt,km,jmw) c c c adv_f4n = 4th order advective flux c common /mw/ adv_f4n(imt,km,1:jmw-1,nt) c c c c these workspace arrays are recalculated for each component of the c equations so do not have to be moved as the MW moves northward. c c adv_fe = advective flux across the eastern face of a cell c adv_fn = advective flux across the northern face of a cell c (removed in most cases and put directly into the c statement functions for speed optimization.) c adv_fb = advective flux across the bottom face of a cell c c diff_fe = diffusive flux across the eastern face of a cell c diff_fn = diffusive flux across the northern face of a cell c diff_fb = diffusive flux across the bottom face of a cell c source = source term c common /mw/ adv_fe(imt,km,2:jmw-1) common /mw/ adv_fn(imt,km,1:jmw-1) common /mw/ adv_fb(imt,0:km,2:jmw-1) c common /mw/ diff_fe(imt,km,2:jmw-1) common /mw/ diff_fn(imt,km,1:jmw-1) common /mw/ diff_fb(imt,0:km,2:jmw-1) common /mw/ diff_fbiso(imt,0:km,2:jmw-1) c common /mw/ source(imt,km,2:jmw-1) c common /mw/ zzi(imt,km,2:jmw-1) c c these grid factors are for optimizations c common /mw/ cstdxtr(imt,jmw), cstdxur(imt,jmw) common /mw/ cstdxt2r(imt,jmw),ah_cstdxur(imt,jmw) common /mw/ csudxtr(imt,jmw), csudxur(imt,jmw) common /mw/ csudxu2r(imt,jmw),am_csudxtr(imt,jmw) c c advmet = coeff for metric advection. c common /advec/ advmet(jmt,2) c c c c smf = surface momentum flux c 1 => zonal wind stress (dynes/cm**2) c 2 => meridional wind stress (dynes/cm**2) c bmf = bottom momentum flux c 1 => zonal bottom drag (dynes/cm**2) c 2 => meridional bottom drag (dynes/cm**2) c stf = surface tracer flux c 1 => surface heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec) c 2 => surface salt flux (grams of salt/cm**2/sec) c stf_turb ! = the turbulent part of surface tracer flux c btf = bottom tracer flux (for consistency but normally zero!) c 1 => bottom heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec) c 2 => bottom salt flux (grams of salt/cm**2/sec) c common /mw/ smf(imt,1:jmw,2), bmf(imt,1:jmw,2) common /mw/ stf(imt,1:jmw,nt), btf(imt,1:jmw,nt) common /mw/ stf_turb(imt,1:jmw,nt) c c c====================== include file "scalar.h" ======================== c c various scalar quantities: c c dtts = time step for density & tracers (in seconds) c dtuv = time step for baroclinic velocity (in seconds) c dtsf = time step for barotropic velocity (in seconds) c c2dtts = 2*dtts c c2dtuv = 2*dtuv c c2dtsf = 2*dtsf c smooth = smoothing coefficient for robert time filter c acor = (>0, 0) = (implicit, explicit) treatment of coriolis c term for internal and external modes. c rho0 = mean density for Boussinesq approximation c rho0r = 1/rho0 c omega = earth`s rotation rate (radians/sec) c radius = earth`s radius (cm) c grav = earth`s gravitational acceleration (cm/sec**2) c cdbot = bottom drag coefficient c ncon = number of passes through convective code in tracer c gcor = time centering for coriolis term c c taux0 = constant zonal windstress (dynes/cm**2) for idealized c equatorial studies c tauy0 = constant meridional windstress (dynes/cm**2) for c idealized equatorial studies c common /scalar/ dtts, dtuv, dtsf, c2dtts, c2dtuv, c2dtsf, acor common /scalar/ smooth common /scalar/ rho0, rho0r, omega, radius, grav, cdbot, gcor common /scalar/ taux0, tauy0 common /scalri/ ncon c c bring in non dimensional constants c c====================== include file "ndcon.h" ======================== c c various non dimensional quantities: c c radian = degrees per radian c pi = something good to eat c common /ndcon/ radian, pi c c !====================== 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 centered leapfrog. ! ! init = (true,false) indicates that this run is a ! (start from initial 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 ! !----------------------------------------------------------------------- ! 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 its diagnostic ! setting a switch = 0.0 causes the diagnostic to be ! done every time step, and setting a switch > 0.0 causes the ! diagnostic to be done repeatedly at the specified interval ! from a specified starting date. ! ! 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 ! ! pressint = number of days between writing pressure field ! ! 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. ! ! rhoint = number of days between writing instantaneous terms ! in the locally referenced potential density equation. ! Relevant only when ! "save_density_terms" ! is enabled. see "iounit.h" for more details. ! all data is at time level "tau" so it is one time step ! before the date shown in "stamp". ! ! diffint = number of days between computing nonconstant ! diffusivity. ! ! diffsnapint = number of days between saving nonconstant ! diffusivity information to netcdf file. ! ! 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. ! ! stabint = number of days between sampling for various stability ! criteria. ! ! snapint = number of days between writing instantaneous samples ! of data from the MOM grid. note: only when "snapshots" ! is enabled. see "iounit.h" for more details. ! all data is at time level "tau" so it is one time step ! before the date shown in "stamp". ! ! timavgint= interval (days) for writing time mean data ! (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) ! ! trajint = number of days between writing particle trajectories. ! particle trajectories are always written on the last ! time step of the run to allow for restarting. ! !----------------------------------------------------------------------- ! outputs from tmngr.F: logical switches !----------------------------------------------------------------------- ! ! rundays = integration time in days (from "runlen") ! ! the following are logical counterparts to the above switches and ! are set 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". ! ! pressts = (false,true) = (don`t, do) write pressure on this time ! step. based on "pressint". ! ! ! 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 diagnosti! ! surface pressure on this time step. based on "dspint" ! ! rhots = (false,true) =(don`t, do) do calculation of diagnostic ! for locally referenced potential density terms. ! Based on "rhoint". ! ! diffts = (false,true) =(don`t, do) do calculation of new ! nonconstant diffusivity. ! Based on "diffint". ! ! diffsnapts = (false,true) =(don`t, do) save fields related ! to nonconstant diffusivity calculation. ! Based on "diffsnapint". ! ! 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" ! ! trajts = (false,true) = (don`t, do) print particle trajectories ! on this time step. based on "trajint" ! ! snapts = (false,true) = (don`t, do) save a snapshot of the data ! on this time step. based on "snapint" ! ! timats = (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" ! ! ! 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 ! and not on a specified interval in days. ! ! 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 ! 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. ! ! author: r. c. pacanowski e-mail=> rcp@gfdl.gov ! ! logical eb, leapfrog, euler1, euler2, forward, eots logical init, first, restrt logical itavg, itmb, itrmb, igyre character *8 rununits ! common /switcc/ rununits common /switcr/ runlen, rundays common /switci/ nmix common /switcl/ eb, leapfrog, euler1, euler2, forward, eots common /switcl/ init, first, restrt common /switcl/ itavg, itmb, itrmb, igyre ! logical eoday, eomon, midmon, eoyear, eorun common /switcl/ eoday, eomon, midmon, eoyear, eorun ! c====================== include file "casimsw.h" ======================== c c holds switching information for assimilation c defaults and namelist settings are handled in setassim.F c it is "included" in switch.h c c rassim = true when within 1/2 time step of a specified interval c else ... false c rassint = interval in days for triggering a new assimilation c logical rassim common /switcr/ rassint common /switcl/ rassim c ! !----------------------------------------------------------------------- ! ! 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 two variables in common. The ! following naming convention is used. ! ! 1) an interval (real) for diagnosti! output (e.g,. glenint) ! 2) a switch (logical) for the interval (e.g., glents ) ! !----------------------------------------------------------------------- ! logical tavgts, glents, trmbts, prxzts common /switcr/ tavgint, glenint, trmbint, prxzint common /switcl/ tavgts, glents, trmbts, prxzts ! logical extts, trajts, exconvts, cmixts common /switcr/ extint, trajint, exconvint, cmixint common /switcl/ extts, trajts, exconvts, cmixts ! logical vmsfts, gyrets, tsits, snapts, fctts common /switcr/ vmsfint, gyreint, tsiint, snapint, fctint common /switcl/ vmsfts, gyrets, tsits, snapts, fctts ! logical tyzts common /switcr/ tyzint common /switcl/ tyzts ! logical stabts, zmbcts, crossts, pressts common /switcr/ stabint, zmbcint, crossint, pressint common /switcl/ stabts, zmbcts, crossts, pressts ! logical osegs, osege common /switcr/ segtim common /switcl/ osegs, osege ! logical rhots common /switcr/ rhoint common /switcl/ rhots ! logical diffts, diffsnapts common /switcr/ diffint, diffsnapint common /switcl/ diffts, diffsnapts ! !----------------------------------------------------------------------- ! ! 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 four 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 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. ! !----------------------------------------------------------------------- ! logical xbtts, dspts, tmbts, timavgts logical xbtperts, dspperts, tmbperts, timavgperts common /switcr/ xbtint, dspint, tmbint, timavgint common /switcr/ xbtper, dspper, tmbper, timavgper common /switcl/ xbtts, dspts, tmbts, timavgts common /switcl/ xbtperts, dspperts, tmbperts, timavgperts !====================== include file "taskrows.h" ====================== ! ! max_tasks = maximum number of processors ! num_processors = requested number of processors ! pn = this processor number (from 1 to num_processors) ! ! jstask3(n) = southernmost buffer jrow for processor "n" ! jetask3(n) = northernmost buffer jrow for processor "n" ! jscomp3(n) = southernmost computed jrow for processor "n" ! jecomp3(n) = northernmost computed jrow for processor "n" ! num_loads3(n)= number of times processor "n`s" MW is loaded ! window_size3(n)= the size of the memory window for processor "n" ! ! jstask = southernmost buffer jrow for this processor (pn) ! jetask = northernmost buffer jrow for this processor (pn) ! jscomp = southernmost computed jrow for this processor (pn) ! jecomp = northernmost computed jrow for this processor (pn) ! num_loads= number of times processor "pn`s" MW is loaded ! window_size= the size of the memory window for this processor (pn) ! It is set as follows: ! if defined max_window then ! jecomp3-jscomp3+1 + 2*jbuf ! for second, fourth, and sixth order windows ! on processor "pn". The window_size must be the ! same on all processors. ! else ! 3 (for second order windows) ! 5 (for fourth order windows) ! 7 (for sixth order windows) ! regardless of the number of processors ! endif ! parameter (max_tasks = 2048) integer pn, window_size3, window_size common /taskrows/ pn common /taskrows/ num_processors common /taskrows/ jstask3(max_tasks), jetask3(max_tasks) common /taskrows/ jscomp3(max_tasks), jecomp3(max_tasks) common /taskrows/ num_loads3(max_tasks) common /taskrows/ jstask, jetask, num_loads, jscomp, jecomp common /taskrows/ window_size3(max_tasks), window_size !====================== include file "tmngr.h" ========================= ! ! time manager variables ! ! time manager inputs: ! ! refrun = T ==> logical switches referenced to the start of run ! refinit = T ==> logical switches referenced to initial condition ! time: year0, month0, day0, hour0, min0, sec0 ! refuser = T ==> logical switches referenced to user specified ! time: ryear, rmonth, rday, rhour, rmin, rsec ! ! 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 ! ic_stamp= 32 character date and time for initial conditions ! stop_stamp = 32 character date and time for last time step ! ! 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 ! !----------------------------------------------------------------------- ! ! logical refrun, refinit, refuser, eqmon, eqyear ! integer year0, month0, day0, hour0, min0, sec0 integer ryear, rmonth, rday, rhour, rmin, rsec character *32 pstamp, stamp, ic_stamp, stop_stamp, xstamp ! common /tmngrl/ refrun, refinit, refuser, eqmon, eqyear ! common /tmngrii/ itt0, itt, irstdy, msrsdy ! common /tmngrr/ dayoyr, relyr, prelyr, relday ! common /tmngrc/ stamp, pstamp, ic_stamp, stop_stamp, xstamp ! common /tmngrii/ year0, month0, day0, hour0, min0, sec0 common /tmngrii/ ryear, rmonth, rday, rhour, rmin, rsec ! character*10 dayname character*12 monname common /tmngrc/ dayname(7), monname(12) ! real yrlen common /tmngrr/ yrlen ! character(len=80) :: timeunits common /tmngrc/ timeunits save cps, itt1st c logical goto_2nd_pass c c----------------------------------------------------------------------- c integrate equations for one time step c----------------------------------------------------------------------- c c----------------------------------------------------------------------- c update timestep counter, set corresponding model time, and set c time dependent logical switches which determine program flow. c----------------------------------------------------------------------- c itt = itt + 1 call tmngr (dtts) call set_time_switches c write(stdout,'(a,i7,1x,a32)') ' ts=',itt,stamp call flush(stdout) c c----------------------------------------------------------------------- c read system clock to calculate compute_time/grid_point/timestep c (neglecting overhead of "setocn") c----------------------------------------------------------------------- c if (first) then itt1st = itt cps = cpusec() endif c c----------------------------------------------------------------------- c initialize diagnostic variables c----------------------------------------------------------------------- c call diag_initial (jscomp, jecomp) c c----------------------------------------------------------------------- c adjust various quantities for leapfrog/mixing timesteps c c leapfrog----------> h(tau+1) = h(tau-1) + 2dt*F(tau) c c forward-----------> tau-1 <= tau c h(tau+1) = h(tau-1) + dt*F(tau) c c euler backward: tau-1 <= tau c euler1----------> h(tau` ) = h(tau-1) + dt*F(tau) c euler2----------> h(tau+1) = h(tau-1) + dt*F(tau`) c----------------------------------------------------------------------- c if (leapfrog) then c c normal leapfrog time step c euler1 = .false. euler2 = .false. forward = .false. eots = .true. c c2dtts = c2*dtts c2dtuv = c2*dtuv c2dtsf = c2*dtsf else c c mixing time step (forward step or euler backward step) c if (eb) then euler1 = .true. euler2 = .false. forward = .false. eots = .false. else euler1 = .false. euler2 = .false. forward = .true. eots = .true. endif c c2dtts = dtts c2dtuv = dtuv c2dtsf = dtsf c endif c if (.not.leapfrog) then do jrow=jstask,jetask do i=1,imt ubarm1(i,jrow,1) = ubar(i,jrow,1) ubarm1(i,jrow,2) = ubar(i,jrow,2) enddo enddo endif c c c----------------------------------------------------------------------- c set time centering "gcor" for coriolis term c----------------------------------------------------------------------- c if (acor .eq. c0) then gcor = c1 elseif (acor .ne. c0) then gcor = c0 endif c c c----------------------------------------------------------------------- c update pointers to tau-1, tau, & tau+1 data on disk. c for latitude rows they point to latdisk c for 2D fields they point to records on kflds c----------------------------------------------------------------------- c taum1disk = mod(itt+1,2) + 1 taudisk = mod(itt ,2) + 1 taup1disk = taum1disk c c----------------------------------------------------------------------- c update pointers (indices) to tau-1, tau, & tau+1 data in the MW c----------------------------------------------------------------------- c c c time levels are being held constant in time. c c c======================================================================= c c SOLVE THE BAROCLINIC AND TRACER EQUATIONS c c Latitude rows are solved using a MW (memory window) approach. c Refer to the manual for details. c c======================================================================= c goto_2nd_pass = .true. do while (goto_2nd_pass) c c c----------------------------------------------------------------------- c Load the processor`s memory window in preparation for solving c tracer and baroclinic equations. (See the manual) c----------------------------------------------------------------------- c do mw=1,num_loads c ncrows = window_size-2*jbuf joff = jscomp - jbuf - 1 + (mw-1)*ncrows js = (jstask - joff)*(1/mw) + (1 + num_rows_to_copy)*(1-1/mw) je = min(window_size,jetask-joff) c is = 2 ie = imt - 1 c c set rows within the MW for solving prognostic equations c (rows are clipped at jrow=2 and jrow=jmt-1) c js_calc = 1+jbuf je_calc = ncrows+jbuf c c----------------------------------------------------------------------- c load prognostic and related variables into the MW c----------------------------------------------------------------------- c call loadmw (joff, js, je, is, ie, latdisk(taum1disk) &, latdisk(taudisk)) c c----------------------------------------------------------------------- c add correction from assimilation c----------------------------------------------------------------------- c ! call add_correction (joff, js_calc, je_calc) call add_correction (joff, js, je) c----------------------------------------------------------------------- c compute pressure gradients c----------------------------------------------------------------------- c call pressure_gradients (joff, js, je) c c----------------------------------------------------------------------- c calculate advection velocities for T cells and U cells c----------------------------------------------------------------------- c call adv_vel (joff, js, je, is, ie) c c----------------------------------------------------------------------- c set horizontal mixing coefficients for momentum and tracers c----------------------------------------------------------------------- c call hmixc (joff, js, je, is, ie) c c c----------------------------------------------------------------------- c calculate isopycnal diffusion tensor components (and c 1 advective velocities) for use with tracers c----------------------------------------------------------------------- c call isopyc (joff, js, je, is, ie) c c----------------------------------------------------------------------- c set vertical boundary conditions for momentum and tracers c----------------------------------------------------------------------- c call setvbc (joff, js, je, is, ie) c c----------------------------------------------------------------------- c set vertical mixing coefficients for momentum and tracers c----------------------------------------------------------------------- c call vmixc (joff, js, je, is, ie) c c c c c c----------------------------------------------------------------------- c compute tracers and internal mode velocities c----------------------------------------------------------------------- c call tracer (joff, js_calc, je_calc, is, ie) call baroclinic (joff, js_calc, je_calc, is, ie) c c----------------------------------------------------------------------- c extract first guess for assimilation c----------------------------------------------------------------------- c if (nassim .ge. 1 .and. nassim .le. nassdur) then ! call extract_guess (joff, js_calc, je_calc) call extract_guess (joff, js, je) endif c c----------------------------------------------------------------------- c calculate diagnostics c----------------------------------------------------------------------- c call diag_mw_inside (joff, js_calc, je_calc, is, ie) c c c----------------------------------------------------------------------- c robert time filter for velocity and tracers c----------------------------------------------------------------------- c call r_time_filter (joff, js_calc, je_calc, is, ie &, latdisk(taudisk)) c c----------------------------------------------------------------------- c write prognostic variables from the MW to disk "tau+1" c----------------------------------------------------------------------- c call putmw (joff, js_calc, je_calc, latdisk(taup1disk)) c enddo c c----------------------------------------------------------------------- c update boundary latitudes when using more than one processor c----------------------------------------------------------------------- c ! update tau and tau+1 disk call mpp_update_domains( ramdisk3d, mom_ramdisk_domain ) c c======================================================================= c c SOLVE THE BAROTROPIC EQUATION c c======================================================================= c call barotropic (cori, itt, jstask, jetask) c c----------------------------------------------------------------------- c if this is the 1st pass of an euler backward timestep, set the c disk pointers so the proper time levels are read on the 2nd pass c and go back to do the 2nd pass. c----------------------------------------------------------------------- c if (euler1) then eots = .true. euler1 = .false. euler2 = .true. ntaux = taum1disk taum1disk = taudisk taudisk = taup1disk taup1disk = ntaux else goto_2nd_pass = .false. endif c enddo c c c----------------------------------------------------------------------- c output all remaining diagnostics c----------------------------------------------------------------------- c call diag_mw_outside(jscomp, jecomp) c c----------------------------------------------------------------------- c if it`s the last timestep then clean things up otherwise return c----------------------------------------------------------------------- c if (eorun) then c c calculate cp_time/grid_point/timestep (neglecting setocn) c ncalls = itt - itt1st + 1 cpe = cpusec() write (stdout,'(/,1x, a, f9.3, a,1x,a,1pg14.7,a)') & 'Timing Summary: ', cpe-cps, ' cp seconds,' &, '(',(cpe-cps)/(ncalls*imt*jmt*km) &, ' cp seconds per grid point per time step)' c c write a restart c if (restrt) call write_restart (jstask, jetask) c write (stdout,'(1x,a)') 'MOMdone' call flush(stdout) call oclose (kflds) call oclose (latdisk(1)) call oclose (latdisk(2)) endif c c return end function cpusec () c c user supplied function to measure cpu time in seconds c c real(kind=4) :: second cpusec = second() return end