[Dart-dev] DART/branches Revision: 11963

[dart at ucar.edu]

nancy at ucar.edu
2017-09-29 17:00:21 -0600 (Fri, 29 Sep 2017)
681
change the handling for surface obs - interpolate from
the z face array, not the z center array, which changes
depending on the height of level 2.   also define surface
obs as having a scale height of -log(1) instead of computing
two values which may be almost but not quite exactly 1.

if the vertical location of an obs is at level 1, return
a vertical location of lower=1, upper=1 and fraction = 0.
this will work correctly with both 2d and 3d fields and
will not generate an out-of-bounds array access for 2d
fields (which have no level 2).  

support putting surface variables into the state vector
to compute forward operators for temperature, winds,
and specific humidity.




Modified: DART/branches/rma_trunk/models/mpas_atm/model_mod.f90
===================================================================
--- DART/branches/rma_trunk/models/mpas_atm/model_mod.f90	2017-09-29 22:56:59 UTC (rev 11962)
+++ DART/branches/rma_trunk/models/mpas_atm/model_mod.f90	2017-09-29 23:00:21 UTC (rev 11963)
@@ -60,6 +60,10 @@
                              get_name_for_quantity,   &
                              QTY_SURFACE_ELEVATION,  &
                              QTY_SURFACE_PRESSURE,   &
+                             QTY_10M_U_WIND_COMPONENT, &
+                             QTY_10M_V_WIND_COMPONENT, &
+                             QTY_2M_TEMPERATURE,  &
+                             QTY_2M_SPECIFIC_HUMIDITY,  &
                              QTY_VERTICAL_VELOCITY,  &
                              QTY_POTENTIAL_TEMPERATURE, &
                              QTY_EDGE_NORMAL_SPEED,  &
@@ -1046,6 +1050,14 @@
          goodkind = .true.
       case (QTY_PRECIPITABLE_WATER) 
          goodkind = .true.
+      case (QTY_10M_U_WIND_COMPONENT)
+         goodkind = .true.
+      case (QTY_10M_V_WIND_COMPONENT)
+         goodkind = .true.
+      case (QTY_2M_TEMPERATURE)
+         goodkind = .true.
+      case (QTY_2M_SPECIFIC_HUMIDITY)
+         goodkind = .true.
       case (QTY_U_WIND_COMPONENT,QTY_V_WIND_COMPONENT)
          ! if the reconstructed winds at the cell centers aren't there,
          ! we can use the edge normal winds, if the user allows it.
@@ -1154,21 +1166,16 @@
       print *, 'model_interpolate: SH ', istatus, expected_obs, trim(locstring)
 
 else if (obs_kind == QTY_SURFACE_ELEVATION) then
-   do e = 1, ens_size
-      location_tmp(e) = set_location(llv(1),llv(2),1.0_r8,VERTISLEVEL)
-   enddo
-   ! why do you have to call vert_convert for surface?
-   call convert_vert_distrib(state_handle, ens_size, location_tmp, QTY_SURFACE_ELEVATION, VERTISHEIGHT, istatus)
-   where (istatus /= 0) expected_obs = missing_r8   ! FIXME: this might not be necessary
-   if ( all(istatus /= 0 ) ) goto 100
 
-   do e = 1, ens_size
-      if (istatus(e) /= 0) then
-         expected_obs(e) = MISSING_R8
-      else
-         expected_obs(e) = query_location(location_tmp(e), 'VLOC')
-      endif
-   enddo
+   call compute_elevation_with_barycentric(location, expected_obs(1), istatus(1))
+   istatus(2:ens_size) = istatus(1)
+   if (istatus(1) /= 0) then
+      expected_obs = missing_r8
+      goto 100
+   endif
+
+   expected_obs(2:ens_size) = expected_obs(1)
+
    if (debug > 10) &
       print *, 'model_interpolate: SURFACE_ELEVATION', istatus, expected_obs, trim(locstring)
 
@@ -1181,16 +1188,6 @@
    where (istatus /= 0) expected_obs = missing_r8   ! FIXME: this might not be necessary
    if ( all(istatus /= 0) ) goto 100
 
-else if (obs_kind == QTY_SURFACE_PRESSURE) then
-   tvars(1) = ivar
-   location_tmp(1) = set_location(llv(1),llv(2),1.0_r8,VERTISSURFACE)
-   call compute_scalar_with_barycentric(state_handle, ens_size, location_tmp(1), 1, tvars, values, istatus)
-   expected_obs = values(1, :)
-   where (istatus /= 0) expected_obs = missing_r8   ! FIXME: this might not be necessary
-   if ( all(istatus /= 0) ) goto 100
-   if (debug > 10) &
-      print *, 'model_interpolate: SURFACE_PRESSURE', istatus, expected_obs, trim(locstring)
-  
 else
    ! direct interpolation, kind is in the state vector
 
@@ -1230,6 +1227,7 @@
 ! returning consistent values and rc codes, both these tests can
 ! be removed for speed.  FIXME.
 do e = 1, ens_size
+
    if ((istatus(e) /= 0 .and. expected_obs(e) /= MISSING_R8) .or. &
        (istatus(e) == 0 .and. expected_obs(e) == MISSING_R8)) then
 
@@ -4247,20 +4245,35 @@
    enddo
 
 else if(is_vertical(location, "SURFACE")) then
-   new_location(1) = set_location(llv(1), llv(2), 1.0_r8, VERTISLEVEL)
 
-   ! Need to get base offsets for the potential temperature, density, and water
-   ! vapor mixing fields in the state vector
-   ivars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE)
-   ivars(2) = get_progvar_index_from_kind(QTY_DENSITY)
-   ivars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO)
+   ivars(1) = get_progvar_index_from_kind(QTY_SURFACE_PRESSURE)
 
-   call compute_scalar_with_barycentric (state_handle, ens_size, new_location(1), 3, ivars, values, istatus)
-   if ( all(istatus/= 0) ) return