[Dart-dev] [5903] DART/branches/development: The COSMOS observation converters pass my preliminary tests.

nancy at ucar.edu nancy at ucar.edu
Thu Oct 18 16:29:49 MDT 2012


Revision: 5903
Author:   thoar
Date:     2012-10-18 16:29:48 -0600 (Thu, 18 Oct 2012)
Log Message:
-----------
The COSMOS observation converters pass my preliminary tests.

Modified Paths:
--------------
    DART/branches/development/obs_def/obs_def_COSMOS_mod.f90

Added Paths:
-----------
    DART/branches/development/observations/COSMOS/
    DART/branches/development/observations/COSMOS/COSMOS_level2_to_obs.f90
    DART/branches/development/observations/COSMOS/COSMOS_to_obs.f90
    DART/branches/development/observations/COSMOS/COSMOS_to_obs.nml
    DART/branches/development/observations/COSMOS/data/
    DART/branches/development/observations/COSMOS/data/COSMIC_parlist.netCDF.txt
    DART/branches/development/observations/COSMOS/data/COSMIC_parlist_station.txt
    DART/branches/development/observations/COSMOS/data/COSMOS_SantaRita_2010.txt
    DART/branches/development/observations/COSMOS/work/
    DART/branches/development/observations/COSMOS/work/input.nml
    DART/branches/development/observations/COSMOS/work/mkmf_COSMOS_level2_to_obs
    DART/branches/development/observations/COSMOS/work/mkmf_COSMOS_to_obs
    DART/branches/development/observations/COSMOS/work/mkmf_obs_sequence_tool
    DART/branches/development/observations/COSMOS/work/mkmf_preprocess
    DART/branches/development/observations/COSMOS/work/path_names_COSMOS_level2_to_obs
    DART/branches/development/observations/COSMOS/work/path_names_COSMOS_to_obs
    DART/branches/development/observations/COSMOS/work/path_names_obs_sequence_tool
    DART/branches/development/observations/COSMOS/work/path_names_preprocess
    DART/branches/development/observations/COSMOS/work/quickbuild.csh

-------------- next part --------------
Modified: DART/branches/development/obs_def/obs_def_COSMOS_mod.f90
===================================================================
--- DART/branches/development/obs_def/obs_def_COSMOS_mod.f90	2012-10-18 22:24:56 UTC (rev 5902)
+++ DART/branches/development/obs_def/obs_def_COSMOS_mod.f90	2012-10-18 22:29:48 UTC (rev 5903)
@@ -2,7 +2,29 @@
 ! provided by UCAR, "as is", without charge, subject to all terms of use at
 ! http://www.image.ucar.edu/DAReS/DART/DART_download
 
-! Created by Rafael Rosolem (09/30/2011) for COSMOS based on file by Ave Arellano
+!----------------------------------------------------------------------
+! This module provides support for observations from COSMOS.
+! Each COSMOS installation has soil properties that have been measured
+! and are needed by the COSMIC algorithm that converts moisture to
+! counts of neutron intensity. These properties constitute additional
+! metadata for each observation. The routines in this module read and
+! write that metadata and provide the 'get_expected_neutron_intensity()'
+! routine that is the 'observation operator'
+!
+!  OBS            1
+!           -1           2          -1
+! obdef
+! loc3d
+!      4.188790204786391        0.6986026390077337         137.7093918137252      3
+! kind
+!            1
+! cosmic
+!    <an array of 4 parameters>
+!    <an array of 4 parameters>
+!    <neutron_intensity_key>
+!      0          0
+!    4.0000000000000000
+!----------------------------------------------------------------------
 
 ! BEGIN DART PREPROCESS KIND LIST
 ! COSMOS_NEUTRON_INTENSITY,    KIND_NEUTRON_INTENSITY
@@ -10,14 +32,12 @@
 
 
 ! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
-!   use obs_def_COSMOS_mod, only : read_neutron_intensity, &
-!                                 write_neutron_intensity, &
-!                          get_expected_neutron_intensity, &
-!                           set_obs_def_neutron_intensity, &
-!                           interactive_neutron_intensity
+!   use obs_def_COSMOS_mod, only : read_cosmos_metadata, &
+!                                 write_cosmos_metadata, &
+!                           interactive_cosmos_metadata, &
+!                          get_expected_neutron_intensity
 ! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
 
-! TJH Questions for Nancy: AFAICT none of this should be public ... 
 
 ! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
 !      case(COSMOS_NEUTRON_INTENSITY)
@@ -26,29 +46,23 @@
 
 
 ! BEGIN DART PREPROCESS READ_OBS_DEF
-!      case(COSMOS_NEUTRON_INTENSITY)
-!         call read_neutron_intensity(obs_def%key, ifile, fform)
+!   case(COSMOS_NEUTRON_INTENSITY)
+!      call read_cosmos_metadata(obs_def%key, key, ifile, fform)
 ! END DART PREPROCESS READ_OBS_DEF
 
 
 ! BEGIN DART PREPROCESS WRITE_OBS_DEF
-!      case(COSMOS_NEUTRON_INTENSITY)
-!         call write_neutron_intensity(obs_def%key, ifile, fform)
+!   case(COSMOS_NEUTRON_INTENSITY)
+!      call write_cosmos_metadata(obs_def%key, ifile, fform)
 ! END DART PREPROCESS WRITE_OBS_DEF
 
 
 ! BEGIN DART PREPROCESS INTERACTIVE_OBS_DEF
-!      case(COSMOS_NEUTRON_INTENSITY)
-!         call interactive_neutron_intensity(obs_def%key)
+!   case(COSMOS_NEUTRON_INTENSITY)
+!      call interactive_cosmos_metadata(obs_def%key)
 ! END DART PREPROCESS INTERACTIVE_OBS_DEF
 
 
