[Dart-dev] [3940] DART/trunk/models: Initial import of the POP ocean model V 2.0

nancy at ucar.edu nancy at ucar.edu
Mon Jun 22 14:21:55 MDT 2009


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090622/28b5e4f2/attachment-0001.html 
-------------- next part --------------
Added: DART/trunk/models/POP/README
===================================================================
--- DART/trunk/models/POP/README	                        (rev 0)
+++ DART/trunk/models/POP/README	2009-06-22 20:21:54 UTC (rev 3940)
@@ -0,0 +1,7 @@
+This is a development trunk.
+This is a work in progress and may not even work ... AT ALL.
+Do not use at your own risk ... DO NOT USE!!!
+
+When this model is supported, this file will reflect that.
+
+Tim

Added: DART/trunk/models/POP/dart_pop_mod.f90
===================================================================
--- DART/trunk/models/POP/dart_pop_mod.f90	                        (rev 0)
+++ DART/trunk/models/POP/dart_pop_mod.f90	2009-06-22 20:21:54 UTC (rev 3940)
@@ -0,0 +1,317 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2009, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+module dart_pop_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 dart_pop_mod


Property changes on: DART/trunk/models/POP/dart_pop_mod.f90
___________________________________________________________________
Name: svn:keywords
   + Date Revision Author HeadURL Id
Name: svn:eol-style
   + native

Added: DART/trunk/models/POP/dart_to_pop.f90
===================================================================
--- DART/trunk/models/POP/dart_to_pop.f90	                        (rev 0)
+++ DART/trunk/models/POP/dart_to_pop.f90	2009-06-22 20:21:54 UTC (rev 3940)
@@ -0,0 +1,139 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2009, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+ 
+program dart_to_pop
+
+!----------------------------------------------------------------------
+! purpose: interface between DART and the POP model
+!
+! method: Read DART state vector (in file 'assim_model_state_ic') and 
+!         write out POP "snapshot" files.
+!
+! author: Tim Hoar 5Apr08
+!
+!----------------------------------------------------------------------
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+use        types_mod, only : r4, r8, SECPERDAY
+use    utilities_mod, only : E_ERR, E_WARN, E_MSG, error_handler, open_file, &
+                             initialize_utilities, finalize_utilities, &
+                             find_namelist_in_file, check_namelist_read, &
+                             logfileunit
+use        model_mod, only : sv_to_snapshot_files, static_init_model, &
+                             get_model_size, DARTtime_to_POPtime, &
+                             get_model_time_step, write_data_namelistfile, &
+                             set_model_end_time
+use  assim_model_mod, only : open_restart_read, aread_state_restart, close_restart
+use time_manager_mod, only : time_type, get_time, print_time, print_date, &
+                             operator(-)
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+   source   = "$URL$", &
+   revision = "$Revision$", &
+   revdate  = "$Date$"
+
+character (len = 128) :: file_in  = 'assim_model_state_ic'
+
+!------------------------------------------------------------------
+! The time manager namelist variables
+! some types/etc come from   <POPsource>/pkg/cal/cal.h
+! some useful insight from   cal_set.F, cal_convdate.F
+!
+! startDate_1 (integer) yyyymmdd   "start date of the integration"
+! startDate_2 (integer) hhmmss
+!------------------------------------------------------------------
+
+character(len=9) :: TheCalendar = 'gregorian'
+integer          :: startDate_1 = 19530101
+integer          :: startDate_2 =          60000
+logical          :: calendarDumps = .false.
+
+NAMELIST /CAL_NML/ TheCalendar, startDate_1, startDate_2, calendarDumps
+
+!----------------------------------------------------------------------
+
+integer               :: iunit, io, x_size
+integer               :: secs, days
+type(time_type)       :: model_time, adv_to_time
+type(time_type)       :: model_timestep, offset
+real(r8), allocatable :: statevector(:)
+
+!----------------------------------------------------------------------
+
+call initialize_utilities('dart_to_pop')
+call static_init_model()
+
+! POP calendar information. The namelist is already read in 
+! static_init_model(), so no further bulletproofing is needed here.
+call find_namelist_in_file("data.cal", "CAL_NML", iunit)
+read(iunit, nml = CAL_NML, iostat = io)
+call check_namelist_read(iunit, io, "CAL_NML")
+
+x_size = get_model_size()
+allocate(statevector(x_size))
+
+!----------------------------------------------------------------------
+! Reads the valid time, the state, and the target time.
+!----------------------------------------------------------------------
+
+iunit = open_restart_read(file_in)
+call aread_state_restart(model_time, statevector, iunit, adv_to_time)
+call close_restart(iunit)
+
+!----------------------------------------------------------------------
+! update the CAL_NML variables so we can rewrite that namelist.
+! Ultimately, we want to keep data:&PARM03:startTime = 0.,
+!----------------------------------------------------------------------
+
+call DARTtime_to_POPtime(  model_time, startDate_1, startDate_2)
+call sv_to_snapshot_files(statevector, startDate_1, startDate_2)
+
+iunit = open_file('data.cal.DART',form='formatted',action='rewind')
+write(iunit, nml=CAL_NML)
+close(iunit)
+
+!----------------------------------------------------------------------
+! convert the adv_to_time to the appropriate number of seconds.
+!----------------------------------------------------------------------
+
+model_timestep = get_model_time_step()
+offset         = adv_to_time - model_time
+
+call set_model_end_time(offset)
+call write_data_namelistfile()
+
+!----------------------------------------------------------------------
+! Log what we think we're doing, and exit.
+!----------------------------------------------------------------------
+
+call get_time(offset, secs, days)
+
+call print_date( model_time,'dart_to_pop:dart model date')
+call print_date(adv_to_time,'dart_to_pop:advance_to date')
+call print_time( model_time,'dart_to_pop:dart model time')
+call print_time(adv_to_time,'dart_to_pop:advance_to time')
+call print_time(     offset,'dart_to_pop:a distance of')
+write(    *      ,'(''dart_to_pop:PARM03   endTime '',i,'' seconds'')') &
+                   (secs + days*SECPERDAY)
+
+call print_date( model_time,'dart_to_pop:dart model date',logfileunit)
+call print_date(adv_to_time,'dart_to_pop:advance_to date',logfileunit)
+call print_time( model_time,'dart_to_pop:dart model time',logfileunit)
+call print_time(adv_to_time,'dart_to_pop:advance_to time',logfileunit)
+call print_time(     offset,'dart_to_pop:  a distance of',logfileunit)
+write(logfileunit,'(''dart_to_pop:PARM03   endTime '',i,'' seconds'')') &
+                   (secs + days*SECPERDAY)
+
+call finalize_utilities()
+
+end program dart_to_pop


