[Dart-dev] [6104] DART/branches/development: The forward observation operator for flux tower observations can now
nancy at ucar.edu
nancy at ucar.edu
Fri May 3 13:32:21 MDT 2013
Revision: 6104
Author: thoar
Date: 2013-05-03 13:32:20 -0600 (Fri, 03 May 2013)
Log Message:
-----------
The forward observation operator for flux tower observations can now
read the CLM history files for global/regional AND unstructured grids.
The unstructured grid scenario is only supported for a single gridcell.
It can be extended by a more comprehensive 'which gridcell is closest
to the observation' logic. At present I can think of nothing more than
brute force for unstructured grids.
This observation operator is slow enough already. It opens all the
netCDF files ONCE (and never closes them, actually), query, and
reads a scalar for each observation. These files are small, if it
gets to the point where speed is needed, a new strategy will be needed.
Modified Paths:
--------------
DART/branches/development/models/clm/work/input.nml
DART/branches/development/obs_def/obs_def_tower_mod.f90
Added Paths:
-----------
DART/branches/development/obs_def/obs_def_tower_mod.nml
-------------- next part --------------
Modified: DART/branches/development/models/clm/work/input.nml
===================================================================
--- DART/branches/development/models/clm/work/input.nml 2013-05-03 17:23:42 UTC (rev 6103)
+++ DART/branches/development/models/clm/work/input.nml 2013-05-03 19:32:20 UTC (rev 6104)
@@ -122,14 +122,13 @@
input_obs_def_mod_file = '../../../obs_def/DEFAULT_obs_def_mod.F90',
output_obs_def_mod_file = '../../../obs_def/obs_def_mod.f90',
input_files = '../../../obs_def/obs_def_tower_mod.f90',
- '../../../obs_def/obs_def_COSMOS_mod.f90',
+ '../../../obs_def/obs_def_COSMOS_mod.f90'
/
# casename will get overwritten in the assimilate.csh script.
&obs_def_tower_nml
- casename = '../dummy',
- verbose = .true.,
+ casename = '../clm_dart',
debug = .true.
/
@@ -281,16 +280,15 @@
rat_cri = 3.0,
input_qc_threshold = 1,
Nregions = 4,
- lonlim1 = 0.0, 0.0, 0.0, 290,
- lonlim2 = 360.0, 360.0, 360.0, 380,
- latlim1 = 20.0, 45.0, 0.0, 20,
- latlim2 = 75.0, 75.0, 75.0, 50,
- reg_names = '20N-75N', '45N-75N', '0-75N', 'North Atlantic',
+ lonlim1 = 0.0, 0.0, 0.0, 180.0,
+ lonlim2 = 360.0, 360.0, 180.0, 360.0,
+ latlim1 = -90.0, 0.0, -90.0, -90.0,
+ latlim2 = 0.0, 90.0, 90.0, 90.0,
+ reg_names = 'South', 'North', 'East', 'West',
+ hlevel = 0.0, 0.05, 0.1, 0.2, 0.50, 1.0, 2.0, 5.0, 10.0, 20.0, 40.0,
print_mismatched_locs = .false.,
create_rank_histogram = .true.,
verbose = .true.,
- hlevel = 10.0, 20.0, 30.0, 40.0, 100.0,
- 200.0, 500.0, 1000.0, 2000.0,
/
Modified: DART/branches/development/obs_def/obs_def_tower_mod.f90
===================================================================
--- DART/branches/development/obs_def/obs_def_tower_mod.f90 2013-05-03 17:23:42 UTC (rev 6103)
+++ DART/branches/development/obs_def/obs_def_tower_mod.f90 2013-05-03 19:32:20 UTC (rev 6104)
@@ -26,18 +26,18 @@
!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
-! use obs_def_tower_mod, only : get_scalar_from_3Dhistory
+! use obs_def_tower_mod, only : get_scalar_from_history
! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
! case(TOWER_LATENT_HEAT_FLUX)
-! call get_scalar_from_3Dhistory('EFLX_LH_TOT_R', state_time, ens_index, location, obs_time, obs_val, istatus)
+! call get_scalar_from_history('EFLX_LH_TOT_R', state_time, ens_index, location, obs_time, obs_val, istatus)
! case(TOWER_SENSIBLE_HEAT_FLUX)
-! call get_scalar_from_3Dhistory('FSH', state_time, ens_index, location, obs_time, obs_val, istatus)
+! call get_scalar_from_history('FSH', state_time, ens_index, location, obs_time, obs_val, istatus)
! case(TOWER_NETC_ECO_EXCHANGE)
-! call get_scalar_from_3Dhistory('NEP', state_time, ens_index, location, obs_time, obs_val, istatus)
+! call get_scalar_from_history('NEP', state_time, ens_index, location, obs_time, obs_val, istatus)
! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
!-----------------------------------------------------------------------------
@@ -66,6 +66,14 @@
! BEGIN DART PREPROCESS MODULE CODE
module obs_def_tower_mod
+! This is the forward observation operator code for CLM for flux tower observations.
+! The DART model state for CLM generally does not have all the pieces necessary
+! to apply the forward observation operator directly, so this code gets what it
+! needs from a CLM history file. These history files have become complicated now
+! that CLM is trying to support unstructured grids. Sometimes the variables of
+! interest are shaped NEP(time, lat, lon), sometimes NEP(time, lndgrid).
+! 'single column' runs may appear as either lat=lon=1 or lndgrid=1
+!
! <next few lines under version control, do not edit>
! $URL$
! $Id$
@@ -73,7 +81,8 @@
! $Date$
use types_mod, only : r4, r8, digits12, MISSING_R8, PI, deg2rad
-use location_mod, only : location_type, get_location
+use location_mod, only : location_type, get_location, get_dist, &
+ set_location, VERTISUNDEF
use time_manager_mod, only : time_type, get_date, set_date, print_date, print_time, &
get_time, set_time, operator(-), operator(/=)
use utilities_mod, only : register_module, E_ERR, E_MSG, error_handler, &
@@ -87,7 +96,7 @@
implicit none
private
-public :: get_scalar_from_3Dhistory
+public :: get_scalar_from_history
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -96,22 +105,23 @@
revdate = "$Date$"
logical :: module_initialized = .false.
+logical :: unstructured = .false.
character(len=129) :: string1, string2, string3
integer :: nlon, nlat, ntime, ens_size
type(time_type) :: initialization_time
-real(r8) :: edgeNorth,edgeEast,edgeSouth,edgeWest
character(len=129), allocatable, dimension(:) :: fname
integer, allocatable, dimension(:) :: ncid
-real(r8), allocatable, dimension(:) :: lon, lat
+real(r8), allocatable, dimension(:) :: lon, lat, area
real(digits12), allocatable, dimension(:) :: rtime
+real(r8), parameter :: RAD2KM = 40030.0_r8/(2.0_r8 * PI) ! (mean radius of earth ~6371km)
+
! namelist items
-character(len=256) :: casename = 'clm_tim'
-logical :: verbose = .false.
+character(len=256) :: casename = 'clm_dart'
logical :: debug = .false.
-namelist /obs_def_tower_nml/ casename, verbose, debug
+namelist /obs_def_tower_nml/ casename, debug
contains
@@ -128,7 +138,7 @@
! that have the observations, etc.
integer :: iunit, io, i
-integer :: dimid, varid
+integer :: varid
integer :: year, month, day, hour, minute, second, leftover
integer, allocatable, dimension(:) :: yyyymmdd,sssss
type(time_type) :: tower_time
@@ -139,7 +149,7 @@
string1 = 'model time does not match initialization time'
string2 = 'model time does not match initialization time'
string3 = 'model time does not match initialization time'
- call error_handler(E_ERR, 'initialize_routine', string1, &
+ call error_handler(E_ERR, 'obs_def_tower.initialize_routine', string1, &
source, revision, revdate, text2=string2,text3=string3)
endif
return
@@ -157,13 +167,13 @@
read(iunit, nml = obs_def_tower_nml, iostat = io)
call check_namelist_read(iunit, io, "obs_def_tower_nml")
-! Record the namelist values used for the run ...
+! Record the namelist values used for the run ...
if (do_nml_file()) write(nmlfileunit, nml=obs_def_tower_nml)
if (do_nml_term()) write( * , nml=obs_def_tower_nml)
! Need to know what day we are trying to assimilate.
! The model stops at midnight, we want all the observations for THE PREVIOUS DAY.
-! The CLM h0 files contain everything from 00:00 to 23:30 for the date in the filename.
+! The CLM h1 files contain everything from 00:00 to 23:30 for the date in the filename.
! The data for [23:30 -> 00:00] get put in the file for the next day.
tower_time = model_time - set_time(0,1)
@@ -172,115 +182,122 @@
! Figure out how many files (i.e. ensemble size) and construct their names.
! The CLM h0 files are constructed such that the midnight that starts the
-! day is IN the file. The last time in the file is 23:30 ...
+! day is IN the file. The last time in the file is 23:30 ...
+! In a perfect model scenario (i.e. 1 instance) CESM constructs filenames
+! without the instance counter in them. We check for that first. If it
+! exists we know the filename and ensemble size.
+
100 format (A,'.clm2_',I4.4,'.h1.',I4.4,'-',I2.2,'-',I2.2,'-',I5.5,'.nc')
+110 format (A,'.clm2' ,'.h1.',I4.4,'-',I2.2,'-',I2.2,'-',I5.5,'.nc')
+write(string1,110) trim(casename),year,month,day,second
+
+if( file_exist(string1) ) then
+ ! We know we are in a perfect model scenario
+ ens_size = 1
+ if (debug .and. do_output()) write(*,*)'Ensemble size is believed to be ',ens_size
+ goto 200
+endif
+
ens_size = 0
ENSEMBLESIZE : do i = 1,200
write(string1,100) trim(casename),i,year,month,day,second
if( file_exist(string1) ) then
- if(verbose .and. do_output()) write(*,*)'observation file "',trim(string1),'" exists.'
+ if(debug .and. do_output()) &
+ write(*,*)'observation file "',trim(string1),'" exists.'
ens_size = ens_size + 1
else
- if(verbose .and. do_output()) write(*,*)'WARNING observation file "',trim(string1),'" does not exist.'
+ if(debug .and. do_output()) &
+ write(*,*)'WARNING observation file "',trim(string1),'" does not exist.'
exit ENSEMBLESIZE
endif
+
enddo ENSEMBLESIZE
if (ens_size < 2) then
write(string1,100) trim(casename),1,year,month,day,second
write(string2,*)'cannot find files to use for observation operator.'
- write(string3,*)'trying files with names like "',trim(string1),'"'
- call error_handler(E_ERR, 'initialize_routine', string2, &
- source, revision, revdate, text2=string3)
+ write(string3,*)'trying files with names like "',trim(string1),'" -or-'
+ write(string1,110) trim(casename),year,month,day,second
+ call error_handler(E_ERR, 'obs_def_tower.initialize_routine', string2, &
+ source, revision, revdate, text2=string3,text3='"'//trim(string1)//'"')
elseif (ens_size >= 200) then
write(string2,*)'ensemble size (',ens_size,') is unnaturally large.'
write(string3,*)'trying files with names like "',trim(string1),'"'
- call error_handler(E_ERR, 'initialize_routine', string2, &
+ call error_handler(E_ERR, 'obs_def_tower.initialize_routine', string2, &
source, revision, revdate, text2=string3)
else
- if (verbose .and. do_output()) write(*,*)'Ensemble size is believed to be ',ens_size
+ if (debug .and. do_output()) then
+ write(string1,*)'Ensemble size is believed to be ',ens_size
+ call error_handler(E_MSG, 'obs_def_tower.initialize_routine', string1, &
+ source, revision, revdate )
+ endif
endif
+200 continue
+
allocate(fname(ens_size),ncid(ens_size))
ncid = 0
ENSEMBLE : do i = 1,ens_size
- write(fname(i),100) trim(casename),i,year,month,day,second
- call nc_check(nf90_open(trim(fname(i)), nf90_nowrite, ncid(i)), &
- 'initialize_routine','open '//trim(fname(i)))
+ if (ens_size == 1) then
+ write(fname(i),110) trim(casename),year,month,day,second
+ call nc_check(nf90_open(trim(fname(i)), nf90_nowrite, ncid(i)), &
+ 'obs_def_tower.initialize_routine','open '//trim(fname(i)))
+ exit ENSEMBLE
+ else
+ write(fname(i),100) trim(casename),i,year,month,day,second
+ call nc_check(nf90_open(trim(fname(i)), nf90_nowrite, ncid(i)), &
+ 'obs_def_tower.initialize_routine','open '//trim(fname(i)))
+ endif
enddo ENSEMBLE
i = 1
-! Harvest information from the first observation file.
! FIXME All other files will be opened to make sure they have the same dimensions.
-call nc_check(nf90_inq_dimid(ncid(i), 'lon', dimid), &
- 'initialize_routine','inq_dimid lon '//trim(fname(i)))
-call nc_check(nf90_inquire_dimension(ncid(i), dimid, len=nlon), &
- 'initialize_routine','inquire_dimension lon '//trim(fname(i)))
+call GetDimensions(ncid(i), fname(i)) ! determines values for nlon, nlat, ntime
-call nc_check(nf90_inq_dimid(ncid(i), 'lat', dimid), &
- 'initialize_routine','inq_dimid lat '//trim(fname(i)))
-call nc_check(nf90_inquire_dimension(ncid(i), dimid, len=nlat), &
- 'initialize_routine','inquire_dimension lat '//trim(fname(i)))
+allocate(lon(nlon), lat(nlat), rtime(ntime), yyyymmdd(ntime), sssss(ntime))
-call nc_check(nf90_inq_dimid(ncid(i), 'time', dimid), &
- 'initialize_routine','inq_dimid time '//trim(fname(i)))
-call nc_check(nf90_inquire_dimension(ncid(i), dimid, len=ntime), &
- 'initialize_routine','inquire_dimension lat '//trim(fname(i)))
+call nc_check(nf90_inq_varid(ncid(i), 'lon', varid), &
+ 'obs_def_tower.initialize_routine','inq_varid lon '//trim(fname(i)))
+call nc_check(nf90_get_var( ncid(i), varid, lon), &
+ 'obs_def_tower.initialize_routine','get_var lon'//trim(fname(i)))
-allocate(lon(nlon),lat(nlat))
-allocate(rtime(ntime),yyyymmdd(ntime),sssss(ntime))
+call nc_check(nf90_inq_varid(ncid(i), 'lat', varid), &
+ 'obs_def_tower.initialize_routine','inq_varid lat '//trim(fname(i)))
+call nc_check(nf90_get_var( ncid(i), varid, lat), &
+ 'obs_def_tower.initialize_routine','get_var lat'//trim(fname(i)))
-call nc_check(nf90_inq_varid(ncid(i), 'lon', varid), &
- 'initialize_routine','inq_varid lon '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, lon), 'initialize_routine', 'get_var lon')
-
-call nc_check(nf90_inq_varid(ncid(i), 'lat', varid), &
- 'initialize_routine','inq_varid lat '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, lat), 'initialize_routine', 'get_var lat')
-
call nc_check(nf90_inq_varid(ncid(i), 'mcdate', varid), &
- 'initialize_routine','inq_varid mcdate '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, yyyymmdd), 'initialize_routine', 'get_var yyyymmdd')
+ 'obs_def_tower.initialize_routine','inq_varid mcdate '//trim(fname(i)))
+call nc_check(nf90_get_var( ncid(i), varid, yyyymmdd), &
+ 'obs_def_tower.initialize_routine','get_var yyyymmdd'//trim(fname(i)))
-call nc_check(nf90_inq_varid(ncid(i), 'mcsec', varid), &
- 'initialize_routine','inq_varid mcsec '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, sssss), 'initialize_routine', 'get_var sssss')
+call nc_check(nf90_inq_varid(ncid(i), 'mcsec', varid), &
+ 'obs_def_tower.initialize_routine','inq_varid mcsec '//trim(fname(i)))
+call nc_check(nf90_get_var( ncid(i), varid, sssss), &
+ 'obs_def_tower.initialize_routine','get_var sssss'//trim(fname(i)))
-! Determine the geographic boundaries of the contents of the history file.
+if ( (nlon == 1) .and. (nlat ==1) ) then
+ allocate(area(nlon))
+ call nc_check(nf90_inq_varid(ncid(i), 'area', varid), &
+ 'obs_def_tower.initialize_routine','inq_varid area '//trim(fname(i)))
+ call nc_check(nf90_get_var(ncid(i), varid, area), &
+ 'obs_def_tower.initialize_routine', 'get_var area'//trim(fname(i)))
+ if (debug .and. do_output()) write(*,*)'obs_def_tower area',area(nlon)
+endif
-call nc_check(nf90_inq_varid(ncid(i), 'edgen', varid), &
- 'initialize_routine','inq_varid edgen '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, edgeNorth), 'initialize_routine', 'get_var edgeNorth')
-
-call nc_check(nf90_inq_varid(ncid(i), 'edgee', varid), &
- 'initialize_routine','inq_varid edgee '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, edgeEast), 'initialize_routine', 'get_var edgeEast')
-
-call nc_check(nf90_inq_varid(ncid(i), 'edges', varid), &
- 'initialize_routine','inq_varid edges '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, edgeSouth), 'initialize_routine', 'get_var edgeSouth')
-
-call nc_check(nf90_inq_varid(ncid(i), 'edgew', varid), &
- 'initialize_routine','inq_varid edgew '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, edgeWest), 'initialize_routine', 'get_var edgeWest')
-
-! call nc_check(nf90_inq_varid(ncid(i), 'time', varid), &
-! 'initialize_routine','inq_varid time '//trim(fname(i)))
-! call nc_check(nf90_get_var(ncid(i), varid, time), 'initialize_routine', 'get_var time')
-
! Convert time in file to a time compatible with the observation sequence file.
do i = 1,ntime
@@ -319,8 +336,98 @@
end subroutine initialize_module
+!======================================================================
-subroutine get_scalar_from_3Dhistory(varstring, state_time, ens_index, location, obs_time, obs_val, istatus )
+
+subroutine GetDimensions(ncid, fname)
+
+! Harvest information from the first observation file.
+! The SingleColumMode files have
+! float lat(lndgrid) ;
+! float lon(lndgrid) ;
+! float area(lndgrid) ;
+! while the 2D files have
+! float lat(lat) ;
+! float lon(lon) ;
+! float area(lat, lon) ;
+
+integer, intent(in) :: ncid
+character(len=*), intent(in) :: fname
+
+! integer, intent(out) :: nlon, nlat, ntime ... module variables
+
+integer, dimension(NF90_MAX_VAR_DIMS) :: londimids, latdimids
+integer :: lonvarid, lonndims
+integer :: latvarid, latndims
+integer :: dimid
+
+call nc_check(nf90_inq_varid(ncid, 'lon', lonvarid), &
+ 'obs_def_tower.GetDimensions','inq_varid lon '//trim(fname))
+call nc_check(nf90_inq_varid(ncid, 'lat', latvarid), &
+ 'obs_def_tower.GetDimensions','inq_varid lat '//trim(fname))
+
+call nc_check(nf90_inquire_variable(ncid, lonvarid, ndims=lonndims, dimids=londimids),&
+ 'obs_def_tower.GetDimensions','inquire lon '//trim(fname))
+call nc_check(nf90_inquire_variable(ncid, latvarid, ndims=latndims, dimids=latdimids),&
+ 'obs_def_tower.GetDimensions','inquire lat '//trim(fname))
+
+if ( (lonndims /= 1) .or. (latndims /= 1) ) then
+ write(string1,*) 'Require "lon" and "lat" variables to be 1D. They are ', &
+ lonndims, latndims
+ call error_handler(E_ERR,'obs_def_tower.GetDimensions',string1,source,revision,revdate)
+endif
+
+call nc_check(nf90_inquire_dimension(ncid, londimids(1), len=nlon), &
+ 'obs_def_tower.GetDimensions','inquire_dimension lon '//trim(fname))
+call nc_check(nf90_inquire_dimension(ncid, latdimids(1), len=nlat), &
+ 'obs_def_tower.GetDimensions','inquire_dimension lat '//trim(fname))
+
+call nc_check(nf90_inq_dimid(ncid, 'time', dimid), &
+ 'obs_def_tower.GetDimensions','inq_dimid time '//trim(fname))
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=ntime), &
+ 'obs_def_tower.GetDimensions','inquire_dimension time '//trim(fname))
+
+if ( nf90_inq_dimid(ncid, 'lndgrid', dimid) == NF90_NOERR) then
+ unstructured = .true.
+ if ((nlon /= 1) .or. (nlat /= 1)) then
+ string1 = 'unstructured grids with more than a single gridcell are not supported.'
+ call error_handler(E_ERR,'obs_def_tower.GetDimensions',string1,source,revision,revdate)
+ endif
+endif
+
+end subroutine GetDimensions
+
+
+!======================================================================
+
+
+subroutine get_scalar_from_history(varstring, state_time, ens_index, location, &
+ obs_time, obs_val, istatus)
+
+character(len=*), intent(in) :: varstring
+type(time_type), intent(in) :: state_time
+integer, intent(in) :: ens_index
+type(location_type), intent(in) :: location
+type(time_type), intent(in) :: obs_time
+real(r8), intent(out) :: obs_val
+integer, intent(out) :: istatus
+
+if ( .not. module_initialized ) call initialize_module(state_time)
+
+if ( unstructured ) then
+ call get_scalar_from_2Dhistory(varstring,ens_index,location,obs_time,obs_val,istatus)
+else
+ call get_scalar_from_3Dhistory(varstring,ens_index,location,obs_time,obs_val,istatus)
+endif
+
+end subroutine get_scalar_from_history
+
+
+!======================================================================
+
+
+subroutine get_scalar_from_3Dhistory(varstring, ens_index, location, obs_time, &
+ obs_val, istatus)
! the routine must return values for:
! obs_val -- the computed forward operator value
! istatus -- return code: 0=ok, > 0 is error, < 0 reserved for system use
@@ -335,7 +442,6 @@
! NEP:missing_value = 1.e+36f ;
character(len=*), intent(in) :: varstring
-type(time_type), intent(in) :: state_time
integer, intent(in) :: ens_index
type(location_type), intent(in) :: location
type(time_type), intent(in) :: obs_time
@@ -349,92 +455,123 @@
integer :: gridloni, gridlatj, timei
integer :: varid, xtype, ndims, natts, dimlen
integer :: io1, io2, second, day
-real(r8) :: loc_lon, loc_lat
+real(r8) :: loc_lon, loc_lat, radius, distance
real(r4), dimension(1) :: hyperslab
real(r4) :: spvalR4
real(r8) :: scale_factor, add_offset
real(digits12) :: otime
-character(len=20) :: strshort
+character(len=NF90_MAX_NAME+20) :: strshort
-if ( .not. module_initialized ) call initialize_module(state_time)
+type(location_type) :: gridloc
obs_val = MISSING_R8
istatus = 1
+!----------------------------------------------------------------------
! if observation is outside region encompassed in the history file - fail
loc = get_location(location) ! loc is in DEGREES
loc_lon = loc(1)
loc_lat = loc(2)
-if ( .not. is_longitude_between(loc_lon, edgeWest, edgeEast, doradians=.FALSE.)) return
-if ((loc_lat < edgeSouth) .or. (loc_lat > edgeNorth)) return
+if ( (nlon==1) .and. (nlat==1) ) then
+ ! Defining the region if running in an unstructured grid is tricky.
+ ! Have lat, lon, and the area of the gridcell which we assume to be basically square.
+ ! The square root of the area defines the length of the edge of the gridcell.
+ ! Half the hypotenuse defines the radius of a circle. Any ob within
+ ! that radius is close enough.
+
+ gridloc = set_location(lon(1),lat(1), 0.0_r8, VERTISUNDEF)
+ distance = get_dist(gridloc, location, no_vert = .TRUE.) * RAD2KM ! planet earth
+ radius = sqrt(2.0_r8 * area(1))/2.0_r8
+
+ if (debug .and. do_output()) then
+ write(string1,*)' observation lon, lat is ',loc_lon, loc_lat
+ write(string2,*)'gridcell lon, lat is ',lon(1),lat(1)
+ write(string3,*)'area,radius is ',area(1),radius,' distance ',distance
+ call error_handler(E_MSG, 'obs_def_tower.get_scalar_from_3Dhistory', &
+ string1, source, revision, revdate, text2=string2, text3=string3)
+ endif
+
+ if ( distance > radius ) return
+
+else
+ if ( .not. is_longitude_between(loc_lon, lon(1), lon(nlon), doradians=.FALSE.)) return
+ if ((loc_lat < lat(1)) .or. (loc_lat > lat(nlat))) return
+endif
+
+!----------------------------------------------------------------------
! Now that we know the observation operator is possible, continue ...
write(strshort,'(''ens_index '',i4,1x,A)')ens_index,trim(varstring)
if (ens_index > ens_size) then
- write(string1,*)'believed to have ',ens_size,'ensemble members for observation operator.'
+ write(string1,*)'Known to have ',ens_size,'ensemble members for observation operator.'
write(string2,*)'asking to use operator for ensemble member ',ens_index
- call error_handler(E_ERR, 'get_scalar_from_3Dhistory', &
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_3Dhistory', &
string1, source, revision, revdate, text2=string2)
endif
+!----------------------------------------------------------------------
! bombproofing ... make sure the netcdf file is open.
-write(*,*)'ncid(',ens_index,') is ',ncid(ens_index)
call nc_check(nf90_inquire(ncid(ens_index)), &
- 'get_scalar_from_3Dhistory', 'inquire '//trim(strshort))
+ 'obs_def_tower.get_scalar_from_3Dhistory', 'inquire '//trim(strshort))
! bombproofing ... make sure the variable is the shape and size we expect
call nc_check(nf90_inq_varid(ncid(ens_index), trim(varstring), varid), &
- 'get_scalar_from_3Dhistory', 'inq_varid '//trim(strshort))
+ 'obs_def_tower.get_scalar_from_3Dhistory', 'inq_varid '//trim(strshort))
call nc_check(nf90_inquire_variable(ncid(ens_index), varid, xtype=xtype, ndims=ndims, &
- dimids=dimids, natts=natts),'get_scalar_from_3Dhistory','inquire variable '//trim(strshort))
+ dimids=dimids, natts=natts), &
+ 'obs_def_tower.get_scalar_from_3Dhistory','inquire variable '//trim(strshort))
if (ndims /= 3) then
write(string1,*)trim(varstring),' is supposed to have 3 dimensions, it has',ndims
- call error_handler(E_ERR, 'get_scalar_from_3Dhistory', &
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_3Dhistory', &
string1, source, revision, revdate)
endif
! If the variable is not a NF90_FLOAT, then the assumptions for processing
! the missing_value, _FillValue, etc., may not be correct.
if (xtype /= NF90_FLOAT) then
- write(string1,*)trim(varstring),' is supposed to be a 32 bit real. xtype = ',NF90_FLOAT,' it is ',xtype
- call error_handler(E_ERR, 'get_scalar_from_3Dhistory', &
+ write(string1,*)trim(varstring),' is supposed to be a 32 bit real. xtype = ', &
+ NF90_FLOAT,' it is ',xtype
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_3Dhistory', &
string1, source, revision, revdate)
endif
! Dimension 1 is longitude
call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(1), len=dimlen), &
- 'get_scalar_from_3Dhistory', 'inquire_dimension 1 '//trim(strshort))
+ 'obs_def_tower.get_scalar_from_3Dhistory', 'inquire_dimension 1 '//trim(strshort))
if (dimlen /= nlon) then
write(string1,*)'LON has length',nlon,trim(varstring),' has ',dimlen,'longitudes.'
- call error_handler(E_ERR, 'get_scalar_from_3Dhistory', &
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_3Dhistory', &
string1, source, revision, revdate)
endif
! Dimension 2 is latitude
call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(2), len=dimlen), &
- 'get_scalar_from_3Dhistory', 'inquire_dimension 2 '//trim(strshort))
+ 'obs_def_tower.get_scalar_from_3Dhistory', 'inquire_dimension 2 '//trim(strshort))
if (dimlen /= nlat) then
write(string1,*)'LAT has length',nlat,trim(varstring),' has ',dimlen,'latitudes.'
- call error_handler(E_ERR, 'get_scalar_from_3Dhistory', &
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_3Dhistory', &
string1, source, revision, revdate)
endif
! Dimension 3 is time
call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(3), len=dimlen), &
- 'get_scalar_from_3Dhistory', 'inquire_dimension 3'//trim(strshort))
+ 'obs_def_tower.get_scalar_from_3Dhistory', 'inquire_dimension 3'//trim(strshort))
if (dimlen /= ntime) then
write(string1,*)'TIME has length',ntime,trim(varstring),' has ',dimlen,'times.'
- call error_handler(E_ERR, 'get_scalar_from_3Dhistory', &
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_3Dhistory', &
string1, source, revision, revdate)
endif
-! Find the grid cell and timestep of interest
+!----------------------------------------------------------------------
+! Find the grid cell and timestep of interest
+! FIXME ... since the history file contents are for the previous 30m,
+! perhaps the closest time is not the best approximation.
! Get the individual locations values
call get_time(obs_time, second, day)
@@ -449,34 +586,36 @@
timei = timeinds(1)
if (debug .and. do_output()) then
- write(*,*)'get_scalar_from_3Dhistory:targetlon, lon, lon index is ', &
+ write(*,*)'obs_def_tower.get_scalar_from_3Dhistory:targetlon, lon, lon index is ', &
loc_lon,lon(gridloni),gridloni
- write(*,*)'get_scalar_from_3Dhistory:targetlat, lat, lat index is ', &
+ write(*,*)'obs_def_tower.get_scalar_from_3Dhistory:targetlat, lat, lat index is ', &
loc_lat,lat(gridlatj),gridlatj
- write(*,*)'get_scalar_from_3Dhistory: targetT, T, T index is ', &
+ write(*,*)'obs_def_tower.get_scalar_from_3Dhistory: targetT, T, T index is ', &
otime,rtime(timei),timei
endif
if ( abs(otime - rtime(timei)) > 30*60 ) then
if (debug .and. do_output()) then
- write(*,*)'get_scalar_from_3Dhistory: no close time ... skipping observation'
- call print_time(obs_time,'get_scalar_from_3Dhistory:observation time')
- call print_date(obs_time,'get_scalar_from_3Dhistory:observation date')
+ write(*,*)'obs_def_tower.get_scalar_from_3Dhistory: no close time ... skipping observation'
+ call print_time(obs_time,'obs_def_tower.get_scalar_from_3Dhistory:observation time')
+ call print_date(obs_time,'obs_def_tower.get_scalar_from_3Dhistory:observation date')
endif
istatus = 2
return
endif
+!----------------------------------------------------------------------
! Grab exactly the scalar we want.
ncstart = (/ gridloni, gridlatj, timei /)
nccount = (/ 1, 1, 1 /)
-call nc_check(nf90_get_var(ncid(ens_index), varid, hyperslab, start=ncstart, count=nccount), &
- 'get_scalar_from_3Dhistory', 'get_var')
+call nc_check(nf90_get_var(ncid(ens_index),varid,hyperslab,start=ncstart,count=nccount), &
+ 'obs_def_tower.get_scalar_from_3Dhistory', 'get_var')
obs_val = hyperslab(1)
+!----------------------------------------------------------------------
! Apply any netCDF attributes ...
io1 = nf90_get_att(ncid(ens_index), varid, '_FillValue' , spvalR4)
@@ -501,6 +640,265 @@
end subroutine get_scalar_from_3Dhistory
+!======================================================================
+
+
+subroutine get_scalar_from_2Dhistory(varstring, ens_index, location, obs_time, &
+ obs_val, istatus)
+! the routine must return values for:
+! obs_val -- the computed forward operator value
+! istatus -- return code: 0=ok, > 0 is error, < 0 reserved for system use
+!
+! The requirement is that the history file variable is a 2D variable shaped similarly:
+!
+! float NEP(time, lndgrid) ;
+! NEP:long_name = "net ecosystem production, blah, blah, blah" ;
+! NEP:units = "gC/m^2/s" ;
+! NEP:cell_methods = "time: mean" ;
+! NEP:_FillValue = 1.e+36f ;
+! NEP:missing_value = 1.e+36f ;
+!
+! Just because it is 2D does not mean it is a single column,
+! although single columns are all that is really supported right now.
+
+character(len=*), intent(in) :: varstring
+integer, intent(in) :: ens_index
+type(location_type), intent(in) :: location
+type(time_type), intent(in) :: obs_time
+real(r8), intent(out) :: obs_val
+integer, intent(out) :: istatus
+
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimids
+real(r8), dimension(3) :: loc
+integer, dimension(2) :: ncstart, nccount
+integer, dimension(1) :: timeinds
+integer :: gridij, timei
+integer :: varid, xtype, ndims, natts, dimlen
+integer :: io1, io2, second, day
+real(r8) :: loc_lon, loc_lat, radius, distance
+real(r4), dimension(1) :: hyperslab
+real(r4) :: spvalR4
+real(r8) :: scale_factor, add_offset
+real(digits12) :: otime
+character(len=NF90_MAX_NAME+20) :: strshort
+
+type(location_type) :: gridloc
+
+obs_val = MISSING_R8
+istatus = 1
+
+!----------------------------------------------------------------------
+! if observation is outside region encompassed in the history file - fail
+loc = get_location(location) ! loc is in DEGREES
+loc_lon = loc(1)
+loc_lat = loc(2)
+
+! Defining the region if running in an unstructured grid is tricky.
+! We have lat, lon, and the area of the gridcell which we assume to be basically square.
+! The square root of the area defines the length of the edge of the gridcell.
+! Half the hypotenuse defines the radius of a circle. Any ob within
+! that radius is close enough.
+
+! TJH FIXME This does not work with unstructured grid.
+! latinds = minloc(abs(lat - loc_lat)) ! these return 'arrays' ...
+! loninds = minloc(abs(lon - loc_lon)) ! these return 'arrays' ...
+! gridij = the closest location
+
+! The "1" in the following reflect the fact that only a single gridcell
+! is currently supported in the unstructured grid configuration.
+! (see GetDimensions())
+
+gridij = 1
+gridloc = set_location(lon(gridij),lat(gridij), 0.0_r8, VERTISUNDEF)
+distance = get_dist(gridloc, location, no_vert = .TRUE.) * RAD2KM
+radius = sqrt(2.0_r8 * area(gridij))/2.0_r8
+
+if (debug .and. do_output()) then
+ write(string1,*)' observation lon, lat is ',loc_lon, loc_lat
+ write(string2,*)'gridcell lon, lat is ',lon(gridij),lat(gridij)
+ write(string3,*)'area,radius is ',area(gridij),radius,' distance ',distance
+ call error_handler(E_MSG, 'obs_def_tower.get_scalar_from_2Dhistory', &
+ string1, source, revision, revdate, text2=string2, text3=string3)
+endif
+
+if ( distance > radius ) return
+
+!----------------------------------------------------------------------
+! Now that we know the observation operator is possible, continue ...
+
+write(strshort,'(''ens_index '',i4,1x,A)')ens_index,trim(varstring)
+
+if (ens_index > ens_size) then
+ write(string1,*)'believed to have ',ens_size,'ensemble members for observation operator.'
+ write(string2,*)'asking to use operator for ensemble member ',ens_index
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_2Dhistory', &
+ string1, source, revision, revdate, text2=string2)
+endif
+
+!----------------------------------------------------------------------
+! bombproofing ... make sure the netcdf file is open.
+
+call nc_check(nf90_inquire(ncid(ens_index)), &
+ 'obs_def_tower.get_scalar_from_2Dhistory', 'inquire '//trim(strshort))
+
+! bombproofing ... make sure the variable is the shape and size we expect
+
+call nc_check(nf90_inq_varid(ncid(ens_index), trim(varstring), varid), &
+ 'obs_def_tower.get_scalar_from_2Dhistory', 'inq_varid '//trim(strshort))
+call nc_check(nf90_inquire_variable(ncid(ens_index), varid, xtype=xtype, ndims=ndims, &
+ dimids=dimids, natts=natts), &
+ 'obs_def_tower.get_scalar_from_2Dhistory','inquire variable '//trim(strshort))
+
+if (ndims /= 2) then
+ write(string1,*)trim(varstring),' is supposed to have 2 dimensions, it has',ndims
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_2Dhistory', &
+ string1, source, revision, revdate)
+endif
+
+! If the variable is not a NF90_FLOAT, then the assumptions for processing
+! the missing_value, _FillValue, etc., may not be correct.
+if (xtype /= NF90_FLOAT) then
+ write(string1,*)trim(varstring),' is supposed to be a 32 bit real. xtype = ', &
+ NF90_FLOAT,' it is ',xtype
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_2Dhistory', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 1 is spatial
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(1), len=dimlen), &
+ 'obs_def_tower.get_scalar_from_2Dhistory', 'inquire_dimension 1 '//trim(strshort))
+if (dimlen /= nlon) then
+ write(string1,*)'LON has length',nlon,trim(varstring),' has ',dimlen,'longitudes.'
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_2Dhistory', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 2 is time
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(2), len=dimlen), &
+ 'obs_def_tower.get_scalar_from_2Dhistory', 'inquire_dimension 2 '//trim(strshort))
+if (dimlen /= ntime) then
+ write(string1,*)'TIME has length',ntime,trim(varstring),' has ',dimlen,'times.'
+ call error_handler(E_ERR, 'obs_def_tower.get_scalar_from_2Dhistory', &
+ string1, source, revision, revdate)
+endif
+
+!----------------------------------------------------------------------
+! Find the timestep of interest
+! FIXME ... since the history file contents are for the previous 30m,
+! perhaps the closest time is not the best approximation.
+
+call get_time(obs_time, second, day)
+otime = real(day,digits12) + real(second,digits12)/86400.0_digits12
+timeinds = minloc(abs(rtime - otime)) ! these return 'arrays' ...
+timei = timeinds(1)
+
+if (debug .and. do_output()) then
+ write(*,*)'obs_def_tower.get_scalar_from_2Dhistory: targetT, T, T index is ', &
+ otime,rtime(timei),timei
+endif
+
+if ( abs(otime - rtime(timei)) > 30*60 ) then
+ if (debug .and. do_output()) then
+ write(*,*)'obs_def_tower.get_scalar_from_2Dhistory: no close time ... skipping observation'
+ call print_time(obs_time,'obs_def_tower.get_scalar_from_2Dhistory:observation time')
+ call print_date(obs_time,'obs_def_tower.get_scalar_from_2Dhistory:observation date')
+ endif
+ istatus = 2
+ return
+endif
+
+!----------------------------------------------------------------------
+! Grab exactly the scalar we want.
+
+ncstart = (/ gridij, timei /)
+nccount = (/ 1, 1 /)
+
+call nc_check(nf90_get_var(ncid(ens_index),varid,hyperslab,start=ncstart,count=nccount), &
+ 'obs_def_tower.get_scalar_from_2Dhistory', 'get_var')
+
+obs_val = hyperslab(1)
+
+!----------------------------------------------------------------------
+! Apply any netCDF attributes ...
+
+io1 = nf90_get_att(ncid(ens_index), varid, '_FillValue' , spvalR4)
+if ((io1 == NF90_NOERR) .and. (hyperslab(1) == spvalR4)) obs_val = MISSING_R8
+
+io2 = nf90_get_att(ncid(ens_index), varid, 'missing_value' , spvalR4)
+if ((io2 == NF90_NOERR) .and. (hyperslab(1) == spvalR4)) obs_val = MISSING_R8
+
+io1 = nf90_get_att(ncid(ens_index), varid, 'scale_factor', scale_factor)
+io2 = nf90_get_att(ncid(ens_index), varid, 'add_offset' , add_offset)
+
+if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val * scale_factor + add_offset
+elseif (io1 == NF90_NOERR) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val * scale_factor
+elseif (io2 == NF90_NOERR) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val + add_offset
+endif
+
+if (obs_val /= MISSING_R8) istatus = 0
+
+end subroutine get_scalar_from_2Dhistory
+
+
+!======================================================================
+
+
+subroutine test_block
+! Defining the region if running in a single column is tricky.
+! We have lat, lon, and the area of the gridcell which we assume to be basically square.
+! The square root of the area defines the length of the edge of the gridcell which
+! can then be interpreted as the diameter of a circle. Any observation within this
+! distance is close enough.
+
+real(r8) :: x1, x2, y1, y2, d1, d2, d3, d4, radius, distance
+type(location_type) :: gridloc, testloc
+
+gridloc = set_location(0.0_r8, 0.0_r8, 0.0_r8, VERTISUNDEF)
+testloc = set_location(1.0_r8, 0.0_r8, 0.0_r8, VERTISUNDEF)
+distance = get_dist(gridloc, testloc, no_vert = .TRUE.) * RAD2KM
+write(*,*)'TJH DEBUG: 1 degree at the equator has distance ',distance,' km.'
+
+x1 = 286.8750 ! gridcell longitude is 287.50
+x2 = 288.1250
+y1 = 42.40837669372559 ! gridcell latitude is 42.8795814514160
+y2 = 43.35078620910645
+
+! compute the distance along the top of the grid cell
+gridloc = set_location(x1, y1, 0.0_r8, VERTISUNDEF)
+testloc = set_location(x2, y1, 0.0_r8, VERTISUNDEF)
+d1 = get_dist(gridloc, testloc, no_vert = .TRUE.) * RAD2KM
+gridloc = set_location(x2, y2, 0.0_r8, VERTISUNDEF)
+d2 = get_dist(gridloc, testloc, no_vert = .TRUE.) * RAD2KM
+testloc = set_location(x1, y2, 0.0_r8, VERTISUNDEF)
+d3 = get_dist(gridloc, testloc, no_vert = .TRUE.) * RAD2KM
+gridloc = set_location(x1, y1, 0.0_r8, VERTISUNDEF)
+d4 = get_dist(gridloc, testloc, no_vert = .TRUE.) * RAD2KM
+
+write(*,*)
+write(*,*)'lengths are ',d1,d2,d3,d4
+write(*,*)'rectangular area is ',((d1+d3)/2.0_r8)*((d2+d4)/2.0_r8)
+
+d1 = sqrt(area(1))
+d2 = sqrt(2.0_r8 * d1**2)
+d3 = d2 / 2.0_r8
+radius = sqrt(2.0_r8 * area(1))/2.0_r8
+write(*,*)'length of one side of a square is x = sqrt(area) = ',d1
+write(*,*)'diagonal is sqrt(x**2 + x**2) = ',d2
+write(*,*)'radius is sqrt(x**2 + x**2)/2.0 = ',d3
+write(*,*)'-or- sqrt(area + area)/2.0 = ',radius
+
+write(*,*)'TJH DEBUG: Radius,area is ',distance,PI*distance**2, &
+ ' gridcell area is ',area(1)
+
+end subroutine test_block
+
+
+!======================================================================
+
+
end module obs_def_tower_mod
! END DART PREPROCESS MODULE CODE
Added: DART/branches/development/obs_def/obs_def_tower_mod.nml
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list