-! BEGIN DART PREPROCESS SET_OBS_DEF_NEUTRON_INTENSITY
-!      case(COSMOS_NEUTRON_INTENSITY)
-!         call set_obs_def_neutron_intensity(obs_def%key)
-! END DART PREPROCESS SET_OBS_DEF_NEUTRON_INTENSITY
-
-
 ! BEGIN DART PREPROCESS MODULE CODE
 module obs_def_COSMOS_mod
 
@@ -58,112 +72,339 @@
 ! $Revision$
 ! $Date$
 
-use        types_mod, only : r8, PI
+use        types_mod, only : r8, PI, metadatalength, MISSING_R8
 use    utilities_mod, only : register_module, error_handler, E_ERR, E_WARN, E_MSG, &
-                             logfileunit, get_unit, open_file, close_file
+                             logfileunit, get_unit, open_file, close_file, nc_check, &
+                             file_exist, ascii_file_format
 use     location_mod, only : location_type, set_location, get_location, &
                              vert_is_height,   VERTISHEIGHT,            &
-                             vert_is_level,    VERTISLEVEL
+                             vert_is_level,    VERTISLEVEL,             &
+                             set_location_missing
 use        model_mod, only : model_interpolate
 use     obs_kind_mod, only : KIND_GEOPOTENTIAL_HEIGHT, KIND_SOIL_MOISTURE
 
+use typesizes
+use netcdf
+
 implicit none
 private
 
-public :: set_obs_def_neutron_intensity,  &
-          get_expected_neutron_intensity, &
-          interactive_neutron_intensity,  &
-          read_neutron_intensity,         &
-          write_neutron_intensity
+public ::            set_cosmos_metadata, &
+                     get_cosmos_metadata, &
+                    read_cosmos_metadata, &
+                   write_cosmos_metadata, &
+             interactive_cosmos_metadata, &
+          get_expected_neutron_intensity
 
-integer, parameter :: nlayers        = 4 ! Number of soil layers used in Noah model
-integer            :: num_neutron_intensity = 0 ! observation counter
+integer :: num_neutron_intensity = 0 ! observation counter ... useful length of metadata arrays
 
 ! version controlled file description for error handling, do not edit
 character(len=128), parameter :: &
    source   = "$URL$", &
-   revision = "$Revision $", &
+   revision = "$Revision$", &
    revdate  = "$Date$"
 
-character(len=256) :: string1
+character(len=256) :: string1, string2
 logical, save      :: module_initialized = .false.
 
+! Metadata for COSMOS observations.
+! There are soil parameters for each site that must be added to each
+! observation in the sequence. Also COSMIC parameters ...
+
+type site_metadata
+   private
+!  character(len=metadatalength) :: sitename
+!  type(location_type) :: location
+!  real(r8)            :: latitude
+!  real(r8)            :: longitude
+!  real(r8)            :: elevation
+   real(r8)            :: bd         ! Dry Soil Bulk Density [  g / cm^3]
+   real(r8)            :: lattwat    ! Lattice Water Content [M^3 /  M^3]
+   real(r8)            :: N          ! High Energy Neutron Intensity
+   real(r8)            :: alpha      ! Ratio of Fast Neutron Creation Factor
+   real(r8)            :: L1         ! High Energy   Soil Attenuation Length
+   real(r8)            :: L2         ! High Energy  Water Attenuation Length
+   real(r8)            :: L3         ! Fast Neutron  Soil Attenuation Length
+   real(r8)            :: L4         ! Fast Neutron Water Attenuation Length
+end type site_metadata
+
+type(site_metadata), allocatable, dimension(:) :: observation_metadata
+character(len=6), parameter :: COSMOSSTRING = 'cosmic'
+
+logical :: debug = .FALSE.
+integer :: MAXcosmoskey = 24*366  ! one year of hourly data - to start
+integer ::    cosmoskey = 0
+
 !----------------------------------------------------------------------------
 contains
 !----------------------------------------------------------------------------
 
 
+
   subroutine initialize_module
 !----------------------------------------------------------------------------
 ! subroutine initialize_module
 !
-! a DART tradition
+
+integer :: i
+
 call register_module(source, revision, revdate)
 
 module_initialized = .true.
 
+allocate(observation_metadata(MAXcosmoskey))
+
+do i = 1,MAXcosmoskey
+   observation_metadata(i)%bd       = MISSING_R8
+   observation_metadata(i)%lattwat  = MISSING_R8
+   observation_metadata(i)%N        = MISSING_R8
+   observation_metadata(i)%alpha    = MISSING_R8
+   observation_metadata(i)%L1       = MISSING_R8
+   observation_metadata(i)%L2       = MISSING_R8
+   observation_metadata(i)%L3       = MISSING_R8
+   observation_metadata(i)%L4       = MISSING_R8
+!  observation_metadata(i)%location = set_location_missing()
+enddo
+
 end subroutine initialize_module
 
 
- subroutine read_neutron_intensity(key, ifile, fform)
+
+ subroutine set_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
 !----------------------------------------------------------------------
-!subroutine read_neutron_intensity(key, ifile, fform)
+!subroutine set_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
 !
