Improvements for the Primitive Equation model

current suggestions

1. DOCUMENTATION
   • Need references like these.
   • need a users' guide.

2. Two-way nesting:
   • Need to make sure that certain auxiliary variables (e.g. D0START for forcing) are consistent between nested domains.
   • Consider these related notes in Grids, PE_mask, Cond_Topo and PE_initial.
   • Consider making nesting fault tolerant (i.e. when one domain fails, other domains continue as best as their able).
   • Cjl is interested in writing interpolated boundary data to auxiliary arrays, to defer decisions on their use. Especially the transport streamfunction / surface pressure p00.
   • CJL also had this to say.
   • Try to rewrite to avoid using Shapiro filter storage.
   • land masks -- The sending of data from under the mask is protected against, PE can still receive data under mask when receiving data from smaller domain. Play it save and block this.
   • Try to allow data passing only at certain time steps. Could be useful for large domain to pass data every time step, but small domain with smaller time step (1/3) to pass only every third time step.
   • Instead of directly interpolating boundary values from larger grid, try:
     • Constraining the averages of the fine boundary values.
     • Interpolating values in a band (1 or 2 grid points) around outside of fine grid. Use radiation conditions which make use of information on both sides.
   • Examine the following papers.
     • Perkins, A., et al.  (1997)  J. Geophys. Res., 102, 3483-3500.
     • Ginis, I., A. Richardson & L. Rothstein.  (1998)  Mon Weather Rev., 126, 1054-1079.

3. Check all netCDF files to make sure that they're of the correct type. (e.g. initialization & assimilation files should come from PE_inintial not OAG). Provide understandable messages when this isn't the case.
   • While we're on the subject of input checking, make sure that all needed files are present & correct at initialization. e.g. During 1-way nesting the file extraction.dat is checked for existence but the PE_inintial files it references are not check for existence until they're first used. They would also need to be checked for sizing (i.e. their grids fit in the allocated space), compatibility (e.g. rigid lid streamfunction vs. free surface), and type (e.g. initialization versus assimilation).

4. Islands:
   • Islands basically work in sor scheme. A problem occurs when the "island boxes" overlap. Revise the resetting to only work on points corresponding to that box.
   • Need to get islands working in conjugate gradient scheme.

5. Changes to tracer balance codes:
   • Instead of computing the total change as the sum of the parts, use ta-tb and compare to the sum.
   • Need to include effect of Shapiro filter. This entails:
     1. an extra term added to the tracer balance storage in the defnrg file to accommodate the Shapiro filter changes
     2. the fields before and after the Shapiro filter would have to be preserved, and their difference taken
     3. the total change written out in diag.F would have to be read, corrected by the difference computed in (2) and rewritten
   • Double check on the use of land masks in these balances.

6. In drifter codes:
   • Don't seem to be recording the initial drifter position in the track netCDF file. Seems to start at first computed time step. Check into this.
   • Need to rework the velocity codes to use sensible extrapolation under the land masks.
   • Revise track_bound.F to recognize locations under the mask as out-of-bounds.
   • Revise initialization to flag and remove out-of-bounds initial locations from working list. Don't record them in the track output netCDF file (saving space).
   • Density drifters jump around a lot. Instead of forcing float to hit a target density, allow vertical velocity based upon the local buoyancy.
     • A vertical acceleration is more natural. Figure out how to integrate this to get a vertical velocity for a Runge-Kutta scheme.
   • Density drifters pick their target density based on their initial location. Might want to be able to specify a target value.
   • In 2-way nesting, it would be nice if a drifter could move "between" domains.

7. In vertical diffusivity codes:
   • Currently replace background diffusivity with wind mixing and replace both with unstable values. At each step, should use a max test to avoid replacing with a smaller value.
   • The explicit mixing option in the tracers uses the variable kmt(i). The problem is: it isn't readily apparent how kmt(i) is initialized. Is this one of those cases in which the cheating oget is hiding the details?

8. Bring back a restart:
   • Allow for several time levels.
   • Make sure all of our additions are included
   • Be sure boundary relaxation maintains smooth fields

9. In defout.F: A working array is filled with the 2D positions at every requested output level. The problem -- the user can request more levels than exist in model (level numbers > km are interpolated to that depth in meters). This can lead to unpredictable problems with writing beyond array bounds.
   • Probable correction: Write each level separately. Don't load into volume array and dump in one shot. Check to see that this solves all possible similar problems.

10. In clinic.F & relax.F, the internal documentation suggests that the variable ztd should remain zero under the mask (except for islands). The problem: the boundary condition routines and the nesting don't ensure this.
    • Probable correction: Re-write those routines to pay attention to the isz & iez indexing.

11. Remove CPP option gridold. If all vertical specification in current PE_initial output isn't available, construct the necessary fields from what data is available.

12. CPP option linear_physics does NOT work with biology. The problem is the background profile matrix is always sized tinit(mprof,3) but loops 220 and 230 in tracer.F assume a size of tinit(mprof,nt+1).

