[Dart-dev] [4470] DART/trunk/models/NCOMMAS/model_mod.f90: Has Lou' s model_interpolate code, all the latest/greatest from Ted,

nancy at ucar.edu nancy at ucar.edu
Fri Aug 6 16:10:32 MDT 2010


Revision: 4470
Author:   thoar
Date:     2010-08-06 16:10:32 -0600 (Fri, 06 Aug 2010)
Log Message:
-----------
Has Lou's model_interpolate code, all the latest/greatest from Ted,
all the POP stuff pruned out - ncommas_to_dart and dart_to_ncommas run
without crashing - might even handle restart files with a TIME coordinate
that is not a singleton ... 

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

-------------- next part --------------
Modified: DART/trunk/models/NCOMMAS/model_mod.f90
===================================================================
--- DART/trunk/models/NCOMMAS/model_mod.f90	2010-08-05 22:35:44 UTC (rev 4469)
+++ DART/trunk/models/NCOMMAS/model_mod.f90	2010-08-06 22:10:32 UTC (rev 4470)
@@ -13,7 +13,8 @@
 ! This is the interface between the ncommas model and DART.
 
 ! Modules that are absolutely required for use are listed
-use        types_mod, only : r4, r8, SECPERDAY, MISSING_R8, rad2deg, PI
+use        types_mod, only : r4, r8, digits12, SECPERDAY, MISSING_R8,          &
+                             rad2deg, deg2rad, PI
 use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
                              print_time, print_date, set_calendar_type,        &
                              operator(*),  operator(+), operator(-),           &
@@ -74,10 +75,16 @@
 
 ! generally useful routines for various support purposes.
 ! the interfaces here can be changed as appropriate.
-public :: get_gridsize, restart_file_to_sv, sv_to_restart_file, &
-          get_ncommas_restart_filename, get_base_time, get_state_time
 
+public :: get_gridsize,                 &
+          restart_file_to_sv,           &
+          sv_to_restart_file,           &
+          get_ncommas_restart_filename, &
+          get_base_time,                &
+          get_state_time
+
 ! version controlled file description for error handling, do not edit
+
 character(len=128), parameter :: &
    source   = '$URL$', &
    revision = '$Revision$', &
@@ -87,14 +94,16 @@
 logical, save :: module_initialized = .false.
 
 ! Storage for a random sequence for perturbing a single initial state
+
 type(random_seq_type) :: random_seq
 
 ! things which can/should be in the model_nml
-integer  :: assimilation_period_days = 1
-integer  :: assimilation_period_seconds = 0
-real(r8) :: model_perturbation_amplitude = 0.2
-logical  :: output_state_vector = .true.
-integer  :: debug = 0   ! turn up for more and more debug messages
+
+integer            :: assimilation_period_days = 0
+integer            :: assimilation_period_seconds = 60
+real(r8)           :: model_perturbation_amplitude = 0.2
+logical            :: output_state_vector = .true.
+integer            :: debug = 0   ! turn up for more and more debug messages
 character(len=32)  :: calendar
 character(len=256) :: ncommas_restart_filename = 'ncommas_restart.nc'
 
@@ -159,17 +168,29 @@
 
 type(progvartype), dimension(max_state_variables) :: progvar
 
-! Grid parameters - the values will be read from a
-! ncommas restart file. 
+! little stuff here and there..
 
+real(digits12), parameter :: rearth=1000.0_digits12 * 6367.0_digits12 ! radius of earth (m)
+
+! Grid parameters - the values will be read from an ncommas restart file. 
 ! Each spatial dimension has a staggered counterpart.
+
 integer :: nxc=-1, nyc=-1, nzc=-1    ! scalar grid positions
 integer :: nxe=-1, nye=-1, nze=-1    ! staggered grid positions
 
+! Global storage for the grid's reference lat/lon
+
+real(r8) :: ref_lat, ref_lon
+real(r8) :: xg_pos, yg_pos, hgt_offset
+
 ! locations of cell centers (C) and edges (E) for each axis.
+
+real(r8), allocatable :: XC(:), XE(:)
+real(r8), allocatable :: YC(:), YE(:)
 real(r8), allocatable :: ZC(:), ZE(:)
 
 ! These arrays store the longitude and latitude of the lower left corner
+
 real(r8), allocatable :: ULAT(:,:), ULON(:,:)  ! XE,YC
 real(r8), allocatable :: VLAT(:,:), VLON(:,:)  ! XC,YE
 real(r8), allocatable :: WLAT(:,:), WLON(:,:)  ! XC,YC
@@ -181,6 +202,7 @@
 
 ! set this to true if you want to print out the current time
 ! after each N observations are processed, for benchmarking.
+
 logical :: print_timestamps = .false.
 integer :: print_every_Nth  = 10000
 
@@ -207,60 +229,13 @@
       MODULE PROCEDURE get_state_time_fname
 END INTERFACE
 
