[Dart-dev] DART/branches Revision: 12465

[dart at ucar.edu]

thoar at ucar.edu
2018-03-25 21:24:38 -0600 (Sun, 25 Mar 2018)
124
get_grid_vertval() works.
Must check compute_gridcell for pft-based quantities.
Have not checked the snow handling at all.




Modified: DART/branches/cesm_clm/models/clm/model_mod.f90
===================================================================
--- DART/branches/cesm_clm/models/clm/model_mod.f90	2018-03-26 00:22:05 UTC (rev 12464)
+++ DART/branches/cesm_clm/models/clm/model_mod.f90	2018-03-26 03:24:38 UTC (rev 12465)
@@ -47,7 +47,7 @@
                              get_close_state, convert_vertical_obs,           &
                              convert_vertical_state
 
-use    utilities_mod, only : register_module, error_handler,                   &
+use    utilities_mod, only : register_module, error_handler, E_ALLMSG,         &
                              E_ERR, E_WARN, E_MSG, logfileunit, get_unit,      &
                              do_output, to_upper,                              &
                              find_namelist_in_file, check_namelist_read,       &
@@ -1970,10 +1970,10 @@
 llat      = loc_array(2)
 lheight   = loc_array(3)
 
-call write_location(0,location,string2)
+call write_location(0,location,charstring=string2)
 if ((debug > 6) .and. do_output()) then
-   write (string3,*) llon, llat, lheight
-   call error_handler(E_MSG, routine, 'requesting interpolation at ', &
+   write (string3,'(3(2x,f18.12))') llon, llat, lheight
+   call error_handler(E_MSG, routine, 'requesting interpolation at', &
                       text2=string2, text3=string3)
 endif
 
@@ -2207,8 +2207,8 @@
 
 !------------------------------------------------------------------
 !> Calculate the expected vertical value for the gridcell.
-!> Each gridcell value is an area-weighted value of an unknown number of
-!> column-based quantities.
+!> Each gridcell value is an area-weighted value of an unknown number
+!> of column-based quantities.
 
 subroutine get_grid_vertval(state_handle, ens_size, location, qty_index, interp_val, istatus)
 
@@ -2232,17 +2232,16 @@
 real(r8) :: topwght, botwght
 real(r8), dimension(1) :: loninds,latinds
 
-real(r8), allocatable, dimension(:, :) :: above, below
-real(r8), allocatable, dimension(:, :) :: area_above
-real(r8), allocatable, dimension(:, :) :: area_below
-integer :: counter, counter_above, counter_below
+real(r8), allocatable :: above(:,:), area_above(:,:)
+real(r8), allocatable :: below(:,:), area_below(:,:)
+integer,  allocatable :: counter_above(:), counter_below(:)
+integer :: counter
+integer :: levelabove, levelbelow
 integer :: imem
 real(r8) :: state(ens_size)
 character(len=obstypelength) :: varstring
+logical :: matched
 
-call error_handler(E_MSG,routine,'currently being tested ...', &
-                    source, revision, revdate)
-
 if ( .not. module_initialized ) call static_init_model
 
 ! Let's assume failure.  Set return val to missing, then the code can
@@ -2252,7 +2251,7 @@
 ! make any error codes set here be in the 20s
 
 interp_val = MISSING_R8  ! the DART bad value flag
-istatus    = 99          ! unknown error
+istatus    = 0           ! presume that it will work
 
 loc        = get_location(location)  ! loc is in DEGREES
 loc_lon    = loc(1)
@@ -2259,6 +2258,8 @@
 loc_lat    = loc(2)
 loc_lev    = loc(3)
 
+!>@todo might need to relax this if we get observations IN the snow layers.
+!> what happens with canopy?
 if ( loc_lev < 0.0_r8 ) then
    write(string1,*)'Cannot support above-ground vertical interpolation.'
    write(string2,*)'requested a value at a depth of ',loc_lev
@@ -2303,29 +2304,35 @@
 ! The above-ground 'depths' are calculated from ZISNO and are negative.
 ! The 'depths' are all positive numbers, increasingly positive is deeper.
 ! The variables currently supported use the subsurface definitions in
-! the module variable LEVNGRND.
+! the module variable LEVNGRND -- number of layers is nlevgrnd
 
 if (loc_lev  <= LEVGRND(1)) then  ! the top level is so close to the surface
    depthabove = LEVGRND(1)        ! just use the top level
    depthbelow = LEVGRND(1)
+   levelabove = 1
+   levelbelow = 1
 elseif (loc_lev >= maxval(LEVGRND)) then  ! at depth, however ... do we
    depthabove    = maxval(LEVGRND)        ! fail or just use the deepest
    depthbelow    = maxval(LEVGRND)        ! I am using the deepest.
+   levelabove    = nlevgrnd
+   levelbelow    = nlevgrnd
 else
-
-   LAYERS : do indexi = 2,size(LEVGRND)
-      if (loc_lev < LEVGRND(indexi)) then