-integer,          intent(out)          :: key
+! Fill the module storage metadata for a particular observation.
+
+integer,  intent(out) :: key
+real(r8), intent(in)  :: bd, lattwat, N, alpha, L1, L2, L3, L4
+
+if ( .not. module_initialized ) call initialize_module
+
+cosmoskey = cosmoskey + 1  ! increase module storage used counter
+
+! Make sure the new key is within the length of the metadata arrays.
+call key_out_of_range(cosmoskey,'set_cosmos_metadata')
+
+key = cosmoskey ! now that we know its legal
+
+observation_metadata(key)%bd      = bd
+observation_metadata(key)%lattwat = lattwat
+observation_metadata(key)%N       = N
+observation_metadata(key)%alpha   = alpha
+observation_metadata(key)%L1      = L1
+observation_metadata(key)%L2      = L2
+observation_metadata(key)%L3      = L3
+observation_metadata(key)%L4      = L4
+
+end subroutine set_cosmos_metadata
+
+
+
+ subroutine get_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
+!----------------------------------------------------------------------
+!subroutine get_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
+!
+! Query the metadata in module storage for a particular observation.
+! This can be useful for post-processing routines, etc.
+
+integer,  intent(in)  :: key
+real(r8), intent(out) :: bd, lattwat, N, alpha, L1, L2, L3, L4
+
+if ( .not. module_initialized ) call initialize_module
+
+! Make sure the desired key is within the length of the metadata arrays.
+call key_out_of_range(key,'get_cosmos_metadata')
+
+bd       = observation_metadata(key)%bd
+lattwat  = observation_metadata(key)%lattwat
+N        = observation_metadata(key)%N
+alpha    = observation_metadata(key)%alpha
+L1       = observation_metadata(key)%L1
+L2       = observation_metadata(key)%L2
+L3       = observation_metadata(key)%L3
+L4       = observation_metadata(key)%L4
+
+end subroutine get_cosmos_metadata
+
+
+
+ subroutine read_cosmos_metadata(key,       obsID, ifile, fform)
+!----------------------------------------------------------------------
+!subroutine read_cosmos_metadata(obs_def%key, key, ifile, fform)
+!
+! This routine reads the metadata for neutron intensity observations.
+!
+integer,          intent(out)          :: key    ! index into local metadata
+integer,          intent(in)           :: obsID
 integer,          intent(in)           :: ifile
 character(len=*), intent(in), optional :: fform
 
 ! temp variables
-integer           :: keyin
-integer           :: counts1
-character(len=32) :: fileformat
+logical           :: is_asciifile
+integer           :: ierr
+character(len=6)  :: header
+integer           :: oldkey
+real(r8)          :: bd, lattwat, N, alpha, L1, L2, L3, L4
 
 if ( .not. module_initialized ) call initialize_module
 
-fileformat = "ascii"   ! supply default
-if(present(fform)) fileformat = trim(adjustl(fform))
+is_asciifile = ascii_file_format(fform)
 
-SELECT CASE (fileformat)
-   CASE ("unf", "UNF", "unformatted", "UNFORMATTED")
-!      read(ifile) keyin
-   CASE DEFAULT
-!      read(ifile, *) keyin
-END SELECT
+write(string2,*)'observation #',obsID
 
-counts1 = counts1 + 1
-key     = counts1
-call set_obs_def_neutron_intensity(key)
+if ( is_asciifile ) then
+   read(ifile, *, iostat=ierr) header
+   call check_iostat(ierr,'read_cosmos_metadata','header',string2)
+   if (trim(header) /= trim(COSMOSSTRING)) then
+       write(string1,*)"Expected neutron_intensity header ["//COSMOSSTRING//"] in input file, got ["//header//"]"
+       call error_handler(E_ERR, 'read_cosmos_metadata', string1, source, revision, revdate, text2=string2)
+   endif
+   read(ifile, *, iostat=ierr) bd, lattwat, N, alpha
+   call check_iostat(ierr,'read_cosmos_metadata','bd -> alpha',string2)
+   read(ifile, *, iostat=ierr) L1, L2, L3, L4
+   call check_iostat(ierr,'read_cosmos_metadata','L1 -> L4',string2)
+   read(ifile, *, iostat=ierr) oldkey
+   call check_iostat(ierr,'read_cosmos_metadata','oldkey',string2)
+else
+   read(ifile, iostat=ierr) header
+   call  check_iostat(ierr,'read_cosmos_metadata','header',string2)
+   if (trim(header) /= trim(COSMOSSTRING)) then
+       write(string1,*)"Expected neutron_intensity header ["//COSMOSSTRING//"] in input file, got ["//header//"]"
+       call error_handler(E_ERR, 'read_cosmos_metadata', string1, source, revision, revdate, text2=string2)
+   endif
+   read(ifile, iostat=ierr) bd, lattwat, N, alpha
+   call  check_iostat(ierr,'read_cosmos_metadata','bd -> alpha',string2)
+   read(ifile, iostat=ierr) L1, L2, L3, L4
+   call  check_iostat(ierr,'read_cosmos_metadata','L1 -> L4',string2)
+   read(ifile, iostat=ierr) oldkey
+   call  check_iostat(ierr,'read_cosmos_metadata','oldkey',string2)
+endif
 
-end subroutine read_neutron_intensity
+! The oldkey is thrown away.
 
+! Store the metadata in module storage and record the new length of the metadata arrays.
+call set_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
 
+! The new 'key' is returned.
 
- subroutine write_neutron_intensity(key, ifile, fform)
+end subroutine read_cosmos_metadata
+
+
+
+ subroutine write_cosmos_metadata(key, ifile, fform)
 !----------------------------------------------------------------------
-!subroutine write_neutron_intensity(key, ifile, fform)
+!subroutine write_cosmos_metadata(key, ifile, fform)
+!
+! writes the metadata for neutron intensity observations.
 
-integer,          intent(in)           :: key
-integer,          intent(in)           :: ifile
-character(len=*), intent(in), optional :: fform
+integer,           intent(in)           :: key
+integer,           intent(in)           :: ifile
+character(len=*),  intent(in), optional :: fform
 
