[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