program pom08 ! ! !----------------------------------------------------------------------! ! For recent changes search for !lyo:_20080415: ! ! ! ! Main changes: ! ! ! ! (1). George Mellor found bugs in subr.profq (also in pom2k); use ! ! the new one in this version. ! ! (2). nsmolar=1 case (for wet-and-dry runs with nwad=1) is ! ! incomplete, and should NOT be used ! ! ! ! ... l.oey (Apr/15/2008) ! ! ! !----------------------------------------------------------------------! ! This POM code has wetting and drying (WAD) capability ! ! ! ! Set nwad=1 for WAD,set =0 to turn off WAD ! ! Set BOTH nwad=0 & nsmolar=0 to recover the original non-WAD POM ! ! ! ! For details see: ! ! ! ! O2005 = Oey [2005; OceanModelling, 9, 133-150] ! ! O2006 = Oey [2006; OceanModelling, 13, 176-195] ! ! O2007 = Oey, Ezer et al [2007; Ocean Dynamics, 57, 205-221] ! ! ! ! or, http://www.aos.princeton.edu/WWWPUBLIC/PROFS/ ! ! also, .../WWWPUBLIC/PROFS/pom98_wad_release_descriptions.html ! ! ! ! For changes that I made, search for the following keyword: ! ! ! !lyo:!wad: or simply "wad:" ! !tne:!wad: (changes taken from T.Ezer's pom98_wad seamount case)! ! ! ! Three (3) files are needed to run the model test cases: ! ! ! ! pom08.f (code), pom08.c (common blocks etc) & pom08.n (netcdf) ! ! ! ! see the runscript runpom08* provided with these files. ! ! ! ! "DOUBLE PRECISION (i.e. REAL*8)" option can be used to ! ! compile, with or without WAD; e.g. w/Intel ifort, use -r8: ! ! ! ! ifort pom08.f -o a.out -r8 ! ! ! ! the -r8 option is the same as specifying -real_size 64 or ! ! -auto_double. However, the code works with or without -r8 ! ! ! ! A link to the netCDF library also needs to be specified. For this! ! and other details, see the runscript: ! ! ! ! runpom08* ! ! ! ! NOTE: for nsmolar=1, the PARAMETER (IM=???,JM=???) in subroutines! ! wadadvt2d & wadsmoladif must match that specified in pom08.c ! ! This is done using "include 'grid'" ! ! ! ! Send bugs and/or improvements to: lyo@princeton.edu ! ! ! ! ... l.oey (May/02/2007) ! ! (Jul/24,30/2007) ! ! (Aug/12,18/2007) ! ! ! ! Support from the Minerals Management Service (MMS) for making ! ! the development of WAD in POM possible is acknowledged. ! ! ! !----------------------------------------------------------------------! ! ! ! ! clyo:Notes: ! c This version incorporates WAD into pom2k.f, and also corrects ! c the bug found by Charles Tang, Glen Carter and Alain Caya (and ! c corrected in the 2006-05-03 version of pom2k.f). ! c Search for "clyo:"for all the changes I have made. ! c ! c ... lyo (Jun/14/2006) ! c ! C ********************************************************************** C * * C * The last code change as rcorded in pom2k.change was on * C * * C * 2006-05-03 * C * (adding IC from file) * C * * C * FUNCTION : This is a version of the three dimensional, time * C * dependent, primitive equation, ocean model * C * developed by Alan Blumberg and George Mellor with * C * subsequent contributions by Leo Oey, Steve Brenner * C * and others. It is now called the Princeton Ocean * C * Model. Two references are: * C * * C * Blumberg, A.F. and G.L. Mellor; Diagnostic and * C * prognostic numerical circulation studies of the * C * South Atlantic Bight, J. Geophys. Res. 88, * C * 4579-4592, 1983. * C * * C * Blumberg, A.F. and G.L. Mellor; A description of a * C * three-dimensional coastal ocean circulation model,* C * Three-Dimensional Coastal Ocean Models, Coastal * C * and Estuarine Sciences, 4, N.S. Heaps, ed., * C * American Geophysical Union, 1-16, 1987. * C * * C * In subroutine profq the model makes use of the * C * turbulence closure sub-model described in: * C * * C * Mellor, G.L. and T. Yamada; Development of a * C * turbulence closure model for geophysical fluid * C * problems, Rev. Geophys. Space Phys., 20, No. 4, * C * 851-875, 1982. * C * (note recent profq that includes breaking waves) * C * * C * A user's guide is available: * C * * C * Mellor, G.L.; User's guide for a three-dimensional, * C * primitive equation, numerical ocean model. * C * Princeton University Report, 1998. * C * * C * In October 2001, the source code underwent * C * revision by John Hunter of the University of * C * Tasmania. Major aspects of the revision were: * C * * C * (1) The revision was based on pom98 updated to * C * 12/9/2001. * C * (2) Declaration of all variables. * C * (3) Rationalisation of the input of all constants. * C * (4) Modifications to the "printer" output. * C * (5) Output to a netCDF file. * C * (6) Inclusion of surface freshwater flux. * C * (7) Inclusion of atmospheric pressure. * C * (8) Inclusion of an additional problem to check (6) * C * and (7), above. * C * (9) Inclusion of option for Smolarkiewicz * C * advection scheme. * C * * C * This revised version is functionally almost * C * equivalent to pom98. The output to device 6 from * C * the "seamount" problem should be almost the same, * C * any differences being due to minor format changes * C * and improvements in rounding. * C * * C * This revision was helped by the following people: * C * Tal Ezer, Peter Holloway, George Mellor, Rich * C * Signell, Ian Webster, Brian Williams and Emma Young.* C * * C ********************************************************************** C * * C * GENERAL NOTES * C * * C * 1. All units are S.I. (M.K.S.) unless otherwise * C * stated. NOTE that time is in days from the start * C * of the run. * C * * C * 2. "b", and "f" refers to backward, * C * central and forward time levels. * C * * C * 3. NetCDF output may be used. In order to omit/use * C * netCDF, comment/uncomment all statements * C * carrying the comment "*netCDF*" at the end of * C * the line (or set netcdf_file='nonetcdf') * C * * C * 4. NetCDF is version 3. An attempt has been made to * C * conform to the NetCDF Climate and Forecast (CF) * C * Metadata Conventions, but this may not yet be * C * complete (see: * C * * C * http://www.cgd.ucar.edu/cms/eaton/cf-metadata/index.html) * C * * C * 5. In order to use netCDF, the program should be * C * compiled with the appropriate library. For * C * example, if using g77, you may need to type: * C * * C * g77 -o pom2k pom2k.f /usr/lib/libnetcdf.a * C * * C * You should also have the "include" file of * C * netCDF subroutines (pom2k.n). * C * * C * 6. In order to use netCDF, you may need to change * C * the name of the "include" file in the statement: * C * * C * include '/usr/include/netcdf.inc' * C * * C * in subroutine write_netcdf * C * * C ********************************************************************** C * * C * SOFTWARE LICENSING * C * * C * This program is free software; you can redistribute * C * it and/or modify it under the terms of the GNU * C * General Public License as published by the Free * C * Software Foundation, either Version 2 of the * C * license, or (at your option) any later version. * C * * C * This program is distributed in the hope that it * C * will be useful, but without any warranty; without * C * even the implied warranty of merchantability or * C * fitness for a particular purpose. See the GNU * C * General Public License for more details. * C * * C * A copy of the GNU General Public License is * C * available at http://www.gnu.org/copyleft/gpl.html * C * or by writing to the Free Software Foundation, Inc.,* C * 59 Temple Place - Suite 330, Boston, MA 02111, USA. * C * * C ********************************************************************** C implicit none C !lyo:!wad:Add WAD variables in pom08.c include 'pom08.c' C C New declarations plus ispi,isp2i: C real aam_init,atot real cbcmax,cbcmin,darea real days,dte2,dvol real eaver real horcon real ispi,isp2i real period,prtd1,prtd2 real saver,smoth,sw,swtch real taver,time0 real vamax,vtot,tsalt real z0b real tatm,satm real wadsmoth !lyo:!wad: real cflmin !lyo:_20080415: integer io(100),jo(100),ko(100) integer i,iend,iext,imax,ispadv,isplit,iswtch integer j,jmax integer k integer nadv,nbct,nbcs,nitera,nread integer iproblem integer nsmolar !lyo:!wad: logical lramp character*120 netcdf_file C C*********************************************************************** C C source should agree with source_c in pom08.c and source_n in C pom08.n. C source='pom08 2008-04-18' C if(source.ne.source_c) then write(6,7) 7 format(/'Incompatible versions of program and include files ', $ '..... program terminated; in pom08.f'/) stop endif C C*********************************************************************** C small=1.e-9 ! Small value C pi=atan(1.e0)*4.e0 ! PI C C*********************************************************************** C C Input of filenames and constants: C C NOTE that the array sizes im, jm and kb should be set in C pom08.c or the 'grid' file created by runpom08 C C----------------------------------------------------------------------- C title='Run 1 ' ! run's title C C----------------------------------------------------------------------- C netcdf_file='pom08.nc' ! netCDF output file c netcdf_file='nonetcdf' ! disable netCDF output C C----------------------------------------------------------------------- C C Problem number: C C iproblem problem initialisation C type subroutine C C 1 seamount seamount C C 2 conservation box C box C C 3 IC from file file2ic C C 4n WAD prob#n, n=1,2 or 3 !lyo:!wad: C iproblem=41 C ! ! !----------------------------------------------------------------------! ! ! ! WAD Control Parameters: ! ! ! !lyo:!wad:set nwad=1 for WAD run, =0 to recover non-WAD run with min ! ! depths set at say 10m (defined in routine called depending on ! ! value of iproblem) ! ! ! nwad=1 ! ! !lyo:!wad:set nsmolar=1 if water depth is to be solved by Smolarkiewicz! ! See O2005 or O2006. ! ! NOTE: nwad=1 w/nsmolar=1 may require finer grid AND even a larger! ! hc, below, say hc=0.1 instead of default hc=0.05 ! ! ! nsmolar=0 ! ! !lyo:!wad:Define hhi0, later hhi will be set to it: ! ! ! ! hhi = HIghest water [above MSL] ever expected, this should ! ! be the observed_highest + say 1 meter to make sure; ! ! here, MSL = Mean Sea Level ! ! ! Choose hhi0=20 (below) assuming that the highest tide or surge ! ! is not expected to rise above 20m wrt MSL ! ! ! hhi0=20.e0 ! ! !lyo:!wad:Define hc in meters: ! ! ! ! hc = Thinnest fluid layer (i.e. water depth) below which cell ! ! is considered dry; default is 0.05m but can be changed in ! ! "params" below; used only if nwad=1 ! ! ! hc=0.05e0 ! ! !lyo:!wad:Set wadsmoth.gt.0 to (cosmetically) smooth out isolated, ! ! thin-fluid wet cells. Smoothing is done if the wet cell is thin, ! ! thickness <= hc*(1+wadsmoth); recommended non-zero wadsmoth is ! ! ~0.05 but it should not be >~0.1 ! wadsmoth=0.e0 ! ! !----------------------------------------------------------------------! ! ! C C----------------------------------------------------------------------- C C mode description C C 2 2-D calculation (bottom stress calculated in advave) C C 3 3-D calculation (bottom stress calculated in profu,v) C C 4 3-D calculation with t and s held fixed C mode=3 C C----------------------------------------------------------------------- C C Advection scheme: C C nadv Advection scheme C C 1 Centred scheme, as originally provide in POM C 2 Smolarkiewicz iterative upstream scheme, based on C subroutines provided by Gianmaria Sannino and Vincenzo C Artale C nadv=1 C C----------------------------------------------------------------------- C C Constants for Smolarkiewicz iterative upstream scheme. C C Number of iterations. This should be in the range 1 - 4. 1 is C standard upstream differencing; 3 adds 50% CPU time to POM: C nitera=3 !=2 !lyo:!wad: C C Smoothing parameter. This should preferably be 1, but 0 < sw < 1 C gives smoother solutions with less overshoot when nitera > 1: C sw=1.e0 !=0.5e0 !lyo:!wad: C C----------------------------------------------------------------------- C C Index to indicate whether run to start from restart file C (nread=0: no restart input file; nread=1: restart input file): C nread=0 C C----------------------------------------------------------------------- C C External (2-D) time step (secs.) according to CFL: ! ! !----------------------------------------------------------------------! ! ! !lyo:!wad:For 3-d runs w/rivers (strong stratification) ! ! strong currents can (and physically they should) develop ! ! in thin film of water over nearly dry cells, so dte and ! ! isplit (below) may need to be smaller. For dx~dy~1km, ! ! I have used values as small as DTE=3.0 and ISPLIT=5. ! ! ! ! Note#1: Since alpha=0 when nwad=1 (see below), the corresponding ! ! DTE needs to be smaller than the DTE used for nwad=0. ! ! Note#2: ISPLIT should also be smaller than its value for no WAD; ! ! see Appendix A of Oey [2006; Ocean Modell. 13, 176-195]. ! ! ! ! For example, for the seamount problem that Tal Ezer has tested ! ! (im,jm=65,49), he found that DTE=60 & ISPLIT=30 worked for ! ! nwad=0 but blew up for nwad=1. However, by decreasing ISPLIT ! ! =5->10, the model with nwad=1 works for DTE=20->60, except that ! ! isolated cells w/thin fluid layers (~5cm) appear, surrounded by ! ! cells which are mostly dry; this probably is to be expected for ! ! the relatively large dx&dy >2~4km used. Using smaller DTE=15 ! ! & ISPLIT=5, or DTE=6 & ISPLIT=10 (as used below) can remove the ! ! isolated thin-fluid cells. ! !----------------------------------------------------------------------! C dte=6.e0 C C----------------------------------------------------------------------- C C / C (dti/dte; dimensionless): C isplit=10 !tne:!wad:dte=6.0 & isplit=30 work w/nwad=0 & hmax=200m C C----------------------------------------------------------------------- C C Date and time of start of initial run of model in format (i.e. C UDUNITS convention) C C YYYY-MM-DD HH:MM:SS <+/->HH:MM C C where "<+/->HH:MM" is the time zone (positive eastwards from C Coordinated Universal Time). NOTE that the climatological time C axis (i.e. beginning of year zero, which does not exist in the C real-world calendar) has been used here. Insert your own date C and time as required: C time_start='2000-01-01 00:00:00 +00:00' C C----------------------------------------------------------------------- C days=0.025e0 ! run duration in days C C----------------------------------------------------------------------- C prtd1=0.0125e0 ! Initial print interval (days) C C----------------------------------------------------------------------- C prtd2=1.e0 ! Final print interval (days) C C----------------------------------------------------------------------- C swtch=1000.e0 ! Time to switch from prtd1 to prtd2 C C----------------------------------------------------------------------- C iskp=1 ! Printout skip interval in i C C----------------------------------------------------------------------- C jskp=1 ! Printout skip interval in j C C----------------------------------------------------------------------- C C Logical for inertial ramp (.true. if inertial ramp to be applied C to wind stress and baroclinic forcing, otherwise .false.) C lramp=.false. C C----------------------------------------------------------------------- C C Reference density (recommended values: 1025 for seawater, C 1000 for freswater; S.I. units): C rhoref=1025.e0 C C----------------------------------------------------------------------- C tbias=0.e0 ! Temperature bias (deg. C) C C----------------------------------------------------------------------- C sbias=0.e0 ! Salinity bias C C----------------------------------------------------------------------- C grav=9.806e0 ! gravity constant (S.I. units) C C----------------------------------------------------------------------- C kappa=0.4e0 ! von Karman's constant C C----------------------------------------------------------------------- C z0b=.01e0 ! Bottom roughness (metres) C C----------------------------------------------------------------------- C zsh=.01e0 ! Bottom Log-layer shift (metres) !lyo:!wad: C C----------------------------------------------------------------------- C cbcmin=.0025e0 ! Minimum bottom friction coeff. C C----------------------------------------------------------------------- C cbcmax=1.e0 ! Maximum bottom friction coeff. C C----------------------------------------------------------------------- C !lyo:!wad:I would use horcon=0.1 (un-related to WAD) horcon=0.2e0 ! Smagorinsky diffusivity coeff. C C----------------------------------------------------------------------- C C Inverse horizontal turbulent Prandtl number C (ah/am; dimensionless): C C NOTE that tprni=0.e0 yields zero horizontal diffusivity! C tprni=.2e0 C C----------------------------------------------------------------------- C C Background viscosity used in subroutines profq, proft, profu and C profv (S.I. units): C umol=2.e-5 C C----------------------------------------------------------------------- C C Maximum depth used in radiation boundary condition in subroutine C bcond (metres): C hmax=4500.e0 C C----------------------------------------------------------------------- C C Maximum magnitude of vaf (used in check that essentially tests C for CFL violation): C vmaxl=100.e0 C C----------------------------------------------------------------------- C C Maximum allowable value of: C C / C C for two adjacent cells (dimensionless). This is used in subroutine C slpmax. If >= 1, then slpmax is not applied: C slmax=2.e0 C C----------------------------------------------------------------------- C C Integers defining the number of logarithmic layers at the C surface and bottom (used by subroutine depth). The number of C logarithmic layers are kl1-2 at the surface and kb-kl2-1 C at the bottom. For no log portions, set kl1=2 and kl2=kb-1: C kl1=6 kl2=kb-2 C C----------------------------------------------------------------------- C C Water type, used in subroutine proft. C C ntp Jerlov water type C C 1 i C 2 ia C 3 ib C 4 ii C 5 iii C ntp=2 C C----------------------------------------------------------------------- C C Surface temperature boundary condition, used in subroutine proft: C C nbct prescribed prescribed short wave C temperature flux penetration C C 1 no yes no C 2 no yes yes C 3 yes no no C 4 yes no yes C nbct=1 C C----------------------------------------------------------------------- C C Surface salinity boundary condition, used in subroutine proft: C C nbcs prescribed prescribed C salinity flux C C 1 no yes C 3 yes no C C NOTE that only 1 and 3 are allowed for salinity. C nbcs=1 C C----------------------------------------------------------------------- C C Step interval during which external (2-D) mode advective terms are C not updated (dimensionless): C ispadv=5 C C----------------------------------------------------------------------- C C Constant in temporal filter used to prevent solution splitting C (dimensionless): C smoth=0.10e0 C C----------------------------------------------------------------------- C C Weight used for surface slope term in external (2-D) dynamic C equation (a value of alph0 = 0.e0 is perfectly acceptable, but the C value, alph0=.225e0 permits a longer time step): C alph0=0.225e0 !lyo:!wad:use alph0 instead of alpha - defined later C C----------------------------------------------------------------------- C C Initial value of aam: C aam_init=500.e0 ! ! !----------------------------------------------------------------------! !lyo:!wad:Tidal amplitude for iproblem=41 w/ or w/o WAD: ! ! ! !lyo:_20080415: tidamp=0.e0 if (iproblem .eq. 41) tidamp=9.e0 !----------------------------------------------------------------------! ! ! C C End of input of constants C*********************************************************************** C C --- Above are the default parameters, alternatively one can C --- use parameters from a file created by runscript runpom2k_pow_wad !clyo:wad: C include 'params' C clyo:wad:beg: c Overwrite some input constants: see "params" above in runpom08 clyo:wad:end: c C*********************************************************************** C C Calculate some constants: C dti=dte*float(isplit) dte2=dte*2 dti2=dti*2 ! ! !----------------------------------------------------------------------! !lyo:!wad: ! ! ! ! Dry-cell T/S relaxation using leapfrog-trapezoidal: ! dtrat=dti/(0.5*86400.) !1/2 day relaxation time-scale cwetrlx1=(1.-dtrat)/(1.+dtrat); cwetrlx2=2.*dtrat/(1.+dtrat) ! ! alpha=alph0*float(1-nwad) !lyo:!wad: must=0 for WAD ! ! hhi=hhi0*float(nwad) ! ! nsmolar=nsmolar*nwad ! ! ! Check nsmolar: !lyo:_20080415: ! ! ! if( (nsmolar.ne.0).and.(mode.ne.2) )then write(6,'('' Stopped; incorrect inputs of '')') write(6,'('' nsmolar = '',i10)') nsmolar write(6,'('' mode = '',i10)') mode write(6,'('' For nsmolar ne 0, the present version works '')') write(6,'('' only if mode = 2. So either run in barotropic '')') write(6,'('' mode = 2, or set nsmolar=0 and resubmit '')') stop endif ! ! ! Check iproblem & nwad compatibility: ! ! ! if( (iproblem.ne.41).and.(nwad.eq.1) )then !lyo:_20080415: write(6,'('' iproblem = '',i10)') iproblem write(6,'('' nwad = '',i10)') nwad write(6,'('' iproblem must = 41 for nwad = 1 '')') !lyo:_20080415: write(6,'('' Stopped: iproblem & nwad are not compatible '')') stop endif ! ! !----------------------------------------------------------------------! ! ! c iend=max0(nint(days*24.e0*3600.e0/dti),2) iprint=nint(prtd1*24.e0*3600.e0/dti) iswtch=nint(swtch*24.e0*3600.e0/dti) C ispi=1.e0/float(isplit) isp2i=1.e0/(2.e0*float(isplit)) C C----------------------------------------------------------------------- C C Print initial summary: C write(6,'(/,'' source = '',a40)') source write(6,'('' title = '',a40/)') title write(6,'('' iproblem = '',i10)') iproblem write(6,'('' nwad = '',i10)') nwad !lyo:!wad: write(6,'('' nsmolar = '',i10)') nsmolar !lyo:!wad: write(6,'('' hhi = '',f10.3)') hhi !lyo:!wad: write(6,'('' hc = '',f10.3)') hc !lyo:!wad: write(6,'('' wadsmoth = '',f10.3)') wadsmoth !lyo:!wad: write(6,'('' mode = '',i10)') mode write(6,'('' nadv = '',i10)') nadv write(6,'('' nitera = '',i10)') nitera write(6,'('' sw = '',f10.4)') sw write(6,'('' nread = '',i10)') nread write(6,'('' dte = '',f10.2)') dte write(6,'('' dti = '',f10.1)') dti write(6,'('' isplit = '',i10)') isplit write(6,'('' time_start = '',a26)') time_start write(6,'('' days = '',f10.4)') days write(6,'('' iend = '',i10)') iend write(6,'('' prtd1 = '',f10.4)') prtd1 write(6,'('' iprint = '',i10)') iprint write(6,'('' prtd2 = '',f10.4)') prtd2 write(6,'('' swtch = '',f10.2)') swtch write(6,'('' iswtch = '',i10)') iswtch write(6,'('' iskp, jskp = '',i5'','',i5)') iskp,jskp write(6,'('' lramp = '',l10)') lramp write(6,'('' rhoref = '',f10.3)') rhoref write(6,'('' tbias = '',f10.3)') tbias write(6,'('' sbias = '',f10.3)') sbias write(6,'('' grav = '',f10.4)') grav write(6,'('' kappa = '',f10.4)') kappa write(6,'('' z0b = '',f10.6)') z0b write(6,'('' zsh = '',f10.6)') zsh !lyo:!wad: write(6,'('' cbcmin = '',f10.6)') cbcmin write(6,'('' cbcmax = '',f10.6)') cbcmax write(6,'('' horcon = '',f10.3)') horcon write(6,'('' tprni = '',f10.4)') tprni write(6,'('' umol = '',f10.4)') umol write(6,'('' hmax = '',f10.2)') hmax write(6,'('' vmaxl = '',f10.4)') vmaxl write(6,'('' slmax = '',f10.4)') slmax write(6,'('' kl1, kl2 = '',i5,'','',i5)') kl1,kl2 write(6,'('' ntp = '',i10)') ntp write(6,'('' nbct = '',i10)') nbct write(6,'('' nbcs = '',i10)') nbcs write(6,'('' ispadv = '',i10)') ispadv write(6,'('' smoth = '',f10.4)') smoth write(6,'('' alpha = '',f10.4)') alpha write(6,'('' tidamp = '',f10.4)') tidamp C C----------------------------------------------------------------------- C C Initialise boundary arrays: C do i=1,im vabn(i)=0.e0 vabs(i)=0.e0 eln(i)=0.e0 !lyo:!wad:note:keep=0.e0 inst. of -hhi els(i)=0.e0 ! .. redefined later in wadseamount etc. do k=1,kb vbn(i,k)=0.e0 vbs(i,k)=0.e0 tbn(i,k)=0.e0 tbs(i,k)=0.e0 sbn(i,k)=0.e0 sbs(i,k)=0.e0 end do end do C do j=1,jm uabe(j)=0.e0 uabw(j)=0.e0 ele(j)=0.e0 !lyo:!wad:note:keep=0.e0 inst. of -hhi elw(j)=0.e0 ! .. redefined later in wadseamount etc. do k=1,kb ube(j,k)=0.e0 ubw(j,k)=0.e0 tbe(j,k)=0.e0 tbw(j,k)=0.e0 sbe(j,k)=0.e0 sbw(j,k)=0.e0 end do end do C C----------------------------------------------------------------------- C C Initialise 2-D and 3-D arrays for safety (this may be overwritten C later): C do j=1,jm do i=1,im uab(i,j)=0.e0 vab(i,j)=0.e0 elb(i,j)=0.e0 !lyo:!wad:note:keep=0.e0 inst. of -hhi etb(i,j)=0.e0 ! .. redefined later in wadseamount etc. e_atmos(i,j)=0.e0 !lyo:!wad:note:=0.e0 is correct vfluxb(i,j)=0.e0 vfluxf(i,j)=0.e0 wusurf(i,j)=0.e0 wvsurf(i,j)=0.e0 wtsurf(i,j)=0.e0 wssurf(i,j)=0.e0 swrad(i,j)=0.e0 drx2d(i,j)=0.e0 dry2d(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im ub(i,j,k)=0.e0 vb(i,j,k)=0.e0 end do end do end do C C----------------------------------------------------------------------- C C Set up sigma layers: C if(iproblem.ne.3) call depth C C----------------------------------------------------------------------- C C Read in grid data, and initial and lateral boundary conditions: C if(iproblem.eq.1) then call seamount else if(iproblem.eq.2) then call box else if(iproblem.eq.3) then call file2ic else if(iproblem.eq.41) then !lyo:!wad: call wadseamount else write(6,8) 8 format(/' Invalid value of iproblem ..... program terminated'/) stop endif ! !lyo:_20080415: ! Make sure that wetmask = fsm for non-WAD runs: if (nwad.eq.0) then wetmask(:,:)=fsm(:,:) endif ! C C Inertial period for temporal filter: C !lyo:_20080415: if (cor(im/2,jm/2).ne.0.0) then period=(2.e0*pi)/abs(cor(im/2,jm/2))/86400.e0 else period=1.0 !set to 1day endif C C Initialise time: C time0=0.e0 time=0.e0 C C Initial conditions: C C NOTE that lateral thermodynamic boundary conditions are often set C equal to the initial conditions and are held constant thereafter. C Users can of course create variable boundary conditions. C do i=1,im do j=1,jm ua(i,j)=uab(i,j) va(i,j)=vab(i,j) el(i,j)=elb(i,j) !lyo:!wad:note:already defined et(i,j)=etb(i,j) ! in wadseamount w/hhi etc etf(i,j)=et(i,j) d(i,j)=h(i,j)+el(i,j) dt(i,j)=h(i,j)+et(i,j) w(i,j,1)=vfluxf(i,j) end do end do C do k=1,kb do j=1,jm do i=1,im l(i,j,k)=0.1*dt(i,j) q2b(i,j,k)=small q2lb(i,j,k)=l(i,j,k)*q2b(i,j,k) kh(i,j,k)=l(i,j,k)*sqrt(q2b(i,j,k)) km(i,j,k)=kh(i,j,k) kq(i,j,k)=kh(i,j,k) aam(i,j,k)=aam_init end do end do end do if (horcon.le.0.0) aam(:,:,:)=-horcon !lyo:_20080415: C do k=1,kbm1 do i=1,im do j=1,jm q2(i,j,k)=q2b(i,j,k) q2l(i,j,k)=q2lb(i,j,k) t(i,j,k)=tb(i,j,k) s(i,j,k)=sb(i,j,k) u(i,j,k)=ub(i,j,k) v(i,j,k)=vb(i,j,k) end do end do end do C call dens(s,t,rho) C call baropg C do k=1,kbm1 do j=1,jm do i=1,im drx2d(i,j)=drx2d(i,j)+drhox(i,j,k)*dz(k) dry2d(i,j)=dry2d(i,j)+drhoy(i,j,k)*dz(k) end do end do end do C !lyo:!wad:cbc_change:begins:-------------------------------------------! ! ! !----------------------------------------------------------------------! ! Calculate bottom friction coefficient: ! ! ... see O2006, Appendix A; also Mellor,JPO,2002:Osc.Blayer... ! !----------------------------------------------------------------------! ! ! do j=1,jm do i=1,im cbc(i,j)=cbcmin ! cbc(i,j)=(kappa/log((1.e0+zz(kbm1))*h(i,j)/z0b))**2 cbc(i,j)=(kappa/log((zsh+(1.e0+zz(kbm1))*d(i,j))/z0b))**2 !lyo:!wad:lyo's originals:---------------------------------------------! !vers.1: cbc(i,j)=(kappa/log((zsh+(zz(kbm1)-z(kb))*d(i,j))/z0b))**2 ! !vers.2: cbc(i,j)=max(cbcmin, ! ! $ (kappa/log((zsh+(zz(kbm1)-z(kb))*d(i,j))/z0b))**2) ! !lyo:!wad:lyo's originals:---------------------------------------------! ! cbc(i,j)=max(cbcmin,cbc(i,j)) C C If the following is invoked, then it is probable that the wrong C choice of z0b or vertical spacing has been made: C cbc(i,j)=min(cbcmax,cbc(i,j)) end do end do !lyo:!wad:cbc_change:ends:---------------------------------------------! C C Calculate external (2-D) CFL time step: C !lyo:_20080415: cflmin=1.e10 do j=1,jm do i=1,im tps(i,j)=0.5e0/sqrt(1.e0/dx(i,j)**2+1.e0/dy(i,j)**2) $ /sqrt(grav*(h(i,j)+small))*fsm(i,j) if (fsm(i,j).ne.0.0) cflmin=min(cflmin,tps(i,j)) !lyo:_20080415: end do end do C C----------------------------------------------------------------------- C C The following data are needed for a seamless restart. if nread=1, C data had been created by a previous run (see write(71) at end of C this program). nread=0 denotes a first time run. C if(nread.eq.1) $ read(70) time0, $ wubot,wvbot,aam2d,ua,uab,va,vab,el,elb,et,etb,egb, $ utb,vtb,u,ub,w,v,vb,t,tb,s,sb,rho, $ adx2d,ady2d,advua,advva, $ km,kh,kq,l,q2,q2b,aam,q2l,q2lb $ ,wetmask,wmarsh !lyo:!wad: C do j=1,jm do i=1,im d(i,j)=h(i,j)+el(i,j) dt(i,j)=h(i,j)+et(i,j) end do end do C ! ! !----------------------------------------------------------------------! !lyo:!wad:Update cbc: ! ! ! cbc(:,:)=cbcmin !lyo:_20080415: if (mode.ne.2) then !lyo:_20080415: do j=1,jm do i=1,im ! cbc(i,j)=cbcmin !lyo:_20080415: cbc(i,j)=(kappa/log((zsh+(1.e0+zz(kbm1))*d(i,j))/z0b))**2 cbc(i,j)=max(cbcmin,cbc(i,j)) cbc(i,j)=min(cbcmax,cbc(i,j)) end do end do endif ! ! !----------------------------------------------------------------------! ! time=time0 C C----------------------------------------------------------------------- C C Print geometry and other initial fields (select statements as C desired): C call prxy('grid increment in x, dx ', $ time,dx ,im,iskp,jm,jskp,0.e0) C call prxy('grid increment in y, dy ', $ time,dy ,im,iskp,jm,jskp,0.e0) C call prxy('Easting of elevation points, east_e ', $ time,east_e ,im,iskp,jm,jskp,0.e0) C call prxy('Northing of elevation points, north_e ', $ time,north_e,im,iskp,jm,jskp,0.e0) C call prxy('Easting of cell corners, east_c ', $ time,east_c ,im,iskp,jm,jskp,0.e0) C call prxy('Northing of cell corners, north_c ', $ time,north_c,im,iskp,jm,jskp,0.e0) C call prxy('Rotation angle of x-axis wrt. east, rot ', $ time,rot,im,iskp,jm,jskp,0.e0) C call prxy('Undisturbed water depth, h ', $ time,h ,im,iskp,jm,jskp,1.e1) C call prxy('Land mask, fsm ', !lyo:!wad: $ time,fsm,im,iskp,jm,jskp,1.e0) C call prxy('Wet & dry mask, wetmask ', !lyo:!wad: $ time,wetmask,im,iskp,jm,jskp,1.e0) C call prxy('u-velocity mask, dum ', $ time,dum,im,iskp,jm,jskp,1.e0) C call prxy('v-velocity mask, dvm ', $ time,dvm,im,iskp,jm,jskp,1.e0) C write(6,'('' Min CFL Time-step, cflmin = '',f10.4)') cflmin !lyo:_20080415: call prxy('External (2-D) CFL time step, tps ', $ time,tps,im,iskp,jm,jskp,1.e0) C call prxy('Bottom friction coefficient, cbc ', !lyo:!wad: $ time,cbc,im,iskp,jm,jskp,0.e0) C C Set sections for output: C ko(1)=1 ko(2)=kb/2 ko(3)=kb-1 C call prxyz('Horizontally-averaged rho, rmean ', $ time,rmean,im,iskp,jm,jskp,kb,ko,3,1.e-5) C C Set sections for output: C jo(1)=1 jo(2)=jm/2 jo(3)=jm-1 C call prxz('Horizontally-averaged rho, rmean ', $ time,rmean,im,iskp,jm,kb,jo,3,1.e-5,dt,zz) C C Set sections for output: C io(1)=1 io(2)=im/2 io(3)=im-1 C call pryz('Horizontally-averaged rho, rmean ', $ time,rmean,im,jm,jskp,kb,io,3,1.e-5,dt,zz) C C----------------------------------------------------------------------- C C Initial conditions: C C Select print statements in printall as desired: C call printall C C----------------------------------------------------------------------- C C Initialise netCDF output and output initial set of data: C if(netcdf_file.ne.'nonetcdf') then call write_netcdf(netcdf_file,1) ! *netCDF* call write_netcdf(netcdf_file,2) ! *netCDF* endif C C----------------------------------------------------------------------- C do 9000 iint=1,iend ! Begin internal (3-D) mode C time=dti*float(iint)/86400.e0+time0 C if(lramp) then ramp=time/period if(ramp.gt.1.e0) ramp=1.e0 else ramp=1.e0 endif C C write(6,2) mode,iint,time C 2 format(' mode,iint,time =',2i5,f9.2) C C----------------------------------------------------------------------- C C Set time dependent, surface and lateral boundary conditions. C The latter will be used in subroutine bcond. Users may C wish to create a subroutine to supply wusurf, wvsurf, wtsurf, C wssurf, swrad and vflux. C C Introduce simple wind stress. Value is negative for westerly or C southerly winds. The following wind stress has been tapered C along the boundary to suppress numerically induced oscilations C near the boundary (Jamart and Ozer, J.G.R., 91, 10621-10631). C To make a healthy surface Ekman layer, it would be well to set C kl1=9. C do j=2,jmm1 do i=2,imm1 c if(iproblem.ne.3) then ! constant wind read in file2ic c c wusurf(i,j)=ramp*(1.e-4*cos(pi*(j-1)/jmm1)) wusurf(i,j)=1.00*(1.e-4*cos(pi*(j-1)/jmm1)) $ *.25e0*(dvm(i,j+1)+dvm(i-1,j+1) $ +dvm(i-1,j)+dvm(i,j)) C --- no wind ---- c wusurf(i,j)=0.e0 !lyo:_20080415: wvsurf(i,j)=0.e0 endif e_atmos(i,j)=0.e0 vfluxf(i,j)=0.e0 C C Set w(i,j,1)=vflux(i,j).ne.0 if one wishes non-zero flow across C the sea surface. See calculation of elf(i,j) below and subroutines C vertvl, advt1 (or advt2). If w(1,j,1)=0, and, additionally, there C is no net flow across lateral boundaries, the basin volume will be C constant; if also vflux(i,j).ne.0, then, for example, the average C salinity will change and, unrealistically, so will total salt. C w(i,j,1)=vfluxf(i,j) C C Set wtsurf to the sensible heat, the latent heat (which involves C only the evaporative component of vflux) and the long wave C radiation: C wtsurf(i,j)=0.e0 C C Set swrad to the short wave radiation: C swrad(i,j)=0.e0 C C To account for change in temperature of flow crossing the sea C surface (generally quite small compared to latent heat effect) C tatm=t(i,j,1)+tbias ! an approximation wtsurf(i,j)=wtsurf(i,j)+vfluxf(i,j)*(tatm-t(i,j,1)-tbias) C C Set the salinity of water vapor/precipitation which enters/leaves C the atmosphere (or e.g., an ice cover) C satm=0.e0 wssurf(i,j)= vfluxf(i,j)*(satm-s(i,j,1)-sbias) C end do end do C clyo: ! call powdriver(iprint,nread,z0b,cbcmin,iend/iprint,fsm) c C----------------------------------------------------------------------- C C Set lateral viscosity: C C If mode=2 then initial values of aam2d are used. If one wishes C to use Smagorinsky lateral viscosity and diffusion for an C external (2-D) mode calculation, then appropiate code can be C adapted from that below and installed just before the end of the C "if(mode.eq.2)" loop in subroutine advave. C C Calculate Smagorinsky lateral viscosity: C C ( hor visc = horcon*dx*dy*sqrt((du/dx)**2+(dv/dy)**2 C +.5*(du/dy+dv/dx)**2) ) C if(mode.ne.2) then call advct(a,c,ee) call baropg C if (horcon.gt.0.0) then !lyo:_20080415: do k=1,kbm1 do j=2,jmm1 do i=2,imm1 aam(i,j,k)=horcon*dx(i,j)*dy(i,j) $ *sqrt( ((u(i+1,j,k)-u(i,j,k))/dx(i,j))**2 $ +((v(i,j+1,k)-v(i,j,k))/dy(i,j))**2 $ +.5e0*(.25e0*(u(i,j+1,k)+u(i+1,j+1,k) $ -u(i,j-1,k)-u(i+1,j-1,k)) $ /dy(i,j) $ +.25e0*(v(i+1,j,k)+v(i+1,j+1,k) $ -v(i-1,j,k)-v(i-1,j+1,k)) $ /dx(i,j)) **2) end do end do end do endif !if (horcon.gt.0.0) then !lyo:_20080415: C C Form vertical averages of 3-D fields for use in external (2-D) C mode: C do j=1,jm do i=1,im adx2d(i,j)=0.e0 ady2d(i,j)=0.e0 drx2d(i,j)=0.e0 dry2d(i,j)=0.e0 aam2d(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im adx2d(i,j)=adx2d(i,j)+advx(i,j,k)*dz(k) ady2d(i,j)=ady2d(i,j)+advy(i,j,k)*dz(k) drx2d(i,j)=drx2d(i,j)+drhox(i,j,k)*dz(k) dry2d(i,j)=dry2d(i,j)+drhoy(i,j,k)*dz(k) aam2d(i,j)=aam2d(i,j)+aam(i,j,k)*dz(k) end do end do end do C call advave(tps) C do j=1,jm do i=1,im adx2d(i,j)=adx2d(i,j)-advua(i,j) ady2d(i,j)=ady2d(i,j)-advva(i,j) end do end do C endif C clyo:beg:changes by Tang et al. do j=1,jm do i=1,im egf(i,j)=el(i,j)*ispi end do end do C do j=1,jm do i=2,im utf(i,j)=ua(i,j)*(d(i,j)+d(i-1,j))*isp2i end do end do do j=2,jm do i=1,im vtf(i,j)=va(i,j)*(d(i,j)+d(i,j-1))*isp2i end do end do clyo:end: C C----------------------------------------------------------------------- C clyomoving:decide to skip wriv & wmarsh for now; wriv in particular c needs to make consistent with vflux do 8000 iext=1,isplit ! Begin external (2-D) mode C C write(6,3) iext,time C 3 format(' iext,time =',i5,f9.2) C ! ! !----------------------------------------------------------------------! !lyo:!wad:Water depth can be solved by Smolarkiewicz to prevent oscil. ! ! ! ! nsmolar=1 case (for wet-and-dry runs with nwad=1) is ! ! incomplete, and should NOT be used ! ! ! if (nsmolar.eq.1 .and. mode.eq.2) then !lyo:_20080415: dd0(:,:)=h(:,:)+ el(:,:)*fsm(:,:) ddb(:,:)=h(:,:)+elb(:,:)*fsm(:,:) tps(:,:)=0.0 !aam2d(:,:) call wadadvt2d(ddb,dd0,ddf,nitera,sw,dx,dy,ua,va,art,aru,arv, 1 fsm,dum,dvm,tps,dte2) elf(:,:)=ddf(:,:)-h(:,:) else ! do j=2,jm do i=2,im fluxua(i,j)=.25e0*(d(i,j)+d(i-1,j)) $ *(dy(i,j)+dy(i-1,j))*ua(i,j) fluxva(i,j)=.25e0*(d(i,j)+d(i,j-1)) $ *(dx(i,j)+dx(i,j-1))*va(i,j) end do end do C C NOTE addition of surface freshwater flux, w(i,j,1)=vflux, compared C with pom98.f. See also modifications to subroutine vertvl. C do j=2,jmm1 do i=2,imm1 elf(i,j)=elb(i,j) $ +dte2*(-(fluxua(i+1,j)-fluxua(i,j) $ +fluxva(i,j+1)-fluxva(i,j))/art(i,j) $ -vfluxf(i,j)) end do end do ! ! endif ! ! !----------------------------------------------------------------------! C call bcond(1) ! ! !----------------------------------------------------------------------! !lyo:!wad:Check wet/dry on ELF: ! ! ! if (nwad.eq.1) then do j=1,jm; do i=1,im DWET(I,J)=H(I,J)+ELF(I,J)*fsm(i,j) enddo; enddo do j=1,jm; do i=1,im wetmask(i,j)=fsm(i,j) if (dwet(i,j).le.hc) wetmask(i,j)=0.0 enddo; enddo if (wadsmoth.gt.0.0) then ! Smooth isolated wet spots: ! do j=1,jm; do i=1,im if( (wetmask(i,j) .gt. $ (wetmask(i-1,j)+wetmask(i+1,j)+wetmask(i,j-1)+wetmask(i,j+1))) $ .and. (dwet(i,j).le.hc*(1.e0+wadsmoth)) ) then ! $ .and. (dwet(i,j).le.hc*1.01) ) then wetmask(i,j)=0.0 ! elf(i,j)=elb(i,j); dwet(i,j)=h(i,j)+elf(i,j)*fsm(i,j) endif enddo; enddo endif endif ! ! !----------------------------------------------------------------------! ! ! if(mod(iext,ispadv).eq.0) call advave(tps) C do j=2,jmm1 do i=2,im uaf(i,j)=adx2d(i,j)+advua(i,j) $ -aru(i,j)*.25e0 $ *(cor(i,j)*d(i,j)*(va(i,j+1)+va(i,j)) $ +cor(i-1,j)*d(i-1,j)*(va(i-1,j+1)+va(i-1,j))) $ +.25e0*grav*(dy(i,j)+dy(i-1,j)) $ *(d(i,j)+d(i-1,j)) $ *((1.e0-2.e0*alpha) $ *(el(i,j)-el(i-1,j)) $ +alpha*(elb(i,j)-elb(i-1,j) $ +elf(i,j)-elf(i-1,j)) $ +e_atmos(i,j)-e_atmos(i-1,j)) $ +drx2d(i,j)+aru(i,j)*(wusurf(i,j)-wubot(i,j)) end do end do C do j=2,jmm1 do i=2,im uaf(i,j)=((h(i,j)+elb(i,j)+h(i-1,j)+elb(i-1,j)) $ *aru(i,j)*uab(i,j) $ -4.e0*dte*uaf(i,j)) $ /((h(i,j)+elf(i,j)+h(i-1,j)+elf(i-1,j)) $ *aru(i,j)) end do end do C do j=2,jm do i=2,imm1 vaf(i,j)=ady2d(i,j)+advva(i,j) $ +arv(i,j)*.25e0 $ *(cor(i,j)*d(i,j)*(ua(i+1,j)+ua(i,j)) $ +cor(i,j-1)*d(i,j-1)*(ua(i+1,j-1)+ua(i,j-1))) $ +.25e0*grav*(dx(i,j)+dx(i,j-1)) $ *(d(i,j)+d(i,j-1)) $ *((1.e0-2.e0*alpha)*(el(i,j)-el(i,j-1)) $ +alpha*(elb(i,j)-elb(i,j-1) $ +elf(i,j)-elf(i,j-1)) $ +e_atmos(i,j)-e_atmos(i,j-1)) $ +dry2d(i,j)+arv(i,j)*(wvsurf(i,j)-wvbot(i,j)) end do end do C do j=2,jm do i=2,imm1 vaf(i,j)=((h(i,j)+elb(i,j)+h(i,j-1)+elb(i,j-1)) $ *vab(i,j)*arv(i,j) $ -4.e0*dte*vaf(i,j)) $ /((h(i,j)+elf(i,j)+h(i,j-1)+elf(i,j-1)) $ *arv(i,j)) end do end do C call bcond(2) C if(iext.eq.(isplit-2))then do j=1,jm do i=1,im etf(i,j)=.25e0*smoth*elf(i,j) end do end do C else if(iext.eq.(isplit-1)) then C do j=1,jm do i=1,im etf(i,j)=etf(i,j)+.5e0*(1.-.5e0*smoth)*elf(i,j) end do end do C else if(iext.eq.isplit) then C do j=1,jm do i=1,im etf(i,j)=(etf(i,j)+.5e0*elf(i,j))*fsm(i,j) end do end do C endif ! ! !----------------------------------------------------------------------! !lyo:!wad:Check wet/dry on UAF & VAF: ! ! ! if (nwad.eq.1) then do j=1,jm; do i=2,im if (0.5*(dwet(i,j)+dwet(i-1,j)).le.hc) then uaf(i,j)=0.0 else !Set flux=0 if coming from dry cell: if (wetmask(i-1,j).eq.0.0 .and. uaf(i,j).gt.0.0) uaf(i,j)=0.0 if (wetmask(i,j) .eq.0.0 .and. uaf(i,j).lt.0.0) uaf(i,j)=0.0 endif enddo; enddo ! do j=2,jm; do i=1,im if (0.5*(dwet(i,j)+dwet(i,j-1)).le.hc) then vaf(i,j)=0.0 else !Set flux=0 if coming from dry cell: if (wetmask(i,j-1).eq.0.0 .and. vaf(i,j).gt.0.0) vaf(i,j)=0.0 if (wetmask(i,j) .eq.0.0 .and. vaf(i,j).lt.0.0) vaf(i,j)=0.0 endif enddo; enddo endif ! ! !----------------------------------------------------------------------! C C Stop if velocity condition violated (generally due to CFL C criterion not being satisfied): C vamax=0.e0 C do j=1,jm do i=1,im if(abs(vaf(i,j)).ge.vamax) then vamax=abs(vaf(i,j)) imax=i jmax=j endif end do end do C if(vamax.le.vmaxl) then C C Apply filter to remove time split and reset time sequence: C do j=1,jm do i=1,im ua(i,j)=ua(i,j) $ +.5e0*smoth*(uab(i,j)-2.e0*ua(i,j)+uaf(i,j)) va(i,j)=va(i,j) $ +.5e0*smoth*(vab(i,j)-2.e0*va(i,j)+vaf(i,j)) el(i,j)=el(i,j) $ +.5e0*smoth*(elb(i,j)-2.e0*el(i,j)+elf(i,j)) elb(i,j)=el(i,j) el(i,j)=elf(i,j) d(i,j)=h(i,j)+el(i,j) uab(i,j)=ua(i,j) ua(i,j)=uaf(i,j) vab(i,j)=va(i,j) va(i,j)=vaf(i,j) end do end do ! ! !----------------------------------------------------------------------! !lyo:!wad:Update cbc: ! ! ! if (mode.ne.2) then !lyo:_20080415: do j=1,jm do i=1,im cbc(i,j)=cbcmin cbc(i,j)=(kappa/log((zsh+(1.e0+zz(kbm1))*d(i,j))/z0b))**2 cbc(i,j)=max(cbcmin,cbc(i,j)) cbc(i,j)=min(cbcmax,cbc(i,j)) end do end do endif ! ! !----------------------------------------------------------------------! ! if(iext.ne.isplit) then clyo:beg:changes by Tang et al. do j=1,jm do i=1,im egf(i,j)=egf(i,j)+el(i,j)*ispi end do end do do j=1,jm do i=2,im utf(i,j)=utf(i,j)+ua(i,j)*(d(i,j)+d(i-1,j))*isp2i end do end do do j=2,jm do i=1,im vtf(i,j)=vtf(i,j)+va(i,j)*(d(i,j)+d(i,j-1))*isp2i end do end do clyo:end: endif C endif C 8000 continue ! End of external (2-D) mode C C----------------------------------------------------------------------- C if(vamax.le.vmaxl) then C C Continue with internal (3-D) mode calculation: C if((iint.ne.1.or.time0.ne.0.e0).and.mode.ne.2) then C C Adjust u(z) and v(z) such that depth average of (u,v) = (ua,va): C do j=1,jm do i=1,im tps(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tps(i,j)=tps(i,j)+u(i,j,k)*dz(k) end do end do end do C do k=1,kbm1 do j=1,jm do i=2,im u(i,j,k)=(u(i,j,k)-tps(i,j))+ $ (utb(i,j)+utf(i,j))/(dt(i,j)+dt(i-1,j)) end do end do end do C do j=1,jm do i=1,im tps(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tps(i,j)=tps(i,j)+v(i,j,k)*dz(k) end do end do end do C do k=1,kbm1 do j=2,jm do i=1,im v(i,j,k)=(v(i,j,k)-tps(i,j))+ $ (vtb(i,j)+vtf(i,j))/(dt(i,j)+dt(i,j-1)) end do end do end do C C vertvl calculates w from u, v, dt (h+et), etf and etb: C call vertvl(a,c) call bcond(5) C C do k=1,kb do j=1,jm do i=1,im uf(i,j,k)=0.e0 vf(i,j,k)=0.e0 end do end do end do C C Calculate q2f and q2lf using uf, vf, a and c as temporary C variables: C call advq(q2b,q2,uf,a,c) call advq(q2lb,q2l,vf,a,c) call profq(a,c,tps,dtef) call bcond(6) C do k=1,kb do j=1,jm do i=1,im q2(i,j,k)=q2(i,j,k) $ +.5e0*smoth*(uf(i,j,k)+q2b(i,j,k) $ -2.e0*q2(i,j,k)) q2l(i,j,k)=q2l(i,j,k) $ +.5e0*smoth*(vf(i,j,k)+q2lb(i,j,k) $ -2.e0*q2l(i,j,k)) q2b(i,j,k)=q2(i,j,k) q2(i,j,k)=uf(i,j,k) q2lb(i,j,k)=q2l(i,j,k) q2l(i,j,k)=vf(i,j,k) end do end do end do C C Calculate tf and sf using uf, vf, a and c as temporary variables: C if(mode.ne.4) then ! ! !----------------------------------------------------------------------! !lyo:!wad:The followings prepare Oey's [1996] and Oey & Chen's [1992] ! ! implementations of river and temperature surface conditions ! ! i.e. using wriv -- but not implemented yet; they are not ! ! directly for WAD, but ssurfwet is required later for wad ! ! ! do j=1,jm do i=1,im ssurfwet(i,j)=SSURF(i,j)*(1.-WRIV(i,j)/( WRIV(i,j)+ $ 1.E-28) ) enddo enddo ! ! C if(nadv.eq.1) then C call advt1(tb,t,tclim,uf,a,c) call advt1(sb,s,sclim,vf,a,c) C else if(nadv.eq.2) then C call advt2(tb,t,tclim,uf,a,c,nitera,sw) call advt2(sb,s,sclim,vf,a,c,nitera,sw) C else C write(6,9) 9 format(/'Invalid value for nadv ..... ', $ 'program terminated'/) stop C endif C call proft(uf,wtsurf,tsurf,nbct,tps) call proft(vf,wssurf,ssurf,nbcs,tps) call bcond(4) ! ! !----------------------------------------------------------------------! !lyo:!wad:Prepare dry-cell T/S "initial conditions" for next time-step ! ! ! if (nwad.eq.1) then do k=1,kb-1; do j=1,jm; do i=1,im if (wetmask(i,j).eq.0.0) then uf(i,j,k)=(cwetrlx1*tb(i,j,k)+cwetrlx2*tsurf(i,j) $ *fsm(i,j)) vf(i,j,k)=(cwetrlx1*sb(i,j,k)+cwetrlx2*ssurfwet(i,j) $ *fsm(i,j)) endif enddo; enddo; enddo endif !----------------------------------------------------------------------! C do k=1,kb do j=1,jm do i=1,im t(i,j,k)=t(i,j,k) $ +.5e0*smoth*(uf(i,j,k)+tb(i,j,k) $ -2.e0*t(i,j,k)) s(i,j,k)=s(i,j,k) $ +.5e0*smoth*(vf(i,j,k)+sb(i,j,k) $ -2.e0*s(i,j,k)) tb(i,j,k)=t(i,j,k) t(i,j,k)=uf(i,j,k) sb(i,j,k)=s(i,j,k) s(i,j,k)=vf(i,j,k) end do end do end do C call dens(s,t,rho) C endif C C Calculate uf and vf: C call advu call advv call profu call profv call bcond(3) ! ! !----------------------------------------------------------------------! !lyo:!wad:Check wet/dry on UF & VF: ! ! ! if (nwad.eq.1) then do k=1,kb-1 do j=1,jm; do i=2,im if (wetmask(i,j)*wetmask(i-1,j).eq.0.0)uf(i,j,k)=0.0 enddo; enddo do j=2,jm; do i=1,im if (wetmask(i,j)*wetmask(i,j-1).eq.0.0)vf(i,j,k)=0.0 enddo; enddo enddo endif !----------------------------------------------------------------------! C do j=1,jm do i=1,im tps(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tps(i,j)=tps(i,j) $ +(uf(i,j,k)+ub(i,j,k)-2.e0*u(i,j,k))*dz(k) end do end do end do C do k=1,kbm1 do j=1,jm do i=1,im u(i,j,k)=u(i,j,k) $ +.5e0*smoth*(uf(i,j,k)+ub(i,j,k) $ -2.e0*u(i,j,k)-tps(i,j)) end do end do end do C do j=1,jm do i=1,im tps(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tps(i,j)=tps(i,j) $ +(vf(i,j,k)+vb(i,j,k)-2.e0*v(i,j,k))*dz(k) end do end do end do C do k=1,kbm1 do j=1,jm do i=1,im v(i,j,k)=v(i,j,k) $ +.5e0*smoth*(vf(i,j,k)+vb(i,j,k) $ -2.e0*v(i,j,k)-tps(i,j)) end do end do end do C do k=1,kb do j=1,jm do i=1,im ub(i,j,k)=u(i,j,k) u(i,j,k)=uf(i,j,k) vb(i,j,k)=v(i,j,k) v(i,j,k)=vf(i,j,k) end do end do end do C endif C do j=1,jm do i=1,im egb(i,j)=egf(i,j) etb(i,j)=et(i,j) et(i,j)=etf(i,j) dt(i,j)=h(i,j)+et(i,j) utb(i,j)=utf(i,j) vtb(i,j)=vtf(i,j) vfluxb(i,j)=vfluxf(i,j) end do end do C endif C C----------------------------------------------------------------------- C C Beginning of print section: C if(iint.ge.iswtch) iprint=nint(prtd2*24.e0*3600.e0/dti) C if(mod(iint,iprint).eq.0.or.vamax.gt.vmaxl) then C write(6,4) time,iint,iext,iprint 4 format(/ $ '**************************************************', $ '**************************************************', $ '*************************'// $ ' time =',f9.4,', iint =',i8,', iext =',i8,', iprint =',i8,//) C C Select print statements in printall as desired: C call printall C call wadout !lyo:!wad: ! vtot=0.e0 atot=0.e0 taver=0.e0 saver=0.e0 eaver=0.e0 do k=1,kbm1 do j=1,jm do i=1,im darea=dx(i,j)*dy(i,j)*wetmask(i,j) !lyo:!wad: dvol=darea*dt(i,j)*dz(k) vtot=vtot+dvol taver=taver+tb(i,j,k)*dvol saver=saver+sb(i,j,k)*dvol end do end do end do C do j=1,jm do i=1,im darea=dx(i,j)*dy(i,j)*wetmask(i,j) !lyo:!wad: atot=atot+darea eaver=eaver+et(i,j)*darea end do end do C taver=taver/vtot saver=saver/vtot eaver=eaver/atot tsalt=(saver+sbias)*vtot C write(6,5) vtot,atot,eaver,taver,saver,tsalt 5 format('vtot = ',e16.7,' atot = ',e16.7, $ ' eaver =',e16.7/'taver =',e16.7, $ ' saver =',e16.7,' tsalt =',e16.7) C C Write netCDF output: C if(netcdf_file.ne.'nonetcdf') then call write_netcdf(netcdf_file,2) ! *netCDF* endif C if(vamax.gt.vmaxl) then C write(6,4) time,iint,iext,iprint C call printall C write(6,6) vamax,imax,jmax 6 format(/////////////////// $ '************************************************'/ $ '************ abnormal job end ******************'/ $ '************* user terminated ******************'/ $ '************************************************'/ $ ' vamax =',e12.3,' imax,jmax =',2i5) C C Close netCDF file: C if(netcdf_file.ne.'nonetcdf') then call write_netcdf(netcdf_file,3) ! *netCDF* endif C stop C endif C endif C C End of print section C C----------------------------------------------------------------------- C 9000 continue ! End of internal (3-D) mode C C----------------------------------------------------------------------- C write(6,4) time,iint,iext,iprint ! call printall !lyo:_20080415:final printing ! C C Set levels for output: C ko(1)=1 ko(2)=2 ko(3)=kb/2 ko(4)=kb-1 ko(5)=kb C call prxyz('Vertical velocity, w ', $ time,w ,im,iskp,jm,jskp,kb,ko,5,-1.e0) C C call prxyz('Turbulent kinetic energy x 2, q2 ', C $ time,q2 ,im,iskp,jm,jskp,kb,ko,5,-1.e0) C C Save this data for a seamless restart: C clyo: ! call powsave c write(71) time, $ wubot,wvbot,aam2d,ua,uab,va,vab,el,elb,et,etb,egb, $ utb,vtb,u,ub,w,v,vb,t,tb,s,sb,rho,adx2d,ady2d,advua,advva, $ km,kh,kq,l,q2,q2b,aam,q2l,q2lb $ ,wetmask,wmarsh !lyo:!wad: C C Close netCDF file: C if(netcdf_file.ne.'nonetcdf') then call write_netcdf(netcdf_file,3) ! *netCDF* endif C !lyo:!wad:print final successful completion: write(6,10) time 10 format(/2x,'JOB SUCCESSFULLY COMPLT.; time = ',1P1e13.5,' days') ! stop C end C C End of main program C C----------------------------------------------------------------------- C subroutine advave(curv2d) C ********************************************************************** C * * C * FUNCTION : Calculates horizontal advection and diffusion. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real curv2d(im,jm) integer i,j C C u-advection and diffusion: C C Advective fluxes: C do j=1,jm do i=1,im advua(i,j)=0.e0 end do end do C do j=2,jm do i=2,imm1 fluxua(i,j)=.125e0*((d(i+1,j)+d(i,j))*ua(i+1,j) $ +(d(i,j)+d(i-1,j))*ua(i,j)) $ *(ua(i+1,j)+ua(i,j)) end do end do C do j=2,jm do i=2,im fluxva(i,j)=.125e0*((d(i,j)+d(i,j-1))*va(i,j) $ +(d(i-1,j)+d(i-1,j-1))*va(i-1,j)) $ *(ua(i,j)+ua(i,j-1)) end do end do C C Add viscous fluxes: C do j=2,jm do i=2,imm1 fluxua(i,j)=fluxua(i,j) $ -d(i,j)*2.e0*aam2d(i,j)*(uab(i+1,j)-uab(i,j)) $ /dx(i,j) end do end do C do j=2,jm do i=2,im tps(i,j)=.25e0*(d(i,j)+d(i-1,j)+d(i,j-1)+d(i-1,j-1)) $ *(aam2d(i,j)+aam2d(i,j-1) $ +aam2d(i-1,j)+aam2d(i-1,j-1)) $ *((uab(i,j)-uab(i,j-1)) $ /(dy(i,j)+dy(i-1,j)+dy(i,j-1)+dy(i-1,j-1)) $ +(vab(i,j)-vab(i-1,j)) $ /(dx(i,j)+dx(i-1,j)+dx(i,j-1)+dx(i-1,j-1))) fluxua(i,j)=fluxua(i,j)*dy(i,j) fluxva(i,j)=(fluxva(i,j)-tps(i,j))*.25e0 $ *(dx(i,j)+dx(i-1,j)+dx(i,j-1)+dx(i-1,j-1)) end do end do C do j=2,jmm1 do i=2,imm1 advua(i,j)=fluxua(i,j)-fluxua(i-1,j) $ +fluxva(i,j+1)-fluxva(i,j) end do end do C C u-advection and diffusion: C do j=1,jm do i=1,im advva(i,j)=0.e0 end do end do C C Advective fluxes: C do j=2,jm do i=2,im fluxua(i,j)=.125e0*((d(i,j)+d(i-1,j))*ua(i,j) $ +(d(i,j-1)+d(i-1,j-1))*ua(i,j-1)) $ *(va(i-1,j)+va(i,j)) end do end do C do j=2,jmm1 do i=2,im fluxva(i,j)=.125e0*((d(i,j+1)+d(i,j))*va(i,j+1) $ +(d(i,j)+d(i,j-1))*va(i,j)) $ *(va(i,j+1)+va(i,j)) end do end do C C Add viscous fluxes: C do j=2,jmm1 do i=2,im fluxva(i,j)=fluxva(i,j) $ -d(i,j)*2.e0*aam2d(i,j)*(vab(i,j+1)-vab(i,j)) $ /dy(i,j) end do end do C do j=2,jm do i=2,im fluxva(i,j)=fluxva(i,j)*dx(i,j) fluxua(i,j)=(fluxua(i,j)-tps(i,j))*.25e0 $ *(dy(i,j)+dy(i-1,j)+dy(i,j-1)+dy(i-1,j-1)) end do end do C do j=2,jmm1 do i=2,imm1 advva(i,j)=fluxua(i+1,j)-fluxua(i,j) $ +fluxva(i,j)-fluxva(i,j-1) end do end do C if(mode.eq.2) then C do j=2,jmm1 do i=2,imm1 wubot(i,j)=-0.5e0*(cbc(i,j)+cbc(i-1,j)) $ *sqrt(uab(i,j)**2 $ +(.25e0*(vab(i,j)+vab(i,j+1) $ +vab(i-1,j)+vab(i-1,j+1)))**2) $ *uab(i,j) end do end do C do j=2,jmm1 do i=2,imm1 wvbot(i,j)=-0.5e0*(cbc(i,j)+cbc(i,j-1)) $ *sqrt(vab(i,j)**2 $ +(.25e0*(uab(i,j)+uab(i+1,j) $ +uab(i,j-1)+uab(i+1,j-1)))**2) $ *vab(i,j) end do end do C do j=2,jmm1 do i=2,imm1 curv2d(i,j)=.25e0 $ *((va(i,j+1)+va(i,j))*(dy(i+1,j)-dy(i-1,j)) $ -(ua(i+1,j)+ua(i,j))*(dx(i,j+1)-dx(i,j-1))) $ /(dx(i,j)*dy(i,j)) end do end do C do j=2,jmm1 do i=3,imm1 advua(i,j)=advua(i,j)-aru(i,j)*.25e0 $ *(curv2d(i,j)*d(i,j) $ *(va(i,j+1)+va(i,j)) $ +curv2d(i-1,j)*d(i-1,j) $ *(va(i-1,j+1)+va(i-1,j))) end do end do C do j=3,jmm1 do i=2,imm1 advva(i,j)=advva(i,j)+arv(i,j)*.25e0 $ *(curv2d(i,j)*d(i,j) $ *(ua(i+1,j)+ua(i,j)) $ +curv2d(i,j-1)*d(i,j-1) $ *(ua(i+1,j-1)+ua(i,j-1))) end do end do C endif C return C end C subroutine advct(xflux,yflux,curv) C ********************************************************************** C * * C * FUNCTION : Calculates the horizontal portions of momentum * C * advection well in advance of their use in advu and * C * advv so that their vertical integrals (created in * C * the main program) may be used in the external (2-D) * C * mode calculation. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real xflux(im,jm,kb),yflux(im,jm,kb) real curv(im,jm,kb) real dtaam integer i,j,k C do k=1,kb do j=1,jm do i=1,im curv(i,j,k)=0.e0 advx(i,j,k)=0.e0 xflux(i,j,k)=0.e0 yflux(i,j,k)=0.e0 end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 curv(i,j,k)=.25e0*((v(i,j+1,k)+v(i,j,k)) $ *(dy(i+1,j)-dy(i-1,j)) $ -(u(i+1,j,k)+u(i,j,k)) $ *(dx(i,j+1)-dx(i,j-1))) $ /(dx(i,j)*dy(i,j)) end do end do end do C C Calculate x-component of velocity advection: C C Calculate horizontal advective fluxes: C do k=1,kbm1 do j=1,jm do i=2,imm1 xflux(i,j,k)=.125e0*((dt(i+1,j)+dt(i,j))*u(i+1,j,k) $ +(dt(i,j)+dt(i-1,j))*u(i,j,k)) $ *(u(i+1,j,k)+u(i,j,k)) end do end do end do C do k=1,kbm1 do j=2,jm do i=2,im yflux(i,j,k)=.125e0*((dt(i,j)+dt(i,j-1))*v(i,j,k) $ +(dt(i-1,j)+dt(i-1,j-1))*v(i-1,j,k)) $ *(u(i,j,k)+u(i,j-1,k)) end do end do end do C C Add horizontal diffusive fluxes: C do k=1,kbm1 do j=2,jm do i=2,imm1 xflux(i,j,k)=xflux(i,j,k) $ -dt(i,j)*aam(i,j,k)*2.e0 $ *(ub(i+1,j,k)-ub(i,j,k))/dx(i,j) dtaam=.25e0*(dt(i,j)+dt(i-1,j)+dt(i,j-1)+dt(i-1,j-1)) $ *(aam(i,j,k)+aam(i-1,j,k) $ +aam(i,j-1,k)+aam(i-1,j-1,k)) yflux(i,j,k)=yflux(i,j,k) $ -dtaam*((ub(i,j,k)-ub(i,j-1,k)) $ /(dy(i,j)+dy(i-1,j) $ +dy(i,j-1)+dy(i-1,j-1)) $ +(vb(i,j,k)-vb(i-1,j,k)) $ /(dx(i,j)+dx(i-1,j) $ +dx(i,j-1)+dx(i-1,j-1))) C xflux(i,j,k)=dy(i,j)*xflux(i,j,k) yflux(i,j,k)=.25e0*(dx(i,j)+dx(i-1,j) $ +dx(i,j-1)+dx(i-1,j-1))*yflux(i,j,k) end do end do end do C C Do horizontal advection: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 advx(i,j,k)=xflux(i,j,k)-xflux(i-1,j,k) $ +yflux(i,j+1,k)-yflux(i,j,k) end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=3,imm1 advx(i,j,k)=advx(i,j,k) $ -aru(i,j)*.25e0 $ *(curv(i,j,k)*dt(i,j) $ *(v(i,j+1,k)+v(i,j,k)) $ +curv(i-1,j,k)*dt(i-1,j) $ *(v(i-1,j+1,k)+v(i-1,j,k))) end do end do end do C C----------------------------------------------------------------------- C do k=1,kb do j=1,jm do i=1,im advy(i,j,k)=0.e0 xflux(i,j,k)=0.e0 yflux(i,j,k)=0.e0 end do end do end do C C Calculate y-component of velocity advection: C C Calculate horizontal advective fluxes: C do k=1,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=.125e0*((dt(i,j)+dt(i-1,j))*u(i,j,k) $ +(dt(i,j-1)+dt(i-1,j-1))*u(i,j-1,k)) $ *(v(i,j,k)+v(i-1,j,k)) end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=1,im yflux(i,j,k)=.125e0*((dt(i,j+1)+dt(i,j))*v(i,j+1,k) $ +(dt(i,j)+dt(i,j-1))*v(i,j,k)) $ *(v(i,j+1,k)+v(i,j,k)) end do end do end do C C Add horizontal diffusive fluxes: C do k=1,kbm1 do j=2,jmm1 do i=2,im dtaam=.25e0*(dt(i,j)+dt(i-1,j)+dt(i,j-1)+dt(i-1,j-1)) $ *(aam(i,j,k)+aam(i-1,j,k) $ +aam(i,j-1,k)+aam(i-1,j-1,k)) xflux(i,j,k)=xflux(i,j,k) $ -dtaam*((ub(i,j,k)-ub(i,j-1,k)) $ /(dy(i,j)+dy(i-1,j) $ +dy(i,j-1)+dy(i-1,j-1)) $ +(vb(i,j,k)-vb(i-1,j,k)) $ /(dx(i,j)+dx(i-1,j) $ +dx(i,j-1)+dx(i-1,j-1))) yflux(i,j,k)=yflux(i,j,k) $ -dt(i,j)*aam(i,j,k)*2.e0 $ *(vb(i,j+1,k)-vb(i,j,k))/dy(i,j) C xflux(i,j,k)=.25e0*(dy(i,j)+dy(i-1,j) $ +dy(i,j-1)+dy(i-1,j-1))*xflux(i,j,k) yflux(i,j,k)=dx(i,j)*yflux(i,j,k) end do end do end do C C Do horizontal advection: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 advy(i,j,k)=xflux(i+1,j,k)-xflux(i,j,k) $ +yflux(i,j,k)-yflux(i,j-1,k) end do end do end do C do k=1,kbm1 do j=3,jmm1 do i=2,imm1 advy(i,j,k)=advy(i,j,k) $ +arv(i,j)*.25e0 $ *(curv(i,j,k)*dt(i,j) $ *(u(i+1,j,k)+u(i,j,k)) $ +curv(i,j-1,k)*dt(i,j-1) $ *(u(i+1,j-1,k)+u(i,j-1,k))) end do end do end do C return C end C subroutine advq(qb,q,qf,xflux,yflux) C ********************************************************************** C * * C * FUNCTION : Calculates horizontal advection and diffusion, and * C * vertical advection for turbulent quantities. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real qb(im,jm,kb),q(im,jm,kb),qf(im,jm,kb) real xflux(im,jm,kb),yflux(im,jm,kb) integer i,j,k C C Do horizontal advection: C do k=2,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=.125e0*(q(i,j,k)+q(i-1,j,k)) $ *(dt(i,j)+dt(i-1,j))*(u(i,j,k)+u(i,j,k-1)) yflux(i,j,k)=.125e0*(q(i,j,k)+q(i,j-1,k)) $ *(dt(i,j)+dt(i,j-1))*(v(i,j,k)+v(i,j,k-1)) end do end do end do C C Do horizontal diffusion: C do k=2,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=xflux(i,j,k) $ -.25e0*(aam(i,j,k)+aam(i-1,j,k) $ +aam(i,j,k-1)+aam(i-1,j,k-1)) $ *(h(i,j)+h(i-1,j)) $ *(qb(i,j,k)-qb(i-1,j,k))*dum(i,j) $ /(dx(i,j)+dx(i-1,j)) yflux(i,j,k)=yflux(i,j,k) $ -.25e0*(aam(i,j,k)+aam(i,j-1,k) $ +aam(i,j,k-1)+aam(i,j-1,k-1)) $ *(h(i,j)+h(i,j-1)) $ *(qb(i,j,k)-qb(i,j-1,k))*dvm(i,j) $ /(dy(i,j)+dy(i,j-1)) xflux(i,j,k)=.5e0*(dy(i,j)+dy(i-1,j))*xflux(i,j,k) yflux(i,j,k)=.5e0*(dx(i,j)+dx(i,j-1))*yflux(i,j,k) end do end do end do C C Do vertical advection, add flux terms, then step forward in time: C do k=2,kbm1 do j=2,jmm1 do i=2,imm1 qf(i,j,k)=(w(i,j,k-1)*q(i,j,k-1)-w(i,j,k+1)*q(i,j,k+1)) $ *art(i,j)/(dz(k)+dz(k-1)) $ +xflux(i+1,j,k)-xflux(i,j,k) $ +yflux(i,j+1,k)-yflux(i,j,k) qf(i,j,k)=((h(i,j)+etb(i,j))*art(i,j) $ *qb(i,j,k)-dti2*qf(i,j,k)) $ /((h(i,j)+etf(i,j))*art(i,j)) end do end do end do C return C end C subroutine advt1(fb,f,fclim,ff,xflux,yflux) C ********************************************************************** C * * C * FUNCTION : Integrates conservative scalar equations. * C * * C * This is centred scheme, as originally provide in * C * POM (previously called advt). * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real fb(im,jm,kb),f(im,jm,kb),fclim(im,jm,kb),ff(im,jm,kb) real xflux(im,jm,kb),yflux(im,jm,kb) integer i,j,k C do j=1,jm do i=1,im f(i,j,kb)=f(i,j,kbm1) fb(i,j,kb)=fb(i,j,kbm1) end do end do C C Do advective fluxes: C do k=1,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=.25e0*((dt(i,j)+dt(i-1,j)) $ *(f(i,j,k)+f(i-1,j,k))*u(i,j,k)) yflux(i,j,k)=.25e0*((dt(i,j)+dt(i,j-1)) $ *(f(i,j,k)+f(i,j-1,k))*v(i,j,k)) end do end do end do C C Add diffusive fluxes: C do k=1,kb do j=1,jm do i=1,im fb(i,j,k)=fb(i,j,k)-fclim(i,j,k) end do end do end do C do k=1,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=xflux(i,j,k) $ -.5e0*(aam(i,j,k)+aam(i-1,j,k)) $ *(h(i,j)+h(i-1,j))*tprni $ *(fb(i,j,k)-fb(i-1,j,k))*dum(i,j) $ /(dx(i,j)+dx(i-1,j)) yflux(i,j,k)=yflux(i,j,k) $ -.5e0*(aam(i,j,k)+aam(i,j-1,k)) $ *(h(i,j)+h(i,j-1))*tprni $ *(fb(i,j,k)-fb(i,j-1,k))*dvm(i,j) $ /(dy(i,j)+dy(i,j-1)) xflux(i,j,k)=.5e0*(dy(i,j)+dy(i-1,j))*xflux(i,j,k) yflux(i,j,k)=.5e0*(dx(i,j)+dx(i,j-1))*yflux(i,j,k) end do end do end do C do k=1,kb do j=1,jm do i=1,im fb(i,j,k)=fb(i,j,k)+fclim(i,j,k) end do end do end do C C Do vertical advection: C do j=2,jmm1 do i=2,imm1 zflux(i,j,1)=f(i,j,1)*w(i,j,1)*art(i,j) zflux(i,j,kb)=0.e0 end do end do C do k=2,kbm1 do j=2,jmm1 do i=2,imm1 zflux(i,j,k)=.5e0*(f(i,j,k-1)+f(i,j,k))*w(i,j,k)*art(i,j) end do end do end do C C Add net horizontal fluxes and then step forward in time: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 ff(i,j,k)=xflux(i+1,j,k)-xflux(i,j,k) $ +yflux(i,j+1,k)-yflux(i,j,k) $ +(zflux(i,j,k)-zflux(i,j,k+1))/dz(k) C ff(i,j,k)=(fb(i,j,k)*(h(i,j)+etb(i,j))*art(i,j) $ -dti2*ff(i,j,k)) $ /((h(i,j)+etf(i,j))*art(i,j)) end do end do end do C return C end C subroutine advt2(fb,f,fclim,ff,xflux,yflux,nitera,sw) C ********************************************************************** C * * C * FUNCTION : Integrates conservative scalar equations. * C * * C * This is a first-order upstream scheme, which * C * reduces implicit diffusion using the Smolarkiewicz * C * iterative upstream scheme with an antidiffusive * C * velocity. * C * * C * It is based on the subroutines of Gianmaria Sannino * C * (Inter-university Computing Consortium, Rome, Italy)* C * and Vincenzo Artale (Italian National Agency for * C * New Technology and Environment, Rome, Italy), * C * downloaded from the POM FTP site on 1 Nov. 2001. * C * The calculations have been simplified by removing * C * the shock switch option. It should be noted that * C * this implementation does not include cross-terms * C * which are in the original formulation. * C * * C * fb,f,fclim,ff . as used in subroutine advt1 * C * xflux,yflux ... working arrays used to save memory * C * nitera ........ number of iterations. This should * C * be in the range 1 - 4. 1 is * C * standard upstream differencing; * C * 3 adds 50% CPU time to POM. * C * sw ............ smoothing parameter. This should * C * preferably be 1, but 0 < sw < 1 * C * gives smoother solutions with less * C * overshoot when nitera > 1. * C * * C * Reference: * C * * C * Smolarkiewicz, P.K.; A fully multidimensional * C * positive definite advection transport algorithm * C * with small implicit diffusion, Journal of * C * Computational Physics, 54, 325-362, 1984. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real fb(im,jm,kb),f(im,jm,kb),fclim(im,jm,kb),ff(im,jm,kb) real xflux(im,jm,kb),yflux(im,jm,kb) real sw integer nitera real fbmem(im,jm,kb),eta(im,jm) real xmassflux(im,jm,kb),ymassflux(im,jm,kb),zwflux(im,jm,kb) integer i,j,k,itera C C Calculate horizontal mass fluxes: C do k=1,kb do j=1,jm do i=1,im xmassflux(i,j,k)=0.e0 ymassflux(i,j,k)=0.e0 end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=2,im xmassflux(i,j,k)=0.25e0*(dy(i-1,j)+dy(i,j)) $ *(dt(i-1,j)+dt(i,j))*u(i,j,k) end do end do C do j=2,jm do i=2,imm1 ymassflux(i,j,k)=0.25e0*(dx(i,j-1)+dx(i,j)) $ *(dt(i,j-1)+dt(i,j))*v(i,j,k) end do end do end do C do j=1,jm do i=1,im fb(i,j,kb)=fb(i,j,kbm1) eta(i,j)=etb(i,j) end do end do C do k=1,kb do j=1,jm do i=1,im zwflux(i,j,k)=w(i,j,k) fbmem(i,j,k)=fb(i,j,k) end do end do end do C C Start Smolarkiewicz scheme: C do itera=1,nitera C C Upwind advection scheme: C do k=1,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=0.5e0 $ *((xmassflux(i,j,k)+abs(xmassflux(i,j,k))) $ *fbmem(i-1,j,k)+ $ (xmassflux(i,j,k)-abs(xmassflux(i,j,k))) $ *fbmem(i,j,k)) C yflux(i,j,k)=0.5e0 $ *((ymassflux(i,j,k)+abs(ymassflux(i,j,k))) $ *fbmem(i,j-1,k)+ $ (ymassflux(i,j,k)-abs(ymassflux(i,j,k))) $ *fbmem(i,j,k)) end do end do end do C do j=2,jmm1 do i=2,imm1 zflux(i,j,1)=0.e0 if(itera.eq.1) zflux(i,j,1)=w(i,j,1)*f(i,j,1)*art(i,j) zflux(i,j,kb)=0.e0 end do end do C do k=2,kbm1 do j=2,jmm1 do i=2,imm1 zflux(i,j,k)=0.5e0 $ *((zwflux(i,j,k)+abs(zwflux(i,j,k))) $ *fbmem(i,j,k)+ $ (zwflux(i,j,k)-abs(zwflux(i,j,k))) $ *fbmem(i,j,k-1)) zflux(i,j,k)=zflux(i,j,k)*art(i,j) end do end do end do C C Add net advective fluxes and step forward in time: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 ff(i,j,k)=xflux(i+1,j,k)-xflux(i,j,k) $ +yflux(i,j+1,k)-yflux(i,j,k) $ +(zflux(i,j,k)-zflux(i,j,k+1))/dz(k) ff(i,j,k)=(fbmem(i,j,k)*(h(i,j)+eta(i,j))*art(i,j) $ -dti2*ff(i,j,k))/((h(i,j)+etf(i,j))*art(i,j)) end do end do end do C C Calculate antidiffusion velocity: C call smol_adif(xmassflux,ymassflux,zwflux,ff,sw) C do j=1,jm do i=1,im eta(i,j)=etf(i,j) do k=1,kb fbmem(i,j,k)=ff(i,j,k) end do end do end do C C End of Smolarkiewicz scheme C end do C C Add horizontal diffusive fluxes: C do k=1,kb do j=1,jm do i=1,im fb(i,j,k)=fb(i,j,k)-fclim(i,j,k) end do end do end do C do k=1,kbm1 do j=2,jm do i=2,im xmassflux(i,j,k)=0.5e0*(aam(i,j,k)+aam(i-1,j,k)) ymassflux(i,j,k)=0.5e0*(aam(i,j,k)+aam(i,j-1,k)) end do end do end do C do k=1,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=-xmassflux(i,j,k)*(h(i,j)+h(i-1,j))*tprni $ *(fb(i,j,k)-fb(i-1,j,k))*dum(i,j) $ *(dy(i,j)+dy(i-1,j))*0.5e0/(dx(i,j)+dx(i-1,j)) yflux(i,j,k)=-ymassflux(i,j,k)*(h(i,j)+h(i,j-1))*tprni $ *(fb(i,j,k)-fb(i,j-1,k))*dvm(i,j) $ *(dx(i,j)+dx(i,j-1))*0.5e0/(dy(i,j)+dy(i,j-1)) end do end do end do C do k=1,kb do j=1,jm do i=1,im fb(i,j,k)=fb(i,j,k)+fclim(i,j,k) end do end do end do C C Add net horizontal fluxes and step forward in time: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 ff(i,j,k)=ff(i,j,k)-dti2*(xflux(i+1,j,k)-xflux(i,j,k) $ +yflux(i,j+1,k)-yflux(i,j,k)) $ /((h(i,j)+etf(i,j))*art(i,j)) end do end do end do C return C end C subroutine advu C ********************************************************************** C * * C * ROUTINE NAME: advu * C * * C * FUNCTION : Does horizontal and vertical advection of * C * u-momentum, and includes coriolis, surface slope * C * and baroclinic terms. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C integer i,j,k C C Do vertical advection: C do k=1,kb do j=1,jm do i=1,im uf(i,j,k)=0.e0 end do end do end do C do k=2,kbm1 do j=1,jm do i=2,im uf(i,j,k)=.25e0*(w(i,j,k)+w(i-1,j,k)) $ *(u(i,j,k)+u(i,j,k-1)) end do end do end do C C Combine horizontal and vertical advection with coriolis, surface C slope and baroclinic terms: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 uf(i,j,k)=advx(i,j,k) $ +(uf(i,j,k)-uf(i,j,k+1))*aru(i,j)/dz(k) $ -aru(i,j)*.25e0 $ *(cor(i,j)*dt(i,j) $ *(v(i,j+1,k)+v(i,j,k)) $ +cor(i-1,j)*dt(i-1,j) $ *(v(i-1,j+1,k)+v(i-1,j,k))) $ +grav*.125e0*(dt(i,j)+dt(i-1,j)) $ *(egf(i,j)-egf(i-1,j)+egb(i,j)-egb(i-1,j) $ +(e_atmos(i,j)-e_atmos(i-1,j))*2.e0) $ *(dy(i,j)+dy(i-1,j)) $ +drhox(i,j,k) end do end do end do C C Step forward in time: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 uf(i,j,k)=((h(i,j)+etb(i,j)+h(i-1,j)+etb(i-1,j)) $ *aru(i,j)*ub(i,j,k) $ -2.e0*dti2*uf(i,j,k)) $ /((h(i,j)+etf(i,j)+h(i-1,j)+etf(i-1,j)) $ *aru(i,j)) end do end do end do C return C end C subroutine advv C ********************************************************************** C * * C * FUNCTION : Does horizontal and vertical advection of * C * v-momentum, and includes coriolis, surface slope * C * and baroclinic terms. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C integer i,j,k C C Do vertical advection: C do k=1,kb do j=1,jm do i=1,im vf(i,j,k)=0.e0 end do end do end do C do k=2,kbm1 do j=2,jm do i=1,im vf(i,j,k)=.25e0*(w(i,j,k)+w(i,j-1,k)) $ *(v(i,j,k)+v(i,j,k-1)) end do end do end do C C Combine horizontal and vertical advection with coriolis, surface C slope and baroclinic terms: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 vf(i,j,k)=advy(i,j,k) $ +(vf(i,j,k)-vf(i,j,k+1))*arv(i,j)/dz(k) $ +arv(i,j)*.25e0 $ *(cor(i,j)*dt(i,j) $ *(u(i+1,j,k)+u(i,j,k)) $ +cor(i,j-1)*dt(i,j-1) $ *(u(i+1,j-1,k)+u(i,j-1,k))) $ +grav*.125e0*(dt(i,j)+dt(i,j-1)) $ *(egf(i,j)-egf(i,j-1)+egb(i,j)-egb(i,j-1) $ +(e_atmos(i,j)-e_atmos(i,j-1))*2.e0) $ *(dx(i,j)+dx(i,j-1)) $ +drhoy(i,j,k) end do end do end do C C Step forward in time: C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 vf(i,j,k)=((h(i,j)+etb(i,j)+h(i,j-1)+etb(i,j-1)) $ *arv(i,j)*vb(i,j,k) $ -2.e0*dti2*vf(i,j,k)) $ /((h(i,j)+etf(i,j)+h(i,j-1)+etf(i,j-1)) $ *arv(i,j)) end do end do end do C return C end C subroutine areas_masks C ********************************************************************** C * * C * FUNCTION : Calculates areas and masks. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C integer i,j C C Calculate areas of "t" and "s" cells: C do j=1,jm do i=1,im art(i,j)=dx(i,j)*dy(i,j) end do end do C C Calculate areas of "u" and "v" cells: C do j=2,jm do i=2,im aru(i,j)=.25e0*(dx(i,j)+dx(i-1,j))*(dy(i,j)+dy(i-1,j)) arv(i,j)=.25e0*(dx(i,j)+dx(i,j-1))*(dy(i,j)+dy(i,j-1)) end do end do C do j=1,jm aru(1,j)=aru(2,j) arv(1,j)=arv(2,j) end do C do i=1,im aru(i,1)=aru(i,2) arv(i,1)=arv(i,2) end do C C Initialise and set up free surface mask if no WAD: !lyo:!wad: C if (nwad.eq.0) then !lyo:!wad: ! ! do j=1,jm ! do i=1,im ! fsm(i,j)=0.e0 ! dum(i,j)=0.e0 ! dvm(i,j)=0.e0 ! if(h(i,j).gt.1.e0) fsm(i,j)=1.e0 ! end do ! end do C C Set up velocity masks: C ! do j=2,jm ! do i=2,im ! dum(i,j)=fsm(i,j)*fsm(i-1,j) ! dvm(i,j)=fsm(i,j)*fsm(i,j-1) ! end do ! end do ! !lyo:_20080415:the above original pom2k version does NOT define dum & dvm ! along i=j=1; use the folliwngs (from subroutine wadh). do j=1,jm; do i=1,im FSM(I,J)=1.; DUM(I,J)=1.; DVM(I,J)=1. IF(h(I,J).LE.1.0)THEN h(I,J)=1.0; FSM(I,J)=0.; DUM(I,J)=0.; DVM(I,J)=0. ENDIF enddo; enddo DO J=1,JM-1; DO I=1,IM IF(FSM(I,J).EQ.0..AND.FSM(I,J+1).NE.0.)DVM(I,J+1)=0. enddo; enddo DO J=1,JM; DO I=1,IM-1 IF(FSM(I,J).EQ.0..AND.FSM(I+1,J).NE.0.)DUM(I+1,J)=0. enddo; enddo ! endif !if (nwad.eq.0) then !lyo:!wad: C return C end C subroutine baropg C ********************************************************************** C * * C * FUNCTION : Calculates baroclinic pressure gradient. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C integer i,j,k C do k=1,kb do j=1,jm do i=1,im rho(i,j,k)=rho(i,j,k)-rmean(i,j,k) end do end do end do C C Calculate x-component of baroclinic pressure gradient: C do j=2,jmm1 do i=2,imm1 drhox(i,j,1)=.5e0*grav*(-zz(1))*(dt(i,j)+dt(i-1,j)) $ *(rho(i,j,1)-rho(i-1,j,1)) end do end do C do k=2,kbm1 do j=2,jmm1 do i=2,imm1 drhox(i,j,k)=drhox(i,j,k-1) $ +grav*.25e0*(zz(k-1)-zz(k)) $ *(dt(i,j)+dt(i-1,j)) $ *(rho(i,j,k)-rho(i-1,j,k) $ +rho(i,j,k-1)-rho(i-1,j,k-1)) $ +grav*.25e0*(zz(k-1)+zz(k)) $ *(dt(i,j)-dt(i-1,j)) $ *(rho(i,j,k)+rho(i-1,j,k) $ -rho(i,j,k-1)-rho(i-1,j,k-1)) end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 drhox(i,j,k)=.25e0*(dt(i,j)+dt(i-1,j)) $ *drhox(i,j,k)*dum(i,j) $ *(dy(i,j)+dy(i-1,j)) end do end do end do C C Calculate y-component of baroclinic pressure gradient: C do j=2,jmm1 do i=2,imm1 drhoy(i,j,1)=.5e0*grav*(-zz(1))*(dt(i,j)+dt(i,j-1)) $ *(rho(i,j,1)-rho(i,j-1,1)) end do end do C do k=2,kbm1 do j=2,jmm1 do i=2,imm1 drhoy(i,j,k)=drhoy(i,j,k-1) $ +grav*.25e0*(zz(k-1)-zz(k)) $ *(dt(i,j)+dt(i,j-1)) $ *(rho(i,j,k)-rho(i,j-1,k) $ +rho(i,j,k-1)-rho(i,j-1,k-1)) $ +grav*.25e0*(zz(k-1)+zz(k)) $ *(dt(i,j)-dt(i,j-1)) $ *(rho(i,j,k)+rho(i,j-1,k) $ -rho(i,j,k-1)-rho(i,j-1,k-1)) end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 drhoy(i,j,k)=.25e0*(dt(i,j)+dt(i,j-1)) $ *drhoy(i,j,k)*dvm(i,j) $ *(dx(i,j)+dx(i,j-1)) end do end do end do C do k=1,kb do j=2,jmm1 do i=2,imm1 drhox(i,j,k)=ramp*drhox(i,j,k) drhoy(i,j,k)=ramp*drhoy(i,j,k) end do end do end do C do k=1,kb do j=1,jm do i=1,im rho(i,j,k)=rho(i,j,k)+rmean(i,j,k) end do end do end do C return C end C subroutine bcond(idx) C ********************************************************************** C * * C * FUNCTION : Applies open boundary conditions. * C * * C * Closed boundary conditions are automatically * C * enabled through specification of the masks, dum, * C * dvm and fsm, in which case the open boundary * C * conditions, included below, will be overwritten. * C * * C * The C-Grid * C * ********** * C * * C * The diagram below and some of the text was provided * C * by D.-S. Ko. It is for the case where u and v are * C * the primary boundary conditions together with t and * C * s (co-located with e). e = el itself is rather * C * unimportant and is substituted from an adjacent * C * interior point. * C * * C * The dotted line (.......) bounds the interior * C * (non-boundary) grid points. In general, only those * C * variables in the interior are computed and * C * variables at open boundary have to be specified. * C * All interpolations are centered in space except * C * those at lateral open boundary where an upstream * C * scheme is usually used. "BU" indictates a line of * C * points where the boundary conditions should be * C * specified. * C * * C * Horizontal locations of e(el), t and s (etc.) are * C * coincident. Unused points are indicated by "NU" and * C * "*". However, for attractive output, adjacent * C * interior values may be filled in at these points. * C * * C * People not acquainted with sigma coordinates have * C * often asked what kind of boundary condition is * C * applied along closed horizontal boundaries. * C * Although the issue is not as important as it might * C * be for z-level grids, a direct answer is "half- * C * slip" which, of course, is between free slip and * C * non-slip. * C * ********** C-------------------------------------------------------------------------------* C | N O R T H * C | * C | 1 2 3 i-1 i i+1 im-2 im-1 im * C-----+-------------------------------------------------------------------------* C | NU BC BC BC BC * C | v v v v v * C | * C |BC> u* e u e u e . . u e u e u e . . u e u e u e +--v--+--v--+--v-- . +--v--+--v--+--v-- . +--v--+--v--+--v- +--v--+--v--+--v-- . +--v--+--v--+--v-- . +--v--+--v--+--v- u* e u e u e . . u e u e u e . . u e u e u e +--v*-+--v*-+--v*- . +--v*-+--v*-+--v*- . +--v*-+--v*-+--v* 0 for integer output, divisor for a * C * given by scala * C * * C ********************************************************************** C implicit none C integer im,jm real a(im,jm) real time,scala integer iskp,jskp character label*(*) real amx,scale integer i,ib,ie,j,jwr,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.lt.0.e0) scale = 1.e0 if (scala.eq.0.e0) then amx=1.e-12 do j=1,jm,jskp do i=1,im,iskp amx=max(abs(a(i,j)),amx) end do end do scale=10.e0**(int(log10(amx)+100.e0)-103) endif if(scala.gt.0.e0) scale=scala C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f9.4,' days multiply all values by ',1pe8.2) C do ib=1,im,cols*iskp C ie=ib+(cols-1)*iskp if(ie.gt.im) ie=im C if(scala.ge.0.e0) then write(6,3) (i,i=ib,ie,iskp) 3 format(/,2x,24i5,/) else write(6,4) (i,i=ib,ie,iskp) 4 format(/,12i10,/) endif C do j=1,jm,jskp jwr=jm+1-j if(scala.ge.0.e0) then write(6,5) jwr,(nint(a(i,jwr)/scale),i=ib,ie,iskp) 5 format(1x,i3,24i5) else write(6,6) jwr,(a(i,jwr),i=ib,ie,iskp) 6 format(1x,i2,12(e10.2)) endif end do C write(6,7) 7 format(//) C end do C return C end C subroutine prxyz(label,time,a,im,iskp,jm,jskp,kb,ko,nko,scala) C ********************************************************************** C * * C * FUNCTION : Writes horizontal layers of a 3-D field with * C * integers or floating point numbers. * C * * C * label ....... label for output * C * time ........ time (days) * C * a(im,jm,kb).. array to be printed * C * iskp ........ skipping interval for i * C * jskp ........ skipping interval for j * C * ko .......... 1-D array of k-indices for output * C * nko ......... number of elements in ko * C * scala ....... < 0 for floating point numbers output * C * 0 for integer output, divisor for a * C * based on magnitudes of |a| values * C * > 0 for integer output, divisor for a * C * given by scala * C * * C * (NOTE that this combines functions of old prxyz and * C * eprxyz) * C * * C ********************************************************************** C implicit none C integer im,jm,kb real a(im,jm,kb) real time,scala integer ko(*) integer iskp,jskp,nko character label*(*) real amx,scale integer i,ib,ie,j,jwr,k,iko,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.lt.0.e0) scale = 1.e0 if (scala.eq.0.e0) then amx=1.e-12 do iko=1,nko k=ko(iko) do j=1,jm,jskp do i=1,im,iskp amx=max(abs(a(i,j,k)),amx) end do end do end do scale=10.e0**(int(log10(amx)+100.e0)-103) endif if(scala.gt.0.e0) scale=scala C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f9.4,' days multiply all values by ',1pe8.2) C do iko=1,nko C k=ko(iko) C write(6,3) k 3 format(3x,/' Layer k = ',i2) C do ib=1,im,cols*iskp C ie=ib+(cols-1)*iskp if(ie.gt.im) ie=im C if(scala.ge.0.e0) then write(6,4) (i,i=ib,ie,iskp) 4 format(/,2x,24i5,/) else write(6,5) (i,i=ib,ie,iskp) 5 format(/,12i10,/) endif C do j=1,jm,jskp jwr=jm+1-j if(scala.ge.0.e0) then write(6,6) jwr,(nint(a(i,jwr,k)/scale),i=ib,ie,iskp) 6 format(1x,i3,24i5) else write(6,7) jwr,(a(i,jwr,k),i=ib,ie,iskp) 7 format(1x,i2,12(e10.2)) endif end do C write(6,8) 8 format(//) C end do C end do C return C end C subroutine prxz(label,time,a,im,iskp,jm,kb,jo,njo,scala,dt,zz) C ********************************************************************** C * * C * FUNCTION : Writes vertical section of a 3-D field, in the * C * x- or i-direction . * C * * C * label ....... label for output * C * time ........ time (days) * C * a(im,jm,kb).. array to be printed * C * iskp ........ skipping interval for i * C * jo .......... 1-D array of j-indices for output * C * njo ......... number of elements in jo * C * scala ....... < 0 for floating point numbers output * C * 0 for integer output, divisor for a * C * based on magnitudes of |a| values * C * > 0 for integer output, divisor for a * C * given by scala * C * dt(im,jm) ... total depth * C * zz(kb) ...... sigma coordinate * C * * C ********************************************************************** C implicit none C integer im,jm,kb real a(im,jm,kb),dt(im,jm),zz(kb) real time,scala integer jo(*) integer iskp,njo character label*(*) real amx,scale integer i,ib,ie,j,k,ijo,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.lt.0.e0) scale = 1.e0 if (scala.eq.0.e0) then amx=1.e-12 do k=1,kb do ijo=1,njo j=jo(ijo) do i=1,im,iskp amx=max(abs(a(i,j,k)),amx) end do end do end do scale=10.e0**(int(log10(amx)+100.e0)-103) endif if(scala.gt.0.e0) scale=scala C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f9.4,' days multiply all values by ',1pe8.2) C do ijo=1,njo C j=jo(ijo) C write(6,3) j 3 format(3x,/' Section j =',i3) C do ib=1,im,cols*iskp C ie=ib+(cols-1)*iskp if(ie.gt.im) ie=im C if(scala.ge.0.e0) then write(6,4) (i,i=ib,ie,iskp) 4 format(/,' i = ',24i5,/) else write(6,5) (i,i=ib,ie,iskp) 5 format(/,' i = ',12i10,/) endif C write(6,6) (nint(dt(i,j)),i=ib,ie,iskp) 6 format(8x,'d =',24i5.0,/,' z or zz') C do k=1,kb if(scala.ge.0.e0) then write(6,7) k,zz(k),(nint(a(i,j,k)/scale),i=ib,ie,iskp) 7 format(1x,i2,2x,f6.3,24i5) else write(6,8) k,zz(k),(a(i,j,k),i=ib,ie,iskp) 8 format(1x,i2,2x,f6.3,12(e10.2)) endif end do C write(6,9) 9 format(//) C end do C end do C return C end C subroutine pryz(label,time,a,im,jm,jskp,kb,io,nio,scala,dt,zz) C ********************************************************************** C * * C * FUNCTION : Writes vertical section of a 3-D field, in the * C * y- or j-direction. * C * * C * label ....... label for output * C * time ........ time (days) * C * a(im,jm,kb).. array to be printed * C * jskp ........ skipping interval for j * C * io .......... 1-D array of i-indices for output * C * nio ......... number of elements in io * C * scala ....... < 0 for floating point numbers output * C * 0 for integer output, divisor for a * C * based on magnitudes of |a| values * C * > 0 for integer output, divisor for a * C * given by scala * C * dt(im,jm) ... total depth * C * zz(kb) ...... sigma coordinate * C * * C ********************************************************************** C implicit none integer im,jm,kb real a(im,jm,kb),dt(im,jm),zz(kb) real time,scala integer io(*) integer jskp,nio character label*(*) real amx,scale integer i,j,jb,je,k,iio,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.lt.0.e0) scale = 1.e0 if (scala.eq.0.e0) then amx=1.e-12 do k=1,kb do j=1,jm,jskp do iio=1,nio i=io(iio) amx=max(abs(a(i,j,k)),amx) end do end do end do scale=10.e0**(int(log10(amx)+100.e0)-103) endif if(scala.gt.0.e0) scale=scala C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f9.4,' days multiply all values by ',1pe8.2) C do iio=1,nio C i=io(iio) C write(6,3) i 3 format(3x,/' Section i =',i3) C do jb=1,jm,cols*jskp C je=jb+(cols-1)*jskp if(je.gt.jm) je=jm C if(scala.ge.0.e0) then write(6,4) (j,j=jb,je,jskp) 4 format(/,' j = ',24i5,/) else write(6,5) (j,j=jb,je,jskp) 5 format(/,' j = ',12i10,/) endif C write(6,6) (nint(dt(i,j)),j=jb,je,jskp) 6 format(8x,'d =',24i5.0,/,' z or zz') C do k=1,kb if(scala.ge.0.e0) then write(6,7) k,zz(k),(nint(a(i,j,k)/scale),j=jb,je,jskp) 7 format(1x,i2,2x,f6.3,24i5) else write(6,8) k,zz(k),(a(i,j,k),j=jb,je,jskp) 8 format(1x,i2,2x,f6.3,12(e10.2)) endif end do C write(6,9) 9 format(//) C end do C end do C return C end C subroutine seamount C ********************************************************************** C * * C * FUNCTION : Sets up for seamount problem. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real delh,delx,elejmid,elwjmid,ra,vel integer i,j,k C C Set delh > 1.0 for an island or delh < 1.0 for a seamount: C delh=0.9e0 C C Grid size: C delx=8000.e0 C C Radius island or seamount: C ra=25000.e0 C C Current velocity: C vel=0.2e0 C C Set up grid dimensions, areas of free surface cells, and C Coriolis parameter: C do j=1,jm do i=1,im C C For constant grid size: C C dx(i,j)=delx C dy(i,j)=delx C C For variable grid size: C dx(i,j)=delx-delx*sin(pi*float(i)/float(im))/2.e0 dy(i,j)=delx-delx*sin(pi*float(j)/float(jm))/2.e0 C cor(i,j)=1.e-4 C end do end do C C Calculate horizontal coordinates of grid points and rotation C angle. C C NOTE that this is introduced solely for the benefit of any post- C processing software, and in order to conform with the requirements C of the NetCDF Climate and Forecast (CF) Metadata Conventions. C C There are four horizontal coordinate systems, denoted by the C subscripts u, v, e and c ("u" is a u-point, "v" is a v-point, C "e" is an elevation point and "c" is a cell corner), as shown C below. In addition, "east_*" is an easting and "north_*" is a C northing. Hence the coordinates of the "u" points are given by C (east_u,north_u). C C Also, if the centre point of the cell shown below is at C (east_e(i,j),north_e(i,j)), then (east_u(i,j),north_u(i,j)) are C the coordinates of the western of the two "u" points, C (east_v(i,j),north_v(i,j)) are the coordinates of the southern of C the two "v" points, and (east_c(i,j),north_c(i,j)) are the C coordinates of the southwestern corner point of the cell. The C southwest corner of the entire grid is at C (east_c(1,1),north_c(1,1)). C C | | C --c------v-------c-- C | | C | | C | | C | | C u e u C | | C | | C | | C | | C --c------v-------c-- C | | C C C NOTE that the following calculation of east_c and north_c only C works properly for a rectangular grid with east and north aligned C with i and j, respectively: C do j=1,jm east_c(1,j)=0.e0 do i=2,im east_c(i,j)=east_c(i-1,j)+dx(i-1,j) end do end do C do i=1,im north_c(i,1)=0.e0 do j=2,jm north_c(i,j)=north_c(i,j-1)+dy(i,j-1) end do end do C C The following works properly for any grid: C C Elevation points: C do j=1,jm-1 do i=1,im-1 east_e(i,j)=(east_c(i,j)+east_c(i+1,j) $ +east_c(i,j+1)+east_c(i+1,j+1))/4.e0 north_e(i,j)=(north_c(i,j)+north_c(i+1,j) $ +north_c(i,j+1)+north_c(i+1,j+1))/4.e0 end do end do C C Extrapolate ends: C do i=1,im-1 east_e(i,jm) $ =((east_c(i,jm)+east_c(i+1,jm))*3.e0 $ -east_c(i,jm-1)-east_c(i+1,jm-1))/4.e0 north_e(i,jm) $ =((north_c(i,jm)+north_c(i+1,jm))*3.e0 $ -north_c(i,jm-1)-north_c(i+1,jm-1))/4.e0 end do C do j=1,jm-1 east_e(im,j) $ =((east_c(im,j)+east_c(im,j+1))*3.e0 $ -east_c(im-1,j)-east_c(im-1,j+1))/4.e0 north_e(im,j) $ =((north_c(im,j)+north_c(im,j+1))*3.e0 $ -north_c(im-1,j)-north_c(im-1,j+1))/4.e0 end do C east_e(im,jm)=east_e(im-1,jm)+east_e(im,jm-1) $ -(east_e(im-2,jm)+east_e(im,jm-2))/2.e0 north_e(im,jm)=north_e(im-1,jm)+north_e(im,jm-1) $ -(north_e(im-2,jm)+north_e(im,jm-2))/2.e0 C C u-points: C do j=1,jm-1 do i=1,im east_u(i,j)=(east_c(i,j)+east_c(i,j+1))/2.e0 north_u(i,j)=(north_c(i,j)+north_c(i,j+1))/2.e0 end do end do C C Extrapolate ends: C do i=1,im east_u(i,jm)=(east_c(i,jm)*3.e0-east_c(i,jm-1))/2.e0 north_u(i,jm)=(north_c(i,jm)*3.e0-north_c(i,jm-1))/2.e0 end do C C v-points: C do j=1,jm do i=1,im-1 east_v(i,j)=(east_c(i,j)+east_c(i+1,j))/2.e0 north_v(i,j)=(north_c(i,j)+north_c(i+1,j))/2.e0 end do end do C C Extrapolate ends: C do j=1,jm east_v(im,j)=(east_c(im,j)*3.e0-east_c(im-1,j))/2.e0 north_v(im,j)=(north_c(im,j)*3.e0-north_c(im-1,j))/2.e0 end do C C rot is the angle (radians, anticlockwise) of the i-axis relative C to east, averaged to a cell centre: C C (NOTE that the following calculation of rot only works properly C for this particular rectangular grid) C do j=1,jm do i=1,im rot(i,j)=0.e0 end do end do C C Define depth: C hmax=4500.e0 !tne:!wad: deep open ocean (orig val - not wad-related) do i=1,im do j=1,jm C h(i,j)=hmax*(1.e0-delh $ *exp(-((east_c(i,j) $ -east_c((im+1)/2,j))**2 $ +(north_c(i,j) $ -north_c(i,(jm+1)/2))**2) $ /ra**2)) if(h(i,j).lt.1.e0) h(i,j)=1.e0 C end do end do C C Close the north and south boundaries to form a channel: C do i=1,im h(i,1)=1.e0 h(i,jm)=1.e0 end do C C Calculate areas and masks: C call areas_masks C C Adjust bottom topography so that cell to cell variations C in h do not exceed parameter slmax: C if(slmax.lt.1.e0) call slpmax C C Set initial conditions: C do k=1,kbm1 do j=1,jm do i=1,im tb(i,j,k)=5.e0+15.e0*exp(zz(k)*h(i,j)/1000.e0)-tbias sb(i,j,k)=35.e0-sbias tclim(i,j,k)=tb(i,j,k) sclim(i,j,k)=sb(i,j,k) ub(i,j,k)=vel*dum(i,j) end do end do end do C C Initialise uab and vab as necessary C (NOTE that these have already been initialised to zero in the C main program): C do j=1,jm do i=1,im uab(i,j)=vel*dum(i,j) end do end do C C Set surface boundary conditions, e_atmos, vflux, wusurf, C wvsurf, wtsurf, wssurf and swrad, as necessary C (NOTE: C 1. These have all been initialised to zero in the main program. C 2. The temperature and salinity of inflowing water must be C defined relative to tbias and sbias.): C do j=1,jm do i=1,im C No conditions necessary for this problem ! !lyo:!wad:!pom2k_bug:tsurf and ssurf were never defined, but should be: tsurf(i,j)=tb(i,j,1) ssurf(i,j)=sb(i,j,1) ! end do end do C C Initialise elb, etb, dt and aam2d: C do j=1,jm do i=1,im elb(i,j)=-e_atmos(i,j) etb(i,j)=-e_atmos(i,j) dt(i,j)=h(i,j)-e_atmos(i,j) aam2d(i,j)=aam(i,j,1) end do end do C ! ! !----------------------------------------------------------------------! !lyo:!wad: Set up pdens before 1st call dens; used also in profq: ! do k=1,kbm1; do j=1,jm; do i=1,im pdens(i,j,k)=grav*rhoref*(-zz(k)*max(h(i,j)-hhi,0.e0))*1.e-5 enddo; enddo; enddo ! ! call dens(sb,tb,rho) ! ! !----------------------------------------------------------------------! C C Generated horizontally averaged density field (in this C application, the initial condition for density is a function C of z (the vertical cartesian coordinate) -- when this is not C so, make sure that rmean has been area averaged BEFORE transfer C to sigma coordinates): C do k=1,kbm1 do j=1,jm do i=1,im rmean(i,j,k)=rho(i,j,k) end do end do end do C C Set lateral boundary conditions, for use in subroutine bcond C (in the seamount problem, the east and west boundaries are open, C while the south and north boundaries are closed through the C specification of the masks fsm, dum and dvm): C rfe=1.e0 rfw=1.e0 rfn=1.e0 rfs=1.e0 C do j=2,jmm1 uabw(j)=uab(2,j) uabe(j)=uab(imm1,j) C C Set geostrophically conditioned elevations at the boundaries: C ele(j)=ele(j-1)-cor(imm1,j)*uab(imm1,j)/grav*dy(imm1,j-1) elw(j)=elw(j-1)-cor(2,j)*uab(2,j)/grav*dy(2,j-1) end do C C Adjust boundary elevations so that they are zero in the middle C of the channel: C elejmid=ele(jmm1/2) elwjmid=elw(jmm1/2) do j=2,jmm1 ele(j)=(ele(j)-elejmid)*fsm(im,j) elw(j)=(elw(j)-elwjmid)*fsm(2,j) end do C C Set thermodynamic boundary conditions (for the seamount C problem, and other possible applications, lateral thermodynamic C boundary conditions are set equal to the initial conditions and C are held constant thereafter - users may, of course, create C variable boundary conditions): C do k=1,kbm1 C do j=1,jm tbe(j,k)=tb(im,j,k) tbw(j,k)=tb(1,j,k) sbe(j,k)=sb(im,j,k) sbw(j,k)=sb(1,j,k) end do C do i=1,im tbn(i,k)=tb(i,jm,k) tbs(i,k)=tb(i,1,k) sbn(i,k)=sb(i,jm,k) sbs(i,k)=sb(i,1,k) end do C end do C return C end C subroutine slpmax C ********************************************************************** C * * C * FUNCTION : Limits the maximum of: * C * * C * / * C * * C * for two adjacent cells. The maximum possible value * C * is unity. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real mean,del integer i,j,loop C do loop=1,10 C C Sweep right: C do j=2,jm-1 C do i=2,im-1 if(fsm(i,j).ne.0.e0.and.fsm(i+1,j).ne.0.e0) then if(abs(h(i+1,j)-h(i,j))/(h(i,j)+h(i+1,j)).ge.slmax) then mean=(h(i+1,j)+h(i,j))/2.e0 del=sign(slmax,h(i+1,j)-h(i,j)) h(i+1,j)=mean*(1.e0+del) h(i,j)=mean*(1.e0-del) endif endif end do C C Sweep left: C do i=im-1,2,-1 if(fsm(i,j).ne.0.e0.and.fsm(i+1,j).ne.0.e0) then if(abs(h(i+1,j)-h(i,j))/(h(i,j)+h(i+1,j)).ge.slmax) then mean=(h(i+1,j)+h(i,j))/2.e0 del=sign(slmax,h(i+1,j)-h(i,j)) h(i+1,j)=mean*(1.e0+del) h(i,j)=mean*(1.e0-del) endif endif end do C end do C C Sweep up: C do i=2,im-1 C do j=2,jm-1 if(fsm(i,j).ne.0.e0.and.fsm(i,j+1).ne.0.e0) then if(abs(h(i,j+1)-h(i,j))/(h(i,j)+h(i,j+1)).ge.slmax) then mean=(h(i,j+1)+h(i,j))/2.e0 del=sign(slmax,h(i,j+1)-h(i,j)) h(i,j+1)=mean*(1.e0+del) h(i,j)=mean*(1.e0-del) endif endif end do C C Sweep down: C do j=jm-1,2,-1 if(fsm(i,j).ne.0.e0.and.fsm(i,j+1).ne.0.e0) then if(abs(h(i,j+1)-h(i,j))/(h(i,j)+h(i,j+1)).ge.slmax) then mean=(h(i,j+1)+h(i,j))/2.e0 del=sign(slmax,h(i,j+1)-h(i,j)) h(i,j+1)=mean*(1.e0+del) h(i,j)=mean*(1.e0-del) endif endif end do C end do C end do C return C end C subroutine smol_adif(xmassflux,ymassflux,zwflux,ff,sw) C ********************************************************************** C * * C * FUNCTION : Calculates the antidiffusive velocity used to * C * reduce the numerical diffusion associated with the * C * upstream differencing scheme. * C * * C * This is based on a subroutine of Gianmaria Sannino * C * (Inter-university Computing Consortium, Rome, Italy)* C * and Vincenzo Artale (Italian National Agency for * C * New Technology and Environment, Rome, Italy), * C * downloaded from the POM FTP site on 1 Nov. 2001. * C * The calculations have been simplified by removing * C * the shock switch option. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real ff(im,jm,kb) real xmassflux(im,jm,kb),ymassflux(im,jm,kb),zwflux(im,jm,kb) real sw real mol,abs_1,abs_2 real value_min,epsilon real udx,u2dt,vdy,v2dt,wdz,w2dt integer i,j,k C parameter (value_min=1.e-9,epsilon=1.0e-14) C C Apply temperature and salinity mask: C do k=1,kb do i=1,im do j=1,jm ff(i,j,k)=ff(i,j,k)*fsm(i,j) end do end do end do C C Recalculate mass fluxes with antidiffusion velocity: C do k=1,kbm1 do j=2,jmm1 do i=2,im if(ff(i,j,k).lt.value_min.or. $ ff(i-1,j,k).lt.value_min) then xmassflux(i,j,k)=0.e0 else udx=abs(xmassflux(i,j,k)) u2dt=dti2*xmassflux(i,j,k)*xmassflux(i,j,k)*2.e0 $ /(aru(i,j)*(dt(i-1,j)+dt(i,j))) mol=(ff(i,j,k)-ff(i-1,j,k)) $ /(ff(i-1,j,k)+ff(i,j,k)+epsilon) xmassflux(i,j,k)=(udx-u2dt)*mol*sw abs_1=abs(udx) abs_2=abs(u2dt) if(abs_1.lt.abs_2) xmassflux(i,j,k)=0.e0 end if end do end do end do C do k=1,kbm1 do j=2,jm do i=2,imm1 if(ff(i,j,k).lt.value_min.or. $ ff(i,j-1,k).lt.value_min) then ymassflux(i,j,k)=0.e0 else vdy=abs(ymassflux(i,j,k)) v2dt=dti2*ymassflux(i,j,k)*ymassflux(i,j,k)*2.e0 $ /(arv(i,j)*(dt(i,j-1)+dt(i,j))) mol=(ff(i,j,k)-ff(i,j-1,k)) $ /(ff(i,j-1,k)+ff(i,j,k)+epsilon) ymassflux(i,j,k)=(vdy-v2dt)*mol*sw abs_1=abs(vdy) abs_2=abs(v2dt) if(abs_1.lt.abs_2) ymassflux(i,j,k)=0.e0 end if end do end do end do C do k=2,kbm1 do j=2,jmm1 do i=2,imm1 if(ff(i,j,k).lt.value_min.or. $ ff(i,j,k-1).lt.value_min) then zwflux(i,j,k)=0.e0 else wdz=abs(zwflux(i,j,k)) w2dt=dti2*zwflux(i,j,k)*zwflux(i,j,k)/(dzz(k-1)*dt(i,j)) mol=(ff(i,j,k-1)-ff(i,j,k)) $ /(ff(i,j,k)+ff(i,j,k-1)+epsilon) zwflux(i,j,k)=(wdz-w2dt)*mol*sw abs_1=abs(wdz) abs_2=abs(w2dt) if(abs_1.lt.abs_2)zwflux(i,j,k)=0.e0 end if end do end do end do C return C end C subroutine vertvl(xflux,yflux) C ********************************************************************** C * * C * FUNCTION : Calculates vertical velocity. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real xflux(im,jm,kb),yflux(im,jm,kb) integer i,j,k C C Reestablish boundary conditions: C do k=1,kbm1 do j=2,jm do i=2,im xflux(i,j,k)=.25e0*(dy(i,j)+dy(i-1,j)) $ *(dt(i,j)+dt(i-1,j))*u(i,j,k) end do end do end do C do k=1,kbm1 do j=2,jm do i=2,im yflux(i,j,k)=.25e0*(dx(i,j)+dx(i,j-1)) $ *(dt(i,j)+dt(i,j-1))*v(i,j,k) end do end do end do C C NOTE that, if one wishes to include freshwater flux, the C surface velocity should be set to vflux(i,j). See also C change made to 2-D volume conservation equation which C calculates elf. C do j=2,jmm1 do i=2,imm1 w(i,j,1)=0.5*(vfluxb(i,j)+vfluxf(i,j)) end do end do C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 w(i,j,k+1)=w(i,j,k) $ +dz(k)*((xflux(i+1,j,k)-xflux(i,j,k) $ +yflux(i,j+1,k)-yflux(i,j,k)) $ /(dx(i,j)*dy(i,j)) $ +(etf(i,j)-etb(i,j))/dti2) end do end do end do C return C end c !lyo:!wad: ! ! ! !----------------------------------------------------------------------! ! ! ! Here are WAD-related subroutines in alphabatical order; ! ! all names begin with "wad" ! ! ! !----------------------------------------------------------------------! ! ! subroutine wadadvt2d(fb,f,ff,nitera,sw,dx,dy,u,v,art,aru,arv, 1 fsm,dum,dvm,aam,dti2) ! ! ! 2-d version of Smolarkiewicz from /home/lyo/pom/pom2k/ ! ! pom2k.f's subroutine advt2: ! ! ! ! This version gets rid of common block, i.e. the routine is now ! ! self-contained, except that im & jm need to be specified below ! ! if "include 'grid'" is NOT used; also, 2-d calculation, i.e, ! ! solve for D: ! ! ! ! d(D)/dt + d(UD)/dx + d(VD)/dy = Diffusion_of_(D) ! ! ! C ********************************************************************** C * * C * FUNCTION : Integrates conservative scalar equations. * C * * C * This is a first-order upstream scheme, which * C * reduces implicit diffusion using the Smolarkiewicz * C * iterative upstream scheme with an antidiffusive * C * velocity. * C * * C * It is based on the subroutines of Gianmaria Sannino * C * (Inter-university Computing Consortium, Rome, Italy)* C * and Vincenzo Artale (Italian National Agency for * C * New Technology and Environment, Rome, Italy), * C * downloaded from the POM FTP site on 1 Nov. 2001. * C * The calculations have been simplified by removing * C * the shock switch option. It should be noted that * C * this implementation does not include cross-terms * C * which are in the original formulation. * C * * C * fb,f,fclim,ff . as used in subroutine advt1 * C * xflux,yflux ... working arrays used to save memory * C * nitera ........ number of iterations. This should * C * be in the range 1 - 4. 1 is * C * standard upstream differencing; * C * 3 adds 50% CPU time to POM. * C * sw ............ smoothing parameter. This should * C * preferably be 1, but 0 < sw < 1 * C * gives smoother solutions with less * C * overshoot when nitera > 1. * C * * C * Reference: * C * * C * Smolarkiewicz, P.K.; A fully multidimensional * C * positive definite advection transport algorithm * C * with small implicit diffusion, Journal of * C * Computational Physics, 54, 325-362, 1984. * C * * C ********************************************************************** C implicit none C real sw,dti2 integer nitera,im,jm,kb !lyo:!wad:kb not required but defined ! in "grid" below c ! PARAMETER (IM=65,JM=49) ! PARAMETER (IM=131,JM=99) include 'grid' c real fb(im,jm),f(im,jm),ff(im,jm) real dx(im,jm),dy(im,jm),u(im,jm),v(im,jm),aam(im,jm) real fsm(im,jm),dum(im,jm),dvm(im,jm) real art(im,jm),aru(im,jm),arv(im,jm) c integer ix,jx c PARAMETER (IX=62,JX=18) real xflux(im,jm),yflux(im,jm) real fbmem(im,jm) real xmassflux(im,jm),ymassflux(im,jm) integer i,j,k,itera,imm1,jmm1 C C Calculate horizontal mass fluxes: C imm1=im-1; jmm1=jm-1 c do j=1,jm do i=1,im xmassflux(i,j)=0.e0 ymassflux(i,j)=0.e0 end do end do C do j=2,jmm1 do i=2,im xmassflux(i,j)=0.5e0*(dy(i-1,j)+dy(i,j))*u(i,j) end do end do C do j=2,jm do i=2,imm1 ymassflux(i,j)=0.5e0*(dx(i,j-1)+dx(i,j))*v(i,j) end do end do C do j=1,jm do i=1,im fbmem(i,j)=fb(i,j) end do end do C C Start Smolarkiewicz scheme: C do itera=1,nitera C C Upwind advection scheme: C do j=2,jm do i=2,im xflux(i,j)=0.5e0 $ *((xmassflux(i,j)+abs(xmassflux(i,j))) $ *fbmem(i-1,j)+ $ (xmassflux(i,j)-abs(xmassflux(i,j))) $ *fbmem(i,j)) C yflux(i,j)=0.5e0 $ *((ymassflux(i,j)+abs(ymassflux(i,j))) $ *fbmem(i,j-1)+ $ (ymassflux(i,j)-abs(ymassflux(i,j))) $ *fbmem(i,j)) end do end do C C Add net advective fluxes and step forward in time: C do j=2,jmm1 do i=2,imm1 ff(i,j)=xflux(i+1,j)-xflux(i,j) $ +yflux(i,j+1)-yflux(i,j) ff(i,j)=(fbmem(i,j)*art(i,j) $ -dti2*ff(i,j))/(art(i,j)) end do end do C C Calculate antidiffusion velocity: C call wadsmoladif(xmassflux,ymassflux,ff,sw,fsm,aru,arv,dti2) c call wadsmoladif(xmassflux,ymassflux,ff,sw,fsm,aru,arv,dti2,im,jm) C do j=1,jm do i=1,im fbmem(i,j)=ff(i,j) end do end do C C End of Smolarkiewicz scheme C end do C C Add horizontal diffusive fluxes: C do j=2,jm do i=2,im xmassflux(i,j)=0.5e0*(aam(i,j)+aam(i-1,j)) ymassflux(i,j)=0.5e0*(aam(i,j)+aam(i,j-1)) end do end do C do j=2,jm do i=2,im xflux(i,j)=-xmassflux(i,j) $ *(fb(i,j)-fb(i-1,j))*dum(i,j) $ *(dy(i,j)+dy(i-1,j))/(dx(i,j)+dx(i-1,j)) yflux(i,j)=-ymassflux(i,j) $ *(fb(i,j)-fb(i,j-1))*dvm(i,j) $ *(dx(i,j)+dx(i,j-1))/(dy(i,j)+dy(i,j-1)) end do end do C C Add net horizontal fluxes and step forward in time: C do j=2,jmm1 do i=2,imm1 ff(i,j)=ff(i,j)-dti2*(xflux(i+1,j)-xflux(i,j) $ +yflux(i,j+1)-yflux(i,j)) $ /(art(i,j)) end do end do C return C end C ! ! !----------------------------------------------------------------------! ! ! subroutine wadout ! ! ! ********************************************************************** ! * * ! * FUNCTION : WAD Outputs (temporary subroutine to be merged with * ! * the rest of the code later) * ! * * C * Edit as approriate. * C * * ! ********************************************************************** ! ! implicit none C integer i,j,k C include 'pom08.c' C C 2-D horizontal fields: C call prxy('wetmask; =0 are land or dry, =1 are wet ', $ time,wetmask,im,iskp,jm,jskp,1.e0) C do j=1,jm; do i=1,im tps(i,j)=fsm(i,j)-wetmask(i,j) enddo; enddo call prxy('fsm-wetmask; WAD: =1 are dry cells ', $ time,tps,im,iskp,jm,jskp,1.e0) C do j=1,jm; do i=1,im tps(i,j)=(elb(i,j)+hhi)*wetmask(i,j) enddo; enddo call prxy('elb+hhi (m; w.r.t MSL) ', $ time,tps,im,iskp,jm,jskp,0.e0) c do j=1,jm; do i=1,im tps(i,j)=d(i,j)*fsm(i,j) enddo; enddo call prxy('Total Water Depth (m) ', $ time,tps,im,iskp,jm,jskp,0.e0) c ctne: -------------- save for matlab plot c write(88,'(2i5,4f8.3)')im-2,jm-2,time*24.,99.9,99.9,99.9 do j=2,jm-1; do i=2,im-1 write(88,'(2i5,4f8.3)')i,j,tps(i,j),uab(i,j),vab(i,j),d(i,j) enddo; enddo c return end ! ! !----------------------------------------------------------------------! ! ! subroutine wadh ! ! ! ********************************************************************** ! * * ! * FUNCTION : Sets up WAD definitions of topography [see Fig.1 of * ! * O2006] * ! * * ! ********************************************************************** ! ! ! Input (through common block) is "h" defined as follows: ! ! ! ! "h" is wrt MSL (Mean Sea Level) at which level h=0, and h<0 for ! ! depth below the MSL, and h>0 for "land" above the MSL (i.e. for ! ! marshes, wetlands, hills & mountains etc.) ! ! ! ! Also, nwad, hhi, slmax ! ! ! ! Outputs: ! ! ! ! h --> (i.e. becomes) h_old_pom + hhi (which can be zero) ! ! i.e. "h" is wrt to some high land-datum (Absolute Land Boundary! ! or ALB), "hhi" meters above the MSL. All cells other than the ! ! ALB cells are either always wet or can potentially become wet. ! ! ! ! Also, fsm, dum, dvm, wetmask, cell areas etc, and hmax ! ! see more detailed definitions inside the subroutine ! ! ! !----------------------------------------------------------------------! implicit none ! ! include 'pom08.c' ! ! integer npos,nneg real hwatmin integer i,j ! ! !----------------------------------------------------------------------! ! Do a rudimentary check on "h:" ! ! ! npos=0; nneg=0 do j=1,jm; do i=1,im if (h(i,j).ge.0.0) then npos=+1 else nneg=-1 endif if (npos*nneg .lt. 0) go to 101 enddo; enddo ! ! write(6,'('' Stopped in subr. wadh; incorrect defn of h:'')') write(6,'('' h is one-sign only; see comments in subr.wadh'',/)') write(6,'('' npos,nneg ='',2i5,/)') npos,nneg call prxy('Undisturbed water depth, h',time,h,im,1,jm,1,0.0) stop ! ! 101 continue ! ! hkeep(:,:)=h(:,:) !Keep original for plots etc ! ! ! HMIN & HC (all +ve) are defined as: ! ! ! ! h = hmin ---> absolute land area [never flooded], FSM=0.0; ! ! h >= hc ---> water depth wrt High water level, FSM=1.0; ! ! ! ! Note that hc is defined in MAIN -- ! ! hmin=0.01; hc=0.05; hco=hc !hco is used to avoid round- hmin=0.01; hco=hc !hco is used to avoid round- hc=hc*1.0001 !off in initial elevation ! hc=hc*1.01 !tne:!wad:- using a larger ! multiplying factor, "1.01" instead of "1.0001", ! ! can help eliminate singular wet spots due to roundoff! ! ! write(6,'('' Outputs from subr. wadh ---------------------'',/)') write(6,'('' hmin,hc,hco ='',1p3e13.4,/)') hmin,hc,hco ! ! ! Reverse sign of "h", and shift reference if hhi.ne.0, see below..! ! (Note that hhi is defined in MAIN) ! h(:,:)=-h(:,:)+hhi ! ! if (nwad.eq.1) then ! ! ! If nwad=1, then hhi.ne.0.0, and line "h(:,:)=-h(:,:)+hhi" above ! ! already shifts the reference level from MSL to ALB, the Absolute ! ! Land Boundary, according to O2006, Fig.1: ! ! ! do j=1,jm; do i=1,im FSM(I,J)=1.; DUM(I,J)=1.; DVM(I,J)=1. IF(h(I,J).LT.0.0)THEN !Define absolute land areas (never flood): ! temporarily set h=-1.0 (some -ve number) h(I,J)=-1.0; FSM(I,J)=0.; DUM(I,J)=0.; DVM(I,J)=0. ENDIF enddo; enddo DO J=1,JM-1; DO I=1,IM IF(FSM(I,J).EQ.0..AND.FSM(I,J+1).NE.0.)DVM(I,J+1)=0. enddo; enddo DO J=1,JM; DO I=1,IM-1 IF(FSM(I,J).EQ.0..AND.FSM(I+1,J).NE.0.)DUM(I+1,J)=0. enddo; enddo ! ! ! The above yields the following definitions for h: ! ! ! ! h=-1, absolute land areas [never flooded], FSM=0.0, always ! ! h>=0, wet but potentially dry areas, FSM also=1.0, always ! ! ! ! Note that slpmax touches fsm=1 points only - wet or potentially ! ! WAD cells; it changes "h" so a wet cell might become dry; so it ! ! must be called here before wetmask is defined: ! ! ! ! Adjust bottom topography so that cell to cell variations ! ! in h do not exceed parameter slmax: ! ! ! if(slmax.lt.1.e0) call slpmax ! ! ! We now define a wetmask, assuming that at t=0, the free surface ! ! is at the MSL (i.e. elevation=0 wrt MSL). If this run does not ! ! begin with time=0 (or if other special free-surface distribution ! ! is desired), then wetmask needs to be redefined or read in e.g. ! ! from a RESTART file, as done later in MAIN if NREAD.ne.0 ! ! ! wetmask(:,:)=fsm(:,:) do j=1,jm; do i=1,im if (h(i,j).lt.hhi) wetmask(i,j)=0.0 enddo; enddo ! ! ! Finalize definitions of "h": ! ! ! ! h=hmin, absolute land areas [never flooded], FSM=0.0, always ! ! h>=hc, wet but potentially dry areas, FSM also=1.0, always ! ! ... but wetmask can be 0 or 1 ! ! ! do j=1,jm; do i=1,im if (h(i,j).lt.0.0) then h(i,j)=hmin else h(i,j)=h(i,j)+hc endif enddo; enddo ! ! call areas_masks !Calc. cells' areas etc but do not alter fsm... !if nwad=1 ! ! else !nwad.eq.0 ! ! DO J=1,JM; DO I=1,IM IF (h(I,J).LE.0.0) h(I,J)=hmin ENDDO; ENDDO ! ! ! Set minimum water depth to hwatmin: ! ! hwatmin=min water depth, must be >1 to match subr.areas_masks' ! ! definition of land, and large enough to avoid grid cells from ! ! becoming dry ! ! ! hwatmin=10.0 DO J=1,JM; DO I=1,IM IF (h(I,J).GT.hmin.and.h(I,J).LT.hwatmin) h(I,J)=hwatmin ENDDO; ENDDO ! ! call areas_masks !Calc. cells' areas etc & fsm if nwad=0 ! ! ! Adjust bottom topography so that cell to cell variations ! ! in h do not exceed parameter slmax: ! ! ! if(slmax.lt.1.e0) call slpmax ! ! wetmask(:,:)=fsm(:,:) !wetmask is same as fsm if nwad=0 ! ! endif !if (nwad.eq.1) then...else... ! ! ! Print to check: ! ! ! CALL PRXY(' Input h ',TIME, hkeep(1,1),IM,1,JM,1,0.0) CALL PRXY(' h after wadh ',TIME, h(1,1),IM,1,JM,1,0.0) CALL PRXY(' FSM ',TIME, FSM(1,1),IM,1,JM,1,1.0) CALL PRXY(' WETMASK ',TIME, WETMASK(1,1),IM,1,JM,1,1.0) tps(:,:)=FSM(:,:)-WETMASK(:,:) CALL PRXY('FSM-WETMASK',TIME,tps(1,1),IM,1,JM,1,1.0) ! ! ! Double-check masks for nwad=0: ! ! ! if (nwad.eq.0) then do j=1,jm; do i=1,im if (fsm(i,j).ne.wetmask(i,j)) then write(6,'('' Stopped, nwad ='',i4)') nwad write(6,'('' fsm.ne.wetmask @ (i,j) ='',2i4,/)') i,j stop endif enddo; enddo endif ! ! return end ! ! !----------------------------------------------------------------------! ! ! subroutine wadseamount C ********************************************************************** C * * C * FUNCTION : Sets up for seamount problem w/wad. * C * * C ********************************************************************** C implicit none C include 'pom08.c' C real delh,delx,elejmid,elwjmid,ra,vel,hland !lyo:!wad:define hland integer nvar !lyo:!wad:define nvar integer i,j,k C C Set delh > 1.0 for an island or delh < 1.0 for a seamount: C delh=1.15e0 !tne:!wad: island (for wad) =0.9e0 for orig.seamount C C Grid size: C delx=8000.e0 c !lyo:!wad:Define variable or uniform grid option: nvar=1 !=1 for variable grid; =0 otherwise !lyo:!wad:Special test case for smaller delx=4km: if ((nvar.eq.0).or.(im.eq.131.and.jm.eq.99)) delx=delx*0.5 C C Radius island or seamount: C ra=50000.e0 !tne:!wad: C C Current velocity: C vel=0.2e0 C C Set up grid dimensions, areas of free surface cells, and C Coriolis parameter: C do j=1,jm do i=1,im C C For constant grid size: C C dx(i,j)=delx C dy(i,j)=delx C C For variable grid size: C !lyo:!wad:variable or uniform grid-- dx(i,j)=delx-float(nvar)*delx*sin(pi*float(i)/float(im))/2.e0 dy(i,j)=delx-float(nvar)*delx*sin(pi*float(j)/float(jm))/2.e0 C cor(i,j)=1.e-4 C end do end do C C Calculate horizontal coordinates of grid points and rotation C angle. C C NOTE that this is introduced solely for the benefit of any post- C processing software, and in order to conform with the requirements C of the NetCDF Climate and Forecast (CF) Metadata Conventions. C C There are four horizontal coordinate systems, denoted by the C subscripts u, v, e and c ("u" is a u-point, "v" is a v-point, C "e" is an elevation point and "c" is a cell corner), as shown C below. In addition, "east_*" is an easting and "north_*" is a C northing. Hence the coordinates of the "u" points are given by C (east_u,north_u). C C Also, if the centre point of the cell shown below is at C (east_e(i,j),north_e(i,j)), then (east_u(i,j),north_u(i,j)) are C the coordinates of the western of the two "u" points, C (east_v(i,j),north_v(i,j)) are the coordinates of the southern of C the two "v" points, and (east_c(i,j),north_c(i,j)) are the C coordinates of the southwestern corner point of the cell. The C southwest corner of the entire grid is at C (east_c(1,1),north_c(1,1)). C C | | C --c------v-------c-- C | | C | | C | | C | | C u e u C | | C | | C | | C | | C --c------v-------c-- C | | C C C NOTE that the following calculation of east_c and north_c only C works properly for a rectangular grid with east and north aligned C with i and j, respectively: C do j=1,jm east_c(1,j)=0.e0 do i=2,im east_c(i,j)=east_c(i-1,j)+dx(i-1,j) end do end do C do i=1,im north_c(i,1)=0.e0 do j=2,jm north_c(i,j)=north_c(i,j-1)+dy(i,j-1) end do end do C C The following works properly for any grid: C C Elevation points: C do j=1,jm-1 do i=1,im-1 east_e(i,j)=(east_c(i,j)+east_c(i+1,j) $ +east_c(i,j+1)+east_c(i+1,j+1))/4.e0 north_e(i,j)=(north_c(i,j)+north_c(i+1,j) $ +north_c(i,j+1)+north_c(i+1,j+1))/4.e0 end do end do C C Extrapolate ends: C do i=1,im-1 east_e(i,jm) $ =((east_c(i,jm)+east_c(i+1,jm))*3.e0 $ -east_c(i,jm-1)-east_c(i+1,jm-1))/4.e0 north_e(i,jm) $ =((north_c(i,jm)+north_c(i+1,jm))*3.e0 $ -north_c(i,jm-1)-north_c(i+1,jm-1))/4.e0 end do C do j=1,jm-1 east_e(im,j) $ =((east_c(im,j)+east_c(im,j+1))*3.e0 $ -east_c(im-1,j)-east_c(im-1,j+1))/4.e0 north_e(im,j) $ =((north_c(im,j)+north_c(im,j+1))*3.e0 $ -north_c(im-1,j)-north_c(im-1,j+1))/4.e0 end do C east_e(im,jm)=east_e(im-1,jm)+east_e(im,jm-1) $ -(east_e(im-2,jm)+east_e(im,jm-2))/2.e0 north_e(im,jm)=north_e(im-1,jm)+north_e(im,jm-1) $ -(north_e(im-2,jm)+north_e(im,jm-2))/2.e0 C C u-points: C do j=1,jm-1 do i=1,im east_u(i,j)=(east_c(i,j)+east_c(i,j+1))/2.e0 north_u(i,j)=(north_c(i,j)+north_c(i,j+1))/2.e0 end do end do C C Extrapolate ends: C do i=1,im east_u(i,jm)=(east_c(i,jm)*3.e0-east_c(i,jm-1))/2.e0 north_u(i,jm)=(north_c(i,jm)*3.e0-north_c(i,jm-1))/2.e0 end do C C v-points: C do j=1,jm do i=1,im-1 east_v(i,j)=(east_c(i,j)+east_c(i+1,j))/2.e0 north_v(i,j)=(north_c(i,j)+north_c(i+1,j))/2.e0 end do end do C C Extrapolate ends: C do j=1,jm east_v(im,j)=(east_c(im,j)*3.e0-east_c(im-1,j))/2.e0 north_v(im,j)=(north_c(im,j)*3.e0-north_c(im-1,j))/2.e0 end do C C rot is the angle (radians, anticlockwise) of the i-axis relative C to east, averaged to a cell centre: C C (NOTE that the following calculation of rot only works properly C for this particular rectangular grid) C do j=1,jm do i=1,im rot(i,j)=0.e0 end do end do C C Define depth: C !tne:!wad:!lyo:!wad: ! Note that unlike original seamount, h<0 is now water below MSL, and ! h>0 is above the MSL and is either land (if h>hland, see below) ! or is potential wet-and-dry region ! hmax=50.e0; hland=5.e0 !note all pnts are potential WAD if we set !hland >= -hmax*(1.e0-delh) (=7.5m); !i.e. "if" below is NOT satisfied ! do i=1,im do j=1,jm C h(i,j)=-hmax*(1.e0-delh $ *exp(-((east_c(i,j) $ -east_c((im+1)/2,j))**2 $ +(north_c(i,j) $ -north_c(i,(jm+1)/2))**2) $ /ra**2)) !lyo:!wad:Make region near top of seamount absolute land region: if(h(i,j).gt.hland) h(i,j)=hhi+1.e0 C end do end do C C Close the north and south boundaries to form a channel: C do i=1,im h(i,1)=hhi+1.e0 !tne:!wad: h(i,jm)=hhi+1.e0 !tne:!wad: end do C C Calculate areas and masks: C !lyo:!wad:Define "h" consistent with WAD and calculate wetmask, cell areas ! and fsm etc. ! call areas_masks call wadh C C Adjust bottom topography so that cell to cell variations C in h do not exceed parameter slmax: C ! if(slmax.lt.1.e0) call slpmax !lyo:wad:now done in wadh above C C Set initial conditions: C do k=1,kbm1 do j=1,jm do i=1,im !lyo:wad: ! tb(i,j,k)=5.e0+15.e0*exp(zz(k)*h(i,j)/1000.e0)-tbias if (h(i,j).gt.hhi) then !lyo:wad:stratified water below MSL: tb(i,j,k)=5.e0+15.e0*exp(zz(k)*(h(i,j)-hhi)/1000.e0) else !lyo:wad:well-mixed water above MSL: tb(i,j,k)=5.e0+15.e0 endif tb(i,j,k)=tb(i,j,k)-tbias ! sb(i,j,k)=35.e0-sbias tclim(i,j,k)=tb(i,j,k) sclim(i,j,k)=sb(i,j,k) ub(i,j,k)=vel*dum(i,j) end do end do end do C C Initialise uab and vab as necessary C (NOTE that these have already been initialised to zero in the C main program): C do j=1,jm do i=1,im uab(i,j)=vel*dum(i,j) end do end do C C Set surface boundary conditions, e_atmos, vflux, wusurf, C wvsurf, wtsurf, wssurf and swrad, as necessary C (NOTE: C 1. These have all been initialised to zero in the main program. C 2. The temperature and salinity of inflowing water must be C defined relative to tbias and sbias.): C do j=1,jm do i=1,im C No conditions necessary for this problem ! !lyo:!wad:!pom2k_bug:tsurf and ssurf were never defined, but should be: tsurf(i,j)=tb(i,j,1) ssurf(i,j)=sb(i,j,1) ! end do end do C C Initialise elb, etb, dt and aam2d: C do j=1,jm do i=1,im elb(i,j)=(-e_atmos(i,j)-hhi)*fsm(i,j) !lyo:!wad: !lyo:!wad:note:The following "if" is not satisfied if nwad=0 since ! ! then fsm=wetmask at all cells ! ! if (fsm(i,j).ne.0.0.and.wetmask(i,j).eq.0.0) ! dry but not land if (fsm(i,j).ne.wetmask(i,j)) ! dry but not land &elb(i,j)=-h(i,j)+hco !note slightly smaller hco