[Dart-dev] [3941] DART/trunk: This provides support for CODAR observations.
nancy at ucar.edu
nancy at ucar.edu
Tue Jun 23 20:50:06 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090623/aec89aa2/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/MITgcm_ocean/create_ocean_obs.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/create_ocean_obs.f90 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/create_ocean_obs.f90 2009-06-24 02:50:06 UTC (rev 3941)
@@ -13,13 +13,13 @@
! 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
+! an assimilation 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 dart_MITocean_mod, only : real_obs_sequence
use utilities_mod, only : initialize_utilities, register_module, &
do_output, logfileunit, &
error_handler, timestamp, E_ERR, E_MSG, &
@@ -45,6 +45,7 @@
integer :: max_num = 800000
character(len = 129) :: fname = 'raw_ocean_obs.txt'
character(len = 129) :: output_name = 'raw_ocean_obs_seq.out'
+logical :: codar = .false.
real(r8) :: lon1 = 0.0_r8, & ! lower longitude bound
lon2 = 360.0_r8, & ! upper longitude bound
@@ -52,7 +53,7 @@
lat2 = 90.0_r8 ! upper latitude bound
namelist /create_ocean_obs_nml/ year, month, day, tot_days, max_num, &
- fname, output_name, lon1, lon2, lat1, lat2
+ fname, output_name, lon1, lon2, lat1, lat2, codar
! ----------------------------------------------------------------------
! start of executable program code
@@ -75,13 +76,9 @@
if (do_output()) write(logfileunit, nml=create_ocean_obs_nml)
if (do_output()) write( * , nml=create_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
-
! The file is read and parsed into a DART observation sequence linked list
seq = real_obs_sequence(fname, year, month, day1, max_num, &
- lon1, lon2, lat1, lat2)
+ lon1, lon2, lat1, lat2, codar)
call write_obs_seq(seq, output_name)
Modified: DART/trunk/models/MITgcm_ocean/create_ocean_obs.html
===================================================================
--- DART/trunk/models/MITgcm_ocean/create_ocean_obs.html 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/create_ocean_obs.html 2009-06-24 02:50:06 UTC (rev 3941)
@@ -150,7 +150,7 @@
<PRE class="indent1">
types_mod
utilities_mod
-ocean_obs_mod
+dart_MITocean_mod
obs_sequence_mod
</PRE>
Modified: DART/trunk/models/MITgcm_ocean/create_ocean_obs.nml
===================================================================
--- DART/trunk/models/MITgcm_ocean/create_ocean_obs.nml 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/create_ocean_obs.nml 2009-06-24 02:50:06 UTC (rev 3941)
@@ -9,5 +9,6 @@
lon1 = 0.0,
lon2 = 360.0,
lat1 = -90.0,
- lat2 = 90.0 /
+ lat2 = 90.0,
+ codar = .false. /
Copied: DART/trunk/models/MITgcm_ocean/dart_MITocean_mod.f90 (from rev 3938, DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90)
===================================================================
--- DART/trunk/models/MITgcm_ocean/dart_MITocean_mod.f90 (rev 0)
+++ DART/trunk/models/MITgcm_ocean/dart_MITocean_mod.f90 2009-06-24 02:50:06 UTC (rev 3941)
@@ -0,0 +1,337 @@
+! 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 dart_MITocean_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+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, &
+ set_obs_def_key
+use time_manager_mod, only : time_type, get_date, set_time, GREGORIAN, &
+ set_date, set_calendar_type, get_time, &
+ print_date, print_time, operator(>=)
+use utilities_mod, only : get_unit, open_file, close_file, file_exist, &
+ register_module, error_handler, &
+ E_ERR, E_MSG, timestamp, is_longitude_between
+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, &
+ get_first_obs, get_last_obs, get_obs_def
+use obs_kind_mod, only : get_obs_kind_index
+use obs_def_ocean_mod, only : set_radial_vel
+
+implicit none
+private
+
+public :: real_obs_sequence
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+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, &
+ lon1, lon2, lat1, lat2, codar)
+!------------------------------------------------------------------------------
+! this function is to prepare data to DART sequence format
+!
+character(len=129), intent(in) :: obsfile
+integer, intent(in) :: year, month, day, max_num
+real(r8), intent(in) :: lon1, lon2, lat1, lat2
+logical, intent(in) :: codar
+
+type(obs_sequence_type) :: real_obs_sequence
+
+
+type(obs_def_type) :: obs_def
+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, id, defkey
+integer :: which_vert, obstype
+
+real (r8) :: lon, lat, vloc, obs_value
+real (r8) :: aqc, var2, lonc, angle
+type(time_type) :: time, pre_time
+
+character(len = 32) :: obs_kind_name
+character(len = 80) :: label
+character(len = 129) :: copy_meta_data, qc_meta_data
+
+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)
+
+obs_num = 0
+iskip = 0
+defkey = 1
+
+! loop over all observations within the file
+!------------------------------------------------------------------------------
+
+obsloop: do
+
+ if (codar) then
+ read(obs_unit,*,end=200) lon, lat, vloc, angle, obs_value, which_vert, var2, aqc, &
+ obs_kind_name, startdate1, startdate2, id
+ else
+ read(obs_unit,*,end=200) lon, lat, vloc, obs_value, which_vert, var2, aqc, &
+ obs_kind_name, startdate1, startdate2
+ endif
+
+ !print*,''
+ !print*,' Observation ', obs_num+1
+ !print*,' lon lat vloc obs_value ',lon, lat, vloc, obs_value
+ !print*,' which_vert var2 aqc ',which_vert, var2, aqc
+ !print*,' obs_kind_name ',obs_kind_name
+ !print*,' date1 date2 ',startdate1, startdate2
+ !if (codar) print*, 'angle id ',angle, id
+
+ ! 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 latitude is not outside valid limits
+ if ((lat > 90.0_r8) .or. (lat < -90.0_r8)) then
+ write(*,*) 'invalid location. lon,lat = ', lon, lat
+ iskip = iskip + 1
+ cycle obsloop
+ endif
+
+ ! reject observations outside the bounding box
+ if(lat < lat1 .or. lat > lat2 .or. &
+ .not. is_longitude_between(lon, lon1, lon2)) then
+ iskip = iskip + 1
+ cycle obsloop
+ endif
+
+ ! and now make sure lon is between 0 and 360, so the
+ ! dart locations module is happy. using modulo is important
+ ! here; plain mod() does not handle negative numbers the same.
+ lon = modulo(lon, 360.0_r8)
+
+
+ ! assign each observation the correct observation type
+ obstype = get_obs_kind_index(obs_kind_name)
+ if(obstype < 1) then
+ print*, 'unknown observation type [',trim(obs_kind_name),'] ... skipping ...'
+ cycle obsloop
+ !else
+ ! print*,trim(obs_kind_name),' is ',obstype
+ 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, obstype, which_vert, seconds, days, &
+ codar, defkey, id, angle)
+
+ if(obs_num == 1) then ! for the first observation
+
+ call insert_obs_in_seq(real_obs_sequence, obs)
+
+ else ! not the first observation
+
+ if(time >= pre_time) then ! same time or later than prev
+
+ call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
+
+ else ! earlier than prev obs, must search from start
+
+ call insert_obs_in_seq(real_obs_sequence, obs)
+
+ endif
+
+ endif
+
+ call copy_obs(prev_obs, obs)
+ pre_time = time
+
+end do obsloop
+
+200 continue
+
+close(obs_unit)
+
+! Print a little summary
+print*, 'obs used = ', obs_num, ' obs skipped = ', iskip
+
+if ( get_first_obs(real_obs_sequence, obs) ) then
+ call get_obs_def(obs, obs_def)
+ pre_time = get_obs_def_time(obs_def)
+ call print_time(pre_time,' first time in sequence is ')
+ call print_date(pre_time,' first date in sequence is ')
+endif
+if( get_last_obs(real_obs_sequence, obs)) then
+ call get_obs_def(obs, obs_def)
+ time = get_obs_def_time(obs_def)
+ call print_time(time,' last time in sequence is ')
+ call print_date(time,' last date in sequence is ')
+endif
+print*, ''
+
+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, &
+ codar, defkey, id, angle)
+!------------------------------------------------------------------------------
+integer, intent(in) :: num_copies, num_qc
+type(obs_type), intent(inout) :: obs
+real(r8), intent(in) :: lon, lat, vloc, obs_value, var2, aqc, angle
+integer, intent(in) :: obs_kind, which_vert, seconds, days, id
+logical, intent(in) :: codar
+integer, intent(inout) :: defkey
+
+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)
+if (codar) then
+ ! this routine increments defkey and stores the private info in
+ ! the ocean module until time to write.
+ call set_radial_vel(defkey, id, angle)
+ call set_obs_def_key(obsdef0, defkey)
+endif
+
+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 dart_MITocean_mod
Modified: DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m
===================================================================
--- DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/matlab/CheckICS.m 2009-06-24 02:50:06 UTC (rev 3941)
@@ -50,6 +50,12 @@
mitbase = '/ptmp/thoar/MITgcm/ics';
fname = sprintf('%s/filter_ics',mitbase);
+ ncopies = 40;
+
+ mitbase = '/ptmp/thoar/MITgcm/obs_gliders';
+ fname = sprintf('%s/filter_restart',mitbase);
+ ncopies = 20;
+
fid = fopen(fname,'rb','ieee-be');
for i = 1:ncopies
@@ -86,7 +92,8 @@
set(gca,'YDir','normal')
title(sprintf('S copy %02d level %03d day %d sec %d\n',i,ilev,T(2),T(1)))
colorbar;
- pause(0.1)
+ display('Pausing, hit any key to continue ...')
+ pause
end
% Plot temperature for each ensemble member
@@ -98,7 +105,8 @@
set(gca,'YDir','normal')
title(sprintf('T copy %02d level %03d day %d sec %d\n',i,ilev,T(2),T(1)))
colorbar;
- pause(0.1)
+ display('Pausing, hit any key to continue ...')
+ pause
end
% Read the U current component for each ensemble member
@@ -110,7 +118,8 @@
set(gca,'YDir','normal')
title(sprintf('U copy %02d level %03d day %d sec %d\n',i,ilev,T(2),T(1)))
colorbar;
- pause(0.1)
+ display('Pausing, hit any key to continue ...')
+ pause
end
% Read the V current component for each ensemble member
@@ -122,7 +131,8 @@
set(gca,'YDir','normal')
title(sprintf('V copy %02d level %03d day %d sec %d\n',i,ilev,T(2),T(1)))
colorbar;
- pause(0.1)
+ display('Pausing, hit any key to continue ...')
+ pause
end
% Read the sea surface height component for each ensemble member
@@ -133,7 +143,8 @@
set(gca,'YDir','normal')
title(sprintf('SSH copy %02d level %03d day %d sec %d\n',i,ilev,T(2),T(1)))
colorbar;
- pause(0.1)
+ display('Pausing, hit any key to continue ...')
+ pause
end
Modified: DART/trunk/models/MITgcm_ocean/matlab/Check_ud.m
===================================================================
--- DART/trunk/models/MITgcm_ocean/matlab/Check_ud.m 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/matlab/Check_ud.m 2009-06-24 02:50:06 UTC (rev 3941)
@@ -2,8 +2,21 @@
%
% fname = 'assim_model_state_ud.0001';
% dsize = [256 225 40];
+% x = Check_ud(fname,dsize);
%
-% x = Check_ud(fname,dsize);
+% fname = '/ptmp/thoar/MITgcm/Run-2/advance_temp0/assim_model_state_ud';
+% dsize = [256 225 40];
+% xi = Check_ud(fname,dsize);
+%
+% fname = '/fs/image/home/thoar/SVN/DART/models/MITgcm_ocean/ibrahim/assim_model_state_ud';
+% dsize = [256 225 40];
+% xg = Check_ud(fname,dsize);
+%
+% Sdiff = xi.S - xg.S; [min(Sdiff(:)) max(Sdiff(:))]
+% Tdiff = xi.T - xg.T; [min(Tdiff(:)) max(Tdiff(:))]
+% Udiff = xi.U - xg.U; [min(Udiff(:)) max(Udiff(:))]
+% Vdiff = xi.V - xg.V; [min(Vdiff(:)) max(Vdiff(:))]
+% Etadiff = xi.Eta - xg.Eta; [min(Etadiff(:)) max(Etadiff(:))]
modelsize = prod(dsize);
@@ -40,4 +53,6 @@
error('second record mismatch')
end
-chunk.fname = fname;
+chunk.fname = fname;
+chunk.seconds = seconds;
+chunk.days = days;
Deleted: DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90
===================================================================
--- DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/ocean_obs_mod.f90 2009-06-24 02:50:06 UTC (rev 3941)
@@ -1,317 +0,0 @@
-! 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$
-! $Id$
-! $Revision$
-! $Date$
-
-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, &
- print_date, print_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, &
- get_first_obs, get_last_obs, get_obs_def
-
-use obs_kind_mod, only : get_obs_kind_index
-
-implicit none
-private
-
-public :: real_obs_sequence
-
-! version controlled file description for error handling, do not edit
-character(len=128), parameter :: &
- source = "$URL$", &
- revision = "$Revision$", &
- revdate = "$Date$"
-
-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, &
- 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
-real(r8), intent(in) :: lon1, lon2, lat1, lat2
-
-type(obs_sequence_type) :: real_obs_sequence
-
-
-type(obs_def_type) :: obs_def
-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 :: which_vert, obstype
-
-real (r8) :: lon, lat, vloc, obs_value
-real (r8) :: aqc, var2, lonc
-type(time_type) :: time, pre_time
-
-character(len = 32) :: obs_kind_name
-character(len = 80) :: label
-character(len = 129) :: copy_meta_data, qc_meta_data
-
-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)
-
-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, &
- obs_kind_name, startdate1, startdate2
-
- !print*,''
- !print*,' Observation ', obs_num+1
- !print*,' lon lat vloc obs_value ',lon, lat, vloc, obs_value
- !print*,' which_vert var2 aqc ',which_vert, var2, aqc
- !print*,' obs_kind_name ',obs_kind_name
- !print*,' date1 date2 ',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
- obstype = get_obs_kind_index(obs_kind_name)
- if(obstype < 1) then
- print*, 'unknown observation type [',trim(obs_kind_name),'] ... skipping ...'
- cycle obsloop
- !else
- ! print*,trim(obs_kind_name),' is ',obstype
- 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, obstype, 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 a little summary
-print*, 'obs used = ', obs_num, ' obs skipped = ', iskip
-
-if ( get_first_obs(real_obs_sequence, obs) ) then
- call get_obs_def(obs, obs_def)
- pre_time = get_obs_def_time(obs_def)
- call print_time(pre_time,' first time in sequence is ')
- call print_date(pre_time,' first date in sequence is ')
-endif
-if( get_last_obs(real_obs_sequence, obs)) then
- call get_obs_def(obs, obs_def)
- time = get_obs_def_time(obs_def)
- call print_time(time,' last time in sequence is ')
- call print_date(time,' last date in sequence is ')
-endif
-print*, ''
-
-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
Modified: DART/trunk/models/MITgcm_ocean/work/path_names_create_ocean_obs
===================================================================
--- DART/trunk/models/MITgcm_ocean/work/path_names_create_ocean_obs 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/models/MITgcm_ocean/work/path_names_create_ocean_obs 2009-06-24 02:50:06 UTC (rev 3941)
@@ -1,6 +1,6 @@
assim_model/assim_model_mod.f90
models/MITgcm_ocean/create_ocean_obs.f90
-models/MITgcm_ocean/ocean_obs_mod.f90
+models/MITgcm_ocean/dart_MITocean_mod.f90
models/MITgcm_ocean/model_mod.f90
common/types_mod.f90
location/threed_sphere/location_mod.f90
Modified: DART/trunk/obs_def/obs_def_ocean_mod.f90
===================================================================
--- DART/trunk/obs_def/obs_def_ocean_mod.f90 2009-06-22 20:21:54 UTC (rev 3940)
+++ DART/trunk/obs_def/obs_def_ocean_mod.f90 2009-06-24 02:50:06 UTC (rev 3941)
@@ -1,8 +1,14 @@
-! Data Assimilation Research Testbed -- DART
+! Data Assimilation Research Testbed -- DART source
! 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$
+! $Id$
+! $Revision$
+! $Date$
+
! BEGIN DART PREPROCESS KIND LIST
!SALINITY, KIND_SALINITY, COMMON_CODE
!TEMPERATURE, KIND_TEMPERATURE, COMMON_CODE
@@ -35,11 +41,666 @@
!SATELLITE_INFRARED_SST, KIND_TEMPERATURE, COMMON_CODE
!SATELLITE_SSH, KIND_SEA_SURFACE_HEIGHT, COMMON_CODE
!SATELLITE_SSS, KIND_SALINITY, COMMON_CODE
+!CODAR_RADIAL_VELOCITY, KIND_VELOCITY
! END DART PREPROCESS KIND LIST
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
+! From Ibrahim - 19 May 2009
+! "Concerning the radar observations: we can use a format very similar to what we
+! have been using for all the other types of observations:
+!
+!LAT LON DEPTH ANGLE OBS_value Z_LOC ERR_obs QC_flag TYPE_obs STARTDATE1 STARTDATE2
+!
+! where here we only add "ANGLE": the angle of the radial observation.
+!
+! Now what we need is to compute the radial velocity R_model from the state vector
+! (which already contains the horizontal velocities U and V) using
+!
+! 'OBS_value' = U*cos(angle) + V*sin(angle)
+!
+! The angles are in degrees (0 from the east) and the observations in m/s.
+! so to compute 'OBS_value' we need to read the observation file to get the ANGLE,
+! and that's basically the only difference with the assimilation of the other
+! types of observations. We also need to introduce the new obs type
+! (CODAR_RADIAL_VELOCITY)."
+!-----------------------------------------------------------------------------
+! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
+! use obs_def_ocean_mod, only : write_radial_vel, read_radial_vel, &
+! interactive_radial_vel, get_expected_radial_vel
+! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
+!-----------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------------
+! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
+! case(CODAR_RADIAL_VELOCITY)
+! call get_expected_radial_vel(state, location, obs_def%key, obs_val, istatus)
+! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
+!-----------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------------
+! BEGIN DART PREPROCESS READ_OBS_DEF
+! case(CODAR_RADIAL_VELOCITY)
+! call read_radial_vel(obs_def%key, ifile, fileformat)
+! END DART PREPROCESS READ_OBS_DEF
+!-----------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------------
+! BEGIN DART PREPROCESS WRITE_OBS_DEF
+! case(CODAR_RADIAL_VELOCITY)
+! call write_radial_vel(obs_def%key, ifile, fileformat)
+! END DART PREPROCESS WRITE_OBS_DEF
+!-----------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------------
+! BEGIN DART PREPROCESS INTERACTIVE_OBS_DEF
+! case(CODAR_RADIAL_VELOCITY)
+! call interactive_radial_vel(obs_def%key)
+! END DART PREPROCESS INTERACTIVE_OBS_DEF
+!-----------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------------
+! BEGIN DART PREPROCESS MODULE CODE
+module obs_def_ocean_mod
+
+use types_mod, only : r8, missing_r8, PI, deg2rad
+use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, &
+ check_namelist_read, find_namelist_in_file, &
+ nmlfileunit, do_output, do_nml_file, do_nml_term
+use location_mod, only : location_type, write_location, read_location, &
+ interactive_location, get_location
+use assim_model_mod, only : interpolate
+use obs_kind_mod, only : KIND_U_CURRENT_COMPONENT, &
+ KIND_V_CURRENT_COMPONENT
+
+implicit none
+private
+
+public :: read_radial_vel, write_radial_vel, interactive_radial_vel, &
+ get_expected_radial_vel, get_obs_def_radial_vel, set_radial_vel
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+logical :: module_initialized = .false.
+
+character(len=129) :: msgstring ! For error message content
+
+! Derived type for radial velocity. Contains auxiliary information stored
+! with each obs of this type; used to compute the forward operator.
+! See more extensive comments in the interactive_radial_vel() routine for
+! expected units, etc. The instrument ID is currently unused, but may be
+! useful for post processing or diagnostics.
+
+type radial_vel_type
+ private
+ integer :: instrument_id ! ID number for data source
+ real(r8) :: beam_angle ! angle of beam
+end type radial_vel_type
+
+! Cumulative index into radial velocity metadata array
+integer :: velkeycount = 0
+
+! namelist items
+integer :: max_radial_vel_obs = 1000000
+
+namelist /obs_def_ocean_mod_nml/ max_radial_vel_obs
+
+! Module global storage for auxiliary obs data, allocated in init routine
+type(radial_vel_type), allocatable :: radial_vel_data(:)
+
+contains
+
+
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! Start of executable routines
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+
+
+subroutine initialize_module
+
+! Called once to set values and allocate space
+
+integer :: iunit, io, rc
+
+! Prevent multiple calls from executing this code more than once.
+if (module_initialized) return
+
+module_initialized = .true.
+
+! Log the version of this source file.
+call register_module(source, revision, revdate)
+
+! Read the namelist entry.
+call find_namelist_in_file("input.nml", "obs_def_ocean_mod_nml", iunit)
+read(iunit, nml = obs_def_ocean_mod_nml, iostat = io)
+call check_namelist_read(iunit, io, "obs_def_ocean_mod_nml")
+
+! Record the namelist values used for the run ...
+if (do_nml_file()) write(nmlfileunit, nml=obs_def_ocean_mod_nml)
+if (do_nml_term()) write( * , nml=obs_def_ocean_mod_nml)
+
+! Allocate space for the auxiliary information associated with each obs
+! This code must be placed after reading the namelist, so the user can
+! increase or decrease the number of obs supported and use more or less
+! memory at run time.
+allocate(radial_vel_data(max_radial_vel_obs), stat = rc)
+if (rc /= 0) then
+ write(msgstring, *) 'initial allocation failed for radial vel obs data,', &
+ 'itemcount = ', max_radial_vel_obs
+ call error_handler(E_ERR,'initialize_module', msgstring, &
+ source, revision, revdate)
+endif
+
+end subroutine initialize_module
+
+
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! Radial velocity section
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+
+
+subroutine read_radial_vel(velkey, ifile, fform)
+
+! Main read subroutine for the radial velocity observation auxiliary data.
+
+integer, intent(out) :: velkey
+integer, intent(in) :: ifile
+character(len=*), intent(in), optional :: fform
+
+character(len=8) :: header
+logical :: is_asciifile
+integer :: instrument_id
+real(r8) :: beam_angle
+integer :: oldkey
+
+! instrument id: Arbitrary number to distinguish different sources of
+! the data. Not used in any of this code, but carried along in case
+! it is useful in postprocessing or diagnostics.
+
+if ( .not. module_initialized ) call initialize_module
+
+is_asciifile = ascii_file_format(fform)
+
+if (is_asciifile) then
+ ! Read the character identifier for verbose formatted output
+ read(ifile, FMT="(a5)") header
+ if(header /= 'CODAR') then
+ call error_handler(E_ERR,'read_radial_vel', &
+ "Expected header 'CODAR' in input file", &
+ source, revision, revdate)
+ endif
+endif
+
+! read_location is a DART library routine that expects an optional string
+! arg, but the other two read routines are local to this module and we can
+! tell them exactly what format to be reading because we already know it.
+instrument_id = read_instrument_id (ifile, is_asciifile)
+beam_angle = read_beam_angle (ifile, is_asciifile)
+
+! Read in the velkey for this particular observation, however, it will
+! be discarded and a new, unique key will be generated in the set routine.
+if (is_asciifile) then
+ read(ifile, *) oldkey
+else
+ read(ifile) oldkey
+endif
+
+! Generate new unique radial velocity observation key, and set the contents
+! of the private defined type.
+call set_radial_vel(velkey, instrument_id, beam_angle )
+
+end subroutine read_radial_vel
+
+
+!----------------------------------------------------------------------
+
+
+subroutine write_radial_vel(velkey, ifile, fform)
+
+! Write radial velocity auxiliary information to the obs_seq file.
+
+integer, intent(in) :: velkey, ifile
+character(len=*), intent(in), optional :: fform
+
+logical :: is_asciifile
+integer :: instrument_id
+real(r8) :: beam_angle
+
+if ( .not. module_initialized ) call initialize_module
+
+! Simple error check on key number before accessing the array
+call velkey_out_of_range(velkey)
+
+is_asciifile = ascii_file_format(fform)
+
+if (is_asciifile) then
+ ! Write the 5 character identifier for verbose formatted output
+ write(ifile, "('CODAR')")
+endif
+
+! Extract the values for this key and call the appropriate write routines.
+instrument_id = radial_vel_data(velkey)%instrument_id
+beam_angle = radial_vel_data(velkey)%beam_angle
+
+! These two write routines are local to this module, and we have already
+! figured out if it is a unformatted/binary file or formatted/ascii, so
+! go ahead and pass that info directly down to the routines.
+call write_instrument_id(ifile, instrument_id, is_asciifile)
+call write_beam_angle(ifile, beam_angle, is_asciifile)
+
+! Write out the velkey used for this run, however this will be discarded
+! when this observation is read in and a new key will be generated.
+! (It may be useful for correlating error messages or identifying particular
+! observations so we are leaving it as part of the aux data.)
+if (is_asciifile) then
+ write(ifile, *) velkey
+else
+ write(ifile) velkey
+endif
+
+end subroutine write_radial_vel
+
+
+!----------------------------------------------------------------------
+
+
+function read_instrument_id(ifile, is_asciiformat)
+
+! Reads instrument id from file that was written by write_instrument_id.
+! See read_radial_vel for additional discussion.
+
+integer, intent(in) :: ifile
+logical, intent(in) :: is_asciiformat
+integer :: read_instrument_id
+
+character(len=5) :: header
+integer :: instrument_id
+
+if ( .not. module_initialized ) call initialize_module
+
+
+if (is_asciiformat) then
+ read(ifile, "(a5)" ) header
+
+ if(header /= 'InsID') then
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list