[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