-!------------------------------------------------
-! These bits are left over from the POP dipole grid.
-!------------------------------------------------
-
-! The regular grid used for dipole interpolation divides the sphere into
-! a set of regularly spaced lon-lat boxes. The number of boxes in
-! longitude and latitude are set by num_reg_x and num_reg_y. Making the
-! number of regular boxes smaller decreases the computation required for
-! doing each interpolation but increases the static storage requirements
-! and the initialization computation (which seems to be pretty small).
-integer, parameter :: num_reg_x = 90, num_reg_y = 90
-
-! The max_reg_list_num controls the size of temporary storage used for
-! initializing the regular grid. Four arrays
-! of size num_reg_x*num_reg_y*max_reg_list_num are needed. The initialization
-! fails and returns an error if max_reg_list_num is too small. A value of
-! 30 is sufficient for the x3 ncommas grid with 180 regular lon and lat boxes 
-! and a value of 80 is sufficient for for the x1 grid.
-integer, parameter :: max_reg_list_num = 80
-
-! The dipole interpolation keeps a list of how many and which dipole quads
-! overlap each regular lon-lat box. The number for the u and t grids are 
-! stored in u_dipole_num and t_dipole_num. The allocatable arrays
-! u_dipole_lon(lat)_list and t_dipole_lon(lat)_list list the longitude 
-! and latitude indices for the overlapping dipole quads. The entry in
-! u_dipole_start and t_dipole_start for a given regular lon-lat box indicates
-! where the list of dipole quads begins in the u_dipole_lon(lat)_list and
-! t_dipole_lon(lat)_list arrays.
-
-integer :: u_dipole_start(num_reg_x, num_reg_y)
-integer :: u_dipole_num  (num_reg_x, num_reg_y) = 0
-integer :: t_dipole_start(num_reg_x, num_reg_y)
-integer :: t_dipole_num  (num_reg_x, num_reg_y) = 0
-integer, allocatable :: u_dipole_lon_list(:), t_dipole_lon_list(:)
-integer, allocatable :: u_dipole_lat_list(:), t_dipole_lat_list(:)
-
-! Need to check for pole quads: for now we are not interpolating in them
-integer :: pole_x, t_pole_y, u_pole_y
-
-
-! Have a global variable saying whether this is dipole or regular lon-lat grid
-! This should be initialized static_init_model. Code to do this is below.
-logical :: dipole_grid
-
 contains
 
-
-
 !==================================================================
 ! All the REQUIRED interfaces come first - just by convention.
 !==================================================================
 
 
-
 function get_model_size()
 !------------------------------------------------------------------
 ! Done - TJH.
@@ -306,44 +281,99 @@
 
 
 
-subroutine get_state_meta_data(index_in, location, var_type)
-!------------------------------------------------------------------
+!############################################################################
 !
-! Given an integer index into the state vector structure, returns the
-! associated location. A second intent(out) optional argument kind
-! can be returned if the model has more than one type of field (for
-! instance temperature and zonal wind component). This interface is
-! required for all filter applications as it is required for computing
-! the distance between observations and state variables.
+!     ##################################################################
+!     ######                                                      ######
+!     ######            SUBROUTINE GET_STATE_META_DATA            ######
+!     ######                                                      ######
+!     ##################################################################
+!
+!     PURPOSE:
+!
+!     Given an integer index into the state vector structure, returns the
+!     associated array indices for lat, lon, and height, as well as the type.
+!
+!############################################################################
+!
+!     Author:  Tim Hoar, Ted Mansell, Lou Wicker
+!
+!     Creation Date:  August 2010
+!     
+!     Variables needed to be stored in the MODEL_MODULE data structure
+!
+!       PROGVARS => Module's data structure
+!
+!############################################################################
 
-integer,             intent(in)  :: index_in
-type(location_type), intent(out) :: location
-integer,             intent(out), optional :: var_type
+subroutine get_state_meta_data(index_in, location, var_type)
 
-real(r8):: lat, lon, height
-integer :: lon_index, lat_index, height_index, local_var
+! Passed variables
 
-call get_state_indices(index_in, lat_index, lon_index, height_index, local_var)
+  integer, intent(in)            :: index_in
+  type(location_type)            :: location
+  integer, optional, intent(out) :: var_type
+  
+! Local variables
 
-if     (is_on_ugrid(local_var)) then
-   lon = ULON(lon_index, lat_index)
-   lat = ULAT(lon_index, lat_index)
-elseif (is_on_vgrid(local_var)) then
-   lon = VLON(lon_index, lat_index)
-   lat = VLAT(lon_index, lat_index)
-else 
-   lon = WLON(lon_index, lat_index)
-   lat = WLAT(lon_index, lat_index)
-endif
+  integer  :: nxp, nyp, nzp, var_index, iloc, jloc, kloc, nf, n
+  integer  :: index, lat_index, lon_index
+  real(r8) :: height
 
-if (debug > 5) print *, 'lon, lat, height = ', lon, lat, height
+  if ( .not. module_initialized ) call static_init_model
+  
+  index = -1
+  nf    = -1
+  
+  DO n = 1,nfields
+    IF( progvar(n)%index1 < index_in <= progvar(n)%indexN ) THEN
+      nf = n
+      index = index_in - progvar(n)%index1
+      EXIT
+    ENDIF
+  ENDDO
+  
+  IF( index == -1 ) THEN
+     write(string1,*) 'Problem, cannot find base_offset, index_in is: ', index_in
+     call error_handler(E_ERR,'get_state_meta_data',string1,source,revision,revdate)
+  ENDIF
+  
+  nxp = progvar(nf)%dimlens(1)
+  nyp = progvar(nf)%dimlens(2)
+  nzp = progvar(nf)%dimlens(3)
 
-location = set_location(lon, lat, height, VERTISHEIGHT)
-
-if (present(var_type)) then
-   var_type = local_var
-endif
-
+  kloc  = 1 + index / (nxp*nyp)
+  index = index - (kloc-1)*nyp*nxp
+  jloc  = 1 + index / nxp
+  index = index - (jloc-1)*nxp
+  iloc  = index
+  
+  lat_index    = jloc
+  lon_index    = iloc
+  height       = zc(kloc)
+  
+  IF (debug > 5) print *, 'lon, lat, height index = ', lon_index, lat_index, kloc
+  
+! Here we assume two things:
+! 1)  the first three arrays in the state variable structure will be U, V, W
+! 2)  everything else in the state variable is at grid-centers
+  
+  IF(progvar(nf)%dart_kind == KIND_U_WIND_COMPONENT) THEN
+    location = set_location(ulon(iloc,jloc), ulat(iloc,jloc), height, VERTISHEIGHT)
+  ELSEIF(progvar(nf)%dart_kind == KIND_V_WIND_COMPONENT) THEN
+    location = set_location(vlon(iloc,jloc), vlat(iloc,jloc), height, VERTISHEIGHT)
+  ELSEIF(progvar(nf)%dart_kind == KIND_VERTICAL_VELOCITY) THEN
+    height   = ze(kloc)
+    location = set_location(wlon(iloc,jloc), wlat(iloc,jloc), height, VERTISHEIGHT)
+  ELSE
+    location = set_location(wlon(iloc,jloc), wlat(iloc,jloc), height, VERTISHEIGHT)
+  ENDIF
+  
+  IF (present(var_type)) THEN
+     var_type = progvar(nf)%dart_kind
+  ENDIF
+  
+return
 end subroutine get_state_meta_data
 
 
