[Dart-dev] [3369] DART/trunk/models/MITgcm_ocean: adding ability to convert synthetic observations to a DART observation sequence file .

nancy at ucar.edu nancy at ucar.edu
Wed May 21 18:20:20 MDT 2008


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080521/e8fc4b5d/attachment.html
-------------- next part --------------
Added: DART/trunk/models/MITgcm_ocean/create_ocean_obs.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/create_ocean_obs.f90	                        (rev 0)
+++ DART/trunk/models/MITgcm_ocean/create_ocean_obs.f90	2008-05-22 00:20:20 UTC (rev 3369)
@@ -0,0 +1,111 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2007, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+program create_ocean_obs
+
+! <next few lines under version control, do not edit>
+! $URL: http://subversion.ucar.edu/DAReS/DART/trunk/ncep_obs/create_ocean_obs.f90 $
+! $Id: create_ocean_obs.f90 3304 2008-04-09 22:20:02Z nancy $
+! $Revision: 3304 $
+! $Date: 2008-04-09 16:20:02 -0600 (Wed, 09 Apr 2008) $
+
+! Initial program to read the raw ocean observations and insert them
+! into an observation sequence. To make things easy ... we will mandate
+! an assimilatin interval of 1 day - so all observations will be
+! redefined to occur at NOON on the day they were observed.
+
+use types_mod,        only : r8, deg2rad, PI
+use obs_sequence_mod, only : obs_sequence_type, write_obs_seq, &
+                             static_init_obs_sequence, destroy_obs_sequence 
+use    ocean_obs_mod, only : real_obs_sequence
+use    utilities_mod, only : initialize_utilities, register_module, &
+                             do_output, logfileunit, &
+                             error_handler, timestamp, E_ERR, E_MSG, &
+                             find_namelist_in_file, check_namelist_read
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+   source   = "$URL: http://subversion.ucar.edu/DAReS/DART/trunk/ncep_obs/create_ocean_obs.f90 $", &
+   revision = "$Revision: 3304 $", &
+   revdate  = "$Date: 2008-04-09 16:20:02 -0600 (Wed, 09 Apr 2008) $"
+
+type(obs_sequence_type) :: seq
+
+character(len = 129) :: output_name
+character(len = 8 ) :: obsdate
+integer :: iunit, io, ii, day1, kkk, kbeg, kend
+
+! ----------------------------------------------------------------------
+! Declare namelist parameters
+! ----------------------------------------------------------------------
+        
+integer :: year = 1996, month =1, day =1, tot_days = 31
+integer :: max_num = 800000
+character(len = 129) :: fname = 'raw_ocean_obs.txt'
+
+logical :: ADCP = .false., DRIFTERS = .false., GLIDERS = .false., &
+           TMI = .false., SSH = .false.
+
+real(r8) :: lon1 =   0.0_r8,  &   !  lower longitude bound
+            lon2 = 360.0_r8,  &   !  upper longitude bound 
+            lat1 = -90.0_r8,  &   !  lower latitude bound
+            lat2 =  90.0_r8       !  upper latitude bound
+
+namelist /ocean_obs_nml/ year, month, day, tot_days, max_num, &
+        fname, ADCP, DRIFTERS, GLIDERS, TMI, SSH, &
+        lon1, lon2, lat1, lat2
+
+! ----------------------------------------------------------------------
+!  ADPUPA: upper-air reports (mostly radiosonde plus few dropsonde, PIBAL)
+!  AIRCFT: Conv. (AIREP, PIREP) and ASDAR aircraft reports
+!  AIRCAR: ACARS sircraft reports
+!  SATEMP: ATOVS retrived temperature
+!  SFCSHP: SURFACE MARINE reports
+!  ADPSFC: SURFACE LAND SYNOPTIC STATION reports
+!  SATWND: Satellite derived wind reports
+! ----------------------------------------------------------------------
+
+! start of executable program code
+
+call initialize_utilities('create_ocean_obs')
+call register_module(source,revision,revdate)
+
+! Initialize the obs_sequence module ...
+
+call static_init_obs_sequence()
+
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "ocean_obs_nml", iunit)
+read(iunit, nml = ocean_obs_nml, iostat = io)
+call check_namelist_read(iunit, io, "ocean_obs_nml")
+
+! Record the namelist values used for the run ...
+call error_handler(E_MSG,'create_ocean_obs','ocean_obs_nml values are',' ',' ',' ')
+if (do_output()) write(logfileunit, nml=ocean_obs_nml)
+if (do_output()) write(     *     , nml=ocean_obs_nml)
+
+lon1 = min(max(lon1,0.0_r8),360.0_r8)
+lon2 = min(max(lon2,0.0_r8),360.0_r8)
+if ( lon1 > lon2 ) lon2 = lon2 + 360.0_r8
+
+  ! define observation filename
+  write(obsdate, '(i4.4,i2.2)') year, month
+
+  seq = real_obs_sequence(fname, year, month, day1, max_num, &
+        ADCP, DRIFTERS, GLIDERS, TMI, SSH, &
+        lon1, lon2, lat1, lat2)
+
+  output_name = 'obs_seq'//obsdate
+  call write_obs_seq(seq, output_name)
+
+  ! release the memory of the seq.
+  call destroy_obs_sequence(seq)
+
+call timestamp(source,revision,revdate,'end') ! close the log file.
+
+end program create_ocean_obs
+