-character(len=32)                      :: fileformat
+logical  :: is_asciifile
+real(r8) :: bd, lattwat, N, alpha, L1, L2, L3, L4
 
 if ( .not. module_initialized ) call initialize_module
 
-fileformat = "ascii"   ! supply default
-if(present(fform)) fileformat = trim(adjustl(fform))
+! given the index into the local metadata arrays - retrieve
+! the metadata for this particular observation.
 
-SELECT CASE (fileformat)
-   CASE ("unf", "UNF", "unformatted", "UNFORMATTED")
-!         write(ifile) key
-   CASE DEFAULT
-!         write(ifile, *) key
-END SELECT
+call get_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
 
-end subroutine write_neutron_intensity
+is_asciifile = ascii_file_format(fform)
 
+if (is_asciifile) then
+   write(ifile, *) trim(COSMOSSTRING)
+   write(ifile, *) bd, lattwat, N, alpha
+   write(ifile, *) L1, L2, L3, L4
+   write(ifile, *) key
+else
+   write(ifile   ) trim(COSMOSSTRING)
+   write(ifile   ) bd, lattwat, N, alpha
+   write(ifile   ) L1, L2, L3, L4
+   write(ifile   ) key
+endif
 
+end subroutine write_cosmos_metadata
 
+
+
+ subroutine interactive_cosmos_metadata(key)
+!----------------------------------------------------------------------
+!subroutine interactive_cosmos_metadata(key)
+!
+integer, intent(out) :: key
+
+real(r8) :: bd, lattwat, N, alpha, L1, L2, L3, L4
+
+if ( .not. module_initialized ) call initialize_module
+
+! Prompt for input for the required metadata
+
+bd      = interactive('"bd"      dry soil bulk density [g/cm^3]'                ,minvalue = 0.0_r8)
+lattwat = interactive('"lattwat" lattice water content [m^3/m^3]'               ,minvalue = 0.0_r8)
+N       = interactive('"N"       high energy neutron intensity [count]'         ,minvalue = 0.0_r8)
+alpha   = interactive('"alpha"   ratio of fast neutron creation factor [Soil to Water]',minvalue=0.0_r8)
+L1      = interactive('"L1"      high energy soil attenuation length [g/cm^2]'  ,minvalue = 0.0_r8)
+L2      = interactive('"L2"      high energy water attenuation length [g/cm^2]' ,minvalue = 0.0_r8)
+L3      = interactive('"L3"      fast neutron soil attenuation length [g/cm^2]' ,minvalue = 0.0_r8)
+L4      = interactive('"L4"      fast neutron water attenuation length [g/cm^2]',minvalue = 0.0_r8)
+
+call set_cosmos_metadata(key, bd, lattwat, N, alpha, L1, L2, L3, L4)
+
+end subroutine interactive_cosmos_metadata
+
+
+
+function interactive(str1,minvalue,maxvalue)
+real(r8)                       :: interactive
+character(len=*),   intent(in) :: str1
+real(r8), optional, intent(in) :: minvalue
+real(r8), optional, intent(in) :: maxvalue
+
+integer :: i
+
+interactive = MISSING_R8
+
+! Prompt with a minimum amount of error checking
+
+if     (present(minvalue) .and. present(maxvalue)) then
+
+   interactive = minvalue - 1.0_r8
+   MINMAXLOOP : do i = 1,10
+      if ((interactive >= minvalue) .and. (interactive <= maxvalue)) exit MINMAXLOOP
+      write(*, *) 'Enter '//str1
+      read( *, *) interactive
+   end do MINMAXLOOP
+
+elseif (present(minvalue)) then
+
+   interactive = minvalue - 1.0_r8
+   MINLOOP : do i=1,10
+      if (interactive >= minvalue) exit MINLOOP
+      write(*, *) 'Enter '//str1
+      read( *, *) interactive
+   end do MINLOOP
+
+elseif (present(maxvalue)) then
+
+   interactive = maxvalue + 1.0_r8
+   MAXLOOP : do i=1,10
+      if (interactive <= maxvalue) exit MAXLOOP
+      write(*, *) 'Enter '//str1
+      read( *, *) interactive
+   end do MAXLOOP
+
+else ! anything goes ... cannot check
+      write(*, *) 'Enter '//str1
+      read( *, *) interactive
+endif
+
+end function interactive
+
+
+
  subroutine get_expected_neutron_intensity(state, location, key, val, istatus)
 !----------------------------------------------------------------------
 !subroutine get_expected_neutron_intensity(state, location, key, val, istatus)
@@ -172,7 +413,7 @@
 
 real(r8),            intent(in)  :: state(:) ! state vector
 type(location_type), intent(in)  :: location ! location of obs
-integer,             intent(in)  :: key      ! obs key
+integer,             intent(in)  :: key      ! key into module metadata
 real(r8),            intent(out) :: val      ! value of obs
 integer,             intent(out) :: istatus  ! status of the calculation
 
@@ -209,9 +450,9 @@
 real(r8) :: vwclat = 0.0_r8 ! Volumetric "lattice" water content (m3/m3)
 real(r8) :: N      = 0.0_r8 ! High energy neutron flux (-)
 real(r8) :: alpha  = 0.0_r8 ! Ratio of Fast Neutron Creation Factor (Soil to Water), alpha (-)
