[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