[Dart-dev] DART/branches Revision: 12999

[dart at ucar.edu]

nancy at ucar.edu
2019-03-04 11:15:03 -0700 (Mon, 04 Mar 2019)
132
add a namelist for overriding surface altimeter check
(for soyoung); comment out print statements that were
left in for debugging.




Modified: DART/branches/recam/models/mpas_atm/model_mod.f90
===================================================================
--- DART/branches/recam/models/mpas_atm/model_mod.f90	2019-03-04 18:04:16 UTC (rev 12998)
+++ DART/branches/recam/models/mpas_atm/model_mod.f90	2019-03-04 18:15:03 UTC (rev 12999)
@@ -222,6 +222,7 @@
 character(len=32)  :: calendar = 'Gregorian'
 real(r8)           :: highest_obs_pressure_mb   = 100.0_r8    ! do not assimilate obs higher than this level.
 real(r8)           :: sfc_elev_max_diff = -1.0_r8    ! do not assimilate if |model - station| height is larger than this [m].
+logical            :: always_assim_surf_altimeters = .false.
 integer            :: xyzdebug = 0
 integer            :: debug = 0   ! turn up for more and more debug messages
 
@@ -998,7 +999,7 @@
 verttype = nint(query_location(location))
 
 if (debug > 10) &
-print *, 'task ', my_task_id(), ' model_interpolate: obs_kind', obs_kind,' at', trim(locstring)
+   print *, 'task ', my_task_id(), ' model_interpolate: obs_kind', obs_kind,' at', trim(locstring)
 
 cellid = find_closest_cell_center(llv(2), llv(1))
 if (debug > 0) print *, ' model_interpolate: cellid', cellid  ! SYHA
@@ -1010,16 +1011,17 @@
 endif
 
 ! Reject obs if the station height is far way from the model terrain.
-! HK is this the same across the ensemble?
+! HK is this the same across the ensemble? yes, obs location is.
 if(is_vertical(location, "SURFACE").and. sfc_elev_max_diff >= 0) then
    if(abs(llv(3) - zGridFace(1,cellid)) > sfc_elev_max_diff) then
       !Soyoung: No threshold for surface altimeter 
-      !if(obs_kind == QTY_SURFACE_PRESSURE .or. obs_kind == QTY_SURFACE_ELEVATION) then
-      !   istatus = 0
-      !else
+      if (always_assim_surf_altimeters .and. &
+          (obs_kind == QTY_SURFACE_PRESSURE .or. obs_kind == QTY_SURFACE_ELEVATION)) then
+         istatus = 0
+      else
          istatus = 12
          goto 100
-      !endif
+      endif
    endif
 endif
 
@@ -4933,9 +4935,14 @@
 vindex = 1
 nedges = nEdgesOnCell(cellid)
 do i=1, nedges
-print *, 'ft: i: ', i
+!print *, 'ft: i: ', i
    edgeid = edgesOnCell(i, cellid)
-print *, 'ft: edgeid: ', edgeid
+!print *, 'ft: edgeid: ', edgeid
+   if (.not. global_grid .and. &
+      (cellsOnEdge(1, edgeid) == 0 .or. cellsOnEdge(2, edgeid) == 0)) then
+      ier = 14
+      return
+   endif
    if (cellsOnEdge(1, edgeid) /= cellid) then
       neighborcells(i) = cellsOnEdge(1, edgeid)
    else
@@ -4942,7 +4949,7 @@
       neighborcells(i) = cellsOnEdge(2, edgeid)
    endif
    verts(i) = verticesOnCell(i, cellid)
-print *, 'ft: verts: ', verts(i), closest_vert
+!print *, 'ft: verts: ', verts(i), closest_vert
    if (verts(i) == closest_vert) vindex = i
    call latlon_to_xyz(latCell(neighborcells(i)), lonCell(neighborcells(i)), &
       xdata(i), ydata(i), zdata(i))
@@ -4954,7 +4961,7 @@
 
 ! and the observation point
 call latlon_to_xyz(lat, lon, r(1), r(2), r(3))
-print *, 'ft: lat/lon: ', lat, lon
+!print *, 'ft: lat/lon: ', lat, lon
 
 if (all(abs(t1-r) < roundoff)) then   ! Located at a grid point (counting roundoff errors)