[Dart-dev] DART/branches Revision: 12062

[dart at ucar.edu]

nancy at ucar.edu
2017-11-07 15:43:56 -0700 (Tue, 07 Nov 2017)
170
start to flesh out the 'find vertial levels' routine.
should support obs on model levels and pressures.
not height yet, and no support for obs above a
cutoff threshold.




Modified: DART/branches/recam/models/cam-fv/new_model_mod.f90
===================================================================
--- DART/branches/recam/models/cam-fv/new_model_mod.f90	2017-11-07 20:46:16 UTC (rev 12061)
+++ DART/branches/recam/models/cam-fv/new_model_mod.f90	2017-11-07 22:43:56 UTC (rev 12062)
@@ -213,7 +213,6 @@
 call set_calendar_type('GREGORIAN')
 
 call read_grid_info(cam_template_filename, grid_data)
-call read_cam_phis_array(cam_phis_filename)
 
 !>@todo do we need to map_qtys here?
 !>@todo do we need to set the model top related stuff here?
@@ -334,6 +333,10 @@
 !> istatus = 4    cannot locate enclosing vertical levels
 !> istatus = 5    cannot retrieve state vector values
 !> istatus = 6    cannot do vertical interpolation
+!> istatus = 7    cannot interpolate in the quad to get the values
+!> istatus = 8    cannot get vertical levels for an obs on model levels
+!> istatus = 9    cannot get vertical levels for an obs on pressure levels
+!> istatus = 10   also cannot get vertical levels for an obs on pressure levels
 !> istatus = X
 !> istatus = 99   unknown error - shouldn't happen
 !>
@@ -443,7 +446,7 @@
 ! and now here potentially we have different results for different
 ! ensemble members.  the things that can vary are dimensioned by ens_size.
 do i=1, 4
-   !^>@todo FIXME build a vertical column to find vertical numbers... # need option for linear or log scale?
+   !>@todo FIXME build a vertical column to find vertical numbers... # need option for linear or log scale?
    call find_vertical_levels(state_handle, ens_size, &
                              four_lons(i), four_lats(i), lon_lat_vert(3), &
                              which_vert, obs_qty, &
@@ -509,7 +512,7 @@
 real(r8), intent(out) :: out_vals(nitems)
 integer,  intent(out) :: my_status(nitems)
 
-! vert_fracts is is 1 is the bottom level and the inverse is the top
+! vert_fracts of 1 is 100% of the bottom level and 0 is 100% of the top level
 out_vals(:) = (botvals(:)* vert_fracts(:)) + (topvals(:) * (1.0_r8-vert_fracts(:)))
 my_status(:) = 0
 
@@ -592,25 +595,198 @@
 real(r8),            intent(out) :: vert_fracts(ens_size)
 integer,             intent(out) :: my_status(ens_size)
 
-!>@todo FIXME does this really need the state vector yet?
-bot_levs(:) = 10
-top_levs(:) = 9
-vert_fracts(:) = 0.5
-my_status(:) = 0
+integer :: bot1, top1, i, nlevels, varid
+integer(i8) :: state_indx
+real(r8) :: fract1
+real(r8) :: surf_pressure(ens_size)
+real(r8) :: pressure_array(grid_data%lev%nsize)
 
+! assume the worst
+bot_levs(:) = MISSING_I
+top_levs(:) = MISSING_I
+vert_fracts(:) = MISSING_R8
+my_status(:) = 98
+
+! number of vertical levels (midlayer points)
+nlevels = grid_data%lev%nsize
+
 select case (which_vert)
+
    case(VERTISPRESSURE)
+      ! construct a pressure column here and find the model levels
+      ! that enclose this value
+      varid = get_varid_from_kind(domain_id, QTY_SURFACE_PRESSURE)
+      state_indx = get_dart_vector_index(lon_index, lat_index, 1, domain_id, varid)
+      surf_pressure(:) = get_state(state_indx, state_handle)
+
+      !>@todo FIXME: should we figure out now or later? how many unique levels we have?
+      !> for now - do the unique culling later so we don't have to carry that count around.
+      do i=1, ens_size
+         call cam_pressure_levels(surf_pressure(i), nlevels, grid_data%hyam, grid_data%hybm, &
+                                  grid_data%P0, pressure_array)
+         call pressure_to_level(nlevels, pressure_array, vert_val, &
+                                bot_levs(i), top_levs(i), vert_fracts(i), my_status(i))
+
+      enddo
+
    case(VERTISHEIGHT)
+      ! construct a height column here and find the model levels
+      ! that enclose this value
+
    case(VERTISLEVEL)
+      ! this routine returns false if the level number is out of range.
+      if (range_set(vert_val, nlevels, bot1, top1, fract1)) then
+         my_status(:) = 8
+         return
+      endif
+
+      ! because we're given a model level as input, all the ensemble
+      ! members have the same outgoing values.
+      bot_levs(:) = bot1
+      top_levs(:) = top1
+      vert_fracts(:) = fract1
+