-real(r8) :: L1     = 0.0_r8 ! High Energy Soil Attenuation Length (g/cm2)
-real(r8) :: L2     = 0.0_r8 ! High Energy Water Attenuation Length (g/cm2)
-real(r8) :: L3     = 0.0_r8 ! Fast Neutron Soil Attenuation Length (g/cm2)
+real(r8) :: L1     = 0.0_r8 ! High Energy Soil   Attenuation Length (g/cm2)
+real(r8) :: L2     = 0.0_r8 ! High Energy Water  Attenuation Length (g/cm2)
+real(r8) :: L3     = 0.0_r8 ! Fast Neutron Soil  Attenuation Length (g/cm2)
 real(r8) :: L4     = 0.0_r8 ! Fast Neutron Water Attenuation Length (g/cm2)
 real(r8) :: zdeg
 real(r8) :: zrad
@@ -253,36 +494,60 @@
 val = 0.0_r8 ! set return value early
 
 !=================================================================================
+! COSMIC: Site specific-parameters come from the observation metadata
+!=================================================================================
+
+! Make sure the desired key is within the length of the metadata arrays.
+call key_out_of_range(key,'get_expected_neutron_intensity')
+
+bd     = observation_metadata(key)%bd
+vwclat = observation_metadata(key)%lattwat
+N      = observation_metadata(key)%N
+alpha  = observation_metadata(key)%alpha
+L1     = observation_metadata(key)%L1
+L2     = observation_metadata(key)%L2
+L3     = observation_metadata(key)%L3
+L4     = observation_metadata(key)%L4
+
+write(*,*)'TJH debug '
+write(*,*)bd, vwclat, N, alpha
+write(*,*)L1,L2,L3,L4
+write(*,*)
+
+!=================================================================================
 ! Determine the number of soil layers and their depths
 ! using only the standard DART interfaces to model
 ! set the locations for each of the model levels
+!=================================================================================
 
 loc_array = get_location(location) ! loc is in DEGREES
 loc_lon   = loc_array(1)
 loc_lat   = loc_array(2)
 
 nlevels = 0
-COUNTLEVELS : do i = 1,maxlayers 
+COUNTLEVELS : do i = 1,maxlayers
    loc = set_location(loc_lon, loc_lat, real(i,r8), VERTISLEVEL)
    call model_interpolate(state,loc,KIND_GEOPOTENTIAL_HEIGHT,obs_val,istatus)
    if (istatus /= 0) exit COUNTLEVELS
    nlevels = nlevels + 1
 enddo COUNTLEVELS
 
-if (nlevels == maxlayers) then
+if ((nlevels == maxlayers) .or. (nlevels == 0)) then
    write(string1,*) 'FAILED to determine number of soil layers in model.'
-   call error_handler(E_ERR,'get_expected_neutron_intensity',string1,source,revision,revdate)
+   if (debug) call error_handler(E_MSG,'get_expected_neutron_intensity',string1,source,revision,revdate)
+   istatus = 1
+   val     = MISSING_R8
+   return 
 else
-   write(*,*)'get_expected_neutron_intensity: we have ',nlevels,' soil layers'
+!   write(*,*)'get_expected_neutron_intensity: we have ',nlevels,' soil layers'
 endif
 
-!=================================================================================
 ! Now actually find the depths at each level
 ! While we're at it, might as well get the soil moisture there, too.
 
 allocate(layerz(nlevels),soil_moisture(nlevels))
 
-FINDLEVELS : do i = 1,nlevels 
+FINDLEVELS : do i = 1,nlevels
    loc = set_location(loc_lon, loc_lat, real(i,r8), VERTISLEVEL)
    call model_interpolate(state,loc,KIND_GEOPOTENTIAL_HEIGHT,layerz(i),istatus)
    ! not checking this error code because it worked just a few lines earlier
@@ -292,11 +557,14 @@
 
    if (istatus /= 0) then
       write(string1,*) 'FAILED to determine soil moisture for layer',i
-      call error_handler(E_ERR,'get_expected_neutron_intensity',string1,source,revision,revdate)
+      if (debug) call error_handler(E_MSG,'get_expected_neutron_intensity',string1,source,revision,revdate)
+      istatus = 2
+      val     = MISSING_R8
+      return
    endif
 
-   soil_moisture(i) = obs_val 
-   
+   soil_moisture(i) = obs_val
+
 enddo FINDLEVELS
 
 !rr: DART soil layers are negative and given in meters while COSMIC needs them to be
@@ -309,24 +577,6 @@
 endif
 
 !=================================================================================
-! COSMIC: Site specific-parameters
-!=================================================================================
-
-!rr: #####
-!rr: I will need to change find a way to read those parameters externally
-!rr: #####
-
-! SANTA RITA SITE
-bd     = 1.4620_r8
-vwclat = 0.0366_r8
-N      = 399.05099359_r8
-alpha  = 0.25985853017_r8
-L1     = 161.98621864_r8
-L2     = 129.14558985_r8
-L3     = 114.04156391_r8
-L4     = 3.8086191933_r8
-
-!=================================================================================
 ! COSMIC: Allocate arrays and initialize variables
 !=================================================================================
 
@@ -379,14 +629,14 @@
 ! COSMIC: Neutron flux calculation
 !=================================================================================
 
-! At some point, you might want to tinker around with non-uniform 
+! At some point, you might want to tinker around with non-uniform
 ! soil layer thicknesses.
 zthick(1) = dz(1) - 0.0_r8 ! Surface layer
 do i = 2,nlyr
    zthick(i) = dz(i) - dz(i-1) ! Remaining layers
 enddo
 
-if ( 2 == 1 ) then ! TJH DEBUG OUTPUT BLOCK 
+if ( 2 == 1 ) then ! TJH DEBUG OUTPUT BLOCK
    iunit = open_file('cosmos_layers.txt',form='formatted',action='write')
    do i = 1,nlyr
       write(iunit,*)dz(i),vwc(i),zthick(i)