Property changes on: DART/trunk/models/POP/dart_to_pop.f90
___________________________________________________________________
Name: svn:keywords
   + Date Revision Author HeadURL Id
Name: svn:eol-style
   + native

Added: DART/trunk/models/POP/model_mod.f90
===================================================================
--- DART/trunk/models/POP/model_mod.f90	                        (rev 0)
+++ DART/trunk/models/POP/model_mod.f90	2009-06-22 20:21:54 UTC (rev 3940)
@@ -0,0 +1,2963 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2009, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+module model_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+! This is the interface between the POP ocean model and DART.
+
+! Modules that are absolutely required for use are listed
+use        types_mod, only : r4, r8, SECPERDAY
+use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time, &
+                             set_calendar_type, GREGORIAN, print_time, print_date, &
+                             operator(*),  operator(+), operator(-), &
+                             operator(>),  operator(<), operator(/), &
+                             operator(/=), operator(<=)
+use     location_mod, only : location_type,      get_close_maxdist_init, &
+                             get_close_obs_init, get_close_obs, set_location, &
+                             VERTISHEIGHT, get_location, vert_is_height, &
+                             vert_is_level, vert_is_surface
+use    utilities_mod, only : register_module, error_handler, E_ERR, E_WARN, E_MSG, &
+                             logfileunit, get_unit, nc_check, do_output, to_upper, &
+                             find_namelist_in_file, check_namelist_read, &
+                             open_file, file_exist, find_textfile_dims, file_to_text
+use     obs_kind_mod, only : KIND_TEMPERATURE, KIND_SALINITY, KIND_U_CURRENT_COMPONENT, &
+                             KIND_V_CURRENT_COMPONENT, KIND_SEA_SURFACE_HEIGHT
+use mpi_utilities_mod, only: my_task_id
+use random_seq_mod,   only : random_seq_type, init_random_seq, random_gaussian
+
+
+implicit none
+private
+
+! these routines must be public and you cannot change
+! the arguments - they will be called *from* the DART code.
+public :: get_model_size,         &
+          adv_1step,              &
+          get_state_meta_data,    &
+          model_interpolate,      &
+          get_model_time_step,    &
+          end_model,              &
+          static_init_model,      &
+          init_time,              &
+          init_conditions,        &
+          nc_write_model_atts,    &
+          nc_write_model_vars,    &
+          pert_model_state,       &
+          get_close_maxdist_init, &
+          get_close_obs_init,     &
+          get_close_obs,          &
+          ens_mean_for_model
+
+! generally useful routines for various support purposes.
+! the interfaces here can be changed as appropriate.
+public :: prog_var_to_vector, vector_to_prog_var, &
+          POP_meta_type, read_meta, write_meta, &
+          read_snapshot, write_snapshot, get_gridsize, &
+          write_data_namelistfile, set_model_end_time, &
+          snapshot_files_to_sv, sv_to_snapshot_files, &
+          timestep_to_DARTtime, DARTtime_to_POPtime, &
+          DARTtime_to_timestepindex 
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+   source   = "$URL$", &
+   revision = "$Revision$", &
+   revdate  = "$Date$"
+
+character(len=129) :: msgstring
+logical, save :: module_initialized = .false.
+
+! Storage for a random sequence for perturbing a single initial state
+type(random_seq_type) :: random_seq
+
+
+!! FIXME: This is horrid ... 'reclen' is compiler-dependent.
+!! IBM XLF  -- item_size_direct_access == 4,8
+!! IFORT    -- needs compile switch "-assume byterecl"
+integer, parameter :: item_size_direct_access = 4
+
+!------------------------------------------------------------------
+!
+! POP namelist section:  we want to share the 'data' namelist file
+! with the model, so we must declare all possible namelist entries to
+! avoid getting an error from a valid namelist file.  Most of these
+! values are unused in this model_mod code; only a few are needed and
+! those are indicated in comments below.
+!
+!------------------------------------------------------------------
+! The time manager namelist variables
+! some types/etc come from   <POPsource>/pkg/cal/cal.h
+! some useful insight from   cal_set.F, cal_convdate.F
+!
+! startDate_1 (integer) yyyymmdd   "start date of the integration"
+! startDate_2 (integer) hhmmss
+!------------------------------------------------------------------
+
+character(len=9) :: TheCalendar = 'gregorian'
+! integration start date follows: yyyymmddhhmmss
+integer          :: startDate_1 = 19530101
+integer          :: startDate_2 =          60000
+logical          :: calendarDumps = .false.
+
+NAMELIST /CAL_NML/ TheCalendar, startDate_1, startDate_2, calendarDumps
+
+! FIXME: these namelists should probably be in a separate file, and only
+! the actual values needed should be made public, so this isn't so messy.
+
+! must match the value in EEPARAMS.h
+integer, parameter :: MAX_LEN_FNAM = 512
+
+! these come from SIZE.h and can vary.  
+! should they be namelist controlled?
+integer, parameter :: max_nx = 1024
+integer, parameter :: max_ny = 1024
+integer, parameter :: max_nz = 512
+integer, parameter :: max_nr = 512
+
+!-- FIXME: i have not been able to find all of these in the
+!-- original source code.  the ones we're using have been
+!-- checked that they are the right type.  the others might
+!-- cause problems since i don't have a namelist file i can
+!-- examine that has the full set of values specified.
+!-- 
+!-- must match lists declared in ini_parms.f
+
+!--   Time stepping parameters variable declarations
+real(r8) :: pickupSuff, &
+      deltaT, deltaTClock, deltaTmom, &
+      deltaTtracer, dTtracerLev(max_nr), deltaTfreesurf, &
+      abEps, alph_AB, beta_AB, &
+      tauCD, rCD, &
+      baseTime, startTime, endTime, chkPtFreq, &
+      dumpFreq, dumpInitAndLast, adjDumpFreq, taveFreq, tave_lastIter, &
+      diagFreq, monitorFreq, adjMonitorFreq, pChkPtFreq, cAdjFreq, &
+      outputTypesInclusive, &
+      tauThetaClimRelax, tauSaltClimRelax, latBandClimRelax, &
+      tauThetaClimRelax3Dim, tauSaltClimRelax3Dim, tauTr1ClimRelax, &
+      periodicExternalForcing, externForcingPeriod, externForcingCycle
+integer :: nIter0, nTimeSteps, nEndIter, momForcingOutAB, tracForcingOutAB
+logical :: forcing_In_AB, &
+      momDissip_In_AB, doAB_onGtGs, &
+      startFromPickupAB2
+
+!--   Gridding parameters variable declarations 
+logical :: usingCartesianGrid, usingCylindricalGrid, &
+           usingSphericalPolarGrid, usingCurvilinearGrid, &
+           deepAtmosphere
+real(r8) :: dxSpacing, dySpacing, delX(max_nx), delY(max_ny), &
+            phiMin, thetaMin, rSphere, &
+            Ro_SeaLevel, delZ(max_nz), delP, delR(max_nr), delRc(max_nr+1), &
+            rkFac, groundAtK1
+character(len=MAX_LEN_FNAM) :: delXFile, delYFile, &
+                      delRFile, delRcFile, &
+                      horizGridFile
+
+!--   Input files variable declarations 
+character(len=MAX_LEN_FNAM) :: &
+      bathyFile, topoFile, shelfIceFile, &
+      hydrogThetaFile, hydrogSaltFile, diffKrFile, &
+      zonalWindFile, meridWindFile, &
+      thetaClimFile, saltClimFile, &
+      surfQfile, surfQnetFile, surfQswFile, EmPmRfile, saltFluxFile, &
+      lambdaThetaFile, lambdaSaltFile, &
+      uVelInitFile, vVelInitFile, pSurfInitFile, &
+      dQdTFile, ploadFile,tCylIn,tCylOut, &
+      eddyTauxFile, eddyTauyFile, &
+      mdsioLocalDir, &
+      the_run_name
+
+!--   Time stepping parameters namelist
+NAMELIST /PARM03/ &
+      nIter0, nTimeSteps, nEndIter, pickupSuff, &
+      deltaT, deltaTClock, deltaTmom, &
+      deltaTtracer, dTtracerLev, deltaTfreesurf, &
+      forcing_In_AB, momForcingOutAB, tracForcingOutAB, &
+      momDissip_In_AB, doAB_onGtGs, &
+      abEps, alph_AB, beta_AB, startFromPickupAB2, &
+      tauCD, rCD, &
+      baseTime, startTime, endTime, chkPtFreq, &
+      dumpFreq, dumpInitAndLast, adjDumpFreq, taveFreq, tave_lastIter, &
+      diagFreq, monitorFreq, adjMonitorFreq, pChkPtFreq, cAdjFreq, &
+      outputTypesInclusive, &
+      tauThetaClimRelax, tauSaltClimRelax, latBandClimRelax, &
+      tauThetaClimRelax3Dim, tauSaltClimRelax3Dim, tauTr1ClimRelax, &
+      periodicExternalForcing, externForcingPeriod, externForcingCycle, &
+      calendarDumps
+
+!--   Gridding parameters namelist
+NAMELIST /PARM04/ &
+      usingCartesianGrid, usingCylindricalGrid, &
+      dxSpacing, dySpacing, delX, delY, delXFile, delYFile, &
+      usingSphericalPolarGrid, phiMin, thetaMin, rSphere, &
+      usingCurvilinearGrid, horizGridFile, deepAtmosphere, &
+      Ro_SeaLevel, delZ, delP, delR, delRc, delRFile, delRcFile, &
+      rkFac, groundAtK1
+
+!--   Input files namelist
+NAMELIST /PARM05/ &
+      bathyFile, topoFile, shelfIceFile, &
+      hydrogThetaFile, hydrogSaltFile, diffKrFile, &
+      zonalWindFile, meridWindFile, &
+      thetaClimFile, saltClimFile, &
+      surfQfile, surfQnetFile, surfQswFile, EmPmRfile, saltFluxFile, &
+      lambdaThetaFile, lambdaSaltFile, &
+      uVelInitFile, vVelInitFile, pSurfInitFile, &
+      dQdTFile, ploadFile,tCylIn,tCylOut, &
+      eddyTauxFile, eddyTauyFile, &
+      mdsioLocalDir, &
+      the_run_name
+
+!------------------------------------------------------------------
+!
+! The DART state vector (control vector) will consist of:  S, T, U, V, Eta
+! (Salinity, Temperature, U velocity, V velocity, Sea Surface Height).
+! S, T are 3D arrays, located at cell centers.  U is staggered in X
+! and V is staggered in Y (meaning the points are located on the cell
+! faces) but the grids are offset by half a cell, so there are actually
+! the same number of points in each grid. 
+! Eta is a 2D field (X,Y only).  The Z direction is downward.
+!
+!------------------------------------------------------------------
+
+integer, parameter :: n3dfields = 4
+integer, parameter :: n2dfields = 1
+integer, parameter :: nfields   = n3dfields + n2dfields
+
+integer, parameter :: S_index   = 1
+integer, parameter :: T_index   = 2
+integer, parameter :: U_index   = 3
+integer, parameter :: V_index   = 4
+integer, parameter :: Eta_index = 5
+
+! (the absoft compiler likes them to all be the same length during declaration)
+! we trim the blanks off before use anyway, so ...
+character(len=128) :: progvarnames(nfields) = (/'S  ','T  ','U  ','V  ','Eta'/)
+
+integer :: start_index(nfields)
+
+! Grid parameters - the values will be read from a
+! standard POP namelist and filled in here.
+
+integer :: Nx=-1, Ny=-1, Nz=-1    ! grid counts for each field
+
+! locations of cell centers (C) and edges (G) for each axis.
+real(r8), allocatable :: XC(:), XG(:), YC(:), YG(:), ZC(:), ZG(:)
+
+! location information - these grids can either be regularly
+! spaced or the spacing along each axis can vary.
+
+!real(r8) :: lat_origin, lon_origin
+!logical  :: regular_lat, regular_lon, regular_depth
+!real(r8) :: delta_lat, delta_lon, delta_depth
+!real(r8), allocatable :: lat_grid(:), lon_grid(:), depth_grid(:)
+
+real(r8)        :: ocean_dynamics_timestep = 900.0_r4
+integer         :: timestepcount = 0
+type(time_type) :: model_time, model_timestep
+
+integer :: model_size    ! the state vector length
+
+! Skeleton of a model_nml that would be in input.nml
+! This is where dart-related model parms could be set.
+logical  :: output_state_vector = .true.
+integer  :: assimilation_period_days = 7
+integer  :: assimilation_period_seconds = 0
+real(r8) :: model_perturbation_amplitude = 0.2
+
+namelist /model_nml/ output_state_vector,         &
+                     assimilation_period_days,    &
+                     assimilation_period_seconds, &
+                     model_perturbation_amplitude
+
+! /pkg/mdsio/mdsio_write_meta.F writes the .meta files 
+type POP_meta_type
+!  private
+   integer :: nDims
+   integer :: dimList(3)
+   character(len=32) :: dataprec
+   integer :: reclen
+   integer :: nrecords
+   integer :: timeStepNumber    ! optional
+end type POP_meta_type
+
+INTERFACE read_snapshot
+      MODULE PROCEDURE read_2d_snapshot
+      MODULE PROCEDURE read_3d_snapshot
+END INTERFACE
+
+INTERFACE write_snapshot
+      MODULE PROCEDURE write_2d_snapshot
+      MODULE PROCEDURE write_3d_snapshot
+END INTERFACE
+
+INTERFACE vector_to_prog_var
+      MODULE PROCEDURE vector_to_2d_prog_var
+      MODULE PROCEDURE vector_to_3d_prog_var
+END INTERFACE
+
+
+contains
+
+!==================================================================
+
+
+
+subroutine static_init_model()
+!------------------------------------------------------------------
+!
+! Called to do one time initialization of the model. In this case,
+! it reads in the grid information and then the model data.
+
+integer :: i, iunit, io
+integer :: ss, dd
+
+! The Plan:
+!
+!   read the standard POP namelist file 'data.cal' for the calendar info
+!
+!   read the standard POP namelist file 'data' for the
+!   time stepping info and the grid info.
+!
+!   open the grid data files to get the actual grid coordinates
+!
+!   Compute the model size.
+!
+!   set the index numbers where the field types change
+!
+!   set the grid location info
+!
+
+! Print module information to log file and stdout.
+call register_module(source, revision, revdate)
+
+! Since this routine calls other routines that could call this routine
+! we'll say we've been initialized pretty dang early.
+module_initialized = .true.
+
+
+! Read the DART namelist for this model
+call find_namelist_in_file("input.nml", "model_nml", iunit)
+read(iunit, nml = model_nml, iostat = io)
+call check_namelist_read(iunit, io, "model_nml")
+
+! Record the namelist values used for the run
+call error_handler(E_MSG,'static_init_model','model_nml values are',' ',' ',' ')
+if (do_output()) write(logfileunit, nml=model_nml)
+if (do_output()) write(     *     , nml=model_nml)
+
+! POP calendar information
+call find_namelist_in_file("data.cal", "CAL_NML", iunit)
+read(iunit, nml = CAL_NML, iostat = io)
+call check_namelist_read(iunit, io, "CAL_NML")
+
+if (index(TheCalendar,'g') > 0 ) then
+   call set_calendar_type(GREGORIAN)
+elseif (index(TheCalendar,'G') > 0 )then
+   call set_calendar_type(GREGORIAN)
+else
+   write(msgstring,*)"namelist data.cal indicates a ",trim(TheCalendar)," calendar."
+   call error_handler(E_MSG,"static_init_model", msgstring, source, revision, revdate)
+   write(msgstring,*)"only have support for Gregorian"
+   call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
+endif
+if (do_output()) write(*,*)'model_mod:namelist cal_NML',startDate_1,startDate_2
+if (do_output()) write(*,nml=CAL_NML)
+
+! Time stepping parameters are in PARM03
+call find_namelist_in_file("data", "PARM03", iunit)
+read(iunit, nml = PARM03, iostat = io)
+call check_namelist_read(iunit, io, "PARM03")
+
+if ((deltaTmom   == deltaTtracer) .and. &
+    (deltaTmom   == deltaTClock ) .and. &
+    (deltaTClock == deltaTtracer)) then
+   ocean_dynamics_timestep = deltaTmom                    ! need a time_type version
+else
+   write(msgstring,*)"namelist PARM03 has deltaTmom /= deltaTtracer /= deltaTClock"
+   call error_handler(E_MSG,"static_init_model", msgstring, source, revision, revdate)
+   write(msgstring,*)"values were ",deltaTmom, deltaTtracer, deltaTClock
+   call error_handler(E_MSG,"static_init_model", msgstring, source, revision, revdate)
+   write(msgstring,*)"At present, DART only supports equal values."
+   call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
+endif
+
+! Define the assimilation period as the model_timestep
+! Ensure model_timestep is multiple of ocean_dynamics_timestep
+
+model_time     = timestep_to_DARTtime(timestepcount)
+model_timestep = set_model_time_step(assimilation_period_seconds, &
+                                     assimilation_period_days,    &
+                                     ocean_dynamics_timestep)
+
+call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
+
+write(msgstring,*)"assimilation period is ",dd," days ",ss," seconds"
+call error_handler(E_MSG,'static_init_model',msgstring,source,revision,revdate)
+if (do_output()) write(logfileunit,*)msgstring
+
+! Grid-related variables are in PARM04
+delX(:) = 0.0_r4
+delY(:) = 0.0_r4
+delZ(:) = 0.0_r4
+call find_namelist_in_file("data", "PARM04", iunit)
+read(iunit, nml = PARM04, iostat = io)
+call check_namelist_read(iunit, io, "PARM04")
+
+! Input datasets are in PARM05
+call find_namelist_in_file("data", "PARM05", iunit)
+read(iunit, nml = PARM05, iostat = io)
+call check_namelist_read(iunit, io, "PARM05")
+
+! The only way I know to compute the number of
+! levels/lats/lons is to set the default value of delZ to 0.0
+! before reading the namelist.  now loop until you get back
+! to zero and that is the end of the list.
+! Not a very satisfying/robust solution ...
+
+Nx = -1
+do i=1, size(delX)
+ if (delX(i) == 0.0_r4) then
+    Nx = i-1
+    exit
+ endif
+enddo
+if (Nx == -1) then
+   write(msgstring,*)'could not figure out number of longitudes from delX in namelist'
+   call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
+endif
+
+Ny = -1
+do i=1, size(delY)
+ if (delY(i) == 0.0_r4) then
+    Ny = i-1
+    exit
+ endif
+enddo
+if (Ny == -1) then
+   write(msgstring,*)'could not figure out number of latitudes from delY in namelist'
+   call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
+endif
+
+Nz = -1
+do i=1, size(delZ)
+ if (delZ(i) == 0.0_r4) then
+    Nz = i-1
+    exit
+ endif
+enddo
+if (Nz == -1) then
+   write(msgstring,*)'could not figure out number of depth levels from delZ in namelist'
+   call error_handler(E_ERR,"static_init_model", msgstring, source, revision, revdate)
+endif
+
+! We know enough to allocate grid variables. 
+
+allocate(XC(Nx), YC(Ny), ZC(Nz))
+allocate(XG(Nx), YG(Ny), ZG(Nz))
+
+! XG (the grid edges) and XC (the grid centroids) must be computed.
+
+XG(1) = thetaMin
+XC(1) = thetaMin + 0.5_r8 * delX(1)
+do i=2, Nx
+ XG(i) = XG(i-1) + delX(i-1)
+ XC(i) = XC(i-1) + 0.5_r8 * delX(i-1) + 0.5_r8 * delX(i) 
+enddo
+
+! YG (the grid edges) and YC (the grid centroids) must be computed.
+
+YG(1) = phiMin
+YC(1) = phiMin + 0.5_r8 * delY(1)
+do i=2, Ny
+ YG(i) = YG(i-1) + delY(i-1)
+ YC(i) = YC(i-1) + 0.5_r8 * delY(i-1) + 0.5_r8 * delY(i) 
+enddo
+
+! the namelist contains a list of thicknesses of each depth level (delZ)
+! ZG (the grid edges) and ZC (the grid centroids) must be computed.
+
+ZG(1) = 0.0_r8
+ZC(1) = -0.5_r8 * delZ(1)
+do i=2, Nz
+ ZG(i) = ZG(i-1) - delZ(i-1)
+ ZC(i) = ZC(i-1) - 0.5_r8 * delZ(i-1) - 0.5_r8 * delZ(i) 
+enddo
+
+! record where in the state vector the data type changes
+! from one type to another, by computing the starting
+! index for each block of data.
+start_index(S_index)   = 1
+start_index(T_index)   = start_index(S_index) + (Nx * Ny * Nz)
+start_index(U_index)   = start_index(T_index) + (Nx * Ny * Nz)
+start_index(V_index)   = start_index(U_index) + (Nx * Ny * Nz)
+start_index(Eta_index) = start_index(V_index) + (Nx * Ny * Nz)
+
+! in spite of the staggering, all grids are the same size
+! and offset by half a grid cell.  4 are 3D and 1 is 2D.
+!  e.g. S,T,U,V = 256 x 225 x 70
+!  e.g. Eta = 256 x 225
+
+if (do_output()) write(logfileunit, *) 'Using grid size : '
+if (do_output()) write(logfileunit, *) '  Nx, Ny, Nz = ', Nx, Ny, Nz
+if (do_output()) write(     *     , *) 'Using grid size : '
+if (do_output()) write(     *     , *) '  Nx, Ny, Nz = ', Nx, Ny, Nz
+!print *, ' 3d field size: ', n3dfields * (Nx * Ny * Nz)
+!print *, ' 2d field size: ', n2dfields * (Nx * Ny)
+model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
+if (do_output()) write(*,*) 'model_size = ', model_size
+
+end subroutine static_init_model
+
+
+
+
+subroutine init_conditions(x)
+!------------------------------------------------------------------
+!
+! Returns a model state vector, x, that is some sort of appropriate
+! initial condition for starting up a long integration of the model.
+! At present, this is only used if the namelist parameter 
+! start_from_restart is set to .false. in the program perfect_model_obs.
+! If this option is not to be used in perfect_model_obs, or if no 
+! synthetic data experiments using perfect_model_obs are planned, 
+! this can be a NULL INTERFACE.
+
+real(r8), intent(out) :: x(:)
+
+if ( .not. module_initialized ) call static_init_model
+ 
+x = 0.0_r8
+
+end subroutine init_conditions
+
+
+
+subroutine adv_1step(x, time)
+!------------------------------------------------------------------
+!
+! Does a single timestep advance of the model. The input value of
+! the vector x is the starting condition and x is updated to reflect
+! the changed state after a timestep. The time argument is intent
+! in and is used for models that need to know the date/time to 
+! compute a timestep, for instance for radiation computations.
+! This interface is only called IF the namelist parameter
+! async is set to 0 in perfect_model_obs or filter -OR- if the 
+! program integrate_model is to be used to advance the model
+! state as a separate executable. If none of these options
+! are used (the model will only be advanced as a separate 
+! model-specific executable), this can be a NULL INTERFACE.
+
+real(r8),        intent(inout) :: x(:)
+type(time_type), intent(in)    :: time
+
+if ( .not. module_initialized ) call static_init_model
+
+if (do_output()) then
+   call print_time(time,'NULL interface adv_1step (no advance) DART time is')
+   call print_time(time,'NULL interface adv_1step (no advance) DART time is',logfileunit)
+endif
+
+end subroutine adv_1step
+
+
+
+function get_model_size()
+!------------------------------------------------------------------
+!
+! Returns the size of the model as an integer. Required for all
+! applications.
+
+integer :: get_model_size
+
+if ( .not. module_initialized ) call static_init_model
+
+get_model_size = model_size
+
+end function get_model_size
+
+
+
+subroutine init_time(time)
+!------------------------------------------------------------------
+!
+! Companion interface to init_conditions. Returns a time that is somehow 
+! appropriate for starting up a long integration of the model.
+! At present, this is only used if the namelist parameter 
+! start_from_restart is set to .false. in the program perfect_model_obs.
+! If this option is not to be used in perfect_model_obs, or if no 
+! synthetic data experiments using perfect_model_obs are planned, 
+! this can be a NULL INTERFACE.
+
+type(time_type), intent(out) :: time
+
+if ( .not. module_initialized ) call static_init_model
+
+! for now, just set to 0
+time = set_time(0,0)
+
+end subroutine init_time
+
+
+
+subroutine model_interpolate(x, location, obs_type, interp_val, istatus)
+!=======================================================================
+!
+
+real(r8),            intent(in) :: x(:)
+type(location_type), intent(in) :: location
+integer,             intent(in) :: obs_type
+real(r8),           intent(out) :: interp_val
+integer,            intent(out) :: istatus
+
+! Model interpolate will interpolate any state variable (S, T, U, V, Eta) to
+! the given location given a state vector. The type of the variable being
+! interpolated is obs_type since normally this is used to find the expected
+! value of an observation at some location. The interpolated value is 
+! returned in interp_val and istatus is 0 for success.
+
+! Local storage
+real(r8)       :: loc_array(3), llon, llat, lheight
+integer        :: base_offset, offset, ind
+integer        :: hgt_bot, hgt_top
+real(r8)       :: hgt_fract
+real(r8)       :: top_val, bot_val
+integer        :: hstatus
+
+if ( .not. module_initialized ) call static_init_model
+

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list