@@ -374,6 +404,7 @@
 ! the dart_ncommas_mod module.
 
 ! Local variables - all the important ones have module scope
+
 integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
 character(len=NF90_MAX_NAME)          :: varname
 character(len=paramname_length)       :: kind_string
@@ -382,6 +413,7 @@
 integer :: ss, dd
 integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID
 logical :: shapeok
+integer :: dimtime
 
 if ( module_initialized ) return ! only need to do this once.
 
@@ -431,11 +463,15 @@
 allocate(ULAT(nxe,nyc), ULON(nxe,nyc))
 allocate(VLAT(nxc,nye), VLON(nxc,nye))
 allocate(WLAT(nxc,nyc), WLON(nxc,nyc))
+allocate(  XC(  nxc  ),   XE(  nxe  ))
+allocate(  YC(  nyc  ),   YE(  nye  ))
 allocate(  ZC(  nzc  ),   ZE(  nze  ))
 
-call get_grid(nxc, nxe, nyc, nye, nzc, nze, &
-              ULAT, ULON, VLAT, VLON, WLAT, WLON, ZC, ZE)
-
+CALL GET_GRID(nxc, nxe, nyc, nye, nzc, nze,       &
+              XC, XE, YC, YE, ZC, ZE,             &
+              ULAT, ULON, VLAT, VLON, WLAT, WLON, &
+              xg_pos, yg_pos, ref_lat, ref_lon, hgt_offset)
+              
 !---------------------------------------------------------------
 ! Compile the list of ncommas variables to use in the creation
 ! of the DART state vector. Required to determine model_size.
@@ -503,6 +539,7 @@
    ! need to skip it. 
    varsize = 1
    dimlen  = 1
+   dimtime = 0
    progvar(ivar)%numdims = 0
    DimensionLoop : do i = 1,numdims
 
@@ -530,7 +567,7 @@
       if (trim(progvar(ivar)%storder) == 'xyz3d') shapeok = .true.
    endif
    if ( .not. shapeok ) then
-      write(string1,*)'unable to handle storage order of '//trim(progvar(ivar)%storder)
+      write(string1,*)'unable to handle storage order of '//trim(progvar(ivar)%storder)//' numdims,timedim = ',progvar(ivar)%numdims,dimtime
       call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate, &
                                               text2=string2)
    endif
@@ -586,9 +623,6 @@
 
 allocate( ens_mean(model_size) )
 
-! Initialize the interpolation routines
-call init_interp()
-
 end subroutine static_init_model
 
 
@@ -1383,7 +1417,7 @@
 if (.not. horiz_dist_only) then
 !  if (base_which /= wrf%dom(1)%vert_coord) then
 !     call vert_interpolate(ens_mean, base_obs_loc, base_obs_kind, istatus1)
-!  elseif (base_array(3) == missing_r8) then
+!  elseif (base_array(3) == MISSING_R8) then
 !     istatus1 = 1
 !  endif
 endif
@@ -1418,7 +1452,7 @@
       ! or vert_interpolate returned error (istatus2=1)
       local_obs_array = get_location(local_obs_loc)
       if (( (.not. horiz_dist_only)             .and. &
-            (local_obs_array(3) == missing_r8)) .or.  &
+            (local_obs_array(3) == MISSING_R8)) .or.  &
             (istatus2 == 1)                   ) then
             dist(k) = 1.0e9
       else
@@ -1453,17 +1487,22 @@
 
 
 
-subroutine get_gridsize(num_x, num_y, num_z)
- integer, intent(out) :: num_x, num_y, num_z
+subroutine get_gridsize(num_xc, num_xe, num_yc, num_ye, num_zc, num_ze )
+ integer, intent(out) :: num_xc, num_yc, num_zc
+ integer, intent(out) :: num_xe, num_ye, num_ze
 !------------------------------------------------------------------
 ! public utility routine.
 
 if ( .not. module_initialized ) call static_init_model
 
- num_x = nxc
- num_y = nyc
- num_z = nzc
+ num_xc = nxc
+ num_yc = nyc
+ num_zc = nzc
 
+ num_xe = nxe
+ num_ye = nye
+ num_ze = nze
+
 end subroutine get_gridsize
 
 
@@ -1681,10 +1720,10 @@
 real(r8), allocatable, dimension(:,:,:)     :: data_3d_array
 real(r8), allocatable, dimension(:,:,:,:)   :: data_4d_array
 
-integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, mystart, mycount
 character(len=NF90_MAX_NAME) :: varname 
-integer :: VarID, numdims, dimlen
-integer :: ncid
+integer :: VarID, ncNdims, dimlen
+integer :: ncid, TimeDimID, TimeDimLength
 character(len=256) :: myerrorstring 
 
 if ( .not. module_initialized ) call static_init_model