@@ -431,7 +681,7 @@
 
    ! This second loop needs to be done for the distribution of angles for fast neutron release
    ! the intent is to loop from 0 to 89.5 by 0.5 degrees - or similar.
-   ! Because Fortran loop indices are integers, we have to divide the indices by 10 - you get the idea.  
+   ! Because Fortran loop indices are integers, we have to divide the indices by 10 - you get the idea.
 
    do angle=0,maxangle,angledz
       zdeg     = real(angle,r8)/10.0_r8   ! 0.0  0.5  1.0  1.5 ...
@@ -452,52 +702,61 @@
 
 enddo
 
-!=================================================================================
-
-! ... and finally calculated what the neutron intensity (which is basically totflux)
-val = totflux
-
 deallocate(dz, zthick, vwc, hiflux, fastpot, &
            h2oeffdens, idegrad, fastflux, isoimass, iwatmass)
 
-! assume all is well for now
-istatus = 0
+!=================================================================================
+! ... and finally set the return the neutron intensity (i.e. totflux)
 
+val     = totflux
+istatus = 0        ! assume all is well if we get this far
+
 return
 end subroutine get_expected_neutron_intensity
 
 
 
- subroutine interactive_neutron_intensity(key)
+
+
+ subroutine check_iostat(istat, routine, varname, msgstring)
 !----------------------------------------------------------------------
-!subroutine interactive_neutron_intensity(key)
 !
-integer, intent(out) :: key
 
-if ( .not. module_initialized ) call initialize_module
+integer,          intent(in) :: istat
+character(len=*), intent(in) :: routine
+character(len=*), intent(in) :: varname
+character(len=*), intent(in) :: msgstring
 
-! Increment the index
-num_neutron_intensity = num_neutron_intensity + 1
-key = num_neutron_intensity
+if ( istat /= 0 ) then
+   write(string1,*)'istat should be 0 but is ',istat,' for '//varname
+   call error_handler(E_ERR, routine, string1, source, revision, revdate, text2=msgstring)
+end if
 
-! Otherwise, prompt for input for the three required beasts
-write(*, *) 'Creating an interactive_COSMOS observation'
+end subroutine check_iostat
 
-end subroutine interactive_neutron_intensity
 
 
-
- subroutine set_obs_def_neutron_intensity(key)
+subroutine key_out_of_range(key, routine)
 !----------------------------------------------------------------------
-! Allows passing of obs_def special information
+! Make sure we are addrssing within the metadata arrays
 
-integer, intent(in) :: key
+integer,          intent(in) :: key
+character(len=*), intent(in) :: routine
 
-if ( .not. module_initialized ) call initialize_module
+! fine -- no problem.
+if (key <= MAXcosmoskey) return
 
-end subroutine set_obs_def_neutron_intensity
+! Bad news.  Tell the user.
+write(string1, *) 'key (',key,') exceeds Nmax_neutron_intensity (', MAXcosmoskey,')'
+write(string2, *) 'Increase MAXcosmoskey in namelist and rerun.'
+call error_handler(E_ERR,routine,string1,source,revision,revdate,text2=string2)
 
+! FIXME ... reallocate and get on with it.
 
+end subroutine key_out_of_range
+
+
+
 end module obs_def_COSMOS_mod
 
 ! END DART PREPROCESS MODULE CODE