13. PE Biology
    • Put a hierarchy so accidental multiple biological CPP options have a default setting.
      • Probably should also check the number of options chosen by run time so an appropriate error message can be printed.
      • This suggests a hierarchy everywhere BUT where the CPP options are printed. Then can have the error message immediately following the bio-option list.
    • in bioDuse.in, check the description of the biopos switch against that in /home/ross/duse/work/HOPS/PE/Src/bioDuse.in
    • Look for Duse database routine in ~/work/HOPS/PE/Src/dbsbprmDuse.F or the underlying RCS directory. Decide what, if any, long-term plans exist for it.
    • A copy of the biological model Carlos give to Katya was in ~afmisops/Katya/Model/PE_DandF.

14. Assimilation:
    • Change the assimilation procedure for the transport streamfunction. Either:
      • Assimilate differences in transport streamfunction. Let:

                   p     = transport streamfunction
                   P     = horizontal average of p over sea points
                   pa,Pa = similar quantities for assimilation field
              

        assimilate (pa-Pa) into (p-P) instead of pa into p.
      • Assimilate barotropic velocities. Back-out transport streamfunction.
    • Check setting of PE assimilation weights. They may be renormalized. This would be bad for fields which are expected to either be entirely assimilated or entirely ignored.

15. Pressure Gradient Calculation:
    • If keep pressure bias option, rewrite it to be done automatically. Don't require 2 PE runs.
    • Look into reformulating the pressure gradient in terms of the difference of the pressure on the box faces. Might give better conservation properties.
    • Examine the following papers.
      • Song, Y.T.  (1998)  A general pressure gradient formulation for ocean models. Part I: Scheme design and diagnostic analysis.  MON WEATHER REV, 126(12), 3213-3230.
      • Song, Y.T. & D.G. Wright  (1998)  A general pressure gradient formulation for ocean models. Part II: Energy, momentum, and bottom torque consistency.  MON WEATHER REV, 126(12), 3231-3247.
      • Sheng, J., R.J. Greatbach & D. Wright  (2001)  Improving predictive skill of ocean circulation models by adjusting currents towards climatology.  J. Geophys. Res., (in press).
        • The adjusting is done by modifying the density the pressure gradient feels with climatological values.

16. Take implicit Coriolis out from the transport relaxation operator. Should be simply a matter of "untwisting" the implicit Coriolis matrix before forming the barotropic vorticity tendency.

17. Laplacian viscosity/diffusion:
    • Currently introduces crap starting at bottom level (i.e. negative temperatures)
    • How about metric diffusion terms aren't consistent with integral formulation of the divergence of a gradient.

18. Point tracer sources (CPP option pttrcsrc):
    • As I recall, the coding was rather anemic. Check out & improve.
    • Need a sample input file & description in Readme.pe.

19. infld.F:
    • From Jeff: Pat - I think there's a bug in infld.F:

      cjad          if(imax.ne.imt) then
                if(ntrc.ne.nt) then
                  write(stdout,930) 'NT, NTRC: ',nt,ntrc
                  call exitus('INFLD')
      

20. Surface forcing:
    • From Jeff:

      It's getting a whole row of mixed layer depth incorrect (doesn't look like the fill value) at 9.9692e+36... unless the fill value was multiplied by the conversion factor. It's the boundary row at the northern part of the grid, it seems. Also, the units aren't as specified in the file (file says meters, units are cm).

    • Look in Jeff's or Larry's codes for the solar penetration option -Dsolar, it should appear in tsource.F. Decide what to do.

21. For energetics, try using counters to position data in arrays.

22. Either fix bug (RR96) in boundary relaxation or eliminate that option.

23. Make a version of PE that can be compiled in "case sensitive" FORTRAN. Use to check for violations of FLOAT/FLoaT rules.

24. remove restriction that PE model has its arrays dimensioned exactly. This will be especially handy for middle nodes in an extended telescopic nest.

25. Include Julian starting time in PE output files.

26. Surface ice. (Actually this isn't nitty-gritty, but I don't have anyplace else to keep it.)

27. On diagnostic initialization, make sure that NO filtering is done until fields are released.

28. PE model probably cannot be run with a horizontal grid smaller than 5x5. On a 4x4, jmtm2=2 and wrtcdf clearly fails to write all the information necessary. step MAY not be calling osav all the times it needs either.

Blasts from the past

1. (11/94, PFJL) infld.F:
   • "use (i,j,k) or (j,i,k): don't have to flip the data"

2. (11/94, PFJL) hvgrid.F:
   • "see set latitude in radians... Is this OK?"

3. (11/94, PFJL) for input/output netCDF files, make kc a variable. Don't use iflag(8) to pass it.

4. (11/94, PFJL) Change masking logic to preventing calculations on land points, not just wiping out the results on them.

5. (12/94, PFJL) Check xtr_trc, xtr_vel & hvgrid
   • "Error in position of u(1,1) & T(1,1)? rterm & cm2deg have to be saved and are not."
   • May be an error in bilinear interpolation (map_t).

6. (12/94, PFJL) "Replace iflag(9) with CPP option dblsigma"
   • PJH: I believe the PE code uses dblsigma but checks iflag(9) from the input file.

7. (9/95, PFJL) Request that PE model write velocity data from slab jmt-1 into netCDF output slab jmt as well. Now, the PE just leaves fill values in that slab.