[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