@@ -1720,14 +1759,22 @@
 if (do_output()) &
     call print_date(statedate,'date of restart file '//trim(filename))
 
-! Start filling the real Nd arrays from the single 1d state vector.
 
-indx = 1
+! The DART prognostic variables are only defined for a single time.
+! We already checked the assumption that variables are xy2d or xyz3d ...
+! IF the netCDF variable has a TIME dimension, it must be the last dimension,
+! and we need to read the LAST timestep and effectively squeeze out the
+! singleton dimension when we stuff it into the DART state vector. 
 
-! If these arrays have an extra dimension (like TIME), and it is only 1,
-! we might have to make the query code smarter, and might have to set
-! a start and count on the get_var() calls.  If it is more than 1, then
-! we have a problem.
+TimeDimID = FindTimeDimension( ncid )
+
+if ( TimeDimID > 0 ) then
+   call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=TimeDimLength), &
+            'sv_to_restart_file', 'inquire timedimlength '//trim(filename))
+else
+   TimeDimLength = 0
+endif
+
 do ivar=1, nfields
 
    varname = trim(progvar(ivar)%varname)
@@ -1738,70 +1785,96 @@
    call nc_check(nf90_inq_varid(ncid,   varname, VarID), &
             'sv_to_restart_file', 'inq_varid '//trim(myerrorstring))
 
-   call nc_check(nf90_inquire_variable(ncid,VarId,dimids=dimIDs,ndims=numdims), &
+   call nc_check(nf90_inquire_variable(ncid,VarId,dimids=dimIDs,ndims=ncNdims), &
             'sv_to_restart_file', 'inquire '//trim(myerrorstring))
 
-   do i = 1,numdims
+   mystart = 1   ! These are arrays, actually.
+   mycount = 1
+   ! Only checking the shape of the variable - sans TIME
+   DimCheck : do i = 1,progvar(ivar)%numdims
+
       write(string1,'(''inquire dimension'',i2,A)') i,trim(myerrorstring)
       call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
             'sv_to_restart_file', string1)
 
       if ( dimlen /= progvar(ivar)%dimlens(i) ) then
-         write(string1,*) trim(myerrorstring),'dim/dimlen',i,dimlen,'not',progvar(ivar)%dimlens(i)
+         write(string1,*) trim(myerrorstring),' dim/dimlen ',i,dimlen,' not ',progvar(ivar)%dimlens(i)
          call error_handler(E_ERR,'sv_to_restart_file',string1,source,revision,revdate)
       endif
-   enddo
 
-   if (numdims == 1) then
-      ni = progvar(ivar)%dimlens(1)
+      mycount(i) = dimlen
+
+   enddo DimCheck
+
+   where(dimIDs == TimeDimID) mystart = TimeDimLength
+   where(dimIDs == TimeDimID) mycount = 1   ! only the latest one
+
+   if ( debug > 5 ) then
+      write(*,*)'sv_to_restart_file '//trim(varname)//' start is ',mystart(1:ncNdims)
+      write(*,*)'sv_to_restart_file '//trim(varname)//' count is ',mycount(1:ncNdims)
+   endif
+
+   indx = progvar(ivar)%index1
+
+   if (ncNdims == 1) then
+      ni = mycount(1)
       allocate(data_1d_array(ni))
 
-      do i = 1, ni   ! size(data_1d_array,1)
+      do i = 1, ni
          data_1d_array(i) = state_vector(indx) 
          indx = indx + 1
       enddo
 
-      call nc_check(nf90_put_var(ncid, VarID, data_1d_array), &
+      call nc_check(nf90_put_var(ncid, VarID, data_1d_array, &
+        start=mystart(1:ncNdims), count=mycount(1:ncNdims)), &
             'sv_to_restart_file', 'put_var '//trim(varname))
       deallocate(data_1d_array)
-   elseif (numdims == 2) then
-      ni = progvar(ivar)%dimlens(1)
-      nj = progvar(ivar)%dimlens(2)
+
+   elseif (ncNdims == 2) then
+
+      ni = mycount(1)
+      nj = mycount(2)
       allocate(data_2d_array(ni, nj))
 
-      do j = 1, nj   ! size(data_2d_array,2)
-      do i = 1, ni   ! size(data_2d_array,1)
+      do j = 1, nj
+      do i = 1, ni
          data_2d_array(i, j) = state_vector(indx) 
          indx = indx + 1
       enddo
       enddo
 
-      call nc_check(nf90_put_var(ncid, VarID, data_2d_array), &
+      call nc_check(nf90_put_var(ncid, VarID, data_2d_array, &
+        start=mystart(1:ncNdims), count=mycount(1:ncNdims)), &
             'sv_to_restart_file', 'put_var '//trim(varname))
       deallocate(data_2d_array)
-   elseif (numdims == 3) then
-      ni = progvar(ivar)%dimlens(1)
-      nj = progvar(ivar)%dimlens(2)
-      nk = progvar(ivar)%dimlens(3)
+
+   elseif (ncNdims == 3) then
+
+      ni = mycount(1)
+      nj = mycount(2)
+      nk = mycount(3)
       allocate(data_3d_array(ni, nj, nk))
    
-      do k = 1, nk   ! size(data_3d_array,3)
-      do j = 1, nj   ! size(data_3d_array,2)
-      do i = 1, ni   ! size(data_3d_array,1)
+      do k = 1, nk
+      do j = 1, nj
+      do i = 1, ni
          data_3d_array(i, j, k) = state_vector(indx) 
          indx = indx + 1
       enddo
       enddo
       enddo
 
-      call nc_check(nf90_put_var(ncid, VarID, data_3d_array), &
+      call nc_check(nf90_put_var(ncid, VarID, data_3d_array, &
+        start=mystart(1:ncNdims), count=mycount(1:ncNdims)), &
             'sv_to_restart_file', 'put_var '//trim(varname))
       deallocate(data_3d_array)
-   elseif (numdims == 4) then
-      ni = progvar(ivar)%dimlens(1)
-      nj = progvar(ivar)%dimlens(2)
-      nk = progvar(ivar)%dimlens(3)
-      nl = progvar(ivar)%dimlens(4)
+
+   elseif (ncNdims == 4) then
+
+      ni = mycount(1)
+      nj = mycount(2)
+      nk = mycount(3)
+      nl = mycount(4)
       allocate(data_4d_array(ni, nj, nk, nl))
 
       do l = 1, nl   ! size(data_4d_array,3)
@@ -1815,21 +1888,30 @@
       enddo
       enddo
 
-      call nc_check(nf90_put_var(ncid, VarID, data_4d_array), &
+      call nc_check(nf90_put_var(ncid, VarID, data_4d_array, &
+        start=mystart(1:ncNdims), count=mycount(1:ncNdims)), &
             'sv_to_restart_file', 'put_var '//trim(varname))
       deallocate(data_4d_array)
+
    else
-      write(string1, *) 'no support for data array of dimension ', numdims
+      write(string1, *) 'no support for data array of dimension ', ncNdims
       call error_handler(E_ERR,'sv_to_restart_file', string1, &
                         source,revision,revdate)
    endif
 
+   indx = indx - 1
+   if ( indx /= progvar(ivar)%indexN ) then
+      write(string1, *)'Variable '//trim(varname)//' filled wrong.'
+      write(string2, *)'Should have ended at ',progvar(ivar)%indexN,' actually ended at ',indx
+      call error_handler(E_ERR,'sv_to_restart_file', string1, &
+                        source,revision,revdate,text2=string2)
+   endif
+
 enddo
 
 call nc_check(nf90_close(ncid), &
              'sv_to_restart_file','close '//trim(filename))
 
-
 end subroutine sv_to_restart_file
 
 
@@ -1840,977 +1922,255 @@
 
 
 
-subroutine init_interp()
+!############################################################################
+!
+!     ##################################################################
+!     ######                                                      ######
+!     ######            SUBROUTINE MODEL_INTERPOLATE              ######
+!     ######                                                      ######
+!     ##################################################################
+!
+!
+!     PURPOSE:
+!
+!     For a given lat, lon, and height, interpolate the correct state value
+!     to that location for the filter from the NCOMMAS state vectors
+!
+!############################################################################
+!
+!     Author:  Tim Hoar, Ted Mansell, Lou Wicker
+!
+!     Creation Date:  August 2010
+!     
+!     Variables needed to be stored in the MODEL_MODULE data structure
+!
+!       XG_POS   = time dependent X (meters) offset of model from its lat/lon point
+!       YG_POS   = time dependent Y (meters) offset of model from its lat/lon point
+!       XE, XC   = 1D arrays storing the local grid edge and center coords (meters)
+!       YE, YC   = 1D arrays storing the local grid edge and center coords (meters)
+!       ZE, ZC   = 1D arrays storing the local grid edge and center coords (meters)
+!       REF_LAT  = grid reference latitude
+!       REF_LON  = grid reference longitude
+!
+!       ERROR codes:
+!
+!       ISTATUS = 99:  general error in case something terrible goes wrong...
+!       ISTATUS = 15:  dont know what to do with vertical coord supplied
+!       ISTATUS = 16:  Obs_type is not found
+!       ISTATUS = 11:  index from from xi is outside domain WRT XE
+!       ISTATUS = 12:  index from from xi is outside domain WRT XC
+!       ISTATUS = 21:  index from from yi is outside domain WRT YE
+!       ISTATUS = 22:  index from from yi is outside domain WRT YC
+!       ISTATUS = 31:  index from from zi is outside domain WRT ZE
+!       ISTATUS = 32:  index from from zi is outside domain WRT ZC
+!       ISTATUS = 33:  index from location(3) is not bounded by 1 --> nz-1
+!       
+!
+!############################################################################
 
-! Initializes data structures needed for ncommas interpolation for
-! either dipole or irregular grid.
-! This should be called at static_init_model time to avoid 
-! having all this temporary storage in the middle of a run.
 
-integer :: i
-
-! Determine whether this is a irregular lon-lat grid or a dipole.
-! Do this by seeing if the lons have the same values at both
-! the first and last latitude row; this is not the case for dipole. 
-
-dipole_grid = .false.
-do i = 1, nxc
-   if(ulon(i, 1) /= ulon(i, nyc)) then
-      dipole_grid = .true.
-!     call init_dipole_interp()
-      return
-   endif
-enddo
-
-end subroutine init_interp
-
-
-!------------------------------------------------------------
-
-
-subroutine get_reg_box_indices(lon, lat, x_ind, y_ind)
-
-! Given a longitude and latitude in degrees returns the index of the regular
-! lon-lat box that contains the point.
-
-real(r8), intent(in)  :: lon, lat
-integer,  intent(out) :: x_ind, y_ind
-
-call get_reg_lon_box(lon, x_ind)
-call get_reg_lat_box(lat, y_ind)
-
-end subroutine get_reg_box_indices
-
-!------------------------------------------------------------
-
-subroutine get_reg_lon_box(lon, x_ind)
-
-! Determine which regular longitude box a longitude is in.
-
-real(r8), intent(in)  :: lon
-integer,  intent(out) :: x_ind
-
-x_ind = int(num_reg_x * lon / 360.0_r8) + 1
-
-! Watch out for exactly at top; assume all lats and lons in legal range
-if(lon == 360.0_r8) x_ind = num_reg_x
-
-end subroutine get_reg_lon_box
-
-!------------------------------------------------------------
-
-subroutine get_reg_lat_box(lat, y_ind)
-
-! Determine which regular latitude box a latitude is in.
-
-real(r8), intent(in)  :: lat
-integer,  intent(out) :: y_ind
-
-y_ind = int(num_reg_y * (lat + 90.0_r8) / 180.0_r8) + 1
-
-! Watch out for exactly at top; assume all lats and lons in legal range
-if(lat == 90.0_r8)  y_ind = num_reg_y
-
-end subroutine get_reg_lat_box
-
-!------------------------------------------------------------
-
-subroutine reg_box_overlap(x_corners, y_corners, is_pole, reg_lon_ind, reg_lat_ind)
-
-real(r8), intent(in)  :: x_corners(4), y_corners(4)
-logical,  intent(in)  :: is_pole
-integer,  intent(out) :: reg_lon_ind(2), reg_lat_ind(2)
-
-! Find a set of regular lat lon boxes that covers all of the area covered by 
-! a dipole grid qaud whose corners are given by the dimension four x_corners 
-! and y_corners arrays.  The two dimensional arrays reg_lon_ind and reg_lat_ind
-! return the first and last indices of the regular boxes in latitude and
-! longitude respectively. These indices may wraparound for reg_lon_ind.  
-! A special computation is needed for a dipole quad that has the true north 
-! pole in its interior. The logical is_pole is set to true if this is the case.
-! This can only happen for the t grid.  If the longitude boxes overlap 0
-! degrees, the indices returned are adjusted by adding the total number of
-! boxes to the second index (e.g. the indices might be 88 and 93 for a case
-! with 90 longitude boxes).
-
-real(r8) :: lat_min, lat_max, lon_min, lon_max
-integer  :: i
-
-!  A quad containing the pole is fundamentally different
-if(is_pole) then
-   ! Need all longitude boxes
-   reg_lon_ind(1) = 1
-   reg_lon_ind(2) = num_reg_x
-   ! Need to cover from lowest latitude to top box
-   lat_min = minval(y_corners)
-   reg_lat_ind(1) = int(num_reg_y * (lat_min + 90.0_r8) / 180.0_r8) + 1
-   call get_reg_lat_box(lat_min, reg_lat_ind(1))
-   reg_lat_ind(2) = num_reg_y
-else
-   ! All other quads do not contain pole (pole could be on edge but no problem)
-   ! This is specific to the dipole ncommas grids that do not go to the south pole
-   ! Finding the range of latitudes is cake
-   lat_min = minval(y_corners)
-   lat_max = maxval(y_corners)
-
-   ! Figure out the indices of the regular boxes for min and max lats
-   call get_reg_lat_box(lat_min, reg_lat_ind(1))
-   call get_reg_lat_box(lat_max, reg_lat_ind(2))
-
-   ! Lons are much trickier. Need to make sure to wraparound the
-   ! right way. There is no guarantee on direction of lons in the
-   ! high latitude dipole rows.
-   ! All longitudes for non-pole rows have to be within 180 degrees
-   ! of one another. 
-   lon_min = minval(x_corners)
-   lon_max = maxval(x_corners)
-   if((lon_max - lon_min) > 180.0_r8) then
-      ! If the max longitude value is more than 180 
-      ! degrees larger than the min, then there must be wraparound.
-      ! Then, find the smallest value > 180 and the largest < 180 to get range.
-      lon_min = 360.0_r8
-      lon_max = 0.0_r8
-      do i=1, 4
-         if(x_corners(i) > 180.0_r8 .and. x_corners(i) < lon_min) lon_min = x_corners(i)
-         if(x_corners(i) < 180.0_r8 .and. x_corners(i) > lon_max) lon_max = x_corners(i)
-      enddo
-   endif
-
-   ! Get the indices for the extreme longitudes
-   call get_reg_lon_box(lon_min, reg_lon_ind(1))
-   call get_reg_lon_box(lon_max, reg_lon_ind(2))
-
-   ! Watch for wraparound again; make sure that second index is greater than first
-   if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
-endif
-
-end subroutine reg_box_overlap
-
-!------------------------------------------------------------
-
-subroutine get_quad_corners(x, i, j, corners)
-
-real(r8), intent(in)  :: x(:, :)
-integer,  intent(in)  :: i, j
-real(r8), intent(out) :: corners(4)
-
-! Grabs the corners for a given quadrilateral from the global array of lower
-! right corners. Note that corners go counterclockwise around the quad.
-
-integer :: ip1
-
-! Have to worry about wrapping in longitude but not in latitude
-ip1 = i + 1
-if(ip1 > nxc) ip1 = 1
-
-corners(1) = x(i,   j  ) 
-corners(2) = x(ip1, j  )
-corners(3) = x(ip1, j+1)
-corners(4) = x(i,   j+1)
-
-end subroutine get_quad_corners
-
-!------------------------------------------------------------
-
-subroutine update_reg_list(reg_list_num, reg_list_lon, reg_list_lat, &
-   reg_lon_ind, reg_lat_ind, dipole_lon_index, dipole_lat_index)
-
-integer, intent(inout) :: reg_list_num(:, :), reg_list_lon(:, :, :), reg_list_lat(:, :, :)
-integer, intent(inout) :: reg_lon_ind(2), reg_lat_ind(2)
-integer, intent(in)    :: dipole_lon_index, dipole_lat_index
- 
-! Updates the data structure listing dipole quads that are in a given regular box
-integer :: ind_x, index_x, ind_y
-
-! Loop through indices for each possible regular cell
-! Have to watch for wraparound in longitude
-if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
-
-do ind_x = reg_lon_ind(1), reg_lon_ind(2)
-   ! Inside loop, need to go back to wraparound indices to find right box
-   index_x = ind_x
-   if(index_x > num_reg_x) index_x = index_x - num_reg_x
-   
-   do ind_y = reg_lat_ind(1), reg_lat_ind(2)
-      ! Make sure the list storage isn't full
-      if(reg_list_num(index_x, ind_y) >= max_reg_list_num) then
-         write(string1,*) 'max_reg_list_num (',max_reg_list_num,') is too small ... increase'
-         call error_handler(E_ERR, 'update_reg_list', string1, source, revision, revdate)
-      endif
-
-      ! Increment the count
-      reg_list_num(index_x, ind_y) = reg_list_num(index_x, ind_y) + 1
-      ! Store this quad in the list for this regular box
-      reg_list_lon(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lon_index
-      reg_list_lat(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lat_index
-   enddo
-enddo
-
-end subroutine update_reg_list
-
-!------------------------------------------------------------
-
 subroutine model_interpolate(x, location, obs_type, interp_val, istatus)
 
-real(r8),            intent(in) :: x(:)
-type(location_type), intent(in) :: location
-integer,             intent(in) :: obs_type
-real(r8),           intent(out) :: interp_val
-integer,            intent(out) :: istatus
+! Passed variables
 
-! Model interpolate will interpolate any state variable (S, T, U, V, PSURF) to
-! the given location given a state vector. The type of the variable being
-! interpolated is obs_type since normally this is used to find the expected
-! value of an observation at some location. The interpolated value is 
-! returned in interp_val and istatus is 0 for success.
+  real(r8),            intent(in)  :: x(:)
+  type(location_type), intent(in)  :: location
+  integer,             intent(in)  :: obs_type
+  real(r8),            intent(out) :: interp_val
+  integer,             intent(out) :: istatus
 
 ! Local storage
-real(r8)       :: loc_array(3), llon, llat, lheight
-integer        :: base_offset, offset, ind, i, ival
-integer        :: hgt_bot, hgt_top
-real(r8)       :: hgt_fract
-real(r8)       :: top_val, bot_val
-integer        :: hstatus
 
-if ( .not. module_initialized ) call static_init_model
+  real(r8)         :: loc_array(3), llon, llat
+  real(r8)         :: lheight
+  integer          :: base_offset, offset
+  integer          :: nf, i, n0, n1
+  integer          :: iloc, jloc, kloc, nxp, nyp, nzp
+  real(r8)         :: xi, yi, xf, yf, zf, q1, q2, vt, vb
 
+  IF ( .not. module_initialized ) call static_init_model
+
 ! Let's assume failure.  Set return val to missing, then the code can
 ! just set istatus to something indicating why it failed, and return.
 ! If the interpolation is good, the interp_val will be set to the 
 ! good value, and the last line here sets istatus to 0.
 ! make any error codes set here be in the 10s
 
-interp_val = MISSING_R8     ! the DART bad value flag
-istatus = 99                ! unknown error
+  interp_val = MISSING_R8     ! the DART bad value flag
+  istatus = 99                ! unknown error
 
 ! Get the individual locations values
-loc_array = get_location(location)
-llon    = loc_array(1)
-llat    = loc_array(2)
-lheight = loc_array(3)
 
-if (debug > 1) print *, 'requesting interpolation at ', llon, llat, lheight
+  loc_array = get_location(location)
+  llon      = loc_array(1)
+  llat      = loc_array(2)
+  lheight   = loc_array(3)
 
-if( vert_is_height(location) ) then
-   ! Nothing to do 
-elseif ( vert_is_surface(location) ) then
-   ! Nothing to do 
-elseif (vert_is_level(location)) then
-   ! convert the level index to an actual height 
-   ind = nint(loc_array(3))
-   if ( (ind < 1) .or. (ind > size(zc)) ) then 
-      istatus = 11
-      return
-   else
-      lheight = zc(ind)
-   endif
-else   ! if pressure or undefined, we don't know what to do
-   istatus = 17
-   return
-endif
+  IF (debug > 1) print *, 'requesting interpolation at ', llon, llat, lheight
 
-! Do horizontal interpolations for the appropriate levels
-! Find the basic offset of this field
+  IF( vert_is_height(location) ) THEN        ! Nothing to do 
+  
+  ELSEIF ( vert_is_surface(location) ) THEN  ! Nothing to do
+ 
+  ELSEIF (vert_is_level(location)) THEN      ! convert the level index to an actual height
+                                             ! This code is wrong (I think) if variable == "W"
+     kloc = nint(loc_array(3))
+     IF( (kloc < 1) .or. (kloc > size(zc)) ) THEN 
+        istatus = 33
+        return
+     ELSE
+        lheight = zc(kloc)
+     ENDIF
+  ELSE   ! if we don't know what to do
+     istatus = 15
+     return
+  ENDIF
 
-ival = -1
-do i = i, nfields
-   if (progvar(ival)%dart_kind == obs_type) exit
-enddo
+! Find what field this is....
 
+  nf = -1
+  DO i = 1,nfields
+     IF (progvar(i)%dart_kind == obs_type) THEN
+       nf = i
+       exit
+     ENDIF
+  ENDDO
+
 ! Not a legal type for interpolation, return istatus error
-if (ival < 0) then
-   istatus = 15
-   return
-endif
 
-base_offset = progvar(ival)%index1
+  IF (nf < 0) then
+     istatus = 16
+     return
+  ENDIF
 
-if (debug > 1) print *, 'base offset now ', base_offset
+  base_offset = progvar(nf)%index1
 
-! surface variables are simpler
-if( vert_is_surface(location) ) then
-   call lon_lat_interpolate(x(base_offset:), llon, llat, obs_type, 1, interp_val, istatus)
-   return
-endif
+  IF (debug > 1) print *, 'base offset now ', base_offset
 
-! Get the bounding vertical levels and the fraction between bottom and top
-call height_bounds(lheight, nzc, zc, hgt_bot, hgt_top, hgt_fract, hstatus)
-if(hstatus /= 0) then
-   istatus = 12
-   return
-endif
+! Not implemented...
 
-! Find the base location for the bottom height and interpolate horizontally 
-!  on this level.  Do bottom first in case it is below the ocean floor; can
-!  avoid the second horizontal interpolation.
-offset = base_offset + (hgt_bot - 1) * nxc * nyc
-if (debug > 1) print *, 'relative bot height offset = ', offset - base_offset
-if (debug > 1) print *, 'absolute bot height offset = ', offset
-call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_bot, bot_val, istatus)
-! Failed istatus from interpolate means give up
-if(istatus /= 0) return
+  IF( vert_is_surface(location) ) THEN
+     istatus = 15
+     return
+  ENDIF
+  
+! From input lat/lon position, find x/y distance from grid lat/lon
+! Assumed lat/lon input is in degrees (.true. flag is used to convert inside routine)
 
-! Find the base location for the top height and interpolate horizontally 
-!  on this level.
-offset = base_offset + (hgt_top - 1) * nxc * nyc
-if (debug > 1) print *, 'relative top height offset = ', offset - base_offset
-if (debug > 1) print *, 'absolute top height offset = ', offset
-call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_top, top_val, istatus)
-! Failed istatus from interpolate means give up
-if(istatus /= 0) return
+  CALL ll_to_xy(xi, yi, 0, ref_lat, ref_lon, llat, llon, .true.)
+  
+  IF( debug > 1 ) THEN
+    print *, 'XMIN / X-LOC from lat/lon / XMAX:  ', xe(1), xi, xe(nxe)
+    print *, 'YMIN / Y-LOC from lat/lon / YMAX:  ', ye(1), yi, ye(nye)
+  ENDIF
+  
+! Using nf, get the true dimensions of the variable
 
+   nxp = progvar(nf)%dimlens(1)
+   nyp = progvar(nf)%dimlens(2)
+   nzp = progvar(nf)%dimlens(3)   
+   
+! Get X (LON) index
 
-! Then weight them by the vertical fraction and return
-interp_val = bot_val + hgt_fract * (top_val - bot_val)
-if (debug > 1) print *, 'model_interp: interp val = ',interp_val
+  IF(obs_type == 1) THEN
+    iloc = find_index(xi, xe, nxe)
+    IF( iloc == -1) THEN  ! Error occurred, observation is outside of the domain
+     istatus = 11
+     return
+    ENDIF
+    xf = (xe(iloc+1) - xi) / (xe(iloc+1) - xe(iloc))
+  ELSE
+    iloc = find_index(xi, xc, nxc)
+    IF( iloc == -1) THEN  ! Error occurred, observation is outside of the domain
+     istatus = 12
+     return
+    ENDIF
+    xf = (xc(iloc+1) - xi) / (xc(iloc+1) - xc(iloc))
+  ENDIF
+  
+! Get Y (LAT) index
 
+   IF(obs_type == 2) THEN
+    jloc = find_index(yi, ye, nye)
+    IF( jloc == -1) THEN  ! Error occurred, observation is outside of the vertical domain
+     istatus = 21
+     return
+    ENDIF
+    yf = (ye(jloc+1) - yi) / (ye(jloc+1) - ye(jloc))
+  ELSE
+    jloc = find_index(yi, yc, nyc) 
+    IF( jloc == -1) THEN  ! Error occurred, observation is outside of the vertical domain
+     istatus = 22
+     return
+    ENDIF
+    yf = (yc(jloc+1) - yi) / (yc(jloc+1) - yc(jloc))
+  ENDIF
+  
+! Get Z (HGT) index
+
+  IF(obs_type == 3) THEN
+    kloc = find_index(lheight, ze, nze)
+    IF( kloc == -1) THEN  ! Error occurred, observation is outside of the vertical domain
+     istatus = 31
+     return
+    ENDIF
+    zf = (ze(kloc+1) - lheight) / (ze(kloc+1) - ze(kloc))
+  ELSE
+    kloc = find_index(lheight, zc, nzc)
+    IF( kloc == -1) THEN  ! Error occurred, observation is outside of the vertical domain
+     istatus = 32
+     return
+    ENDIF
+    zf = (zc(kloc+1) - lheight) / (zc(kloc+1) - zc(kloc))
+  ENDIF
+
+  IF( debug > 1 ) THEN
+    print *, "ILOC:  ", iloc
+    print *, "XI = ", xi
+    print *, "XF = ", xf
+    print *, "JLOC:  ", jloc
+    print *, "YI = ", yi 
+    print *, "YF = ", yf
+    print *, "ZLOC = ", kloc
+    print *, "ZI = ", lheight 
+    print *, "ZF = ", zf
+  ENDIF
+  
+! Find all the corners of the trilinear cube and then form the 1D->2D->3D interpolation
+! Since base_offset is always the first point (e.g., (1,1,1)), then subtract "1" to get that loc
+  
+  n0 = base_offset + (jloc-1)*nxp + (kloc-1)*nxp*nyp + iloc - 1  ! Location of (i,  j,k) ?                         
+  n1 = base_offset + (jloc-1)*nxp + (kloc-1)*nxp*nyp + iloc      ! Location of (i+1,j,k) ?                    
+  q1 = (1.0_r8-xf)*x(n1) + xf*x(n0)                              ! Value along "j,k" grid line
+  
+  n0 = base_offset + (jloc)*nxp + (kloc-1)*nxp*nyp + iloc - 1    ! Location of (i,  j+1,k) ?                         
+  n1 = base_offset + (jloc)*nxp + (kloc-1)*nxp*nyp + iloc        ! Location of (i+1,j+1,k) ?                    
+  q2 = (1.0_r8-xf)*x(n1) + xf*x(n0)                              ! Value along "j+1,k" grid line
+  
+  vb = (1.0_r8-yf)*q2 + yf*q1                                    ! Binlinear value on the bottom plane
+  
+  n0 = base_offset + (jloc-1)*nxp + (kloc)*nxp*nyp + iloc - 1    ! Location of (i,  j,k+1) ?                         
+  n1 = base_offset + (jloc-1)*nxp + (kloc)*nxp*nyp + iloc        ! Location of (i+1,j,k+1) ?                    
+  q1 = (1.0_r8-xf)*x(n1) + xf*x(n0)                              ! Value along "j,k+1" grid line
+  
+  n0 = base_offset + (jloc)*nxp + (kloc)*nxp*nyp + iloc - 1      ! Location of (i,  j+1,k+1) ?                         

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list