[Dart-dev] [5760] DART/branches/development: level4_to_obs: Applied conversion from umol to gC correctly.
nancy at ucar.edu
nancy at ucar.edu
Tue Jun 12 11:24:22 MDT 2012
Revision: 5760
Author: thoar
Date: 2012-06-12 11:24:22 -0600 (Tue, 12 Jun 2012)
Log Message:
-----------
level4_to_obs: Applied conversion from umol to gC correctly.
makedaily: The tower obs filenames are obs_def.daybefore.YYYYMMDD and actually contain
the obs for the day BEFORE YYYYMMDD ... the model stops at YYYYMMDD, but
we need the obs for the day before to match the CLM history files.
[00:00 to 23:59] from day BEFORE goes into tower observation files to match
the CLM history file generation.
The assimilation window must be > 2 days
in order to assimilate the obs at 00Z of the day 'before'.
DEFAULT_obs_kind_mod.F90 has some of the obvious kinds needed for CLM bgc.
obs_def_tower_mod: now supports reading CLM history files to get the NEP 'observation'.
opens the CLM history files for each ensemble member in the init routine.
Requires that the history file be created with something like:
&clm_inparm
hist_empty_htapes = .true.
hist_fincl1 = 'NEP'
hist_nhtfrq = 1,
hist_mfilt = 48
hist_avgflag_pertape = 'A'
/
Needs the model time from the model_mod routine to be able to
create the history file name FOR THE PREVIOUS DAY.
Still need to check the metadata in the other CLM history files
for the obs operator code.
Modified Paths:
--------------
DART/branches/development/models/clm/shell_scripts/makedaily.sh
DART/branches/development/obs_def/obs_def_tower_mod.f90
DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90
DART/branches/development/observations/Ameriflux/level4_to_obs.f90
DART/branches/development/observations/Ameriflux/work/path_names_level4_to_obs
DART/branches/development/observations/Ameriflux/work/path_names_obs_sequence_tool
-------------- next part --------------
Modified: DART/branches/development/models/clm/shell_scripts/makedaily.sh
===================================================================
--- DART/branches/development/models/clm/shell_scripts/makedaily.sh 2012-06-12 17:04:54 UTC (rev 5759)
+++ DART/branches/development/models/clm/shell_scripts/makedaily.sh 2012-06-12 17:24:22 UTC (rev 5760)
@@ -1,8 +1,8 @@
#!/bin/bash
# BLUEFIRE /usr/local/bin/bash
-# split a yearly file into "daily" files which start at 12:01Z
-# the previous day and end at 12:00Z on the day that matches the
+# split a yearly file into "daily" files which start at 00:00Z
+# the previous day and end at 23:59Z on the day that matches the
# day in the filename.
# set the first and last days to be split. can roll over
@@ -26,7 +26,7 @@
# to set the initial day while we are doing the total day calc.
# make sure there is an initial input.nml for advance_time
-# cp -f input.nml.template input.nml
+cp -f ../work/input.nml.template input.nml || exit 1
# these outputs from advance time (with the -g flag) are
# 2 integers: gregorian_day_number seconds
@@ -35,47 +35,44 @@
mon2=`printf %02d $end_month`
day2=`printf %02d $end_day`
end_d=(`echo ${end_year}${mon2}${day2}00 0 -g | ../work/advance_time`)
- echo $end_d
+ echo last day is day $end_d
mon2=`printf %02d $start_month`
day2=`printf %02d $start_day`
start_d=(`echo ${start_year}${mon2}${day2}00 0 -g | ../work/advance_time`)
- echo $start_d
+ echo first day is day $start_d
# these are a string in the format YYYYMMDDHH
# do them here to prime the loop below which first takes them apart.
+prevday=(`echo ${start_year}${mon2}${day2}00 -1d | ../work/advance_time`)
currday=(`echo ${start_year}${mon2}${day2}00 0 | ../work/advance_time`)
nextday=(`echo ${start_year}${mon2}${day2}00 +1d | ../work/advance_time`)
-prevday=(`echo ${start_year}${mon2}${day2}00 -1d | ../work/advance_time`)
# how many total days are going to be split (for the loop counter)
# (pull out the first of the 2 numbers which are output from advance_time)
let totaldays=${end_d}-${start_d}+1
- echo $totaldays
# loop over each day
let d=1
while (( d <= totaldays)) ; do
-echo "subsetting $d of $totaldays ..."
-#echo $currday $nextday
+ echo "subsetting $d of $totaldays ..."
+ #echo $currday $nextday
# parse out the parts from a string which is YYYYMMDDHH
- # both for the current day, yesterday, and tomorrow
+ # for yesterday(previous), today(current), and tomorrow(next)
+ pyear=${prevday:0:4}
+ pmonth=${prevday:4:2}
+ pday=${prevday:6:2}
+
cyear=${currday:0:4}
cmonth=${currday:4:2}
cday=${currday:6:2}
+
nyear=${nextday:0:4}
nmonth=${nextday:4:2}
nday=${nextday:6:2}
- pyear=${prevday:0:4}
- pmonth=${prevday:4:2}
- pday=${prevday:6:2}
-#echo curr $cyear $cmonth $cday
-#echo next $nyear $nmonth $nday
-#echo prev $pyear $pmonth $pday
-
# compute the equivalent gregorian days here.
g=(`echo ${cyear}${cmonth}${cday}00 -1d -g | ../work/advance_time`)
greg0=${g[0]}
@@ -84,20 +81,12 @@
g=(`echo ${cyear}${cmonth}${cday}00 +1d -g | ../work/advance_time`)
greg2=${g[0]}
- echo starting WOD obs for ${cyear}${cmonth}${cday} gregorian= $greg1
+ echo prev $pyear $pmonth $pday which is gregorian $greg0
+ echo curr $cyear $cmonth $cday which is gregorian $greg1
+ echo next $nyear $nmonth $nday which is gregorian $greg2
-# I have one whopping long obs_seq.in containing 3 obs per day for
-# 4000 days (all taken at midnight) for a perfect model experiment.
-# I don't need to try to determine the right input files - TJH
-#
-# # Determine the pertinent input files.
-# # last 12 hrs yesterdays data plus the first 12 hrs of todays
-# if [[ ${cmonth} == '01' && ${cday} == '01' ]] ; then
-# echo "../yearly_files/obs_seq${pyear}.wod" > olist
-# echo "../yearly_files/obs_seq${cyear}.wod" >> olist
-# else
-# echo "../yearly_files/obs_seq${cyear}.wod" > olist
-# fi
+ # I have annual files ...
+ # I'll need to revisit this when I wrap over year boundaries ... TJH
sed -e "s/YYYY/${cyear}/g" \
-e "s/MM/${cmonth}/g" \
@@ -107,8 +96,6 @@
-e "s/GREG1/${greg1}/g" \
-e "s/GREG2/${greg2}/g" < ../work/input.nml.template > input.nml
-# cat input.nml
-
# make sure output dir exists
if [[ ! -d ../${cyear}${cmonth} ]] ; then
mkdir ../${cyear}${cmonth}
@@ -121,11 +108,10 @@
prevday=$currday
currday=$nextday
nextday=(`echo ${nyear}${nmonth}${nday}00 +1d | ../work/advance_time`)
-#echo $currday $nextday $prevday
+ #echo $currday $nextday $prevday
# advance the loop counter
let d=d+1
-#echo d=$d
done
Modified: DART/branches/development/obs_def/obs_def_tower_mod.f90
===================================================================
--- DART/branches/development/obs_def/obs_def_tower_mod.f90 2012-06-12 17:04:54 UTC (rev 5759)
+++ DART/branches/development/obs_def/obs_def_tower_mod.f90 2012-06-12 17:24:22 UTC (rev 5760)
@@ -33,11 +33,11 @@
!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
! case(TOWER_LATENT_HEAT_FLUX)
-! call get_expected_latent_heat_flux(state, state_time, ens_index, location, obs_def%key, obs_val, istatus)
+! call get_expected_latent_heat_flux(state, state_time, ens_index, location, obs_time, obs_key, obs_val, istatus)
! case(TOWER_SENSIBLE_HEAT_FLUX)
-! call get_expected_sensible_heat_flux(state, state_time, ens_index, location, obs_def%key, obs_val, istatus)
+! call get_expected_sensible_heat_flux(state, state_time, ens_index, location, obs_time, obs_key, obs_val, istatus)
! case(TOWER_NETC_ECO_EXCHANGE)
-! call get_expected_net_C_production(state, state_time, ens_index, location, obs_def%key, obs_val, istatus)
+! call get_expected_net_C_production(state, state_time, ens_index, location, obs_time, obs_key, obs_val, istatus)
! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
!-----------------------------------------------------------------------------
@@ -72,13 +72,16 @@
! $Revision$
! $Date$
-use types_mod, only : r4, r8, MISSING_R8, PI, deg2rad
+use types_mod, only : r4, r8, digits12, MISSING_R8, PI, deg2rad
use location_mod, only : location_type, get_location
-use time_manager_mod, only : time_type
+use time_manager_mod, only : time_type, get_date, set_date, print_date, print_time, &
+ get_time, set_time, operator(-)
+use model_mod, only : get_model_time
use utilities_mod, only : register_module, E_ERR, E_MSG, error_handler, &
check_namelist_read, find_namelist_in_file, &
nmlfileunit, do_output, do_nml_file, do_nml_term, &
- nc_check
+ nc_check, file_exist
+
use typesizes
use netcdf
@@ -95,20 +98,22 @@
revision = "$Revision$", &
revdate = "$Date$"
-logical :: module_initialized = .false.
+logical :: module_initialized = .false.
+character(len=129) :: string1, string2, string3
+integer :: nlon, nlat, ntime, ens_size
-character(len=129) :: string1,string2,string3
-
+character(len=129), allocatable, dimension(:) :: fname
integer, allocatable, dimension(:) :: ncid
-character(len=129), allocatable, dimension(:) :: fname
-real(r8), allocatable, dimension(:) :: lon, lat, time
-integer :: nlon, nlat, ntime, ens_size
+real(r8), allocatable, dimension(:) :: lon, lat
+real(digits12), allocatable, dimension(:) :: rtime
+
! namelist items
character(len=256) :: casename = 'clm_tim'
-logical :: debug = .true.
+logical :: verbose = .false.
+logical :: debug = .false.
-namelist /obs_def_tower_nml/ casename, debug
+namelist /obs_def_tower_nml/ casename, verbose, debug
contains
@@ -120,10 +125,14 @@
subroutine initialize_module
-! Called once to set values and allocate space
+! Called once to set values and allocate space, open all the CLM files
+! that have the observations, etc.
integer :: iunit, io, rc, i
integer :: dimid, varid
+integer :: year, month, day, hour, minute, second, leftover
+integer, allocatable, dimension(:) :: yyyymmdd,sssss
+type(time_type) :: model_time
! Prevent multiple calls from executing this code more than once.
if (module_initialized) return
@@ -142,71 +151,148 @@
if (do_nml_file()) write(nmlfileunit, nml=obs_def_tower_nml)
if (do_nml_term()) write( * , nml=obs_def_tower_nml)
-! FIXME: Figure out how many files (i.e. ensemble size)
-ens_size = 4
+! 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 data for [23:30 -> 00:00] get put in the file for the next day.
+model_time = get_model_time()
+model_time = model_time - set_time(0,1)
+call get_date(model_time, year, month, day, hour, minute, second)
+second = second + minute*60 + hour*3600
+
+! 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 ...
+
+100 format (A,'.clm2_',I4.4,'.h0.',I4.4,'-',I2.2,'-',I2.2,'-',I5.5,'.nc')
+
+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.'
+ ens_size = ens_size + 1
+ else
+ if(verbose .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)
+
+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, &
+ source, revision, revdate, text2=string3)
+
+else
+
+ if (verbose .and. do_output()) write(*,*)'Ensemble size is believed to be ',ens_size
+
+endif
+
allocate(fname(ens_size),ncid(ens_size))
+ncid = 0
ENSEMBLE : do i = 1,ens_size
- write(fname(i),'(A,''.clm2_'',I4.4,''.h.nc'')')trim(casename),i
+ 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)))
enddo ENSEMBLE
-ncid = 0
-
i = 1
-if (debug) write(*,*)'Reading observations from ',trim(fname(i))
! Harvest information from the first observation file.
! FIXME All other files will be opened to make sure they have the same dimensions.
-if (ncid(i) == 0) then ! we need to open it
- call nc_check(nf90_open(trim(fname(i)), nf90_nowrite, ncid(i)), &
- 'obs_def_tower_mod:initialize_routine','open '//trim(fname(i)))
-endif
-
call nc_check(nf90_inq_dimid(ncid(i), 'lon', dimid), &
- 'obs_def_tower_mod:initialize_routine','inq_dimid lon '//trim(fname(i)))
+ 'initialize_routine','inq_dimid lon '//trim(fname(i)))
call nc_check(nf90_inquire_dimension(ncid(i), dimid, len=nlon), &
- 'obs_def_tower_mod:initialize_routine','inquire_dimension lon '//trim(fname(i)))
+ 'initialize_routine','inquire_dimension lon '//trim(fname(i)))
call nc_check(nf90_inq_dimid(ncid(i), 'lat', dimid), &
- 'obs_def_tower_mod:initialize_routine','inq_dimid lat '//trim(fname(i)))
+ 'initialize_routine','inq_dimid lat '//trim(fname(i)))
call nc_check(nf90_inquire_dimension(ncid(i), dimid, len=nlat), &
- 'obs_def_tower_mod:initialize_routine','inquire_dimension lat '//trim(fname(i)))
+ 'initialize_routine','inquire_dimension lat '//trim(fname(i)))
call nc_check(nf90_inq_dimid(ncid(i), 'time', dimid), &
- 'obs_def_tower_mod:initialize_routine','inq_dimid time '//trim(fname(i)))
+ 'initialize_routine','inq_dimid time '//trim(fname(i)))
call nc_check(nf90_inquire_dimension(ncid(i), dimid, len=ntime), &
- 'obs_def_tower_mod:initialize_routine','inquire_dimension lat '//trim(fname(i)))
+ 'initialize_routine','inquire_dimension lat '//trim(fname(i)))
-allocate(lon(nlon),lat(nlat),time(ntime))
+allocate(lon(nlon),lat(nlat))
+allocate(rtime(ntime),yyyymmdd(ntime),sssss(ntime))
call nc_check(nf90_inq_varid(ncid(i), 'lon', varid), &
- 'obs_def_tower_mod:initialize_routine','inq_varid lon '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, lon), &
- 'obs_def_tower_mod:initialize_routine', 'get_var')
+ '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), &
- 'obs_def_tower_mod:initialize_routine','inq_varid lat '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, lat), &
- 'obs_def_tower_mod:initialize_routine', 'get_var')
+ '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), 'time', varid), &
- 'obs_def_tower_mod:initialize_routine','inq_varid time '//trim(fname(i)))
-call nc_check(nf90_get_var(ncid(i), varid, time), &
- 'obs_def_tower_mod:initialize_routine', 'get_var')
+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')
-if (debug) write(*,*)'obs_def_tower lon',lon
-if (debug) write(*,*)'obs_def_tower lat',lat
-if (debug) write(*,*)'obs_def_tower time',time
+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), '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
+
+ year = yyyymmdd(i)/10000
+ leftover = yyyymmdd(i) - year*10000
+ month = leftover/100
+ day = leftover - month*100
+
+ hour = sssss(i)/3600
+ leftover = sssss(i) - hour*3600
+ minute = leftover/60
+ second = leftover - minute*60
+
+ model_time = set_date(year, month, day, hour, minute, second)
+ call get_time(model_time, second, day)
+
+ rtime(i) = real(day,digits12) + real(second,digits12)/86400.0_digits12
+
+ if (debug .and. do_output()) then
+ write(*,*)'timestep yyyymmdd sssss',i,yyyymmdd(i),sssss(i)
+ call print_date(model_time,'tower_mod date')
+ call print_time(model_time,'tower_mod time')
+ write(*,*)'tower_mod time as a real ',rtime(i)
+ endif
+
+enddo
+
+if (debug .and. do_output()) write(*,*)'obs_def_tower lon',lon
+if (debug .and. do_output()) write(*,*)'obs_def_tower lat',lat
+
! FIXME
! check all other ensemble member history files to make sure metadata is the same.
+deallocate(yyyymmdd, sssss)
+
end subroutine initialize_module
-subroutine get_expected_latent_heat_flux(state, state_time, ens_index, location, obs_key, obs_val, istatus)
+subroutine get_expected_latent_heat_flux(state, state_time, ens_index, location, obs_time, obs_key, 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
@@ -215,6 +301,7 @@
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
integer, intent(in) :: obs_key
real(r8), intent(out) :: obs_val
integer, intent(out) :: istatus
@@ -232,7 +319,7 @@
-subroutine get_expected_sensible_heat_flux(state, state_time, ens_index, location, obs_key, obs_val, istatus)
+subroutine get_expected_sensible_heat_flux(state, state_time, ens_index, location, obs_time, obs_key, 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
@@ -240,6 +327,7 @@
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
integer, intent(in) :: obs_key
real(r8), intent(out) :: obs_val
integer, intent(out) :: istatus
@@ -256,22 +344,23 @@
end subroutine get_expected_sensible_heat_flux
-subroutine get_expected_net_C_production(state, state_time, ens_index, location, obs_key, obs_val, istatus)
+subroutine get_expected_net_C_production(state, state_time, ens_index, location, obs_time, obs_key, 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
!
-! float NEE(time, lat, lon) ;
-! NEE:long_name = "net ecosystem exchange of carbon, blah, blah, blah" ;
-! NEE:units = "gC/m^2/s" ;
-! NEE:cell_methods = "time: mean" ;
-! NEE:_FillValue = 1.e+36f ;
-! NEE:missing_value = 1.e+36f ;
+! float NEP(time, lat, lon) ;
+! 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 ;
real(r8), intent(in) :: state(:)
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
integer, intent(in) :: obs_key
real(r8), intent(out) :: obs_val
integer, intent(out) :: istatus
@@ -279,20 +368,24 @@
integer, dimension(NF90_MAX_VAR_DIMS) :: dimids
real(r8), dimension(3) :: loc
integer, dimension(3) :: ncstart, nccount
-integer, dimension(1) :: loninds, latinds
-integer :: gridloni, gridlatj
+integer, dimension(1) :: loninds, latinds, timeinds
+integer :: gridloni, gridlatj, timei
integer :: varid, xtype, ndims, natts, dimlen
-integer :: io1, io2
+integer :: io1, io2, second, day
real(r8) :: loc_lon, loc_lat
real(r4), dimension(1) :: hyperslab
real(r4) :: spvalR4
real(r8) :: scale_factor, add_offset
+real(digits12) :: otime
+character(len=20) :: strshort
if ( .not. module_initialized ) call initialize_module
obs_val = MISSING_R8
istatus = 1
+write(strshort,'(''ens_index '',i4)')ens_index
+
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
@@ -302,15 +395,19 @@
! 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_expected_net_C_production', 'inquire '//trim(strshort))
+
! bombproofing ... make sure the variable is the shape and size we expect
-call nc_check(nf90_inq_varid(ncid(ens_index), 'NEE', varid), &
- 'get_expected_net_C_production', 'inq_varid NEE ')
+call nc_check(nf90_inq_varid(ncid(ens_index), 'NEP', varid), &
+ 'get_expected_net_C_production', 'inq_varid NEP '//trim(strshort))
call nc_check(nf90_inquire_variable(ncid(ens_index), varid, xtype=xtype, ndims=ndims, &
- dimids=dimids, natts=natts),'get_expected_net_C_production','inquire variable NEE ')
+ dimids=dimids, natts=natts),'get_expected_net_C_production','inquire variable NEP '//trim(strshort))
if (ndims /= 3) then
- write(string1,*)'NEE is supposed to have 3 dimensions, it has',ndims
+ write(string1,*)'NEP is supposed to have 3 dimensions, it has',ndims
call error_handler(E_ERR, 'get_expected_net_C_production', &
string1, source, revision, revdate)
endif
@@ -318,59 +415,77 @@
! 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,*)'NEE supposed to be a 32 bit real. xtype = ',NF90_FLOAT,' it is ',xtype
+ write(string1,*)'NEP supposed to be a 32 bit real. xtype = ',NF90_FLOAT,' it is ',xtype
call error_handler(E_ERR, 'get_expected_net_C_production', &
string1, source, revision, revdate)
endif
! Dimension 1 is longitude
call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(1), len=dimlen), &
- 'get_expected_net_C_production', 'inquire_dimension NEE 1')
+ 'get_expected_net_C_production', 'inquire_dimension NEP 1'//trim(strshort))
if (dimlen /= nlon) then
- write(string1,*)'LON has length',nlon,'NEE has ',dimlen,'longitudes.'
+ write(string1,*)'LON has length',nlon,'NEP has ',dimlen,'longitudes.'
call error_handler(E_ERR, 'get_expected_net_C_production', &
string1, source, revision, revdate)
endif
! Dimension 2 is latitude
call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(2), len=dimlen), &
- 'get_expected_net_C_production', 'inquire_dimension NEE 2')
+ 'get_expected_net_C_production', 'inquire_dimension NEP 2'//trim(strshort))
if (dimlen /= nlat) then
- write(string1,*)'LAT has length',nlat,'NEE has ',dimlen,'latitudes.'
+ write(string1,*)'LAT has length',nlat,'NEP has ',dimlen,'latitudes.'
call error_handler(E_ERR, 'get_expected_net_C_production', &
string1, source, revision, revdate)
endif
! Dimension 3 is time
call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(3), len=dimlen), &
- 'get_expected_net_C_production', 'inquire_dimension NEE 3')
+ 'get_expected_net_C_production', 'inquire_dimension NEP 3'//trim(strshort))
if (dimlen /= ntime) then
- write(string1,*)'TIME has length',ntime,'NEE has ',dimlen,'times.'
+ write(string1,*)'TIME has length',ntime,'NEP has ',dimlen,'times.'
call error_handler(E_ERR, 'get_expected_net_C_production', &
string1, source, revision, revdate)
endif
-! Find the grid cell of interest
+! Find the grid cell and timestep of interest
! Get the individual locations values
+call get_time(obs_time, second, day)
+otime = real(day,digits12) + real(second,digits12)/86400.0_digits12
loc = get_location(location) ! loc is in DEGREES
loc_lon = loc(1)
loc_lat = loc(2)
+
latinds = minloc(abs(lat - loc_lat)) ! these return 'arrays' ...
loninds = minloc(abs(lon - loc_lon)) ! these return 'arrays' ...
+timeinds = minloc(abs(rtime - otime)) ! these return 'arrays' ...
+
gridlatj = latinds(1)
gridloni = loninds(1)
+timei = timeinds(1)
if (debug .and. do_output()) then
write(*,*)'get_expected_net_C_production:targetlon, lon, lon index is ', &
loc_lon,lon(gridloni),gridloni
write(*,*)'get_expected_net_C_production:targetlat, lat, lat index is ', &
loc_lat,lat(gridlatj),gridlatj
+ write(*,*)'get_expected_net_C_production: targetT, T, T index is ', &
+ otime,rtime(timei),timei
endif
-! For now, just grab the last timestep ...
+if ( abs(otime - rtime(timei)) > 30*60 ) then
+ if (debug .and. do_output()) then
+ write(*,*)'get_expected_net_C_production: no close time ... skipping observation'
+ call print_time(obs_time,'get_expected_net_C_production:observation time')
+ call print_date(obs_time,'get_expected_net_C_production:observation date')
+ endif
+ istatus = 2
+ return
+endif
-ncstart = (/ gridloni, gridlatj, ntime /)
+! 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), &
Modified: DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90
===================================================================
--- DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90 2012-06-12 17:04:54 UTC (rev 5759)
+++ DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90 2012-06-12 17:24:22 UTC (rev 5760)
@@ -231,23 +231,33 @@
KIND_FLASH_RATE_2D = 106
! kinds for CLM - Community Land Model (Tim Hoar)
+! There will be lots of these ... perhaps leave some space ...
integer, parameter, public :: &
KIND_SNOW_THICKNESS = 107, &
KIND_SNOW_WATER = 108, &
KIND_SNOWCOVER_FRAC = 109, &
KIND_LIQUID_WATER = 110, &
KIND_ICE = 111, &
- KIND_LEAF_CARBON = 112, &
- KIND_LEAF_AREA_INDEX = 113, &
- KIND_NET_CARBON_FLUX = 114, &
- KIND_LATENT_HEAT_FLUX = 115, &
- KIND_SENSIBLE_HEAT_FLUX = 116, &
- KIND_RADIATION = 117, &
- KIND_NET_CARBON_PRODUCTION = 118
+ KIND_CARBON = 112, &
+ KIND_SOIL_CARBON = 113, &
+ KIND_ROOT_CARBON = 114, &
+ KIND_STEM_CARBON = 115, &
+ KIND_LEAF_CARBON = 116, &
+ KIND_LEAF_AREA_INDEX = 117, &
+ KIND_NET_CARBON_FLUX = 118, &
+ KIND_LATENT_HEAT_FLUX = 119, &
+ KIND_SENSIBLE_HEAT_FLUX = 120, &
+ KIND_RADIATION = 121, &
+ KIND_NET_CARBON_PRODUCTION = 122, &
+ KIND_NITROGEN = 123, &
+ KIND_SOIL_NITROGEN = 124, &
+ KIND_ROOT_NITROGEN = 125, &
+ KIND_STEM_NITROGEN = 126, &
+ KIND_LEAF_NITROGEN = 127
!! PRIVATE ONLY TO THIS MODULE. see comment below near the max_obs_specific
!! declaration.
-integer, parameter :: max_obs_generic = 118
+integer, parameter :: max_obs_generic = 127
!----------------------------------------------------------------------------
! This list is autogenerated by the 'preprocess' program. To add new
@@ -478,13 +488,22 @@
obs_kind_names(109) = obs_kind_type(KIND_SNOWCOVER_FRAC ,'KIND_SNOWCOVER_FRAC')
obs_kind_names(110) = obs_kind_type(KIND_LIQUID_WATER ,'KIND_LIQUID_WATER')
obs_kind_names(111) = obs_kind_type(KIND_ICE ,'KIND_ICE')
-obs_kind_names(112) = obs_kind_type(KIND_LEAF_CARBON ,'KIND_LEAF_CARBON')
-obs_kind_names(113) = obs_kind_type(KIND_LEAF_AREA_INDEX ,'KIND_LEAF_AREA_INDEX')
-obs_kind_names(114) = obs_kind_type(KIND_NET_CARBON_FLUX ,'KIND_NET_CARBON_FLUX')
-obs_kind_names(115) = obs_kind_type(KIND_LATENT_HEAT_FLUX ,'KIND_LATENT_HEAT_FLUX')
-obs_kind_names(116) = obs_kind_type(KIND_SENSIBLE_HEAT_FLUX ,'KIND_SENSIBLE_HEAT_FLUX')
-obs_kind_names(117) = obs_kind_type(KIND_RADIATION ,'KIND_RADIATION')
-obs_kind_names(118) = obs_kind_type(KIND_NET_CARBON_PRODUCTION ,'KIND_NET_CARBON_PRODUCTION')
+obs_kind_names(112) = obs_kind_type(KIND_CARBON ,'KIND_CARBON')
+obs_kind_names(113) = obs_kind_type(KIND_SOIL_CARBON ,'KIND_SOIL_CARBON')
+obs_kind_names(114) = obs_kind_type(KIND_ROOT_CARBON ,'KIND_ROOT_CARBON')
+obs_kind_names(115) = obs_kind_type(KIND_STEM_CARBON ,'KIND_STEM_CARBON')
+obs_kind_names(116) = obs_kind_type(KIND_LEAF_CARBON ,'KIND_LEAF_CARBON')
+obs_kind_names(117) = obs_kind_type(KIND_LEAF_AREA_INDEX ,'KIND_LEAF_AREA_INDEX')
+obs_kind_names(118) = obs_kind_type(KIND_NET_CARBON_FLUX ,'KIND_NET_CARBON_FLUX')
+obs_kind_names(119) = obs_kind_type(KIND_LATENT_HEAT_FLUX ,'KIND_LATENT_HEAT_FLUX')
+obs_kind_names(120) = obs_kind_type(KIND_SENSIBLE_HEAT_FLUX ,'KIND_SENSIBLE_HEAT_FLUX')
+obs_kind_names(121) = obs_kind_type(KIND_RADIATION ,'KIND_RADIATION')
+obs_kind_names(122) = obs_kind_type(KIND_NET_CARBON_PRODUCTION ,'KIND_NET_CARBON_PRODUCTION')
+obs_kind_names(123) = obs_kind_type(KIND_NITROGEN ,'KIND_NITROGEN')
+obs_kind_names(124) = obs_kind_type(KIND_SOIL_NITROGEN ,'KIND_SOIL_NITROGEN')
+obs_kind_names(125) = obs_kind_type(KIND_ROOT_NITROGEN ,'KIND_ROOT_NITROGEN')
+obs_kind_names(126) = obs_kind_type(KIND_STEM_NITROGEN ,'KIND_STEM_NITROGEN')
+obs_kind_names(127) = obs_kind_type(KIND_LEAF_NITROGEN ,'KIND_LEAF_NITROGEN')
! count here, then output below
@@ -537,8 +556,8 @@
endif
end do
! Falling off the end is an error
- write(err_string, *) trim(assimilate_these_obs_types(i)), &
- ' from obs_kind_nml is not a legal observation kind'
+ write(err_string, *) '"',trim(assimilate_these_obs_types(i)), &
+ '" from obs_kind_nml is not a legal observation kind to assimilate'
call error_handler(E_ERR, 'initialize_module', err_string, source, revision, revdate)
44 continue
end do
@@ -555,8 +574,8 @@
endif
end do
! Falling off the end is an error
- write(err_string, *) trim(evaluate_these_obs_types(i)), &
- ' from obs_kind_nml is not a legal observation kind'
+ write(err_string, *) '"',trim(evaluate_these_obs_types(i)), &
+ '" from obs_kind_nml is not a legal observation kind to evaluate'
call error_handler(E_ERR, 'initialize_module', err_string, source, revision, revdate)
55 continue
end do
Modified: DART/branches/development/observations/Ameriflux/level4_to_obs.f90
===================================================================
--- DART/branches/development/observations/Ameriflux/level4_to_obs.f90 2012-06-12 17:04:54 UTC (rev 5759)
+++ DART/branches/development/observations/Ameriflux/level4_to_obs.f90 2012-06-12 17:24:22 UTC (rev 5760)
@@ -25,11 +25,13 @@
use utilities_mod, only : initialize_utilities, finalize_utilities, &
register_module, error_handler, E_MSG, E_ERR, &
open_file, close_file, do_nml_file, do_nml_term, &
- check_namelist_read, find_namelist_in_file
+ check_namelist_read, find_namelist_in_file, &
+ nmlfileunit
use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, &
set_date, set_time, get_time, print_time, &
- print_date, operator(-), operator(+)
+ print_date, operator(-), operator(+), operator(>), &
+ operator(<), operator(==), operator(<=), operator(>=)
use location_mod, only : VERTISHEIGHT
@@ -77,7 +79,7 @@
!-----------------------------------------------------------------------
character(len=256) :: input_line, string1, string2, string3
-integer :: nmlfileunit, iline, nlines
+integer :: iline, nlines
logical :: file_exist, first_obs
integer :: n, i, oday, osec, rcio, iunit
integer :: num_copies, num_qc, max_obs
@@ -85,7 +87,7 @@
type(obs_sequence_type) :: obs_seq
type(obs_type) :: obs, prev_obs
type(time_type) :: time_obs, prev_time, offset
-integer, parameter :: umol_to_gC = 12.0_r8 * -1000000.0_r8
+real(r8), parameter :: umol_to_gC = 1.0_r8/(1000000.0_r8 * 12.0_r8)
type towerdata
type(time_type) :: time_obs
@@ -143,7 +145,8 @@
! time setup
call set_calendar_type(GREGORIAN)
-offset = set_time(nint(abs(timezoneoffset)*3600.0_r8),0)
+offset = set_time(nint(abs(timezoneoffset)*3600.0_r8),0)
+prev_time = set_time(0, 0)
if (verbose) print *, 'tower located at lat, lon, elev =', latitude, longitude, elevation
if (verbose) print *, 'flux observations taken at =', flux_height,'m'
@@ -204,11 +207,11 @@
if (rcio > 0) then
write (string1,'(''Cannot read (error '',i3,'') line '',i8,'' in '',A)') &
rcio, iline, trim(text_input_file)
- call error_handler(E_ERR,'count_file_lines', string1, source, revision, revdate)
+ call error_handler(E_ERR,'main', string1, source, revision, revdate)
endif
! parse the line into the tower structure (including the observation time)
- call stringparse(input_line,iline)
+ call stringparse(input_line, iline)
if (iline <= 2) then
write(*,*)''
@@ -234,31 +237,34 @@
! make an obs derived type, and then add it to the sequence
! If the QC value is good, use the observation.
! Increasingly larger QC values are more questionable quality data.
+ ! The observation errors are from page 183, Table 7.1(A) in
+ ! Chapter 7 of a book by A.D. Richardson et al. via Andy Fox.
- if (tower%hQC <= maxgoodqc) then
- oerr = 2.5_r8 ! 10 percent of the 'nighttime' values - total guess
+ if (tower%hQC <= maxgoodqc) then ! Sensible Heat Flux [W m-2]
+ oerr = 10.0_r8 + abs(tower%h)*0.22_r8
qc = real(tower%hQC,r8)
call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%h, &
TOWER_LATENT_HEAT_FLUX, oerr, oday, osec, qc, obs)
call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs)
endif
- if (tower%leQC <= maxgoodqc) then
- oerr = 0.3_r8 ! 10 percent of the 'nighttime' values - total guess
+ if (tower%leQC <= maxgoodqc) then ! Latent Heat Flux [W m-2]
+ oerr = 10.0_r8 + abs(tower%le)*0.32_r8
qc = real(tower%leQC,r8)
call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%le, &
TOWER_SENSIBLE_HEAT_FLUX, oerr, oday, osec, qc, obs)
call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs)
endif
- ! A crude estimate of this would be something like 1+(0.25*|flux|) mumol m^2s^1
- ! which will give a range of ~1-5 at most of our sites. -- Andy Fox
- if (tower%neeQC <= maxgoodqc) then
- oerr = 1.0_r8 + abs(tower%NEE)*0.25_r8
- oerr = oerr * umol_to_gC
- tower%NEE = tower%NEE * umol_to_gC
+ if (tower%neeQC <= maxgoodqc) then ! Net Ecosystem Exchange [umol m-2 s-1]
+ if (tower%nee <= 0) then
+ oerr = (2.0_r8 + abs(tower%nee)*0.1_r8) * umol_to_gC
+ else
+ oerr = (2.0_r8 + abs(tower%nee)*0.4_r8) * umol_to_gC
+ endif
+ tower%nee = -tower%nee * umol_to_gC ! to match convention in CLM [gC m-2 s-1]
qc = real(tower%neeQC,r8)
- call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%Nee, &
+ call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%nee, &
TOWER_NETC_ECO_EXCHANGE, oerr, oday, osec, qc, obs)
call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs)
endif
@@ -293,7 +299,7 @@
! vkind - kind of vertical coordinate (pressure, level, etc)
! obsv - observation value
! okind - observation kind
-! oerr - observation error
+! oerr - observation error (in units of standard deviation)
! day - gregorian day
! sec - gregorian second
! qc - quality control value
@@ -462,13 +468,14 @@
subroutine stringparse(str1,linenum)
-! just declare everything as reals and see how it falls out
+! just declare everything as reals and chunk it
character(len=*), intent(in) :: str1
integer , intent(in) :: linenum
real(r8), dimension(34) :: values
-integer :: ihour, imin, isec, seconds
+integer :: iday, ihour, imin, isec, seconds
+type(time_type) :: time0, time1, time2
values = MISSING_R8
@@ -480,7 +487,9 @@
! Stuff what we want into the tower structure
!
-! Convert to 'CLM-friendly' units.
+! Convert to 'CLM-friendly' units AFTER we determine observation error variance.
+! That happens in the main routine.
+!
! NEE_or_fMDS has units [umolCO2 m-2 s-1]
! H_f has units [W m-2]
! LE_f has units [W m-2]
@@ -498,22 +507,65 @@
tower%h = values(tower%hindex )
tower%hQC = nint(values(tower%hQCindex ))
-! put observation time/date into a dart time format
+! decode the time pieces ... two times ...
+! The LAST line of these files is knackered ... and we have to check that
+! if the doy is greater than the ymd ...
+! 12,31,23.5,366.979 N-1
+! 1, 1, 0.0,367.000 N
+iday = int(tower%doy)
ihour = int(tower%hour)
seconds = nint((tower%hour - real(ihour,r8))*3600)
imin = seconds / 60
isec = seconds - imin * 60
+time0 = set_date(year, tower%month, tower%day, ihour, imin, isec)
-tower%time_obs = set_date(year, tower%month, tower%day, ihour, imin, isec)
+isec = ihour*3600 + imin*60 + isec
+time1 = set_date(year,1,1,0,0,0) + set_time(isec, (iday-1))
+time2 = time0 - time1
+call get_time(time2, isec, iday)
+
+if ( iday > 0 ) then
+ ! we need to change the day ...
+
+ tower%time_obs = time1
+
+ if (verbose) then
+ write(string1,*)'converting ',tower%month,tower%day,tower%hour,tower%doy
+ write(string2,*)'the day-of-year indicates we should amend the month/day values.'
+ call error_handler(E_MSG,'stringparse', string1, source, revision, &
+ revdate, text2=string2)
+
+ call print_date(time0, 'stringparse: using ymd date is')
+ call print_date(time1, 'stringparse: using doy date is')
+ call print_time(time0, 'stringparse: using ymd time is')
+ call print_time(time1, 'stringparse: using doy time is')
+ call print_time(time2, 'stringparse: difference is')
+ endif
+else
+
+ tower%time_obs = time0
+
+endif
+
if (timezoneoffset < 0.0_r8) then
tower%time_obs = tower%time_obs - offset
else
tower%time_obs = tower%time_obs + offset
endif
-
+! The QC values can be 'missing' ... in which case the values are too
+
+if (tower%neeQC < 0) tower%neeQC = maxgoodqc + 1000
+if (tower%leQC < 0) tower%leQC = maxgoodqc + 1000
+if (tower%hQC < 0) tower%hQC = maxgoodqc + 1000
+
+! if (tower%neeQC < maxgoodqc) then
+! write(*,*)'nee umol m-2 s-1 ',tower%nee
+! write(*,*)'nee gC m-2 s-1 ',tower%nee*umol_to_gC
+! endif
+
end subroutine stringparse
Modified: DART/branches/development/observations/Ameriflux/work/path_names_level4_to_obs
===================================================================
--- DART/branches/development/observations/Ameriflux/work/path_names_level4_to_obs 2012-06-12 17:04:54 UTC (rev 5759)
+++ DART/branches/development/observations/Ameriflux/work/path_names_level4_to_obs 2012-06-12 17:24:22 UTC (rev 5760)
@@ -4,7 +4,7 @@
obs_kind/obs_kind_mod.f90
obs_def/obs_def_mod.f90
assim_model/assim_model_mod.f90
-models/template/model_mod.f90
+models/clm/model_mod.f90
common/types_mod.f90
random_seq/random_seq_mod.f90
utilities/utilities_mod.f90
Modified: DART/branches/development/observations/Ameriflux/work/path_names_obs_sequence_tool
===================================================================
--- DART/branches/development/observations/Ameriflux/work/path_names_obs_sequence_tool 2012-06-12 17:04:54 UTC (rev 5759)
+++ DART/branches/development/observations/Ameriflux/work/path_names_obs_sequence_tool 2012-06-12 17:24:22 UTC (rev 5760)
@@ -4,7 +4,7 @@
obs_def/obs_def_mod.f90
cov_cutoff/cov_cutoff_mod.f90
assim_model/assim_model_mod.f90
-models/template/model_mod.f90
+models/clm/model_mod.f90
common/types_mod.f90
location/threed_sphere/location_mod.f90
mpi_utilities/null_mpi_utilities_mod.f90
More information about the Dart-dev
mailing list