program ozpom C C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ORIGINAL : REF:JRH:2002-03-22 * C * * C * SOURCE : ozpom.f * 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 * * 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 a major * 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) Removal of all "spaghetti" code and "goto"s. * C * (4) Rationalisation of the input of all constants. * C * (5) Modifications to the "printer" output. * C * (6) Output to a netCDF file. * C * (7) Inclusion of surface freshwater flux. * C * (8) Inclusion of atmospheric pressure. * C * (9) Inclusion of an additional problem to check (7) * C * and (8), above. * C * (10) Inclusion of option for Smolarkiewicz * C * advection scheme. * C * (11) Conversion to "ozpom", creating a fork from * C * Princeton's pom2k (22/3/2002) * 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. * 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 ozpom ozpom.f /usr/lib/libnetcdf.a * C * * C * You should also have the "include" file of * C * netCDF subroutines (ozpom.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 * NOTE ON TREATMENT OF TEMPERATURE AND SALINITY * C * ********************************************* * C * * C * Temperature and salinity are treated in OzPOM as * C * concentrations PER UNIT VOLUME. For example, the * C * variable called "salinity" would appear to be: * C * * C * / * C * * C * which is not the same as the true salinity, which * C * is: * C * * C * x 1000 / * C * * C * The two are related by the density of seawater * C * (rho_sw) divided by 1000, which is in the * C * approximate range 1 to 1.03, and which may vary * C * significantly throughout the model. In OzPOM, "t" * C * and "s" should be regarded as approximations to the * C * true potential temperature and salinity (which is * C * how they are defined in the subroutine dens), in * C * which case: * C * * C * Heat content / unit volume (at the surface) * C * = t x Cp x rho_sw (J / cu. m.) * C * * C * and: * C * * C * Heat flux due to a vertical current, w, * C * (at the surface) * C * = t x Cp x rho_sw x w (J / sq. m. / s.) * C * and: * C * * C * Heat flux derived from OzPOM variable, wtsurf * C * (at the surface) * C * = Cp x rho_sw x wtsurf (J / sq. m. / s.) * C * * C * and: * C * * C * Salt content / unit volume * C * = s x rho_sw / 1000 (kg / cu. m.) * C * * C * and: * C * * C * Salt flux due to a vertical current, w, * C * = s x rho_sw x w / 1000 (kg / sq. m. / s.) * C * * C * and: * C * * C * Salt flux derived from OzPOM variable, wssurf * C * (at the surface) * C * = rho_sw x wssurf / 1000 (kg / sq. m. / s.) * C * * C * In other words, the above are not as you might * C * expect for variables defined PER UNIT VOLUME, for * C * which there would normally be no dependence on * C * rho_sw. This is the result of a Boussinesq-type * C * approximation in the conservation equations for heat* C * and salt. * 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 * * C * REVISION HISTORY * C * * C * Versions are identified by a number and a date. All * C * source files (ozpom.f, ozpom.c and ozpom.n) must * C * have identical version numbers and dates in order * C * to run. You should not change the version number or * C * the date yourself. * C * * C * Version numbers are now in the format X.Y.Z, where * C * an increment of X indicates a major change in the * C * code (no such change is envisaged at present). An * C * increment in Y indicates a change to the Fortran * C * (e.g. a bug fix), while an increment in Z indicates * C * a change to a comment only. The number X.Y.Z * C * increments like a normal decimal number (e.g. 1.2.0 * C * could follow 1.1.3). * C * * C * DATE VERSION CHANGE * C * (YYYY-MM-DD) * C * * C * 2002-03-22 1 Original version * C * * C * 2002-03-26 1.0.1 Correction to comments in subroutines * C * box and seamount regarding relationship * C * between dx, dy and art, and addition of * C * clarifying comments. Minor modifications * C * to other comments. * C * * C * 2002-03-27 1.0.2 Minor modification to comments in main * C * program, and subroutines box and * C * seamount. * C * * C * 2002-04-03 1.0.3 Minor changes to comments regarding * C * setting of Smagorinsky lateral viscosity.* C * Addition of comment to header of * C * subroutine write_netcdf in ozpom.n. * C * * C * 2002-04-05 1.1.0 (a) Correction of error in calculation * C * of p in subroutine dens. This error * C * resulted in a miscalculation of pressure * C * when calculating density, equivalent to * C * the deviation of atmospheric pressure * C * from a mean value. * C * * C * (b) Definition of new 2-D variable * C * e_atmos0, which is the long-term * C * average of e_atmos. This is introduced * C * so as to better address the problem * C * described by Robertson et al. (2001), * C * and to allow a better specification of * C * the depth when initially calculating * C * the bottom friction coefficient, and the * C * external CFL timestep. * C * * C * (c) Addition of an explanatory comment * C * regarding the meaning of the variables * C * e_atmos and e_atmos0. * C * * C * 2002-04-10 1.1.1 (a) Changes to header of "printer" * C * output of temperature and salinity to * C * indicate application of tbias and sbias. * C * * C * (b) Changes to netCDF output of * C * temperature and salinity to indicate * C * application of tbias and sbias. * C * * C * (c) Removal of comments etc. relating to * C * definition of variable s. * C * * C * 2002-04-10 1.2.0 Inclusion of vflux in I/O for restart, * C * so that first output of vflux after a * C * restart is correct. * C * * C * 2002-05-07 1.2.1 Minor correction to comment defining * C * isplit. * C * * C * 2002-05-17 1.3.0 Initialisation of aam2d(i,j) removed * C * from subroutines "box" and "seamount" and* C * placed in main program. It was * C * previously set to aam(i,j,1) in "box" * C * and "seamount", prior to aam(i,j,k) being* C * initialised. This change makes no * C * difference to the results of the * C * "box" and "seamount" test. However, it * C * does affect the initialisation of runs * C * with mode=2, in which case aam2d was * C * used in advave in the external mode * C * calculation, prior to being initialised. * C * * C * 2002-06-04 1.4.0 Inclusion of "printer" output of values * C * of netcdf_file and aam_init in main * C * program. Modification to format of * C * "printer" output of time in main program * C * and in subroutines prxy, prxyz, prxz and * C * pryz. * C * * C * 2002-06-14 1.5.0 Moved statement: * C * art(i,j)=dx(i,j)*dy(i,j) * C * from subroutines "box" and "seamount" to * C * subroutine "areas_masks". * C * * C * 2002-07-12 1.6.0 In subroutine "profq", modified * C * calculation of boygr so that the sound * C * velocity term is based on the "reference"* C * surface elevation, rather than on the * C * actual surface elevation. This is to * C * ensure consistency with the changes to * C * pom98 made on 2001-9-12. Under certain * C * conditions, versions of pom98/pom2k/ozpom* C * between 2001-9-12 and this modification * C * could exhibit enhanced vertical mixing. * C * Also changed two occurrences of "d0" to * C * "e0". * C * * C * 2002-07-15 1.6.1 Minor cosmetic changes. * C * * C * 2002-09-09 1.7.0 In subroutines "advq", "advt1" and * C * "advt2", changed "h" to "dt" in * C * horizontal diffusion terms, for * C * consistency with subroutine "advct". * C * Removed assymetries in many places by the* C * inclusion of brackets, correcting some * C * loop indices and other re-arrangement * C * (such as modifications to the calculation* C * of aam). Added symmetry-testing problem, * C * symm (iproblem = 3 or 4). Problem * C * my_problem is now iproblem = 5. Removed * C * unused parameters in ozpom.c. Added both * C * lineprinter and netCDF output of surface * C * variables e_atmos, vflux, wssurf, wtsurf,* C * wusurf and wvsurf. * C * * C * 2002-11-25 1.8.0 (a) Generated new versions of dens and * C * profq, based on the new density algorithm* C * of McDougall, Jackett, Wright and Feistel* C * (submitted to J. Atmospheric and Oceanic * C * Technology, 2002). Old versions were * C * renamed dens_old and profq_old. The old * C * or new routines may be selected by * C * setting the new variable, idens, to 1 or * C * 2 respectively in the main program (the * C * facility for using the old density * C * algorithms will be removed in later * C * versions of OzPOM). New arrays, zd and * C * dep0 were declared in common to * C * accommodate these new routines. * C * * C * NOTE that, if using the new density * C * routines, the initialisation routine * C * (e.g. my_problem) must contain the * C * statement (in an appropriate loop): * C * * C * dep0(i,j)=h(i,j)-e_atmos0(i,j) * C * * C * after initialisation of h and e_atmos0. * C * * C * (b) Included check that kl1 and kl2 are * C * in the allowable ranges: * C * * C * 2 <= kl1 <= 13 * C * kb-12 <= kl2 <= kb-1 * C * * C * 2003-09-09 1.9.0 (a) Minor shortening of some declaration * C * lines in order to allow conversion of * C * code from "real" to "double precision". * C * * C * (b) Removed redundant declarations of p, * C * sp and tp from subroutine profq. * C * * C * (c) Minor corrections to comments. * C * * C * (d) Removed tps from argument of advave * C * and declared variable curv2d internally. * C * * C * (e) Removed tps from argument of proft, * C * profq and profq_old and declared variable* C * dh internally. * C * * C * (f) Array tps renamed tmp2d to indicate * C * that it is temporary 2-D array. * C * * C * (g) Removed equivalence from profq and * C * profq_old. * C * * C * (h) Declared all 3-D arrays in common. * C * These are: boygr, curv, fbmem, gh, kn, * C * prod, rad, rhtmp, sm, sh, stf, xflux, * C * xmassflux (renamed xmflx), yflux, * C * ymassflux (renamed ymflx),zwflux (renamed* C * zwflx). This increases number of 3-D * C * arrays by 7, but makes it easier to track* C * where memory is used. Total 3-D arrays * C * are now 50. * C * * C * (i) All calculations of dt(i,j) masked by* C * fsm(i,j) in order to ensure that, * C * outside model, dt is always h(i,j) prior * C * to output. * C * * C * (j) Correction to output of wusurf and * C * wvsurf labels in subroutine printall. * C * * C * (k) Defined vfluxb and vfluxf instead of * C * vflux, and modified main program and * C * subroutines box, printall, seamount, symm* C * and vertvl, so that vertical velocity is * C * calculated correctly for time-varying * C * vfluxb and vfluxf. * C * * C * (l) Inclusion of all surface fluxes in * C * restart file, so that initial output * C * after a restart is correct. * C * * C * 2003-09-11 1.10.0 (a) Check for dryness in calculation of * C * bottom friction coefficient, for safety. * C * * C * (b) Check for dryness in calculation of * C * CFL timestep, for safety. * C * * C * (c) Inserted warnings that, over land, * C * e_atmos0 and e_atmos should be less than * C * h. * C * * C * 2003-09-16 1.10.1 Inserted warning comment in header * C * concerning treatment of temperature and * C * salinity, and associated references. * C * * C * 2003-09-18 1.10.2 In subroutine printall, changed order of * C * declarations for compatibility and * C * removed unwanted indentation. * C * * C * 2003-12-02 1.11.0 In subroutine baropg, added surface * C * elevation terms to baroclinic pressure * C * gradient (these terms are only important * C * when the surface slope is high, as * C * in the case of an ice shelf cavity). * C * * C ********************************************************************** C implicit none C include 'ozpom.c' 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 real z0b real tatm,satm real smag_1,smag_2,smag_3,smag_4,smag_a,smag_b C C Declare large enough for number of sections or layers output C to device 6: C integer io(100),jo(100),ko(100) C integer i,iend,iext,imax,ispadv,isplit,iswtch integer j,jmax integer k integer nadv,nbct,nbcs,nitera,nread integer iproblem C logical lramp C character*120 netcdf_file C C*********************************************************************** C C This shows program name, version number and date - DO NOT CHANGE: C source='ozpom version 1.11.0 2003-12-02 ' C C Check versions of program and include files are the same: C if(source.ne.source_c) then C write(6,7) 7 format(/'Incompatible versions of program and include files ', $ '..... program terminated'/) stop C endif C C*********************************************************************** C C Small value: C small=1.e-9 C C Calculate pi: C pi=atan(1.e0)*4.e0 C C*********************************************************************** C C Input of filenames and constants: C C NOTE that the array sizes im, jm and kb should be set in ozpom.c C C----------------------------------------------------------------------- C C Title of this run: C title='Run 1 ' C C----------------------------------------------------------------------- C C Filename of netCDF output file: C netcdf_file='ozpom.nc' C C----------------------------------------------------------------------- C C Density algorithm (1: old; 2: new): C (included for version 1.8.0 and 1.9.* only) C idens=2 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 symmetry box symm C (x-forcing) C C 4 symmetry box symm C (y-forcing) C C 5 your own my_problem C iproblem=1 C C----------------------------------------------------------------------- C C Model type: 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=2 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=0.5e0 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.): C dte=6.e0 C C----------------------------------------------------------------------- C C / C (dti/dte; dimensionless): C isplit=30 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='0000-01-01 00:00:00 +00:00' C C----------------------------------------------------------------------- C C Duration of run (days): C days=0.025e0 C C----------------------------------------------------------------------- C C Initial print interval (days): C prtd1=0.0125e0 C C----------------------------------------------------------------------- C C Final print interval (days): C prtd2=1.e0 C C----------------------------------------------------------------------- C C Time after which print interval changes from prtd1 to C prtd2 (days): C swtch=1000.e0 C C----------------------------------------------------------------------- C C Skipping interval in i for output: C iskp=2 C C----------------------------------------------------------------------- C C Skipping interval in j for output: C jskp=2 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 C Temperature bias (deg. C): C tbias=0.e0 C C----------------------------------------------------------------------- C C Salinity bias: C sbias=0.e0 C C----------------------------------------------------------------------- C C Acceleration due to gravity (S.I. units): C grav=9.806e0 C C----------------------------------------------------------------------- C C von Karman's constant (dimensionless): C kappa=0.4e0 C C----------------------------------------------------------------------- C C Bottom roughness (metres): C z0b=0.01e0 C C----------------------------------------------------------------------- C C Minimum allowable bottom friction coefficient (dimensionless): C cbcmin=0.0025e0 C C----------------------------------------------------------------------- C C Maximum allowable bottom friction coefficient (dimensionless): C cbcmax=1.e0 C C----------------------------------------------------------------------- C C Coefficient of Smagorinsky diffusivity (dimensionless): C horcon=0.2e0 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=0.0e0 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 layers of reduced thickness at the C surface and bottom (used by subroutine depth). The number of C layers of reduced thickness 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 NOTE that 2 <= kl1 <= 13 AND kb-12 <= 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 alpha = 0.e0 is perfectly acceptable, but the C value, alpha=.225e0 permits a longer time step): C alpha=0.225e0 C C----------------------------------------------------------------------- C C Initial value of aam: C aam_init=500.e0 C C End of input of constants C*********************************************************************** C C Calculate some constants: C dti=dte*float(isplit) dte2=dte*2 dti2=dti*2 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,1) source, $ title, $ netcdf_file, $ iproblem, $ mode, $ nadv, $ nitera, $ sw, $ nread, $ dte, $ dti, $ isplit, $ time_start, $ days, $ iend, $ prtd1, $ iprint, $ prtd2, $ swtch, $ iswtch, $ iskp, $ jskp, $ lramp, $ rhoref, $ tbias, $ sbias, $ grav, $ kappa, $ z0b, $ cbcmin, $ cbcmax, $ horcon, $ tprni, $ umol, $ hmax, $ vmaxl, $ slmax, $ kl1,kl2, $ ntp, $ nbct, $ nbcs, $ ispadv, $ smoth, $ alpha, $ aam_init C 1 format(/' source = ',a40/ $ ' title = ',a40// $ ' netcdf_file= ',a120/ $ ' iproblem = ',i10/ $ ' mode = ',i10/ $ ' nadv = ',i10/ $ ' nitera = ',i10/ $ ' sw = ',f10.4/ $ ' nread = ',i10/ $ ' dte = ',f10.2/ $ ' dti = ',f10.1/ $ ' isplit = ',i10/ $ ' time_start = ',a26/ $ ' days = ',f10.4/ $ ' iend = ',i10/ $ ' prtd1 = ',f10.4/ $ ' iprint = ',i10/ $ ' prtd2 = ',f10.4/ $ ' swtch = ',f10.2/ $ ' iswtch = ',i10/ $ ' iskp = ',i10/ $ ' jskp = ',i10/ $ ' lramp = ',l10/ $ ' rhoref = ',f10.3/ $ ' tbias = ',f10.3/ $ ' sbias = ',f10.3/ $ ' grav = ',f10.4/ $ ' kappa = ',f10.4/ $ ' z0b = ',f10.6/ $ ' cbcmin = ',f10.6/ $ ' cbcmax = ',f10.6/ $ ' horcon = ',f10.3/ $ ' tprni = ',f10.4/ $ ' umol = ',f10.6/ $ ' hmax = ',f10.2/ $ ' vmaxl = ',f10.4/ $ ' slmax = ',f10.4/ $ ' kl1,kl2 = ',i5,',',i4/ $ ' ntp = ',i10/ $ ' nbct = ',i10/ $ ' nbcs = ',i10/ $ ' ispadv = ',i10/ $ ' smoth = ',f10.4/ $ ' alpha = ',f10.4/ $ ' aam_init = ',f10.2) C C----------------------------------------------------------------------- C C Initialise boundary arrays: C do i=1,im vabn(i)=0.e0 vabs(i)=0.e0 eln(i)=0.e0 els(i)=0.e0 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 elw(j)=0.e0 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 etb(i,j)=0.e0 e_atmos0(i,j)=0.e0 e_atmos(i,j)=0.e0 vfluxb(i,j)=0.e0 vfluxf(i,j)=0.e0 tflux(i,j)=0.e0 sflux(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(kl1.lt.2 .or.kl1.gt.13.or. $ kl2.lt.kb-12.or.kl2.gt.kb-1) then write(6,10) 10 format(/' Invalid value of kl1 and/or kl2', $ ' ..... program terminated'/) stop endif C 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 symm(3) else if(iproblem.eq.4) then call symm(4) else if(iproblem.eq.5) then call my_problem else write(6,8) 8 format(/' Invalid value of iproblem ..... program terminated'/) stop endif C C Inertial period for temporal filter: C period=(2.e0*pi)/(abs(cor(im/2,jm/2))*86400.e0) 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) et(i,j)=etb(i,j) etf(i,j)=et(i,j) d(i,j)=h(i,j)+el(i,j) dt(i,j)=h(i,j)+et(i,j)*fsm(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.1e0*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 C do j=1,jm do i=1,im aam2d(i,j)=aam_init end do end do 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 if(idens.eq.1) then C C Included for version 1.8.0 and 1.9.* only: C call dens_old(s,t,rho) else call dens(s,t,dep0,e_atmos0,zz,rho) endif 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 C Calculate bottom friction coefficient: C do j=1,jm do i=1,im C C Check for dryness: C if((1.e0+zz(kbm1))*(h(i,j)-e_atmos0(i,j))/z0b.ge.1.e0) then cbc(i,j)=(kappa/log((1.e0+zz(kbm1)) $ *(h(i,j)-e_atmos0(i,j))/z0b))**2 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)) else cbc(i,j)=cbcmax endif C end do end do C C Calculate external (2-D) CFL time step: C do j=1,jm do i=1,im C C Check for dryness: C if(h(i,j)-e_atmos0(i,j).ge.0) then tmp2d(i,j)=fsm(i,j) $ /(2.e0*sqrt((1.e0/dx(i,j)**2+1.e0/dy(i,j)**2) $ *grav*(h(i,j)-e_atmos0(i,j)+small))) else tmp2d(i,j)=0.e0 endif C 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, $ vfluxf,wssurf,wtsurf,wusurf,wvsurf 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)*fsm(i,j) end do end do C 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('Free surface mask, fsm ', $ time,fsm,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 call prxy('External (2-D) CFL time step, tmp2d ', $ time,tmp2d,im,iskp,jm,jskp,1.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 call write_netcdf(netcdf_file,1) ! *netCDF* call write_netcdf(netcdf_file,2) ! *netCDF* C C----------------------------------------------------------------------- C C Begin numerical integration: C do iint=1,iend 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. vflux (i.e. vfluxb or vfluxf) is the C UPWARD flux of water through the sea surface (i.e. in most cases, C evaporation - precipitation). C C NOTE: read "NOTE ON TREATMENT OF TEMPERATURE AND SALINITY" in C program header. 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 wusurf(i,j)=ramp*(1.e-4*cos(pi*(j-1)/jmm1)) C $ *.25e0*(dvm(i,j+1)+dvm(i-1,j+1) C $ +dvm(i-1,j)+dvm(i,j)) C wusurf(i,j)=0.e0 C wvsurf(i,j)=0.e0 vfluxb(i,j)=vfluxf(i,j) C vfluxf(i,j)=0.e0 C C Initially set wtsurf to the sensible heat, the latent heat (which C inolves 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 For completeness, set wssurf to any source of salt NOT associated C with vflux: C wssurf(i,j)=0.e0 C C Include the following, or something similar, if vflux(i,j,) is C non-zero (NOTE that, if tflux and sflux are constant in time, they C may be set in subroutines seamount, box or my_problem): C C An approximation to the actual temperature of the water C associated with vflux (set tatm to something else, if you know C what it is): C C tatm=t(i,j,1)+tbias C C The heat flux associated with vflux: C C tflux(i,j)=(tatm-tbias)*vfluxf(i,j) C C The total heat flux just beneath the sea surface: C (NOTE that this is all included via the vertical diffusive step) C wtsurf(i,j)=wtsurf(i,j)+tflux(i,j) C C An approximation to the actual salinity of the water associated C with vflux (set satm to something else, if you know what it is): C C satm=0.e0 C C The salt flux associated with vflux: C C sflux(i,j)=(satm-sbias)*vfluxf(i,j) C C The total salt flux just beneath the sea surface: C (NOTE that this is all included via the vertical diffusive step) C wssurf(i,j)=wssurf(i,j)+sflux(i,j) C end do end do C C----------------------------------------------------------------------- C if(mode.ne.2) then call advct call baropg C do k=1,kbm1 do j=2,jmm1 do i=2,imm1 C C Calculate Smagorinsky lateral viscosity for 3-D case: 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 ( hor visc = horcon*dx*dy*sqrt((du/dx)**2+(dv/dy)**2 C +.5*(du/dy+dv/dx)**2) ) C C The following aids symmetry of the derivation of aam with respect C to the i- and j-axes, when the velocity field is nearly uniform: C smag_1=((u(i+1,j,k)-u(i,j,k))/dx(i,j))**2 smag_2=((v(i,j+1,k)-v(i,j,k))/dy(i,j))**2 smag_3=.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) smag_4=.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) smag_a=smag_1+smag_2+.5e0*(smag_3+smag_4)**2 smag_b=smag_2+smag_1+.5e0*(smag_4+smag_3)**2 aam(i,j,k)=horcon*art(i,j)*sqrt((smag_a+smag_b)/2.e0) C end do end do end do 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 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 do j=1,jm do i=1,im egf(i,j)=el(i,j)*ispi end do end do C do j=2,jm do i=2,im utf(i,j)=ua(i,j)*(d(i,j)+d(i-1,j))*isp2i vtf(i,j)=va(i,j)*(d(i,j)+d(i,j-1))*isp2i end do end do C C----------------------------------------------------------------------- C C Beginning of external (2-D) mode: C do iext=1,isplit C C write(6,3) iext,time C 3 format(' iext,time =',i5,f9.2) C 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 freshwater flux (as the velocity, vflux), C compared with pom98.f. See also modifications to subroutine C 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 C call bcond(1) if(mod(iext,ispadv).eq.0) call advave 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))) $ *arv(i,j)*vab(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 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 C if(iext.ne.isplit) then do j=2,jm do i=2,im egf(i,j)=egf(i,j)+el(i,j)*ispi utf(i,j)=utf(i,j)+ua(i,j)*(d(i,j)+d(i-1,j))*isp2i vtf(i,j)=vtf(i,j)+va(i,j)*(d(i,j)+d(i,j-1))*isp2i end do end do endif C endif C end do C C 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 tmp2d(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tmp2d(i,j)=tmp2d(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)-tmp2d(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 tmp2d(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tmp2d(i,j)=tmp2d(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)-tmp2d(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 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 as temporary variables: C call advq(q2b,q2,uf) call advq(q2lb,q2l,vf) C if(idens.eq.1) then C C Included for version 1.8.0 and 1.9.* only: C call profq_old else call profq endif C 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 as temporary variables: C if(mode.ne.4) then C if(nadv.eq.1) then C call advt1(tb,t,tclim,uf) call advt1(sb,s,sclim,vf) C else if(nadv.eq.2) then C call advt2(tb,t,tclim,uf,nitera,sw) call advt2(sb,s,sclim,vf,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) call proft(vf,wssurf,ssurf,nbcs) call bcond(4) 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 if(idens.eq.1) then C C Included for version 1.8.0 and 1.9.* only: C call dens_old(s,t,rho) else call dens(s,t,dep0,e_atmos0,zz,rho) endif C endif C C Calculate uf and vf: C call advu call advv call profu call profv call bcond(3) C do j=1,jm do i=1,im tmp2d(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tmp2d(i,j)=tmp2d(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)-tmp2d(i,j)) end do end do end do C do j=1,jm do i=1,im tmp2d(i,j)=0.e0 end do end do C do k=1,kbm1 do j=1,jm do i=1,im tmp2d(i,j)=tmp2d(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)-tmp2d(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)*fsm(i,j) utb(i,j)=utf(i,j) vtb(i,j)=vtf(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 =',f10.4,', iint =',i8,', iext =',i8, $ ', iprint =',i8,//) C C Select print statements in printall as desired: C call printall C 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=art(i,j)*fsm(i,j) 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=art(i,j)*fsm(i,j) atot=atot+darea eaver=eaver+et(i,j)*darea end do end do C taver=taver/vtot saver=saver/vtot eaver=eaver/atot C write(6,5) vtot,atot,eaver,taver,saver 5 format(' vtot = ',e20.8,' atot = ',e20.8/ $ ' eaver =',e20.8,' taver =',e20.8, $ ' saver =',e20.8/) C C Write netCDF output: C call write_netcdf(netcdf_file,2) ! *netCDF* 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 call write_netcdf(netcdf_file,3) ! *netCDF* C stop C endif C endif C C End of print section C C----------------------------------------------------------------------- C end do C C End of numerical integration C C----------------------------------------------------------------------- C write(6,4) time,iint,iext,iprint 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 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, $ vfluxf,wssurf,wtsurf,wusurf,wvsurf C C Close netCDF file: C call write_netcdf(netcdf_file,3) ! *netCDF* C stop C end C C End of main program C C----------------------------------------------------------------------- C subroutine advave C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: advave * C * * C * FUNCTION : Calculates horizontal advection and diffusion. * C * * C ********************************************************************** C implicit none C include 'ozpom.c' C real curv2d(im,jm) C 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-1)+ua(i,j)) 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))*dy(i,j) end do end do C do j=2,jm do i=2,im tmp2d(i,j)=.25e0*((d(i,j)+d(i-1,j))+(d(i,j-1)+d(i-1,j-1))) $ *((aam2d(i,j)+aam2d(i-1,j)) $ +(aam2d(i,j-1)+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,j-1))+(dx(i-1,j)+dx(i-1,j-1)))) fluxva(i,j)=(fluxva(i,j)-tmp2d(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 v-advection and diffusion: C C Advective fluxes: C do j=1,jm do i=1,im advva(i,j)=0.e0 end do end do 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))*dx(i,j) end do end do C do j=2,jm do i=2,im tmp2d(i,j)=.25e0*((d(i,j)+d(i,j-1))+(d(i-1,j)+d(i-1,j-1))) $ *((aam2d(i,j)+aam2d(i,j-1)) $ +(aam2d(i-1,j)+aam2d(i-1,j-1))) $ *((vab(i,j)-vab(i-1,j)) $ /((dx(i,j)+dx(i,j-1))+(dx(i-1,j)+dx(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)))) fluxua(i,j)=(fluxua(i,j)-tmp2d(i,j))*.25e0 $ *((dy(i,j)+dy(i,j-1))+(dy(i-1,j)+dy(i-1,j-1))) end do end do C do j=2,jmm1 do i=2,imm1 advva(i,j)=(fluxva(i,j)-fluxva(i,j-1)) $ +(fluxua(i+1,j)-fluxua(i,j)) 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-1,j)) $ +(vab(i,j+1)+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,j-1)) $ +(uab(i+1,j)+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))) $ /art(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 C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: advct * 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 'ozpom.c' C real dtaam C 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))) $ /art(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,j-1)) $ +(dx(i-1,j)+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,j-1))+(dt(i-1,j)+dt(i-1,j-1))) $ *((aam(i,j,k)+aam(i,j-1,k)) $ +(aam(i-1,j,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,j-1)) $ +(dx(i-1,j)+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,j-1)) $ +(dy(i-1,j)+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) C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: advq * C * * C * FUNCTION : Calculates horizontal advection and diffusion, and * C * vertical advection for turbulent quantities. * C * * C ********************************************************************** C implicit none C include 'ozpom.c' C real qb(im,jm,kb),q(im,jm,kb),qf(im,jm,kb) C 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))) $ *(dt(i,j)+dt(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))) $ *(dt(i,j)+dt(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) C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: advt1 * 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 'ozpom.c' C real fb(im,jm,kb),f(im,jm,kb),fclim(im,jm,kb) real ff(im,jm,kb) C 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)) $ *(dt(i,j)+dt(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)) $ *(dt(i,j)+dt(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)=0.e0 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,nitera,sw) C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: advt2 * 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 * 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 'ozpom.c' C real fb(im,jm,kb),f(im,jm,kb),fclim(im,jm,kb) real ff(im,jm,kb) C real sw C integer nitera C real eta(im,jm) C integer i,j,k,itera C C Calculate horizontal mass fluxes: C do k=1,kb do j=1,jm do i=1,im xmflx(i,j,k)=0.e0 ymflx(i,j,k)=0.e0 end do end do end do C do k=1,kbm1 do j=2,jmm1 do i=2,im xmflx(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 ymflx(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 zwflx(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 $ *((xmflx(i,j,k)+abs(xmflx(i,j,k))) $ *fbmem(i-1,j,k)+ $ (xmflx(i,j,k)-abs(xmflx(i,j,k))) $ *fbmem(i,j,k)) C yflux(i,j,k)=0.5e0 $ *((ymflx(i,j,k)+abs(ymflx(i,j,k))) $ *fbmem(i,j-1,k)+ $ (ymflx(i,j,k)-abs(ymflx(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 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 $ *((zwflx(i,j,k)+abs(zwflx(i,j,k))) $ *fbmem(i,j,k)+ $ (zwflx(i,j,k)-abs(zwflx(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(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 xmflx(i,j,k)=0.5e0*(aam(i,j,k)+aam(i-1,j,k)) ymflx(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)=-xmflx(i,j,k)*(dt(i,j)+dt(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)=-ymflx(i,j,k)*(dt(i,j)+dt(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 * OzPOM SOURCE CODE * 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 'ozpom.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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: advv * 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 'ozpom.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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: areas_masks * C * * C * FUNCTION : Calculates areas and masks. * C * * C ********************************************************************** C implicit none C include 'ozpom.c' C integer i,j C C Calculate areas of free surface 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: C 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 C return C end C subroutine baropg C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: baropg * C * * C * FUNCTION : Calculates baroclinic pressure gradient. * C * * C ********************************************************************** C implicit none C include 'ozpom.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)) $ +(et(i,j)-et(i-1,j))*2.e0) $ *((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)) $ +(et(i,j)-et(i,j-1))*2.e0) $ *((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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: bcond * 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 C real a(im,jm) C real time,scala C integer iskp,jskp C character label*(*) C real amx,scale C integer i,ib,ie,j,jwr,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.eq.0.e0) then amx=0.e0 do j=1,jm,jskp do i=1,im,iskp amx=max(abs(a(i,j)),amx) end do end do if(amx.lt.1.e-20) then scale=1.e0 else scale=10.e0**(int(log10(amx)+100.e0)-103) endif else if(scala.gt.0.e0) then scale=scala else scale=1.e0 endif C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f10.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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: prxyz * 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 C real a(im,jm,kb) C real time,scala C integer ko(*) C integer iskp,jskp,nko C character label*(*) C real amx,scale C 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.eq.0.e0) then amx=0.e0 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 if(amx.lt.1.e-20) then scale=1.e0 else scale=10.e0**(int(log10(amx)+100.e0)-103) endif else if(scala.gt.0.e0) then scale=scala else scale=1.e0 endif C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f10.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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: prxz * 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 C real a(im,jm,kb),dt(im,jm),zz(kb) C real time,scala C integer jo(*) C integer iskp,njo C character label*(*) C real amx,scale C integer i,ib,ie,j,k,ijo,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.eq.0.e0) then amx=0.e0 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 if(amx.lt.1.e-20) then scale=1.e0 else scale=10.e0**(int(log10(amx)+100.e0)-103) endif else if(scala.gt.0.e0) then scale=scala else scale=1.e0 endif C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f10.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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: pryz * 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 C integer im,jm,kb C real a(im,jm,kb),dt(im,jm),zz(kb) C real time,scala C integer io(*) C integer jskp,nio C character label*(*) C real amx,scale C integer i,j,jb,je,k,iio,cols C if(scala.ge.0.e0) then cols=24 else cols=12 endif C if (scala.eq.0.e0) then amx=0.e0 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 if(amx.lt.1.e-20) then scale=1.e0 else scale=10.e0**(int(log10(amx)+100.e0)-103) endif else if(scala.gt.0.e0) then scale=scala else scale=1.e0 endif C write(6,1) label 1 format(1x,a40/) write(6,2) time,scale 2 format(' Time = ',f10.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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: seamount * C * * C * FUNCTION : Sets up seamount problem. * C * * C * NOTE ON DEFINITION OF GRID SIZES * C * * C * dx(i,j) is the distance between the points occupied * C * by u(i,j,k) and u(i+1,j,k). * C * dy(i,j) is the distance between the points occupied * C * by v(i,j,k) and v(i,j+1,k). * C * art(i,j) is the area of the cell with the point * C * occupied by el(i,j) or t(i,j,k) or s(i,j,k) at its * C * centre. * C * * C * Therefore art(i,j)=dx(i,j)*dy(i,j). * C * * C ********************************************************************** C implicit none C include 'ozpom.c' C real delh,delx,elejmid,elwjmid,ra,vel C integer i,j,k C C Check size of problem is correct: C if(im.ne.65.or.jm.ne.49) then write(6,2) 2 format('For this problem, im = 65 and jm = 49:'/ $ 'edit file ozpom.c as appropriate ', $ '..... program terminated'/) stop endif 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 do i=1,im do j=1,jm C h(i,j)=4500.e0*(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_atmos0, e_atmos, vfluxf, tflux, C sflux, wusurf, wvsurf, wtsurf, wssurf and swrad, as necessary C (NOTE: C 1. These have all been initialised to zero in the main program. C 2. e_atmos is the deviation of the atmospheric pressure from a C mean value, expressed in units of "height of seawater". C 3. e_atmos0 is the long-term average of e_atmos. C 4. NOTE that, over land, e_atmos0 and e_atmos should be less C than h. C 5. The temperature and salinity of inflowing water must be C defined relative to tbias and sbias. C 6. Read "NOTE ON TREATMENT OF TEMPERATURE AND SALINITY" in C program header.): C do j=1,jm do i=1,im C C Long-term average total water depth, allowing for inverse C barometer (for use in deriving density for baroclinic pressure C gradient): C dep0(i,j)=h(i,j)-e_atmos0(i,j) end do end do C C Initialise elb, etb and dt: 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)*fsm(i,j) end do end do C if(idens.eq.1) then C C Included for version 1.8.0 and 1.9.* only: C call dens_old(sb,tb,rho) else call dens(sb,tb,dep0,e_atmos0,zz,rho) endif 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 * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: slpmax * 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 'ozpom.c' C real mean,del C 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(ff,sw) C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: smol_adif * 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 'ozpom.c' C real ff(im,jm,kb) C real sw C real mol,abs_1,abs_2 real value_min,epsilon real udx,u2dt,vdy,v2dt,wdz,w2dt C 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 xmflx(i,j,k)=0.e0 else udx=abs(xmflx(i,j,k)) u2dt=dti2*xmflx(i,j,k)*xmflx(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) xmflx(i,j,k)=(udx-u2dt)*mol*sw abs_1=abs(udx) abs_2=abs(u2dt) if(abs_1.lt.abs_2) xmflx(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 ymflx(i,j,k)=0.e0 else vdy=abs(ymflx(i,j,k)) v2dt=dti2*ymflx(i,j,k)*ymflx(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) ymflx(i,j,k)=(vdy-v2dt)*mol*sw abs_1=abs(vdy) abs_2=abs(v2dt) if(abs_1.lt.abs_2) ymflx(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 zwflx(i,j,k)=0.e0 else wdz=abs(zwflx(i,j,k)) w2dt=dti2*zwflx(i,j,k)*zwflx(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) zwflx(i,j,k)=(wdz-w2dt)*mol*sw abs_1=abs(wdz) abs_2=abs(w2dt) if(abs_1.lt.abs_2)zwflx(i,j,k)=0.e0 end if end do end do end do C return C end C subroutine symm(iproblem) C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: symm * C * * C * FUNCTION : Sets up symmetry box problem. * C * * C * This basin is square, has a flat bottom, is * C * surrounded by walls and is initialised with uniform * C * salinity and temperature. It is forced by an * C * atmospheric pressure gradient and wind stress in * C * the either x- or y-directions. Coriolis is zero, so * C * as to ensure a solution that is both mirror- * C * symmetric about the forcing axis and symmetric with * C * respect to interchange of axes. * C * * C * For notes on grid sizes, horizontal coordinates of * C * grid points, and other details see subroutines * C * seamount or box. * C * * C ********************************************************************** C implicit none C include 'ozpom.c' C integer iproblem C real depth,delx,tatm,satm C integer i,j,k C C Check size of problem is correct: C if(im.ne.49.or.jm.ne.49) then write(6,2) 2 format('For this problem, im = 49 and jm = 49:'/ $ 'edit file ozpom.c as appropriate ', $ '..... program terminated'/) stop endif C C Water depth: C depth=4500.e0 C C Grid size: C delx=8000.e0 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 Note that reflectional symmetry is assured here: C dx(i,j)=delx-delx*cos(pi*float(abs(i*2-im-1)) $ /float(im*2))/2.e0 dy(i,j)=delx-delx*cos(pi*float(abs(j*2-jm-1)) $ /float(jm*2))/2.e0 C cor(i,j)=0.e0 C end do end do C C Calculate horizontal coordinates of grid points and rotation C angle. 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 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 do j=1,jm do i=1,im rot(i,j)=0.e0 end do end do C C Define depth: C do i=1,im do j=1,jm C h(i,j)=depth C C Note that reflectional symmetry is assured here: C h(i,j)=depth $ -depth*(cos(pi*float(abs(i*2-im-1))/float(im*2)) $ *cos(pi*float(abs(j*2-jm-1))/float(jm*2))) $ /4.e0 end do end do C C Close the north and south boundaries: C do i=1,im h(i,1)=1.e0 h(i,jm)=1.e0 end do C C Close the east and west boundaries: C do j=1,jm h(1,j)=1.e0 h(im,j)=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 tbias and sbias here for test (tbias and sbias would C normally only be set in the main program): C tbias=10.e0 sbias=20.e0 write(6,1) tbias,sbias 1 format(/' tbias and sbias changed in subroutine symm to:'/ $ 2f10.3//) C C Set initial conditions: C do k=1,kbm1 do j=1,jm do i=1,im tb(i,j,k)=20.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) end do end do end do C C Set surface boundary conditions, e_atmos0, e_atmos, vfluxf, tflux, C sflux, wusurf, wvsurf, wtsurf, wssurf and swrad, as necessary C (NOTE: C 1. These have all been initialised to zero in the main program. C 2. e_atmos is the deviation of the atmospheric pressure from a C mean value, expressed in units of "height of seawater". C 3. e_atmos0 is the long-term average of e_atmos. C 4. NOTE that, over land, e_atmos0 and e_atmos should be less C than h. C 5. The temperature and salinity of inflowing water must be C defined relative to tbias and sbias. C 6. Read "NOTE ON TREATMENT OF TEMPERATURE AND SALINITY" in C program header.): C do j=1,jm do i=1,im C if(iproblem.eq.3) then e_atmos(i,j)=real(i)/25.e0-1.e0 wusurf(i,j)=0.1e0 else e_atmos(i,j)=real(j)/25.e0-1.e0 wvsurf(i,j)=0.1e0 endif C C Ensure atmospheric pressure cannot make water depth go negative: C e_atmos(i,j)=min(e_atmos(i,j),h(i,j)) e_atmos0(i,j)=e_atmos(i,j) C C Long-term average total water depth, allowing for inverse C barometer (for use in deriving density for baroclinic pressure C gradient): C dep0(i,j)=h(i,j)-e_atmos0(i,j) C vfluxf(i,j)=0.e0 C C See main program, just after "Begin numerical integration", for C an explanation of these terms: C tatm=20.e0 tflux(i,j)=(tatm-tbias)*vfluxf(i,j) satm=35.e0 sflux(i,j)=(satm-sbias)*vfluxf(i,j) C end do end do C C Initialise elb, etb and dt: 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)*fsm(i,j) end do end do C if(idens.eq.1) then C C Included for version 1.8.0 and 1.9.* only: C call dens_old(sb,tb,rho) else call dens(sb,tb,dep0,e_atmos0,zz,rho) endif 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 this problem, all lateral boundaries are closed through C the specification of the masks fsm, dum and dvm): C rfe=1.e0 rfw=1.e0 rfn=1.e0 rfs=1.e0 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 vertvl C ********************************************************************** C * * C * OzPOM SOURCE CODE * C * * C * ROUTINE NAME: vertvl * C * * C * FUNCTION : Calculates vertical velocity. * C * * C ********************************************************************** C implicit none C include 'ozpom.c' C 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, in order to include surface freshwater flux, the C surface velocity has been set to vflux (i.e. vfluxb(i,j) and C vfluxf(i,j)). See also change made to 2-D volume conservation C equation which calculates elf. C do j=1,jm do i=1,im w(i,j,1)=(vfluxb(i,j)+vfluxf(i,j))/2.e0 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))) $ /art(i,j) $ +(etf(i,j)-etb(i,j))/dti2) end do end do end do C return C end C include 'ozpom.n' ! *netCDF* C C End of source code C C----------------------------------------------------------------------- C