[Dart-dev] [9602] DART/trunk/models/ROMS/model_mod.f90: Updating some routines in accordance with newer versions from Rutgers.

nancy at ucar.edu nancy at ucar.edu
Tue Jan 26 10:33:54 MST 2016


Revision: 9602
Author:   thoar
Date:     2016-01-26 10:33:53 -0700 (Tue, 26 Jan 2016)
Log Message:
-----------
Updating some routines in accordance with newer versions from Rutgers.
Much work still to be done. Some of the bigger changes in this version are:

The land/sea masks are now available. The actual dimensions for all the coordinate
systems are explicitly available instead of having to use Nx-1, Ny-1, etc. There
are still instances in the code where Nx-1 etc. are being used.
These could be/should be replaced.

get_val() is more sophisticated than an exhaustive search, but still does not
use land/water mask for anything. Now returns a MISSING_R8 if fails. Used to return
indeterminate values.

get_reg_box_indices() is much different but still wrong ... model_mod_check fails when running
the exhaustive test_interpolate() routine.

lon_lat_interpolate() was using variable 'masked' without actually setting it.
Now checks for valid returns from get_val() at each step before continuing.

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-01-26 16:19:53 UTC (rev 9601)
+++ DART/trunk/models/ROMS/model_mod.f90	2016-01-26 17:33:53 UTC (rev 9602)
@@ -148,10 +148,10 @@
 !> Everything needed to describe a variable. Basically all the metadata from
 !> a netCDF file is stored here as well as all the information about where
 !> the variable is stored in the DART state vector.
-!> @todo FIXME ... replace numxi and numeta with simply nx or something.
-!>       since the structure is an array - one for each variable -
-!>       it doesn't make sense to declare a number of horizontal 
-!>       locations of the edges for a variable that is on cell centers ...
+!> @todo FIXME ... do we need numvertical as opposed to ZonHalf ...
+!>
+!> @todo FIXME ... add a field for what kind of mask to use
+!>
 
 type progvartype
    private
@@ -188,12 +188,32 @@
 
 ! nx, ny and nz are the size of the rho grids.
 integer :: Nx = -1, Ny = -1, Nz = -1
+
+integer :: Nxi_rho
+integer :: Nxi_u
+integer :: Nxi_v
+integer :: Nxi_psi
+integer :: Neta_rho
+integer :: Neta_u
+integer :: Neta_v
+integer :: Neta_psi
+integer :: Nxi_vert
+integer :: Neta_vert
+integer :: Ns_rho
+integer :: Ns_w
+
 real(r8), allocatable :: ULAT(:,:), ULON(:,:), &
                          TLAT(:,:), TLON(:,:), &
                          VLAT(:,:), VLON(:,:), &
                            PM(:,:),   PN(:,:), &
                          ANGL(:,:),   HT(:,:), ZC(:,:,:)
 
+integer, parameter :: i2 = SELECTED_INT_KIND(2) ! need something to coerce to NF90_SHORT
+integer(i2), allocatable :: mask_rho(:,:), &
+                            mask_psi(:,:), &
+                              mask_u(:,:), &
+                              mask_v(:,:)
+
 real(r8)        :: ocean_dynamics_timestep = 900.0_r4
 type(time_type) :: model_timestep
 
@@ -401,16 +421,15 @@
 ! Local storage
 real(r8)       :: loc_array(3), llon, llat, lheight
 integer        :: base_offset, top_offset
-real(r8)       :: tp_val
 integer        :: ivar,obs_kind
 integer        :: x_ind, y_ind,lat_bot, lat_top, lon_bot, lon_top
-real(r8)       :: x_corners(4), y_corners(4),p(4)
+real(r8)       :: x_corners(4), y_corners(4), p(4)
 
 if ( .not. module_initialized ) call static_init_model
 
 ! Successful istatus is 0
-interp_val = 0.0_r8
-istatus = 0
+interp_val = MISSING_R8
+istatus = 99
 
 ! Get the individual locations values
 loc_array = get_location(location)
@@ -443,7 +462,7 @@
 ! Find the basic offset of this field
 
 base_offset = progvar(ivar)%index1
-top_offset = progvar(ivar)%indexN
+top_offset  = progvar(ivar)%indexN
 
 !print *, 'base offset now for ',trim(progvar(ivar)%varname),base_offset,top_offset
 
@@ -451,7 +470,7 @@
 !
 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, 1, interp_val, istatus)
    return
 endif
 
@@ -460,22 +479,22 @@
 !2. do vertical interpolation at four corners
 !3. do a spatial interpolation
 
-   if(progvar(ivar)%kind_string == 'KIND_U_CURRENT_COMPONENT') then
-      call get_reg_box_indices(llon, llat, ULON, ULAT, x_ind, y_ind)
-      ! Getting corners for accurate interpolation
-      call get_quad_corners(ULON, x_ind, y_ind, x_corners)
-      call get_quad_corners(ULAT, x_ind, y_ind, y_corners)
-   elseif (progvar(ivar)%kind_string == 'KIND_V_CURRENT_COMPONENT') then
-      call get_reg_box_indices(llon, llat, VLON, VLAT, x_ind, y_ind)
-      ! Getting corners for accurate interpolation
-      call get_quad_corners(VLON, x_ind, y_ind, x_corners)
-      call get_quad_corners(VLAT, x_ind, y_ind, y_corners)
-   else
-      call get_reg_box_indices(llon, llat, TLON, TLAT, x_ind, y_ind)
-      ! Getting corners for accurate interpolation
-      call get_quad_corners(TLON, x_ind, y_ind, x_corners)
-      call get_quad_corners(TLAT, x_ind, y_ind, y_corners)
-   endif
+if(progvar(ivar)%kind_string == 'KIND_U_CURRENT_COMPONENT') then
+   call get_reg_box_indices(llon, llat, ULON, ULAT, x_ind, y_ind)
+   ! Getting corners for accurate interpolation
+   call get_quad_corners(ULON, x_ind, y_ind, x_corners)
+   call get_quad_corners(ULAT, x_ind, y_ind, y_corners)
+elseif (progvar(ivar)%kind_string == 'KIND_V_CURRENT_COMPONENT') then
+   call get_reg_box_indices(llon, llat, VLON, VLAT, x_ind, y_ind)
+   ! Getting corners for accurate interpolation
+   call get_quad_corners(VLON, x_ind, y_ind, x_corners)
+   call get_quad_corners(VLAT, x_ind, y_ind, y_corners)
+else
+   call get_reg_box_indices(llon, llat, TLON, TLAT, x_ind, y_ind)
+   ! Getting corners for accurate interpolation
+   call get_quad_corners(TLON, x_ind, y_ind, x_corners)
+   call get_quad_corners(TLAT, x_ind, y_ind, y_corners)
+endif
 
 lon_bot=x_ind
 lat_bot=y_ind
