[Dart-dev] [10641] DART/trunk/models/ROMS/model_mod.f90: First stab at using a land mask.

[nancy at ucar.edu]

Revision: 10641
Author:   thoar
Date:     2016-08-16 16:51:16 -0600 (Tue, 16 Aug 2016)
Log Message:
-----------
First stab at using a land mask.
There is still the problem that bounds checking dies in get_reg_box_indices()
Since Nancy is very close to a complete replacement, there is no need to debug
this now.

Modified Paths:
--------------
    DART/trunk/models/ROMS/model_mod.f90

-------------- next part --------------
Modified: DART/trunk/models/ROMS/model_mod.f90
===================================================================
--- DART/trunk/models/ROMS/model_mod.f90	2016-08-16 15:46:21 UTC (rev 10640)
+++ DART/trunk/models/ROMS/model_mod.f90	2016-08-16 22:51:16 UTC (rev 10641)
@@ -171,6 +171,7 @@
    integer :: indexN        ! location in dart state vector of last  occurrence
    integer :: dart_kind
    character(len=paramname_length) :: kind_string
+   character(len=paramname_length) :: mask
    logical  :: clamping     ! does variable need to be range-restricted before
    real(r8) :: range(2)     ! being stuffed back into analysis file.
    logical  :: out_of_range_fail  ! is out of range fatal if range-checking?
@@ -182,7 +183,6 @@
 
 type(progvartype), dimension(max_state_variables) :: progvar
 
-
 ! Grid parameters - the values will be read from a
 ! standard ROMS namelist and filled in here.
 
@@ -467,10 +467,10 @@
 !print *, 'base offset now for ',trim(progvar(ivar)%varname),base_offset,top_offset
 
 ! For Sea Surface Height or Pressure don't need the vertical coordinate
-!
+! but the surface layer is the LAST layer.
 if( vert_is_surface(location) ) then
    call lon_lat_interpolate(x(base_offset:top_offset), llon, llat, &
-   obs_type, 1, interp_val, istatus)
+   obs_type, Ns_rho, interp_val, istatus)
    return
 endif
 
@@ -721,6 +721,16 @@
 
    enddo DimensionLoop
 
+   !> @TODO FIXME (check, really) Are there any variables that use the mask_psi
+
+   if (    progvar(ivar)%varname == 'u') then
+           progvar(ivar)%mask    = 'mask_u'
+   elseif (progvar(ivar)%varname == 'v') then
+           progvar(ivar)%mask    = 'mask_v'
+   else
+           progvar(ivar)%mask    = 'mask_rho'
+   endif
+
    ! this call sets: clamping, bounds, and out_of_range_fail in the progvar entry
    call get_variable_bounds(roms_state_bounds, ivar)
 
@@ -3665,10 +3675,15 @@
 tp   = 0.0_r8
 iidd = 0
 
-!>@ todo FIXME ... check ... get_val() can return MISSING_R8 ... is this possible
+! if the lonid, latid is a missing_r8, then the id is a land point
+! and no interpolation is possible. Just return a MISSING_R8 value.
 
 do i=1,Nz
-    tp(i)=get_val(lonid, latid, i, x, var_type)
+    tp(i) = get_val(lonid, latid, i, x, var_type)
+    if (tp(i) == MISSING_R8) then
+       interp_val = MISSING_R8
+       return
+    endif
     zz(i)=ZC(lonid,latid,i)
 enddo
 
@@ -3701,7 +3716,7 @@
 !>
 !> Returns the DART state value for a given lat, lon, and level index.
 !> 'get_val' will be a MISSING value if this is NOT a valid grid location
-!> (e.g. land, or below the ocean floor).
+!> (e.g. land)
 !>
 !> @param get_val the value of the DART state at the gridcell of interest.
 !> @param lon_index the index of the longitude of interest.
@@ -3730,10 +3745,6 @@
 
 ivar = get_progvar_index_from_kind(var_type)
 
-!> @ todo FIXME Implement the land masking. Must determine which mask
-!>        is appropriate for this variable ... extend progvar?
-!> could check dimension names of the variables against the dimension names of the masks
-
 Ndim3 = progvar(ivar)%numvertical
 Ndim2 = progvar(ivar)%numeta
 Ndim1 = progvar(ivar)%numxi
@@ -3742,6 +3753,19 @@
      (  lat_index < 1 .or.   lat_index > Ndim2) .or. &
      (level_index < 1 .or. level_index > Ndim3) ) return
 
+! Checking the mask. 0 is water, 1 is land. nothing else possible
+
+select case ( trim(progvar(ivar)%mask) )
+   case ('mask_u')
+      if (mask_u(  lon_index, lat_index) /= 0) return
+   case ('mask_v')
+      if (mask_v(  lon_index, lat_index) /= 0) return
+   case ('mask_rho')
+      if (mask_rho(lon_index, lat_index) /= 0) return
+   case ('mask_psi')
+      if (mask_psi(lon_index, lat_index) /= 0) return
+end select
+
 ! implicit assumption on packing order into the DART vector
 
 tt = (level_index - 1) * Ndim1 * Ndim2 + &