Added: DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90	                        (rev 0)
+++ DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90	2008-05-22 00:20:20 UTC (rev 3369)
@@ -0,0 +1,322 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2007, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+module ocean_obs_mod
+
+! <next few lines under version control, do not edit>
+! $URL: http://subversion.ucar.edu/DAReS/DART/trunk/ncep_obs/ocean_obs_mod.f90 $
+! $Id: ocean_obs_mod.f90 3349 2008-05-16 19:35:31Z thoar $
+! $Revision: 3349 $
+! $Date: 2008-05-16 13:35:31 -0600 (Fri, 16 May 2008) $
+
+use types_mod,        only : r8, rad2deg, PI
+use obs_def_mod,      only : obs_def_type, get_obs_def_time, read_obs_def, &
+                             write_obs_def, destroy_obs_def, interactive_obs_def, &
+                             copy_obs_def, set_obs_def_time, set_obs_def_kind, &
+                             set_obs_def_error_variance, set_obs_def_location
+use time_manager_mod, only : time_type, get_date, set_time, GREGORIAN, &
+                             set_date, set_calendar_type, get_time, &
+                             operator(==)
+use    utilities_mod, only : get_unit, open_file, close_file, file_exist, &
+                             register_module, error_handler, &
+                             E_ERR, E_MSG, timestamp
+use     location_mod, only : location_type, set_location, VERTISHEIGHT, VERTISSURFACE
+use obs_sequence_mod, only : init_obs_sequence, init_obs, insert_obs_in_seq, &
+                             set_obs_values, set_qc, obs_sequence_type, obs_type, &
+                             copy_obs, set_copy_meta_data, set_qc_meta_data, set_obs_def
+
+use     obs_kind_mod, only : KIND_SALINITY, &
+                             KIND_TEMPERATURE, &
+                             KIND_U_WIND_COMPONENT, &
+                             KIND_V_WIND_COMPONENT, &
+                             KIND_SEA_SURFACE_HEIGHT
+
+use obs_kind_mod, only : SALINITY
+use obs_kind_mod, only : TEMPERATURE
+use obs_kind_mod, only : U_CURRENT_COMPONENT
+use obs_kind_mod, only : V_CURRENT_COMPONENT
+use obs_kind_mod, only : SEA_SURFACE_HEIGHT
+
+implicit none
+private
+
+public :: real_obs_sequence
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+   source   = "$URL: http://subversion.ucar.edu/DAReS/DART/trunk/ncep_obs/ocean_obs_mod.f90 $", &
+   revision = "$Revision: 3349 $", &
+   revdate  = "$Date: 2008-05-16 13:35:31 -0600 (Fri, 16 May 2008) $"
+
+logical, save :: module_initialized = .false.
+
+! set this to true if you want to print out the current time
+! after each N observations are processed, for benchmarking.
+logical :: print_timestamps = .false.
+integer :: print_every_Nth  = 10000
+
+contains
+
+!-------------------------------------------------
+
+function real_obs_sequence (obsfile, year, month, day, max_num, &
+          ADCP, DRIFTERS, GLIDERS, TMI, SSH, &
+          lon1, lon2, lat1, lat2)
+!------------------------------------------------------------------------------
+!  this function is to prepare data to DART sequence format
+!
+character(len=129), intent(in) :: obsfile
+integer,            intent(in) :: year, month, day, max_num
+logical,            intent(in) :: ADCP, DRIFTERS, GLIDERS, TMI, SSH
+real(r8),           intent(in) :: lon1, lon2, lat1, lat2
+
+type(obs_sequence_type) :: real_obs_sequence
+
+
+type(obs_type) :: obs, prev_obs
+integer :: i, num_copies, num_qc
+integer :: days, seconds
+integer :: yy, mn, dd, hh, mm, ss
+integer :: startdate1, startdate2
+integer :: obs_num, calender_type, iskip
+integer :: obs_unit
+integer :: obs_prof, obs_kind, obs_kind_gen, which_vert, obstype
+
+real (r8) :: obs_err, lon, lat, vloc, obs_value
+real (r8) :: aqc, var2, lon2c, lonc
+type(time_type) :: time, pre_time
+
+character(len = 8 ) :: obsdate
+character(len = 80) :: label
+character(len = 129) :: copy_meta_data, qc_meta_data
+character(len = 6 ) :: subset
+logical :: pass
+
+if ( .not. module_initialized ) call initialize_module
+
+num_copies  = 1
+num_qc      = 1
+
+! Initialize an obs_sequence 
+
+call init_obs_sequence(real_obs_sequence, num_copies, num_qc, max_num)
+
+! set meta data of obs_seq
+
+do i = 1, num_copies
+   copy_meta_data = 'observation'  
+   call set_copy_meta_data(real_obs_sequence, i, copy_meta_data)
+end do
+
+do i = 1, num_qc
+   qc_meta_data = 'QC index'
+   call set_qc_meta_data(real_obs_sequence, i, qc_meta_data)
+end do
+
+! Initialize the obs variable
+
+call init_obs(obs, num_copies, num_qc)
+call init_obs(prev_obs, num_copies, num_qc)
+
+! set observation time type
+calender_type = GREGORIAN
+call set_calendar_type(calender_type)
+
+! open observation data file
+
+obs_unit  = get_unit()
+open(unit = obs_unit, file = obsfile, form='formatted', status='old')
+print*, 'input file opened= ', trim(obsfile)
+rewind (obs_unit)
+
+!print*, 'ncep obsdates = ', obsdate
+
+obs_num = 0
+iskip   = 0
+
+!  loop over all observations within the file
+!------------------------------------------------------------------------------
+
+obsloop:  do
+
+   read(obs_unit,*,end=200) lon, lat, vloc, obs_value, which_vert, var2, aqc, &
+                              obstype, startdate1, startdate2
+
+   ! Calculate the DART time from the observation time 
+   yy =     startdate1/10000
+   mn = mod(startdate1/100,100)
+   dd = mod(startdate1    ,100)
+   hh =     startdate2/10000
+   mm = mod(startdate2/100,100)
+   ss = mod(startdate2    ,100)
+   time = set_date(yy,mn,dd,hh,mm,ss)
+   call get_time(time,seconds,days)
+
+   ! verify the location is not outside valid limits
+   if((lon > 360.0_r8) .or. (lon <   0.0_r8) .or.  &
+      (lat >  90.0_r8) .or. (lat < -90.0_r8)) then
+      write(*,*) 'invalid location.  lon,lat = ', lon, lat
+      iskip = iskip + 1
+      cycle obsloop
+   endif
+
+   lonc = lon
+   if (lon2 > 360.0_r8 .and. lon < 180.0_r8) lonc = lon + 360.0_r8
+
+   ! reject observations outside the bounding box
+   if(lat < lat1 .or. lat > lat2 .or. lonc < lon1 .or. lonc > lon2) then
+     iskip = iskip + 1
+     cycle obsloop
+   endif
+
+!   assign each observation the correct observation type
+   if(obstype == 1) then
+     obs_kind_gen = KIND_TEMPERATURE
+     obs_kind     =      TEMPERATURE
+   elseif(obstype == 2) then
+    obs_kind_gen = KIND_SALINITY
+    obs_kind     =      SALINITY
+   elseif(obstype == 3) then
+    obs_kind_gen = KIND_U_WIND_COMPONENT
+    obs_kind     =      U_CURRENT_COMPONENT
+   elseif(obstype == 4) then
+    obs_kind_gen = KIND_V_WIND_COMPONENT
+    obs_kind     =      V_CURRENT_COMPONENT
+   elseif(obstype == 5) then
+    obs_kind_gen = KIND_SEA_SURFACE_HEIGHT
+    obs_kind     =      SEA_SURFACE_HEIGHT
+   else
+      print*, 'unknown observation type ... skipping ...'
+      cycle obsloop
+   endif
+
+   obs_num = obs_num + 1
+
+   ! print a reassuring message after every Nth processed obs.
+   ! if requested, print in the form of a timestamp.  
+   ! the default is just a plain string with the current obs count.
+   if(mod(obs_num, print_every_Nth) == 0) then
+       write(label, *) 'obs count = ', obs_num
+       if (print_timestamps) then
+          call timestamp(string1=label, pos='')
+       else
+          write(*,*) trim(label)
+       endif
+   endif
+   if(obs_num == max_num) then
+      print*, 'Max limit for observation count reached.  Increase value in namelist'
+      stop
+   endif
+   
+!   create the obs_def for this observation, add to sequence
+!------------------------------------------------------------------------------
+   
+   call real_obs(num_copies, num_qc, obs, lon, lat, vloc, obs_value, &
+                 var2, aqc, obs_kind, which_vert, seconds, days)
+   
+   if(obs_num == 1) then ! for the first observation 
+
+     call insert_obs_in_seq(real_obs_sequence, obs)
+     call copy_obs(prev_obs, obs)
+     pre_time = time
+
+   else  !  not the first observation 
+
+     if(time == pre_time) then  ! same time as previous observation
+
+       call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
+       call copy_obs(prev_obs, obs)
+       pre_time = time
+
+    else  ! not the same time
+
+       call insert_obs_in_seq(real_obs_sequence, obs)
+       call copy_obs(prev_obs, obs)
+       pre_time = time
+
+    endif
+
+  endif
+
+end do obsloop
+
+200 continue
+
+close(obs_unit)
+
+print*, 'date ', obsdate, ' obs used = ', obs_num, ' obs skipped = ', iskip
+!print*, 'obs_num= ',obs_num,' skipped= ',obsdate,iskip
+
+end function real_obs_sequence
+
+
+
+subroutine real_obs(num_copies, num_qc, obs, lon, lat, vloc, obs_value, &
+                      var2, aqc, obs_kind, which_vert, seconds, days)
+!------------------------------------------------------------------------------
+integer,        intent(in)    :: num_copies, num_qc
+type(obs_type), intent(inout) :: obs
+real(r8),       intent(in)    :: lon, lat, vloc, obs_value, var2, aqc
+integer,        intent(in)    :: obs_kind, which_vert, seconds, days
+
+integer            :: i
+real(r8)           :: aqc01(1), obs_value01(1)
+type(obs_def_type) :: obsdef0
+
+if ( .not. module_initialized ) call initialize_module
+
+! Does real initialization of an observation type
+
+call real_obs_def(obsdef0, lon, lat, vloc, &
+                    var2, obs_kind, which_vert, seconds, days)
+call set_obs_def(obs, obsdef0)
+
+do i = 1, num_copies
+   obs_value01(1) = obs_value
+   call set_obs_values(obs, obs_value01(1:1) )
+end do
+
+do i = 1, num_qc
+   aqc01(1) = aqc
+   call set_qc(obs, aqc01(1:1))
+end do
+
+end subroutine real_obs
+
+
+
+subroutine real_obs_def(obs_def, lon, lat, vloc, &
+                        var2, obs_kind, which_vert, seconds, days)
+!----------------------------------------------------------------------
+type(obs_def_type), intent(inout) :: obs_def
+real(r8),intent(in) :: lon, lat, vloc, var2
+integer, intent(in) :: obs_kind, which_vert, seconds, days
+
+type(location_type) :: loc0
+
+if ( .not. module_initialized ) call initialize_module
+
+! set obs location
+loc0 = set_location(lon, lat, vloc, which_vert )
+call set_obs_def_location(obs_def, loc0)
+
+! set obs kind
+call set_obs_def_kind(obs_def, obs_kind)
+
+call set_obs_def_time(obs_def, set_time(seconds, days) )
+call set_obs_def_error_variance(obs_def, var2)
+
+end subroutine real_obs_def
+
+
+
+subroutine initialize_module
+!-------------------------------------------------
+call register_module(source, revision, revdate)
+module_initialized = .true.
+end subroutine initialize_module
+
+
+end module ocean_obs_mod