@@ -495,16 +514,23 @@
 
 !write(*,*)'model_mod: obs loc vs. model loc ', llon, llat, TLON(x_ind,y_ind),TLAT(x_ind,y_ind)
 
-call vert_interpolate(x(base_offset:top_offset),lon_bot,lat_bot,obs_kind,lheight,tp_val)
-p(1)=tp_val
-call vert_interpolate(x(base_offset:top_offset),lon_top,lat_bot,obs_kind,lheight,tp_val)
-p(2)=tp_val
-call vert_interpolate(x(base_offset:top_offset),lon_top,lat_top,obs_kind,lheight,tp_val)
-p(3)=tp_val
-call vert_interpolate(x(base_offset:top_offset),lon_bot,lat_top,obs_kind,lheight,tp_val)
-p(4)=tp_val
+! If any of these fail, we can exit (fail) immediately.
+! The interp_val and istatus values are initially set to a failed value, so just use those.
 
+call vert_interpolate(x(base_offset:top_offset),lon_bot,lat_bot,obs_kind,lheight, p(1))
+if (p(1) == MISSING_R8) return
+
+call vert_interpolate(x(base_offset:top_offset),lon_top,lat_bot,obs_kind,lheight, p(2))
+if (p(2) == MISSING_R8) return
+
+call vert_interpolate(x(base_offset:top_offset),lon_top,lat_top,obs_kind,lheight, p(3))
+if (p(3) == MISSING_R8) return
+
+call vert_interpolate(x(base_offset:top_offset),lon_bot,lat_top,obs_kind,lheight, p(4))
+if (p(4) == MISSING_R8) return
+
 call quad_bilinear_interp(llon, llat, x_corners, y_corners, p, interp_val)
+istatus = 0  ! If we get this far, all good.
 
 !write(*,*) 'Ending model interpolate ...'
 
@@ -870,9 +896,9 @@
 ! variables if we parse the state vector into prognostic variables.
 
 ! for the dimensions and coordinate variables
-integer :: nxirhoDimID,nxiuDimID,nxivDimID
-integer :: netarhoDimID,netauDimID,netavDimID
-integer :: nsrhoDimID,nswDimID
+integer :: nxirhoDimID, nxiuDimID, nxivDimID, nxipsiDimID, nxivertDimID
+integer :: netarhoDimID, netauDimID, netavDimID, netapsiDimID, netavertDimID
+integer :: nsrhoDimID, nswDimID
 
 ! for the prognostic variables
 integer :: ivar, VarID
