[Dart-dev] DART/branches Revision: 12331

[dart at ucar.edu]

nancy at ucar.edu
2018-01-16 10:49:53 -0700 (Tue, 16 Jan 2018)
198
fix the scale height (negative log) conversions.
fix up the comments some more on the height conversion.
have johnny help me with better mathematical description
of what's actually going on there.




Modified: DART/branches/recam/models/cam-fv/model_mod.f90
===================================================================
--- DART/branches/recam/models/cam-fv/model_mod.f90	2018-01-16 16:44:36 UTC (rev 12330)
+++ DART/branches/recam/models/cam-fv/model_mod.f90	2018-01-16 17:49:53 UTC (rev 12331)
@@ -693,6 +693,7 @@
 !> istatus = 15   can not get indices from given state vector index
 !> istatus = 16   cannot do vertical interpolation for bottom layer
 !> istatus = 17   cannot do vertical interpolation for top layer
+!> istatus = 98   unknown error - shouldn't happen
 !> istatus = 99   unknown error - shouldn't happen
 !>
 
@@ -2630,18 +2631,31 @@
 !#! 
 
 !-----------------------------------------------------------------------
-! Purpose:
-!   To compute geopotential height using the CCM2 hybrid coordinate
-!   vertical slice.  Since the vertical integration matrix is a
-!   function of latitude and longitude, it is not explicitly
-!   computed as for sigma coordinates.  The integration algorithm
-!   is derived from Boville's mods in the ibm file hybrid 1mods
-!   (6/17/88).  All vertical slice arrays are oriented top to
-!   bottom as in CCM2.  This field is on full model levels (aka
-!   "midpoints") not half levels.
-!
-! Equation references are to "Hybrid Coordinates for CCM1"
-!    https://opensky.ucar.edu/islandora/object/technotes%3A149/datastream/PDF/view
+!> This code is using a finite difference method to evaluate an 
+!> integral to solve the hydrostatic equation. 
+!>
+!> The details are in the reference given below.
+!> Don't change this code until you have read the paper and
+!> understand what they're doing.  The paper uses a matrix
+!> while this code gets away with ignoring 'l' and evaluating
+!> the 'k' vector directly. 
+!>
+!> Equation references are to "Hybrid Coordinates for CCM1"
+!>    https://opensky.ucar.edu/islandora/object/technotes%3A149/datastream/PDF/view
+!>
+!> Here is a comment from the NCL function that does the
+!> same thing for them.
+!>
+!> Purpose:
+!>   To compute geopotential height using the CCM2 hybrid coordinate
+!>   vertical slice.  Since the vertical integration matrix is a
+!>   function of latitude and longitude, it is not explicitly
+!>   computed as for sigma coordinates.  The integration algorithm
+!>   is derived from Boville's mods in the ibm file hybrid 1mods
+!>   (6/17/88).  All vertical slice arrays are oriented top to
+!>   bottom as in CCM2.  This field is on full model levels (aka
+!>   "midpoints") not half levels.
+!>
 
 subroutine build_heights(nlevels,p_surf,h_surf,virtual_temp,height_midpts,height_interf,variable_r)
 
@@ -2675,18 +2689,18 @@
    r_g0_tv(:) = const_r       / g0 * virtual_temp(:)
 endif
 
-! calculate the pressure column midpoints, which is 
-! one less than pm_ln() needs. The pressure at the
-! surface interface is at nlevels+1
+! calculate the log of the pressure column midpoints.
+! items 1:nlevels are the midpoints, but NOTICE THAT
+! the pressure at nlevels+1 is the pressure of the 
+! actual surface interface, not a midpoint!!
 
 call cam_p_col_midpts(p_surf, nlevels, p_mid)
    
-p_mid(nlevels+1) = p_surf * grid_data%hybi%vals(nlevels+1)
+p_mid(nlevels+1) = p_surf * grid_data%hybi%vals(nlevels+1)   ! surface interface
    
-! compute top only if the top interface pressure is nonzero.
-where      (pm_mid >  0.0_r8) 
+where      (p_mid >  0.0_r8) 
    pm_ln = log(pm_ln) 
-else where (pm_mid <= 0.0_r8)
+else where (p_mid <= 0.0_r8)
    pm_ln = 0
 end where
 
@@ -2700,12 +2714,6 @@
 !  write(*, 200) i, pm_ln(i)
 !enddo
 
-! Initialize height_midpts to sum of ground height and thickness of
-! top half-layer terms.
-
-pterm(1)       = 0.0_r8
-pterm(nlevels) = 0.0_r8
-
 !        height_midpts(1)=top  ->  height_midpts(nlevels)=bottom
 ! 
 ! level
@@ -2720,34 +2728,49 @@
 !              ---------------------------------------------------
 ! NL+1/2 -----/|||||||||||||||||||||||||||||||||||||||||||||||||||\-----
 
-! Midpoint layer terms
-!        Eq 3.a.109.2  where l=K,k<K  h(k,l) = 1/2 * ln [  p(k+1) / p(k-1) ]