Added: DART/branches/development/observations/COSMOS/COSMOS_level2_to_obs.f90
===================================================================
--- DART/branches/development/observations/COSMOS/COSMOS_level2_to_obs.f90	                        (rev 0)
+++ DART/branches/development/observations/COSMOS/COSMOS_level2_to_obs.f90	2012-10-18 22:29:48 UTC (rev 5903)
@@ -0,0 +1,604 @@
+! DART software - Copyright 2004 - 2011 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+program COSMOS_level2_to_obs
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! COSMOS_level2_to_obs - reads the COSMOS data as defined in the level 2 
+!     product available from the COSMOS data portal at:
+!     http://cosmos.hwr.arizona.edu/Probes/probemap.php
+!
+!     original 3 May       2012   Tim Hoar       NCAR/IMAGe
+!     modified 6 September 2012   Rafael Rosolem University of Arizona 
+!     modified 10 October  2012   Tim Hoar
+!        * include site-specific metadata in obs sequence 
+!
+! QC flags are defined as: 0 = BEST     (corrected for water vapor)
+!                          1 = OK       (NOT corrected for water vapor)
+!
+! The correction for water vapor is intented to be available 
+! in the level 2 data at some point in the future.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use          types_mod, only : r8, MISSING_R8, metadatalength
+
+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, &
+                               nmlfileunit, file_exist, nc_check, to_upper, &
+                               find_textfile_dims
+
+use   time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, &
+                               set_date, set_time, get_time, print_time, &
+                               print_date, operator(-), operator(+), operator(>), &
+                               operator(<), operator(==), operator(<=), operator(>=)
+
+use       location_mod, only : location_type, set_location, VERTISSURFACE
+
+use   obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, &
+                               static_init_obs_sequence, init_obs, write_obs_seq, &
+                               init_obs_sequence, get_num_obs, set_obs_def, &
+                               set_copy_meta_data, set_qc_meta_data, &
+                               set_obs_values, set_qc
+
+use        obs_def_mod, only : obs_def_type, set_obs_def_kind, &
+                               set_obs_def_location, set_obs_def_time, &
+                               set_obs_def_error_variance, set_obs_def_key
+
+use  obs_utilities_mod, only : add_obs_to_seq
+
+use obs_def_COSMOS_mod, only : set_cosmos_metadata
+
+use       obs_kind_mod, only : COSMOS_NEUTRON_INTENSITY
+
+use typesizes
+use netcdf
+
+implicit none
+
+!-----------------------------------------------------------------------
+! version controlled file description for error handling, do not edit
+!-----------------------------------------------------------------------
+
+character(len=128), parameter :: &
+   source   = "$URL$", &
+   revision = "$Revision$", &
+   revdate  = "$Date$"
+
+!-----------------------------------------------------------------------
+! Namelist with default values
+!-----------------------------------------------------------------------
+
+character(len=256) :: site_metadata_file = 'COSMIC_parlist.nc'
+character(len=128) :: text_input_file = 'corcounts.txt'
+character(len=128) :: obs_out_file    = 'obs_seq.out'
+character(len=128) :: sitename        = 'missing'
+integer            :: year
+real(r8)           :: maxgoodqc       = 3
+logical            :: verbose         = .false.
+
+namelist /COSMOS_to_obs_nml/ site_metadata_file, text_input_file, &
+   obs_out_file, sitename, year, maxgoodqc, verbose
+
+!-----------------------------------------------------------------------
+! globally-scoped variables
+!-----------------------------------------------------------------------
+
+character(len=256)      :: input_line, string1, string2, string3
+integer                 :: iline
+logical                 :: first_obs
+integer                 :: oday, osec, rcio, iunit
+integer                 :: num_copies, num_qc, max_obs
+real(r8)                :: oerr, qc
+type(obs_sequence_type) :: obs_seq
+type(obs_type)          :: obs, prev_obs
+type(time_type)         :: prev_time
+integer, parameter      :: LEVEL2QC = 1
+
+! The cosmosdata type holds all the information from the text file
+! from the COSMOS instrument.
+!
+! "You should use this Level 2 data to get the counts and its standard deviation.
+!  They are the last two columns in the file (respectively, 'CORR' and 'ERR')"
+
+type cosmosdata
+  integer           :: dateindex
+  integer           :: timeindex
+  integer           :: neutronindex
+  integer           :: neutronstdindex
+  character(len=20) :: datestring       = 'YYYY-MM-DD'
+  character(len=20) :: timestring       = 'HH:MM'
+  character(len=20) :: neutronstring    = 'CORR'
+  character(len=20) :: neutronstdstring = 'ERR'
+  type(time_type)   :: time_obs
+  integer  :: year
+  integer  :: month
+  integer  :: day
+  integer  :: hour
+  integer  :: minute
+  real(r8) :: neutron
+  real(r8) :: neutronstd
+  real(r8) :: qc
+end type cosmosdata
+
+type(cosmosdata) :: cosmos ! we only need one of these.
+
+! The site_metadata type holds all the site-specific information
+
+type site_metadata
+   character(len=metadatalength) :: sitename
+   type(location_type) :: location
+   real(r8)            :: latitude
+   real(r8)            :: longitude
+   real(r8)            :: elevation
+   real(r8)            :: bd         ! Dry Soil Bulk Density [  g / cm^3]
+   real(r8)            :: lattwat    ! Lattice Water Content [M^3 /  M^3]
+   real(r8)            :: N          ! High Energy Neutron Intensity
+   real(r8)            :: alpha      ! Ratio of Fast Neutron Creation Factor
+   real(r8)            :: L1         ! High Energy   Soil Attenuation Length
+   real(r8)            :: L2         ! High Energy  Water Attenuation Length
+   real(r8)            :: L3         ! Fast Neutron  Soil Attenuation Length
+   real(r8)            :: L4         ! Fast Neutron Water Attenuation Length
+end type site_metadata
+
+type(site_metadata), allocatable, dimension(:) :: cosmos_metadata
+integer  :: nSites, siteIndx, obsindx
+real(r8) :: bd, lattwat, N, alpha, L1, L2, L3, L4   
+
+!-----------------------------------------------------------------------
+! start of executable code
+!-----------------------------------------------------------------------
+
+call initialize_utilities('COSMOS_level2_to_obs')
+
+! Print module information to log file and stdout.
+call register_module(source, revision, revdate)
+
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "COSMOS_to_obs_nml", iunit)
+read(iunit, nml = COSMOS_to_obs_nml, iostat = rcio)
+call check_namelist_read(iunit, rcio, "COSMOS_to_obs_nml")
+
+! Record the namelist values used for the run ...
+if (do_nml_file()) write(nmlfileunit, nml=COSMOS_to_obs_nml)
+if (do_nml_term()) write(     *     , nml=COSMOS_to_obs_nml)
+
+! time setup
+call set_calendar_type(GREGORIAN)
+prev_time = set_time(0, 0)
+
+! Read the COSMOS metadata/parameters for each site.
+! These will be added to the metadata for neutron intensity observations.
+nSites   = read_cosmos_metadata(site_metadata_file)
+siteIndx = find_site_index(sitename)
+bd       = cosmos_metadata(siteIndx)%bd
+lattwat  = cosmos_metadata(siteIndx)%lattwat
+N        = cosmos_metadata(siteIndx)%N
+alpha    = cosmos_metadata(siteIndx)%alpha
+L1       = cosmos_metadata(siteIndx)%L1
+L2       = cosmos_metadata(siteIndx)%L2
+L3       = cosmos_metadata(siteIndx)%L3
+L4       = cosmos_metadata(siteIndx)%L4
+
+if (verbose) print *, 'COSMOS site located at lat, lon, elev  =', &
+                       cosmos_metadata(siteIndx)%latitude, &
+                       cosmos_metadata(siteIndx)%longitude, &
+                       cosmos_metadata(siteIndx)%elevation
+
+! We need to know the maximum number of observations in the input file.
+! Each line has info for a single observation we want (COSMOS neutron counts).
+! Each observation in this series will have a single
+! observation value, its standard deviation, and a quality control flag.  
+! Initialize two empty observations - one to track location
+! in observation sequence - the other is for the new observation.
+
+call find_textfile_dims(text_input_file, max_obs)
+
+iunit = open_file(text_input_file, 'formatted', 'read')
+if (verbose) print *, 'opened input file ' // trim(text_input_file)
+
+num_copies = 1
+num_qc     = 1
+first_obs  = .true.
+
+call static_init_obs_sequence()
+call init_obs(        obs,      num_copies, num_qc)
+call init_obs(        prev_obs, num_copies, num_qc)
+call init_obs_sequence(obs_seq, num_copies, num_qc, max_obs)
+
+! the first one needs to contain the string 'observation' and the
+! second needs the string 'QC'.
+call set_copy_meta_data(obs_seq, 1, 'observation')
+call set_qc_meta_data(  obs_seq, 1, 'COSMOS QC')
+
+! The first line describes all the fields ... column headers, if you will
+call decode_header(iunit)
+
+obsloop: do iline = 2,max_obs
+
+   ! read in entire text line into a buffer
+   read(iunit,'(A)',iostat=rcio) input_line
+   if (rcio < 0) exit obsloop
+   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,'main', string1, source, revision, revdate)
+   endif
+
+   ! parse the line into the cosmos structure (including the observation time)
+   call stringparse(input_line, iline)
+
+   if (iline <= 2) then
+      write(*,*)''
+      write(*,*)'Check of the first observation: (column,string,value)'
+      write(*,*)cosmos%dateindex, cosmos%datestring, cosmos%year, cosmos%month, cosmos%day
+      write(*,*)cosmos%timeindex, cosmos%timestring, cosmos%hour, cosmos%minute
+      write(*,*)cosmos%neutronindex    , cosmos%neutronstring     , cosmos%neutron
+      write(*,*)cosmos%neutronstdindex , cosmos%neutronstdstring  , cosmos%neutronstd
+      call print_date(cosmos%time_obs, 'observation date is')
+      call print_time(cosmos%time_obs, 'observation time is')
+   end if
+
+   if (verbose) call print_date(cosmos%time_obs, 'obs time is')
+
+   call get_time(cosmos%time_obs, osec, oday)
+
+   ! 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.
+   ! Quality Control Flags by Rafael Rosolem (rosolem at email.arizona.edu)
+   ! QC flags define as: 0 = BEST     (corrected for water vapor)
+   !                     1 = OK       (NOT corrected for water vapor)
+   !                     2 = BAD      (not reliable/consistent)
+   !                     3 = MISSING  (missing data, i.e. -9999)
+
+   if (cosmos%qc <= maxgoodqc) then   ! COSMOS neutron counts
+      oerr = cosmos%neutronstd
+      qc   = real(cosmos%qc,r8)
+
+      call set_cosmos_metadata(obsindx, bd, lattwat, N, alpha, L1, L2, L3, L4) 
+
+      call create_3d_obs(cosmos_metadata(siteIndx)%latitude, &
+                         cosmos_metadata(siteIndx)%longitude, 0.0_r8, VERTISSURFACE, &
+           cosmos%neutron, COSMOS_NEUTRON_INTENSITY, oerr, oday, osec, qc, obsindx, obs)
+
+      call add_obs_to_seq(obs_seq, obs, cosmos%time_obs, prev_obs, prev_time, first_obs)
+   endif
+
+end do obsloop
+
+! if we added any obs to the sequence, write it out to a file now.
+if ( get_num_obs(obs_seq) > 0 ) then
+   if (verbose) print *, 'writing obs_seq, obs_count = ', get_num_obs(obs_seq)
+   call write_obs_seq(obs_seq, obs_out_file)
+endif
+
+! end of main program
+call finalize_utilities()
+
+contains
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!   create_3d_obs - subroutine that is used to create an observation
+!                   type from observation data.  
+!
+!   NOTE: assumes the code is using the threed_sphere locations module, 
+!         that the observation has a single data value and a single
+!         qc value.
+!
+!    lat   - latitude of observation
+!    lon   - longitude of observation
+!    vval  - vertical coordinate
+!    vkind - kind of vertical coordinate (pressure, level, etc)
+!    obsv  - observation value
+!    okind - observation kind
+!    oerr  - observation error
+!    day   - gregorian day
+!    sec   - gregorian second
+!    qc    - quality control value
+!    key   - index to metadata in obs_def_COSMOS_mod arrays
+!    obs   - observation type
+!
+!    extended from the observations/utilities/obs_utilities_mod.f90 v 5601
+!    to support the extra metadata -- TJH
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine create_3d_obs(lat, lon, vval, vkind, obsv, okind, oerr, day, sec, qc, key, obs)
+integer,        intent(in)    :: okind, vkind, day, sec
+real(r8),       intent(in)    :: lat, lon, vval, obsv, oerr, qc
+integer,        intent(in)    :: key
+type(obs_type), intent(inout) :: obs
+
+real(r8)              :: obs_val(1), qc_val(1)
+type(obs_def_type)    :: obs_def
+
+call set_obs_def_location(obs_def, set_location(lon, lat, vval, vkind))
+call set_obs_def_kind(obs_def, okind)
+call set_obs_def_time(obs_def, set_time(sec, day))
+call set_obs_def_error_variance(obs_def, oerr * oerr)
+call set_obs_def_key(obs_def, key)
+call set_obs_def(obs, obs_def)
+
+obs_val(1) = obsv
+call set_obs_values(obs, obs_val)
+qc_val(1)  = qc
+call set_qc(obs, qc_val)
+
+end subroutine create_3d_obs

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list