Added: DART/trunk/models/MITgcm_ocean/work/mkmf_create_ocean_obs
===================================================================
--- DART/trunk/models/MITgcm_ocean/work/mkmf_create_ocean_obs	                        (rev 0)
+++ DART/trunk/models/MITgcm_ocean/work/mkmf_create_ocean_obs	2008-05-22 00:20:20 UTC (rev 3369)
@@ -0,0 +1,15 @@
+#!/bin/csh
+#
+# Data Assimilation Research Testbed -- DART
+# Copyright 2004-2007, Data Assimilation Research Section
+# University Corporation for Atmospheric Research
+# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+#
+# <next few lines under version control, do not edit>
+# $URL: http://subversion.ucar.edu/DAReS/DART/trunk/models/cam/work/mkmf_create_ocean_obs $
+# $Id: mkmf_create_ocean_obs 2691 2007-03-11 18:18:09Z thoar $
+# $Revision: 2691 $
+# $Date: 2007-03-11 12:18:09 -0600 (Sun, 11 Mar 2007) $
+#
+../../../mkmf/mkmf -p create_ocean_obs -t ../../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../.." path_names_create_ocean_obs


Property changes on: DART/trunk/models/MITgcm_ocean/work/mkmf_create_ocean_obs
___________________________________________________________________
Name: svn:executable
   + *

Added: DART/trunk/models/MITgcm_ocean/work/path_names_create_ocean_obs
===================================================================
--- DART/trunk/models/MITgcm_ocean/work/path_names_create_ocean_obs	                        (rev 0)
+++ DART/trunk/models/MITgcm_ocean/work/path_names_create_ocean_obs	2008-05-22 00:20:20 UTC (rev 3369)
@@ -0,0 +1,16 @@
+assim_model/assim_model_mod.f90
+models/MITgcm_ocean/create_ocean_obs.f90
+models/MITgcm_ocean/ocean_obs_mod.f90
+models/MITgcm_ocean/model_mod.f90
+common/types_mod.f90
+location/threed_sphere/location_mod.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+obs_def/obs_def_mod.f90
+obs_kind/obs_kind_mod.f90
+obs_model/obs_model_mod.f90
+obs_sequence/obs_sequence_mod.f90
+time_manager/time_manager_mod.f90
+utilities/utilities_mod.f90
+random_nr/random_nr_mod.f90
+random_seq/random_seq_mod.f90
+ensemble_manager/ensemble_manager_mod.f90


More information about the Dart-dev mailing list