@@ -972,30 +998,42 @@
    ! Define the new dimensions IDs
    !> @todo FIXME we have the actual original dimensions ... why are we not using them here
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_rho',  len = Nx, &
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_rho',  len = Nxi_rho, &
         dimid = nxirhoDimID),'nc_write_model_atts', 'xi_rho def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_rho', len = Ny,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_rho', len = Neta_rho,&
         dimid = netarhoDimID),'nc_write_model_atts', 'eta_rho def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='s_rho',   len = Nz,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='s_rho',   len = Ns_rho,&
         dimid = nsrhoDimID),'nc_write_model_atts', 's_rho def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='s_w',     len = Nz+1,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='s_w',     len = Ns_w,&
         dimid = nswDimID),'nc_write_model_atts','s_w def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_u',    len = Nx-1,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_u',    len = Nxi_u,&
         dimid = nxiuDimID),'nc_write_model_atts', 'xi_u def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_v',    len = Nx,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_v',    len = Nxi_v,&
         dimid = nxivDimID),'nc_write_model_atts', 'xi_v def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_u',   len = Ny,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_u',   len = Neta_u,&
         dimid = netauDimID),'nc_write_model_atts', 'eta_u def_dim '//trim(filename))
 
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_v',   len = Ny-1,&
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_v',   len = Neta_v,&
         dimid = netavDimID),'nc_write_model_atts', 'eta_v def_dim '//trim(filename))
 
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_psi',   len = Nxi_psi,&
+        dimid = nxipsiDimID),'nc_write_model_atts', 'xi_psi def_dim '//trim(filename))
+
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_psi',   len = Neta_psi,&
+        dimid = netapsiDimID),'nc_write_model_atts', 'eta_psi def_dim '//trim(filename))
+
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='xi_vert',   len = Nxi_vert,&
+        dimid = nxivertDimID),'nc_write_model_atts', 'xi_vert def_dim '//trim(filename))
+
+   call nc_check(nf90_def_dim(ncid=ncFileID, name='eta_vert',   len = Neta_vert,&
+        dimid = netavertDimID),'nc_write_model_atts', 'eta_vert def_dim '//trim(filename))
+
    ! Create the (empty) Coordinate Variables and the Attributes
 
    call nc_check(nf90_def_var(ncFileID,name='lon_rho', xtype=nf90_double, &
@@ -1054,6 +1092,38 @@
    call nc_check(nf90_put_att(ncFileID,  VarID, 'units', 'm'), &
                  'nc_write_model_atts', 'z_c units '//trim(filename))
 
+   call nc_check(nf90_def_var(ncFileID,name='mask_rho', xtype=nf90_short, &
+                 dimids=(/ nxirhoDimID, netarhoDimID /), varid=VarID),&
+                 'nc_write_model_atts', 'mask_rho def_var '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'long_name', 'mask on RHO-points'), &
+                 'nc_write_model_atts', 'mask_rho long_name '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'comment', '0==land, 1==water'), &
+                 'nc_write_model_atts', 'mask_rho comment '//trim(filename))
+
+   call nc_check(nf90_def_var(ncFileID,name='mask_psi', xtype=nf90_short, &
+                 dimids=(/ nxipsiDimID, netapsiDimID /), varid=VarID),&
+                 'nc_write_model_atts', 'mask_psi def_var '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'long_name', 'mask on psi-points'), &
+                 'nc_write_model_atts', 'mask_psi long_name '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'comment', '0==land, 1==water'), &
+                 'nc_write_model_atts', 'mask_psi units '//trim(filename))
+
+   call nc_check(nf90_def_var(ncFileID,name='mask_u', xtype=nf90_short, &
+                 dimids=(/ nxiuDimID, netauDimID /), varid=VarID),&
+                 'nc_write_model_atts', 'mask_u def_var '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'long_name', 'mask on U-points'), &
+                 'nc_write_model_atts', 'mask_u long_name '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'comment', '0==land, 1==water'), &
+                 'nc_write_model_atts', 'mask_u units '//trim(filename))
+
+   call nc_check(nf90_def_var(ncFileID,name='mask_v', xtype=nf90_short, &
+                 dimids=(/ nxivDimID, netavDimID /), varid=VarID),&
+                 'nc_write_model_atts', 'mask_v def_var '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'long_name', 'mask on V-points'), &
+                 'nc_write_model_atts', 'mask_v long_name '//trim(filename))
+   call nc_check(nf90_put_att(ncFileID,  VarID, 'comment', '0==land, 1==water'), &
+                 'nc_write_model_atts', 'mask_v units '//trim(filename))
+
    ! Create the (empty) Prognostic Variables and the Attributes
 
    do ivar=1, nfields
@@ -1135,6 +1205,26 @@
    call nc_check(nf90_put_var(ncFileID, VarID, ZC ), &
                 'nc_write_model_atts', 'z_c put_var '//trim(filename))
 
+   call nc_check(NF90_inq_varid(ncFileID, 'mask_rho', VarID), &
+                 'nc_write_model_atts', 'mask_rho inq_varid '//trim(filename))
+   call nc_check(nf90_put_var(ncFileID, VarID, mask_rho ), &
+                'nc_write_model_atts', 'mask_rho put_var '//trim(filename))
+
+   call nc_check(NF90_inq_varid(ncFileID, 'mask_psi', VarID), &
+                 'nc_write_model_atts', 'mask_psi inq_varid '//trim(filename))
+   call nc_check(nf90_put_var(ncFileID, VarID, mask_psi ), &
+                'nc_write_model_atts', 'mask_psi put_var '//trim(filename))
+
+   call nc_check(NF90_inq_varid(ncFileID, 'mask_u', VarID), &
+                 'nc_write_model_atts', 'mask_u inq_varid '//trim(filename))
+   call nc_check(nf90_put_var(ncFileID, VarID, mask_u ), &
+                'nc_write_model_atts', 'mask_u put_var '//trim(filename))
+
+   call nc_check(NF90_inq_varid(ncFileID, 'mask_v', VarID), &
+                 'nc_write_model_atts', 'mask_v inq_varid '//trim(filename))
+   call nc_check(nf90_put_var(ncFileID, VarID, mask_v ), &
+                'nc_write_model_atts', 'mask_v put_var '//trim(filename))
+
    endif
 
 ! Flush the buffer and leave netCDF file open
@@ -1409,6 +1499,7 @@
 ! this only executes the first time since counter
 ! gets incremented after the first use and the value
 ! is saved between calls.
+
 if (counter == 1) counter = counter + (my_task_id() * 1000)
 
 call init_random_seq(random_seq, counter)
@@ -1449,17 +1540,17 @@
 !> @param gc precomputed 'get_close_type' to speed up candidate selection 
 !> @param base_obs_loc location of the observation in question
 !> @param base_obs_kind DART KIND of observation in question
-!> @param obs_loc array of comparison locations
-!> @param obs_kind matching array of KINDs for the comparison locations
-!> @param num_close the number of obs_loc locations that are within the prespecified distance (information contained in 'gc')
-!> @param close_ind the indices of the obs_loc locations that are 'close'
+!> @param locs array of comparison locations
+!> @param loc_kind matching array of KINDs for the comparison locations
+!> @param num_close the number of locs locations that are within the prespecified distance (information contained in 'gc')
+!> @param close_ind the indices of the locs locations that are 'close'
 !> @param dist the distances of each of the close locations.
 !>
 
 subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, &
-                         obs_loc, obs_kind, num_close, close_ind, dist)
+                         locs, loc_kind, num_close, close_ind, dist)
 
-! Note that both base_obs_loc and obs_loc are intent(inout), meaning that these
+! Note that both base_obs_loc and locs are intent(inout), meaning that these
 ! locations are possibly modified here and returned as such to the calling routine.
 ! The calling routine is always filter_assim and these arrays are local arrays
 ! within filter_assim. In other words, these modifications will only matter within
@@ -1468,11 +1559,11 @@
 type(get_close_type), intent(in)    :: gc
 type(location_type),  intent(inout) :: base_obs_loc
 integer,              intent(in)    :: base_obs_kind
-type(location_type),  intent(inout) :: obs_loc(:)
-integer,              intent(in)    :: obs_kind(:)
+type(location_type),  intent(inout) :: locs(:)
+integer,              intent(in)    :: loc_kind(:)
 integer,              intent(out)   :: num_close
 integer,              intent(out)   :: close_ind(:)
-real(r8),             intent(out)   :: dist(:)
+real(r8), OPTIONAL,   intent(out)   :: dist(:)
 
 integer                :: ztypeout
 integer                :: t_ind, istatus1, istatus2, k
@@ -1512,18 +1603,18 @@
 endif
 
 if (istatus1 == 0) then
-   call loc_get_close_obs(gc, base_obs_loc, base_obs_kind, obs_loc, obs_kind, &
+   call loc_get_close_obs(gc, base_obs_loc, base_obs_kind, locs, loc_kind, &
                           num_close, close_ind)
 
     do k = 1, num_close
 
       t_ind = close_ind(k)
-      local_obs_loc   = obs_loc(t_ind)
+      local_obs_loc   = locs(t_ind)
       local_obs_which = nint(query_location(local_obs_loc))
 
       if (.not. horiz_dist_only) then
           if (local_obs_which /= vert_localization_coord) then
-              call vert_convert(ens_mean, local_obs_loc, obs_kind(t_ind), ztypeout, istatus2)
+              call vert_convert(ens_mean, local_obs_loc, loc_kind(t_ind), ztypeout, istatus2)
           else
               istatus2 = 0
           endif
@@ -1537,8 +1628,8 @@
             dist(k) = 1.0e9_r8
       else
 
-       ! if (get_obs_kind_var_type(base_obs_kind) == obs_kind(t_ind)) then
-           dist(k) = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, obs_kind(t_ind))
+       ! if (get_obs_kind_var_type(base_obs_kind) == loc_kind(t_ind)) then
+           dist(k) = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, loc_kind(t_ind))
        ! else
        !    dist(k) = 1.0e9_r8
        ! endif
@@ -2093,46 +2184,30 @@
 !> By reading the dimensions first, we can use them in variable
 !> declarations later - which is faster than using allocatable arrays.
 !>
-!> @todo FIXME Since every stagger dimension is specified, we should explicitly 
-!> use them instead of trying to add or subtract 1 from some nominal value 
-!>
 
 subroutine get_grid_dimensions()
 
 integer  :: ncid
 
-integer :: xi_rho
-integer :: xi_u
-integer :: xi_v
-integer :: xi_psi
-integer :: eta_rho
-integer :: eta_u
-integer :: eta_v
-integer :: eta_psi
-integer :: xi_vert
-integer :: eta_vert
-integer :: s_rho
-integer :: s_w
-
 call nc_check(nf90_open(trim(grid_definition_filename), nf90_nowrite, ncid), &
               'get_grid_dimensions', 'open '//trim(grid_definition_filename))
 
-xi_rho   = get_dimension_length(ncid, 'xi_rho',   grid_definition_filename)
-xi_u     = get_dimension_length(ncid, 'xi_u',     grid_definition_filename)
-xi_v     = get_dimension_length(ncid, 'xi_v',     grid_definition_filename)
-xi_psi   = get_dimension_length(ncid, 'xi_psi',   grid_definition_filename)
-eta_rho  = get_dimension_length(ncid, 'eta_rho',  grid_definition_filename)
-eta_u    = get_dimension_length(ncid, 'eta_u',    grid_definition_filename)
-eta_v    = get_dimension_length(ncid, 'eta_v',    grid_definition_filename)
-eta_psi  = get_dimension_length(ncid, 'eta_psi',  grid_definition_filename)
-xi_vert  = get_dimension_length(ncid, 'xi_vert',  grid_definition_filename)
-eta_vert = get_dimension_length(ncid, 'eta_vert', grid_definition_filename)
-s_rho    = get_dimension_length(ncid, 's_rho',    grid_definition_filename)
-s_w      = get_dimension_length(ncid, 's_w',      grid_definition_filename)
+Nxi_rho   = get_dimension_length(ncid, 'xi_rho',   grid_definition_filename)
+Nxi_u     = get_dimension_length(ncid, 'xi_u',     grid_definition_filename)
+Nxi_v     = get_dimension_length(ncid, 'xi_v',     grid_definition_filename)
+Nxi_psi   = get_dimension_length(ncid, 'xi_psi',   grid_definition_filename)
+Neta_rho  = get_dimension_length(ncid, 'eta_rho',  grid_definition_filename)
+Neta_u    = get_dimension_length(ncid, 'eta_u',    grid_definition_filename)
+Neta_v    = get_dimension_length(ncid, 'eta_v',    grid_definition_filename)
+Neta_psi  = get_dimension_length(ncid, 'eta_psi',  grid_definition_filename)
+Nxi_vert  = get_dimension_length(ncid, 'xi_vert',  grid_definition_filename)
+Neta_vert = get_dimension_length(ncid, 'eta_vert', grid_definition_filename)
+Ns_rho    = get_dimension_length(ncid, 's_rho',    grid_definition_filename)
+Ns_w      = get_dimension_length(ncid, 's_w',      grid_definition_filename)
 
-Nx = xi_rho   ! Setting the 'global' variables ... instead of namelist
-Ny = eta_rho  ! Setting the 'global' variables ... instead of namelist
-Nz = s_rho    ! Setting the 'global' variables ... instead of namelist
+Nx =  Nxi_rho  ! Setting the nominal value of the 'global' variables
+Ny = Neta_rho  ! Setting the nominal value of the 'global' variables
+Nz =   Ns_rho  ! Setting the nominal value of the 'global' variables
 
 call nc_check(nf90_close(ncid), &
               'get_grid_dimensions','close '//trim(grid_definition_filename))
@@ -2152,22 +2227,30 @@
 
 subroutine get_grid()
 
-integer  :: k,ncid, VarID,stat,i,j
+integer  :: k,ncid, VarID, stat, i, j
 
 ! The following 'automatic' arrays are more efficient than allocatable arrays.
 ! This is, in part, why the grid dimensions were determined previously.
 
 real(r8) :: dzt0(Nx,Ny)
-real(r8) :: s_rho(Nz),Cs_r(Nz),SSH(Nx,Ny)
-real(r8) :: x_rho(Nx,Ny), &
-            y_rho(Nx,Ny), &
-            x_u(Nx-1,Ny), &
-            y_u(Nx-1,Ny), &
-            x_v(Nx,Ny-1), &
-            y_v(Nx,Ny-1)
+real(r8) :: s_rho(Ns_rho), Cs_r(Ns_rho), SSH(Nx,Ny)
 
+!>@ todo FIXME these aren't needed all the time.
+!>       Should they be allocatable instead of automatic?
+real(r8) :: x_rho(Nxi_rho,Neta_rho), &
+            y_rho(Nxi_rho,Neta_rho), &
+            x_u(Nxi_u,Neta_u), &
+            y_u(Nxi_u,Neta_u), &
+            x_v(Nxi_v,Neta_v), &
+            y_v(Nxi_v,Neta_v)
+
 logical :: nolatlon = .false.
 
+logical :: mymask(Nx,Ny,Nz)
+
+real(r8), allocatable :: datmat(:,:)
+real(r8), parameter :: all_land = 0.001_r8
+
 if (debug > 1) then
    write(string1,*)'..  NX is ',Nx
    write(string2,*)'NY is ',Ny
@@ -2180,17 +2263,21 @@
 call nc_check(nf90_open(trim(grid_definition_filename), nf90_nowrite, ncid), &
       'get_grid', 'open '//trim(grid_definition_filename))
 
-if (.not. allocated(ULAT)) allocate(ULAT(Nx-1,Ny))
-if (.not. allocated(ULON)) allocate(ULON(Nx-1,Ny))
-if (.not. allocated(TLAT)) allocate(TLAT(Nx,Ny))
-if (.not. allocated(TLON)) allocate(TLON(Nx,Ny))
-if (.not. allocated(VLAT)) allocate(VLAT(Nx,Ny-1))
-if (.not. allocated(VLON)) allocate(VLON(Nx,Ny-1))
-if (.not. allocated(HT))   allocate(  HT(Nx,Ny))
+if (.not. allocated(ULAT)) allocate(ULAT(Nxi_u,Neta_u))
+if (.not. allocated(ULON)) allocate(ULON(Nxi_u,Neta_u))
+if (.not. allocated(VLAT)) allocate(VLAT(Nxi_v,Neta_v))
+if (.not. allocated(VLON)) allocate(VLON(Nxi_v,Neta_v))
+if (.not. allocated(TLAT)) allocate(TLAT(Nxi_rho,Neta_rho))
+if (.not. allocated(TLON)) allocate(TLON(Nxi_rho,Neta_rho))
+if (.not. allocated(HT))   allocate(  HT(Nxi_rho,Neta_rho))
+if (.not. allocated(ANGL)) allocate(ANGL(Nxi_rho,Neta_rho))
+
+if (.not. allocated(ZC))   allocate(  ZC(Nx,Ny,Nz))
+
+!> todo are PM, PN used for anything?
+
 if (.not. allocated(PM))   allocate(  PM(Nx,Ny))
 if (.not. allocated(PN))   allocate(  PN(Nx,Ny))
-if (.not. allocated(ANGL)) allocate(ANGL(Nx,Ny))
-if (.not. allocated(ZC))   allocate(  ZC(Nx,Ny,Nz))
 
 ! Read the variables
 
@@ -2312,6 +2399,62 @@
 call nc_check(nf90_get_var( ncid, VarID, Cs_r), &
       'get_grid', 'get_var Cs_r '//trim(grid_definition_filename))
 
+! Allocate the masks and set them all to water.
+! 1 is water, 0 is land (apparently) There are seemingly no intermediate values.
+! In the ROMS netCDF file, these are typed 'double' (overkill), 
+! but we are implementing as short integer.
+! (Could even do logicals if netCDF (output) supported them.)
+
+if (.not. allocated(mask_rho)) allocate(mask_rho(Nxi_rho, Neta_rho))
+if (.not. allocated(mask_psi)) allocate(mask_psi(Nxi_psi, Neta_psi))
+if (.not. allocated(mask_u))   allocate(mask_u(  Nxi_u  , Neta_u  ))
+if (.not. allocated(mask_v))   allocate(mask_v(  Nxi_v  , Neta_v  ))
+
+mask_rho = 0
+mask_psi = 0
+mask_u   = 0
+mask_v   = 0
+
+! Read mask on RHO-points
+
+allocate(datmat(Nxi_rho, Neta_rho))
+call nc_check(nf90_inq_varid(ncid, 'mask_rho', VarID), &
+      'get_grid', 'inq_varid mask_rho '//trim(grid_definition_filename))
+call nc_check(nf90_get_var( ncid, VarID, datmat), &
+      'get_grid', 'get_var mask_rho '//trim(grid_definition_filename))
+where(datmat > all_land) mask_rho = 1
+deallocate(datmat)
+
+! Read mask on PSI-points
+
+allocate(datmat(Nxi_psi, Neta_psi))
+call nc_check(nf90_inq_varid(ncid, 'mask_psi', VarID), &
+      'get_grid', 'inq_varid mask_psi '//trim(grid_definition_filename))
+call nc_check(nf90_get_var( ncid, VarID, datmat), &
+      'get_grid', 'get_var mask_psi '//trim(grid_definition_filename))
+where(datmat > all_land) mask_psi = 1
+deallocate(datmat)
+
+! Read mask on U-points
+
+allocate(datmat(Nxi_u, Neta_u))
+call nc_check(nf90_inq_varid(ncid, 'mask_u', VarID), &
+      'get_grid', 'inq_varid mask_u '//trim(grid_definition_filename))
+call nc_check(nf90_get_var( ncid, VarID, datmat), &
+      'get_grid', 'get_var mask_u '//trim(grid_definition_filename))
+where(datmat > all_land) mask_u = 1
+deallocate(datmat)
+
+! Read mask on V-points
+
+allocate(datmat(Nxi_v, Neta_v))
+call nc_check(nf90_inq_varid(ncid, 'mask_v', VarID), &
+      'get_grid', 'inq_varid mask_v '//trim(grid_definition_filename))
+call nc_check(nf90_get_var( ncid, VarID, datmat), &
+      'get_grid', 'get_var mask_v '//trim(grid_definition_filename))
+where(datmat > all_land) mask_v = 1
+deallocate(datmat)
+
 call nc_check(nf90_close(ncid), &
              'get_var','close '//trim(grid_definition_filename))
 
@@ -2347,6 +2490,25 @@
      ZC(:,:,k) = -( dzt0(:,:) + SSH(:,:)*(1.0_r8 + dzt0(:,:)/HT(:,:)) )
 end do
 
+if (do_output() .and. debug > 3) then
+    write(string1,*)'    min/max ULON ',minval(ULON), maxval(ULON)
+    write(string2,*)'min/max ULAT ',minval(ULAT), maxval(ULAT)
+    call error_handler(E_MSG,'get_grid',string1, text2=string2)
+
+    write(string1,*)'    min/max VLON ',minval(VLON), maxval(VLON)
+    write(string2,*)'min/max VLAT ',minval(VLAT), maxval(VLAT)
+    call error_handler(E_MSG,'get_grid',string1, text2=string2)
+
+    write(string1,*)'    min/max TLON ',minval(TLON), maxval(TLON)
+    write(string2,*)'min/max TLAT ',minval(TLAT), maxval(TLAT)
+    call error_handler(E_MSG,'get_grid',string1, text2=string2)
+
+    mymask = ZC < 1.0E30 .and. ZC > -1.0E30
+    write(string1,*)'    min/max   ZC ',minval(ZC,mymask), maxval(ZC,mymask)
+    call error_handler(E_MSG,'get_grid',string1)
+
+endif
+
 end subroutine get_grid
 
 
@@ -3475,10 +3637,10 @@
 integer  :: lat_bot, lat_top, lon_bot, lon_top
 integer  :: x_ind, y_ind,ivar
 real(r8) :: p(4), x_corners(4), y_corners(4)
-logical  :: masked
 
 ! Succesful return has istatus of 0
-istatus = 0
+istatus    = 99
+interp_val = MISSING_R8
 
 ivar = get_progvar_index_from_kind(var_type)
 
@@ -3518,32 +3680,25 @@
 
 ! Get the values at the four corners of the box or quad
 ! Corners go around counterclockwise from lower left
-p(1) = get_val(lon_bot, lat_bot, height, x, var_type)
-if(masked) then
-   istatus = 3
-   return
-endif
+! If any one of these fail, go no furhter.
 
-p(2) = get_val(lon_top, lat_bot, height, x, var_type)
-if(masked) then
-   istatus = 3
-   return
-endif
+istatus = 3
+   p(1) = get_val(lon_bot, lat_bot, height, x, var_type)
+if(p(1) == MISSING_R8) return
 
-p(3) = get_val(lon_top, lat_top, height, x, var_type)
-if(masked) then
-   istatus = 3
-   return
-endif
+   p(2) = get_val(lon_top, lat_bot, height, x, var_type)
+if(p(2) == MISSING_R8) return
 
-p(4) = get_val(lon_bot, lat_top, height, x, var_type)
-if(masked) then
-   istatus = 3
-   return
-endif
+   p(3) = get_val(lon_top, lat_top, height, x, var_type)
+if(p(3) == MISSING_R8) return
 
+   p(4) = get_val(lon_bot, lat_top, height, x, var_type)
+if(p(4) == MISSING_R8) return
+
 ! Full bilinear interpolation for quads
+! istatus = 0 is good, whatever the interp_val is, it is ...
 
+istatus = 0
 call quad_bilinear_interp(lon, lat, x_corners, y_corners, p, interp_val)
 
 end subroutine lon_lat_interpolate
@@ -3578,8 +3733,11 @@
 integer    :: i,iidd
 real(r8)   :: tp(Nz), zz(Nz)
 
-tp=0.0
-iidd=0
+tp   = 0.0_r8
+iidd = 0
+
+!>@ todo FIXME ... check ... get_val() can return MISSING_R8 ... is this possible
+
 do i=1,Nz
     tp(i)=get_val(lonid, latid, i, x, var_type)
     zz(i)=ZC(lonid,latid,i)
@@ -3613,14 +3771,14 @@
 !-----------------------------------------------------------------------
 !>
 !> Returns the DART state value for a given lat, lon, and level index.
-!> 'get_val' will be a MISSNG value if this is NOT a valid grid location
+!> 'get_val' will be a MISSING value if this is NOT a valid grid location
 !> (e.g. land, or below the ocean floor).
 !>
 !> @param get_val the value of the DART state at the gridcell of interest.
 !> @param lon_index the index of the longitude of interest.
 !> @param lat_index the index of the latitude of interest.
 !> @param level_index the index of the level of interest.
-!> @param x the DART state vector.
+!> @param x the (contiguous) portion of the DART state vector for the variable of interest
 !> @param var_type the DART KIND of interest.
 !>
 !> @ todo FIXME Johnny may have a better way to do this if everything stays
@@ -3635,25 +3793,33 @@
 real(r8), intent(in) :: x(:)
 real(r8)             :: get_val
 
-integer  ::  ii,jj,kk,tt,ivar
+integer :: tt     ! the (absolute) index into the DART state vector
+integer :: ivar
+integer :: Ndim1, Ndim2, Ndim3
 
+get_val = MISSING_R8 ! guilty until proven otherwise
+
 ivar = get_progvar_index_from_kind(var_type)
 
-tt=0
-do ii=1,progvar(ivar)%numvertical
-   do jj=1,progvar(ivar)%numeta
-      do kk=1,progvar(ivar)%numxi
-         tt=tt+1
-         if( (ii==level_index) .and. (jj==lat_index) .and. (kk==lon_index) ) then
-            get_val=x(tt)
-            return
-         else
-            get_val=MISSING_R8
-         endif
-      enddo
-   enddo
-enddo
+!> @ todo FIXME Implement the land masking. Must determine which mask
+!>        is appropriate for this variable ... extend progvar?
 
+Ndim3 = progvar(ivar)%numvertical
+Ndim2 = progvar(ivar)%numeta
+Ndim1 = progvar(ivar)%numxi
+
+if ( (  lon_index < 1 .or.   lon_index > Ndim1) .or. &
+     (  lat_index < 1 .or.   lat_index > Ndim2) .or. &
+     (level_index < 1 .or. level_index > Ndim3) ) return
+
+! implicit assumption on packing order into the DART vector
+
+tt = (level_index - 1) * Ndim1 * Ndim2 + &
+     (  lat_index - 1) * Ndim1 + &
+        lon_index
+
+if(tt > 0 .and. tt <= size(x)) get_val = x(tt)
+
 end function get_val
 
 
@@ -3682,21 +3848,26 @@
 integer, intent(out) :: z_index
 integer, intent(in)  :: var_type
 
-integer :: ivar, Nx, Ny
+integer :: ivar, numxi, numeta
 
-ivar = get_progvar_index_from_kind(var_type)
-Nx = progvar(ivar)%numxi
-Ny = progvar(ivar)%numeta
+ivar   = get_progvar_index_from_kind(var_type)
+numxi  = progvar(ivar)%numxi
+numeta = progvar(ivar)%numeta
 
 if (progvar(ivar)%kind_string=='KIND_SEA_SURFACE_HEIGHT') then
   z_index = 1
 else
-  z_index = ceiling( float(offset) / (Nx * Ny))
+  z_index = ceiling( float(offset) / (numxi * numeta))
 endif
-y_index = ceiling( float(offset - ((z_index-1)*Nx*Ny)) / Nx)
-x_index =  offset - ((z_index-1)*Nx*Ny) - ((y_index-1)*Nx)
+y_index = ceiling( float(offset - ((z_index-1)*numxi*numeta)) / numxi )
+x_index =  offset - ((z_index-1)*numxi*numeta) - ((y_index-1)*numxi)
 
-if (debug > 3) print *, 'lon, lat, depth index = ', x_index, y_index, z_index
+if (do_output() .and. debug > 3) then
+   write(string1,*) 'checking ',trim(progvar(ivar)%varname), ' offset = ', offset
+   write(string2,*) 'lon, lat, depth index = ', x_index, y_index, z_index
+   call error_handler(E_MSG,'get_state_indices:', string1, &
+              source, revision, revdate, text2=string2)
+endif
 
 end subroutine get_state_indices
 
@@ -3876,6 +4047,17 @@
 !>
 !> @todo FIXME should this really error out or silently fail?
 !> At this point it may silently fail.
+!>
+!> ROMS has a structured horizontal grid.
+!>
+!> We group the structured grid cells into regular boxes, 
+!> e.g., 70 x 40 boxes. 
+!>
+!> (1)calculate distances from a given point to the center cell of each box
+!> (2)after we find the box which is closest to the given point,
+!>    we can calculate distances to the cells within the box.
+!> bjchoi 2014/08/07
+!-------------------------------------------------------------
 
 subroutine get_reg_box_indices(lon, lat, LON_al, LAT_al, x_ind, y_ind)
 
@@ -3886,10 +4068,203 @@
 integer,  intent(out) :: x_ind
 integer,  intent(out) :: y_ind
 
+real(r8),allocatable  :: boxLON(:,:), boxLAT(:,:)
+integer               :: i, j, ii, jj, Nx, Ny, nxbox, nybox
+integer               :: num_cell, loc_box, i_start,i_end,j_start,j_end
+real(r8)              :: lon_dif, lat_dif, tp1, tp2
+real(r8),allocatable  :: rtemp(:)
+integer, allocatable  :: itemp(:),jtemp(:)
+integer               :: t, loc1 
+integer               :: min_location(1) ! stupid parser
+
+! Divide the model grid cells into a group of boxes.
+! For example, divide the model domain into 47 x 70 subdomains (or boxes).
+! Each box contains num_cell * num_cell grid cells.
+num_cell = 20 ! please use even number,e.g., 4,6,8,10,12..., for num_cell  
+Nx = size(LON_al,1)
+Ny = size(LON_al,2)
+nxbox = ceiling( float(Nx) / num_cell )
+nybox = ceiling( float(Ny) / num_cell )
+
+! we define the center cell of each box:
+! boxLON and boxLAT
+allocate( boxLON(nxbox,nybox) )
+allocate( boxLAT(nxbox,nybox) )
+allocate( rtemp(Nx*Ny) )
+allocate( itemp(Nx*Ny) )
+allocate( jtemp(Nx*Ny) )
+
+do ii = 1, nxbox-1 
+  do jj = 1, nybox-1
+     i = (num_cell/2) + (ii-1)*num_cell
+     j = (num_cell/2) + (jj-1)*num_cell
+     boxLON(ii,jj) = LON_al(i,j)
+     boxLAT(ii,jj) = LAT_al(i,j)
+  enddo 
+  jj = nybox
+  i = (num_cell/2) + (ii-1)*num_cell
+  j = (num_cell/2) + (jj-1)*num_cell
+  if (j .GT. Ny) j = Ny
+  boxLON(ii,jj) = LON_al(i,j)
+  boxLAT(ii,jj) = LAT_al(i,j)
+enddo 
+
+ii = nxbox 
+i = (num_cell/2) + (ii-1)*num_cell
+if (i .GT. Nx) i = Nx
+do jj = 1, nybox-1
+   j = (num_cell/2) + (jj-1)*num_cell
+   boxLON(ii,jj) = LON_al(i,j)
+   boxLAT(ii,jj) = LAT_al(i,j)
+enddo 
+jj = nybox
+j = (num_cell/2) + (jj-1)*num_cell
+if (j .GT. Ny) j = Ny
+
+! find the box which the obs data belongs to
+t=0
+do ii=1,nxbox
+ do jj=1,nybox
+    t=t+1
+    lon_dif = (lon - boxLON(ii,jj))*111.41288*cos(lat*DEG2RAD) - &
+                 0.09350 * cos(3.0 * lat*DEG2RAD) + &
+                 0.00012 * cos(5.0 * lat*DEG2RAD)
+    lat_dif = (lat - boxLat(ii,jj))*111.13295
+    rtemp(t) = sqrt(lon_dif*lon_dif+lat_dif*lat_dif)
+    itemp(t) = ii
+    jtemp(t) = jj
+ enddo
+enddo
+
+! find minium distance
+min_location = minloc(rtemp)
+loc_box = min_location(1)
+
+! locate the box which contains the obs data (ii, jj)
+! set the ranges of i and j around the box.
+! the neighboring 4 model grid points are within this range.
+ii = itemp(loc_box)
+jj = jtemp(loc_box)
+i_start = max((ii-1)*num_cell - num_cell/2, 1) ! span 2 cells
+i_end   = min(    ii*num_cell + num_cell/2,Nx) ! span 2 cells
+j_start = max((jj-1)*num_cell - num_cell/2, 1)  
+j_end   = min(    jj*num_cell + num_cell/2,Ny)
+
+! now, calculate the distance between the obs and close candiate points
+t=0
+do i=i_start,i_end 
+ do j=j_start,j_end
+    t=t+1
+    lon_dif = (lon - LON_al(i,j))*111.41288*cos(lat*DEG2RAD) - &
+              0.09350 * cos(3.0 * lat*DEG2RAD) + &
+              0.00012 * cos(5.0 * lat*DEG2RAD)
+    lat_dif = (lat - Lat_al(i,j))*111.13295
+    rtemp(t) = sqrt(lon_dif*lon_dif+lat_dif*lat_dif)
+    itemp(t) = i
+    jtemp(t) = j
+  enddo
+enddo
+
+!find minimum distances
+min_location = minloc(rtemp)
+loc1 = min_location(1)
+
+!> @todo FIXME dies here with bounds checking turned on and you attempt
+!>       to interpolate something outside the domain. This WHOLE ROUTINE
+!>       can be replaced with something similar from POP.
+
+!check point location
+tp1 = checkpoint(LON_al(itemp(loc1)  ,jtemp(loc1)  ), &
+                 LAT_al(itemp(loc1)  ,jtemp(loc1)  ), &
+                 LON_al(itemp(loc1)+1,jtemp(loc1)+1), &
+                 LAT_al(itemp(loc1)+1,jtemp(loc1)+1),lon,lat)
+
+tp2 = checkpoint(LON_al(itemp(loc1)  ,jtemp(loc1)  ), &
+                 LAT_al(itemp(loc1)  ,jtemp(loc1)  ), &
+                 LON_al(itemp(loc1)  ,jtemp(loc1)+1), &
+                 LAT_al(itemp(loc1)  ,jtemp(loc1)+1),lon,lat)
+
+if(     (tp1  > 0.0_r8) .and. (tp2  > 0.0_r8)) then
+      x_ind=itemp(loc1)-1
+      y_ind=jtemp(loc1)
+elseif( (tp1  > 0.0_r8) .and. (tp2  < 0.0_r8)) then
+      x_ind=itemp(loc1)
+      y_ind=jtemp(loc1)
+elseif( (tp1  < 0.0_r8) .and. (tp2  > 0.0_r8)) then
+      x_ind=itemp(loc1)-1
+      y_ind=jtemp(loc1)-1
+elseif( (tp1  < 0.0_r8) .and. (tp2  < 0.0_r8)) then
+      x_ind=itemp(loc1)
+      y_ind=jtemp(loc1)-1
+elseif( (tp1 == 0.0_r8) .and. (tp2  > 0.0_r8)) then
+      x_ind=itemp(loc1)-1
+      y_ind=jtemp(loc1)
+elseif( (tp1 == 0.0_r8) .and. (tp2  < 0.0_r8)) then
+      x_ind=itemp(loc1)
+      y_ind=jtemp(loc1)
+elseif( (tp1  > 0.0_r8) .and. (tp2 == 0.0_r8)) then
+      x_ind=itemp(loc1)
+      y_ind=jtemp(loc1)
+elseif( (tp1  < 0.0_r8) .and. (tp2 == 0.0_r8)) then
+      x_ind=itemp(loc1)-1
+      y_ind=jtemp(loc1)-1
+else
+   write(string1,*)'ERROR in finding matching grid point for'
+   write(string2,*)'longitude ',lon
+   write(string3,*)'latitude  ',lat
+   call error_handler(E_MSG,'get_reg_box_indices:', string1, &
+              source, revision, revdate, text2=string2, text3=string3)
+endif
+
+if (do_output() .and. debug > 0) then
+
+   write(string1,*)'checking' 
+   write(string2,*)'longitude ',lon,' is index ',x_ind
+   write(string3,*)'latitude  ',lat,' is index ',y_ind
+   call error_handler(E_MSG,'get_reg_box_indices:', string1, &
+              source, revision, revdate, text2=string2, text3=string3)
+
+endif
+
+deallocate(boxLON, boxLAT)
+deallocate(rtemp)
+deallocate(itemp)
+deallocate(jtemp)
+
+end subroutine get_reg_box_indices
+                                        
+!-----------------------------------------------------------------------
+!>
+!> Given a longitude and latitude in degrees returns the index of the
+!> regular lon-lat box that contains the point. TJH ... this was the original
+!> version that dies when trying to interpolate something outside the domain.
+!> It dies with a run-time error when you have bounds-checking turned on.
+!>
+!> @param lon the longitude of interest
+!> @param lat the latitude of interest
+!> @param LON_al the matrix of longitudes defining the boxes (grids).
+!> @param LAT_al the matrix of latitudes defining the boxes (grids).
+!> @param x_ind the 'x' (longitude) index of the location of interest.
+!> @param y_ind the 'y' (latitude) index of the location of interest.
+!>
+!> @todo FIXME should this really error out or silently fail?
+!> At this point it may silently fail.
+!>
+
+subroutine get_reg_box_indices_org(lon, lat, LON_al, LAT_al, x_ind, y_ind)
+
+real(r8), intent(in)  :: lon
+real(r8), intent(in)  :: lat
+real(r8), intent(in)  :: LON_al(:,:)
+real(r8), intent(in)  :: LAT_al(:,:)
+integer,  intent(out) :: x_ind
+integer,  intent(out) :: y_ind
+
 real(r8)              :: lon_dif,lat_dif,tp1,tp2
 real(r8), allocatable :: rtemp(:)
 integer,  allocatable :: itemp(:),jtemp(:)
 integer               :: i,j,t,loc1
+integer               :: min_location(1) ! stupid parser
 
 ! tricky for rotated domain
 
@@ -3916,7 +4291,8 @@
 
 !find minimum distances
 
-loc1 = FindMinimum(rtemp, 1, t)
+min_location = minloc(rtemp)
+loc1 = min_location(1)
 
 !check point location
 
@@ -3966,39 +4342,11 @@
 deallocate(itemp)
 deallocate(jtemp)
 
-end subroutine get_reg_box_indices
+end subroutine get_reg_box_indices_org
 
 
 !-----------------------------------------------------------------------
 !>
-!> Find the minimum value 
-!> @todo FIXME .... uh ... totally replace with intrinsic 'minval' routine
-!> y = minval(x(iStart:iEnd))
-
-INTEGER FUNCTION FindMinimum(x, iStart, iEnd)
-IMPLICIT  NONE
-real(r8), INTENT(IN) :: x(1:)
-INTEGER,  INTENT(IN) :: iStart, iEnd
-
-INTEGER                            :: Minimum
-INTEGER                            :: Location
-INTEGER                            :: i
-
-Minimum  = x(iStart)              ! assume the first is the min
-Location = iStart                 ! record its position
-DO i = iStart+1, iEnd             ! start with next elements
-   IF (x(i) < Minimum) THEN       !   if x(i) less than the min?
-      Minimum  = x(i)             !      Yes, a new minimum found
-      Location = i                !      record its position
-   END IF
-END DO
-FindMinimum = Location            ! return the position
-
-END FUNCTION  FindMinimum
-
-
-!-----------------------------------------------------------------------
-!>
 !> check (xp,yp) with line (x1,y1)-->(x2,y2)
 !> if checkpoint>0, (xp,yp) is on the left side of the line
 !> if checkpoint<0, (xp,yp) is on the right side of the line


More information about the Dart-dev mailing list