[Dart-dev] [5676] DART/branches/development/models/template: start to flesh out a full example/template directory.
nancy at ucar.edu
nancy at ucar.edu
Mon Apr 9 11:50:48 MDT 2012
Revision: 5676
Author: nancy
Date: 2012-04-09 11:50:48 -0600 (Mon, 09 Apr 2012)
Log Message:
-----------
start to flesh out a full example/template directory. the original
1d-based, simplistic model_mod needs to remain because it is compiled
into generic tools which require a model_mod to satisfy the subroutine
names but is never actually called. the 'full_model_mod.f90' file
is the start of an example full 3d geophysical model interface file,
along with the normal utilities we support.
Added Paths:
-----------
DART/branches/development/models/template/README
DART/branches/development/models/template/dart_to_model.f90
DART/branches/development/models/template/full_model_mod.f90
DART/branches/development/models/template/model_mod_check.f90
DART/branches/development/models/template/model_to_dart.f90
DART/branches/development/models/template/trans_time.f90
DART/branches/development/models/template/work/mkmf_dart_to_model
DART/branches/development/models/template/work/mkmf_model_mod_check
DART/branches/development/models/template/work/mkmf_model_to_dart
DART/branches/development/models/template/work/mkmf_trans_time
DART/branches/development/models/template/work/path_names_dart_to_model
DART/branches/development/models/template/work/path_names_model_mod_check
DART/branches/development/models/template/work/path_names_model_to_dart
DART/branches/development/models/template/work/path_names_trans_time
Removed Paths:
-------------
DART/branches/development/models/template/utils/
-------------- next part --------------
Copied: DART/branches/development/models/template/README (from rev 5672, DART/branches/mpas/models/template/README)
===================================================================
--- DART/branches/development/models/template/README (rev 0)
+++ DART/branches/development/models/template/README 2012-04-09 17:50:48 UTC (rev 5676)
@@ -0,0 +1,135 @@
+
+Hints for porting a new model to DART:
+
+copy this template directory into a DART/models/xxx
+directory for your new model.
+
+if the coordinate system for the model is 1d, you're ok as-is.
+if model coordinates are 3d, edit the work/path_names_* files
+and change location/oned/* to location/threed_sphere/*
+
+if your model is closer to the simpler examples (e.g. lorenz),
+the existing model_mod.f90 is a good place to start.
+if your model is a full 3d geophysical one (e.g. like cam, pop, etc)
+then rename full_model_mod.f90 to model_mod.f90 and start there.
+
+edit all the work/path_names_* files and change models/template/xxx
+to use the name of the directory for your model.
+
+try ./quickbuild.csh and everything should compile at this point.
+
+the required subroutines are these:
+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
+
+in addition, model_mod can contain subroutines that are used
+for other utility programs and we recommend at least the following
+routines be added to model_mod.f90:
+
+public :: model_file_to_dart_vector, & ! converter
+ dart_vector_to_model_file, & ! converter
+ get_gridsize, & ! called by everyone
+ get_model_filename, & ! called by both (set_model_filename?)
+ get_state_time, & ! model_to_sv, static_init_model
+ set_state_time !(?) ! sv_to_model, trans_time
+
+
+edit the model mod and fill in the routines in this order:
+
+1. static_init_model() - make it read in the grid information
+ and the number of variables that will be in the state vector
+ (fill in the progvar derived type). fill in the model_size
+ variable. as part of this work, fill in the get_gridsize()
+ code.
+
+ after number 1 is done, get_model_size() and
+ get_model_time_step() from the template should be ok as-is.
+
+2. model_file_to_dart_vector() - given a model data file, read in
+ the fields and put them into the 1D DART state vector. make
+ sure the order and values match the progvar array.
+
+3. dart_vector_to_model_file() - do the reverse of the previous step.
+
+4. get_state_meta_data() - given an index number into the state vector
+ return the location and type. the code which loops over the
+ progvar should already do the type, but code to compute what
+ lon, lat, and vertical (for a 3d model) or x location (1d)
+ corresponds to this item must be written.
+
+5. model_interpolate() - given a location (lon/lat/vert in 3d, x in 1d)
+ and a state KIND_xxx kind, return the interpolated value the field
+ has at that location. this is probably one of the routines that
+ will take the most code to write.
+
+6. nc_write_model_atts(), nc_write_model_vars() - when filter runs
+ it calls these routines to output model data into a netcdf diagnostic
+ file which is unrelated to the model data files. it is possible to
+ have the ensemble data just be dumped as a single 1D vector but
+ that makes the files less useful. generally it's most useful to
+ put the grid info and dump each variable from the state vector
+ into a separate netcdf variable. the diagnostic files contain the
+ ensemble mean, the ensemble stddev, the inflation information, and
+ then optionally some or all of the individual ensemble members.
+
+for now, ignore these routines:
+ get_close_maxdist_init()
+ get_close_obs_init()
+ get_close_obs()
+ ens_mean_for_model()
+ end_model()
+
+if you have data in a dart initial condition/restart file, then you
+can ignore these routines:
+ init_time()
+ init_conditions()
+otherwise, have them return an initial time and an initial default
+ensemble state.
+
+if your model is NOT subroutine callable, you can ignore this routine:
+ adv_1step()
+otherwise have it call the interface to your model and add the files
+necessary to build your model to all the work/path_names_* files.
+add the model source to a src/ directory.
+
+if you want to let filter add gaussian noise to a single state vector
+to generate an ensemble, you can ignore this routine
+ pert_model_state()
+otherwise fill in code that does whatever perturbation makes sense
+to have an initial ensemble of states. in some cases that means
+adding a different range of values to each different field in the
+state vector.
+
+at this point you should have enough code to test and run simple
+experiments. the 'model_mod_check' utility program can be used
+during this process to check your implementation.
+
+the general flow is:
+
+ ./model_to_dart - read model data and convert it into a dart state vector file
+ ./create_obs_sequence - make a file with a single observation in it
+ ./perfect_model_obs - should interpolate a value for the obs
+ ./dart_to_model - convert the dart vector back into a model data file
+
+ generate an ensemble of states, or set 'start_from_restart' to .false.
+ run ./filter with the single observation
+ look at the Prior_Diag.nc and Posterior_Diag.nc files
+ diff them with ncdiff: ncdiff Posterior_Diag.nc Prior_Diag.nc Innov.nc
+ plot it, with ncview if possible: ncview Innov.nc
+ the difference between the two is the impact of that single observation
+ see if it's at the right location and if the differences seem reasonable
+
Copied: DART/branches/development/models/template/dart_to_model.f90 (from rev 5672, DART/branches/mpas/models/template/dart_to_model.f90)
===================================================================
--- DART/branches/development/models/template/dart_to_model.f90 (rev 0)
+++ DART/branches/development/models/template/dart_to_model.f90 2012-04-09 17:50:48 UTC (rev 5676)
@@ -0,0 +1,191 @@
+! DART software - Copyright 2004 - 2011 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+program dart_to_model
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!----------------------------------------------------------------------
+! purpose: interface between DART and the model model
+!
+! method: Read DART state vector and overwrite values in a model restart file.
+! If the DART state vector has an 'advance_to_time' present, a
+! file called model_in.DART is created with a time_manager_nml namelist
+! appropriate to advance model to the requested time.
+!
+! The dart_to_model_nml namelist setting for advance_time_present
+! determines whether or not the input file has an 'advance_to_time'.
+! Typically, only temporary files like 'assim_model_state_ic' have
+! an 'advance_to_time'.
+!
+! author: Tim Hoar 25 Jun 09, revised 12 July 2010
+!----------------------------------------------------------------------
+
+use types_mod, only : r8
+use utilities_mod, only : initialize_utilities, timestamp, &
+ find_namelist_in_file, check_namelist_read, &
+ logfileunit, open_file, close_file
+use assim_model_mod, only : open_restart_read, aread_state_restart, close_restart
+use time_manager_mod, only : time_type, print_time, print_date, operator(-), &
+ get_time, get_date
+use model_mod, only : static_init_model, sv_to_restart_file, &
+ get_model_size, get_base_time, get_model_restart_dirname
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+!------------------------------------------------------------------
+! The namelist variables
+!------------------------------------------------------------------
+
+character (len = 128) :: dart_to_model_input_file = 'dart.ic'
+logical :: advance_time_present = .false.
+character(len=256) :: model_restart_dirname = 'model_restartdir'
+
+namelist /dart_to_model_nml/ dart_to_model_input_file, &
+ advance_time_present, &
+ model_restart_dirname
+
+!----------------------------------------------------------------------
+
+integer :: iunit, io, x_size, diff1, diff2
+type(time_type) :: model_time, adv_to_time, base_time
+real(r8), allocatable :: statevector(:)
+logical :: verbose = .FALSE.
+
+!----------------------------------------------------------------------
+
+call initialize_utilities(progname='dart_to_model', output_flag=verbose)
+
+!----------------------------------------------------------------------
+! Call model_mod:static_init_model() which reads the model namelists
+! to set grid sizes, etc.
+!----------------------------------------------------------------------
+
+call static_init_model()
+
+x_size = get_model_size()
+allocate(statevector(x_size))
+
+! Read the namelist to get the input dirname.
+
+call find_namelist_in_file("input.nml", "dart_to_model_nml", iunit)
+read(iunit, nml = dart_to_model_nml, iostat = io)
+call check_namelist_read(iunit, io, "dart_to_model_nml")
+
+write(*,*)
+write(*,*) 'dart_to_model: converting DART file ', "'"//trim(dart_to_model_input_file)//"'"
+write(*,*) 'to model restart files in directory ', "'"//trim(model_restart_dirname)//"'"
+
+!----------------------------------------------------------------------
+! Reads the valid time, the state, and the target time.
+!----------------------------------------------------------------------
+
+iunit = open_restart_read(dart_to_model_input_file)
+
+if ( advance_time_present ) then
+ call aread_state_restart(model_time, statevector, iunit, adv_to_time)
+else
+ call aread_state_restart(model_time, statevector, iunit)
+endif
+call close_restart(iunit)
+
+print *, 'read state vector'
+!----------------------------------------------------------------------
+! update the current model state vector
+! Convey the amount of time to integrate the model ...
+! time_manager_nml: stop_option, stop_count increments
+!----------------------------------------------------------------------
+
+print *, 'calling sv to restart file'
+call sv_to_restart_file(statevector, model_restart_dirname, model_time)
+
+if ( advance_time_present ) then
+ call write_model_time_control(model_time, adv_to_time)
+endif
+
+!----------------------------------------------------------------------
+! Log what we think we're doing, and exit.
+!----------------------------------------------------------------------
+
+call print_date( model_time,'dart_to_model:model model date')
+call print_time( model_time,'dart_to_model:DART model time')
+call print_date( model_time,'dart_to_model:model model date',logfileunit)
+call print_time( model_time,'dart_to_model:DART model time',logfileunit)
+
+if ( advance_time_present ) then
+call print_time(adv_to_time,'dart_to_model:advance_to time')
+call print_date(adv_to_time,'dart_to_model:advance_to date')
+call print_time(adv_to_time,'dart_to_model:advance_to time',logfileunit)
+call print_date(adv_to_time,'dart_to_model:advance_to date',logfileunit)
+endif
+
+! When called with 'end', timestamp will call finalize_utilities()
+call timestamp(string1=source, pos='end')
+
+!======================================================================
+contains
+!======================================================================
+
+subroutine write_model_time_control(model_time, adv_to_time)
+! The idea is to write a text file with the following structure:
+!
+!#TIMESTART
+!2003 year
+!06 month
+!21 day
+!00 hour
+!00 minute
+!00 second
+!
+!#TIMEEND
+!2003 year
+!07 month
+!21 day
+!00 hour
+!00 minute
+!00 second
+!
+
+type(time_type), intent(in) :: model_time, adv_to_time
+integer :: iyear,imonth,iday,ihour,imin,isec
+
+iunit = open_file('DART_model_time_control.txt', action='write')
+write(iunit,*)
+
+call get_date(model_time,iyear,imonth,iday,ihour,imin,isec)
+write(iunit,'(''#TIMESTART'')')
+write(iunit,'(i4.4,10x,''year'' )')iyear
+write(iunit,'(i2.2,12x,''month'' )')imonth
+write(iunit,'(i2.2,12x,''day'' )')iday
+write(iunit,'(i2.2,12x,''hour'' )')ihour
+write(iunit,'(i2.2,12x,''minute'')')imin
+write(iunit,'(i2.2,12x,''second'')')isec
+write(iunit,*)
+
+call get_date(adv_to_time,iyear,imonth,iday,ihour,imin,isec)
+write(iunit,'(''#TIMEEND'')')
+write(iunit,'(i4.4,10x,''year'' )')iyear
+write(iunit,'(i2.2,12x,''month'' )')imonth
+write(iunit,'(i2.2,12x,''day'' )')iday
+write(iunit,'(i2.2,12x,''hour'' )')ihour
+write(iunit,'(i2.2,12x,''minute'')')imin
+write(iunit,'(i2.2,12x,''second'')')isec
+write(iunit,*)
+
+call close_file(iunit)
+end subroutine write_model_time_control
+
+
+
+end program dart_to_model
Copied: DART/branches/development/models/template/full_model_mod.f90 (from rev 5672, DART/branches/mpas/models/template/full_model_mod.f90)
===================================================================
--- DART/branches/development/models/template/full_model_mod.f90 (rev 0)
+++ DART/branches/development/models/template/full_model_mod.f90 2012-04-09 17:50:48 UTC (rev 5676)
@@ -0,0 +1,3202 @@
+! DART software - Copyright 2004 - 2011 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+module model_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+! This is the interface between the model model and DART.
+
+! Modules that are absolutely required for use are listed
+use types_mod, only : r4, r8, digits12, SECPERDAY, MISSING_R8, &
+ rad2deg, deg2rad, PI
+use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
+ print_time, print_date, set_calendar_type, &
+ operator(*), operator(+), operator(-), &
+ operator(>), operator(<), operator(/), &
+ operator(/=), operator(<=)
+
+use location_mod, only : location_type, get_dist, query_location, &
+ get_close_maxdist_init, get_close_type, &
+ set_location, get_location, horiz_dist_only, &
+ vert_is_undef, VERTISUNDEF, &
+ vert_is_surface, VERTISSURFACE, &
+ vert_is_level, VERTISLEVEL, &
+ vert_is_pressure, VERTISPRESSURE, &
+ vert_is_height, VERTISHEIGHT, &
+ get_close_obs_init, loc_get_close_obs => get_close_obs
+
+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, close_file
+
+use obs_kind_mod, only : paramname_length, &
+ get_raw_obs_kind_index, &
+ get_raw_obs_kind_name
+
+use mpi_utilities_mod, only: my_task_id
+
+use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
+
+use dart_model_mod, only: get_model_nLons, get_model_nLats, get_model_nAlts, &
+ get_nSpecies, get_nSpeciesTotal, get_nIons, &
+ get_nSpeciesAll, decode_model_indices
+
+use typesizes
+use netcdf
+
+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, &
+ static_init_model, &
+ end_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 :: get_gridsize, &
+ restart_file_to_statevector, &
+ statevector_to_restart_file, &
+ get_model_restart_dirname, &
+ get_base_time, &
+ get_state_time
+
+! version controlled file description for error handling, do not edit
+
+character(len=128), parameter :: &
+ source = '$URL$', &
+ revision = '$Revision$', &
+ revdate = '$Date$'
+
+character(len=256) :: string1, string2
+logical, save :: module_initialized = .false.
+
+! Storage for a random sequence for perturbing a single initial state
+
+type(random_seq_type) :: random_seq
+
+! things which can/should be in the model_nml
+
+integer :: assimilation_period_days = 0
+integer :: assimilation_period_seconds = 60
+real(r8) :: model_perturbation_amplitude = 0.2
+logical :: output_state_vector = .true.
+integer :: debug = 0 ! turn up for more and more debug messages
+character(len=32) :: calendar = 'Gregorian'
+character(len=256) :: model_restart_dirname = 'model_restartdir'
+
+namelist /model_nml/ &
+ model_restart_dirname, &
+ output_state_vector, &
+ assimilation_period_days, & ! for now, this is the timestep
+ assimilation_period_seconds, &
+ model_perturbation_amplitude,&
+ calendar, &
+ debug
+
+!------------------------------------------------------------------
+!
+! The DART state vector may consist of things like:
+!
+! U long_name = "X-WIND COMPONENT" float U(TIME, ALT, LAT, XE)
+! V long_name = "Y-WIND COMPONENT" float V(TIME, ALT, YE, LON)
+! W long_name = "Z-WIND COMPONENT" float W(TIME, ZE, LAT, LON)
+! TH long_name = "POTENTIAL TEMPERATURE" float TH(TIME, ALT, LAT, LON)
+! DBZ long_name = "RADAR REFLECTIVITY" float DBZ(TIME, ALT, LAT, LON)
+! WZ long_name = "VERTICAL VORTICITY" float WZ(TIME, ALT, LAT, LON)
+! PI long_name = "PERT. EXNER" float PI(TIME, ALT, LAT, LON)
+! QV long_name = "VAPOR MIXING RATIO" float QV(TIME, ALT, LAT, LON)
+! QC long_name = "CLOUD MIXING RATIO" float QC(TIME, ALT, LAT, LON)
+! QR long_name = "RAIN MIXING RATIO" float QR(TIME, ALT, LAT, LON)
+! QI long_name = "ICE MIXING RATIO" float QI(TIME, ALT, LAT, LON)
+! QS long_name = "SNOW MIXING RATIO" float QS(TIME, ALT, LAT, LON)
+! QH long_name = "GRAUPEL MIXING RATIO" float QH(TIME, ALT, LAT, LON)
+!
+! The variables in the model restart file that are used to create the
+! DART state vector are specified in the input.nml:model_vars_nml namelist.
+!
+!------------------------------------------------------------------
+
+integer, parameter :: max_state_variables = 80
+integer, parameter :: num_state_table_columns = 2
+character(len=NF90_MAX_NAME) :: model_state_variables(max_state_variables * num_state_table_columns ) = ' '
+character(len=NF90_MAX_NAME) :: variable_table(max_state_variables, num_state_table_columns )
+
+namelist /model_vars_nml/ model_state_variables
+
+integer :: nfields
+
+! Everything needed to describe a variable
+
+type progvartype
+ private
+ character(len=NF90_MAX_NAME) :: varname ! crazy species name
+ character(len=NF90_MAX_NAME) :: long_name
+ character(len=NF90_MAX_NAME) :: units
+ character(len=NF90_MAX_NAME) :: storder
+ character(len=NF90_MAX_NAME) :: model_varname ! NDensityS, IDensityS, ...
+ integer :: model_dim ! dimension defining species
+ integer :: model_index ! 'iSpecies' or u,v,w ...
+ integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens ! nlons, nlats, nalts [, nspecies]
+ integer :: posdef
+ integer :: numdims
+ integer :: varsize ! prod(dimlens(1:numdims))
+ integer :: index1 ! location in dart state vector of first occurrence
+ integer :: indexN ! location in dart state vector of last occurrence
+ integer :: dart_kind
+ character(len=paramname_length) :: kind_string
+end type progvartype
+
+type(progvartype), dimension(max_state_variables) :: progvar
+
+! These are statically defined in ModSize.f90 ...
+! nAlts is the one and only number of altitudes ... no block-dependence
+! nLons, nLats are the number of lons/lats PER block
+! the number of blocks comes from UAM.in
+
+integer :: nLons, nLats, nAlts
+
+! "... keep in mind that if the model resolution is 5 deg latitude,
+! the model will actually go from -87.5 to 87.5 latitude
+! (even though you specify -90 to 90 in the UAM.in file),
+! since the latitudes/longitudes are at cell centers,
+! while the edges are at the boundaries." -- Aaron Ridley
+
+integer :: NgridLon=-1, NgridLat=-1, NgridAlt=-1 ! scalar grid counts
+integer :: nBlocksLon=-1, nBlocksLat=-1 ! number of blocks along each dim
+real(r8) :: LatStart=MISSING_R8, LatEnd=MISSING_R8, LonStart=MISSING_R8
+integer :: nSpeciesTotal=-1, nSpecies=-1, nIons=-1, nSpeciesAll=-1
+
+! scalar grid positions
+
+real(r8), allocatable :: LON(:) ! longitude centers
+real(r8), allocatable :: LAT(:) ! latitude centers
+real(r8), allocatable :: ALT(:) ! vertical level centers
+
+integer :: model_size ! the state vector length
+type(time_type) :: model_time ! valid time of the model state
+type(time_type) :: model_timestep ! smallest time to adv model
+real(r8), allocatable :: ens_mean(:) ! may be needed for forward ops
+
+! 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
+
+integer, parameter :: nGhost = 2 ! number of ghost cells on all edges
+
+!------------------------------------------------------------------
+! The model restart manager namelist variables
+!------------------------------------------------------------------
+
+character(len= 64) :: ew_boundary_type, ns_boundary_type
+
+INTERFACE vector_to_prog_var
+ MODULE PROCEDURE vector_to_1d_prog_var
+ MODULE PROCEDURE vector_to_2d_prog_var
+ MODULE PROCEDURE vector_to_3d_prog_var
+ MODULE PROCEDURE vector_to_4d_prog_var
+END INTERFACE
+
+INTERFACE get_base_time
+ MODULE PROCEDURE get_base_time_ncid
+ MODULE PROCEDURE get_base_time_fname
+END INTERFACE
+
+INTERFACE get_index_range
+ MODULE PROCEDURE get_index_range_string
+ MODULE PROCEDURE get_index_range_int
+END INTERFACE
+
+contains
+
+!==================================================================
+! All the REQUIRED interfaces come first - just by convention.
+!==================================================================
+
+
+function get_model_size()
+!------------------------------------------------------------------
+! Done - TJH.
+! 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 adv_1step(x, time)
+!------------------------------------------------------------------
+! Done - TJH.
+! 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
+
+! FIXME: put an error handler call here - we cannot advance the model
+! this way and it would be an error if filter called it.
+
+end subroutine adv_1step
+
+
+
+subroutine get_state_meta_data(index_in, location, var_type)
+!------------------------------------------------------------------
+! Done - JLA.
+! given an index into the state vector, return its location and
+! if given, the var kind. despite the name, var_type is a generic
+! kind, like those in obs_kind/obs_kind_mod.f90, starting with KIND_
+
+integer, intent(in) :: index_in
+type(location_type) :: location
+integer, optional, intent(out) :: var_type
+
+! Local variables
+
+integer :: lat_index, lon_index, alt_index
+integer :: n, nf, myindx, remainder, remainder2
+
+if ( .not. module_initialized ) call static_init_model
+
+! Find out which of the 3D fields index_in is part of
+nf = -1
+
+FindIndex : do n = 1,nfields
+ if( (progvar(n)%index1 <= index_in) .and. (index_in <= progvar(n)%indexN) ) then
+ nf = n
+ myindx = index_in - progvar(n)%index1 + 1
+ exit FindIndex
+ endif
+enddo FindIndex
+
+if( myindx == -1 ) then
+ write(string1,*) 'Problem, cannot find base_offset, index_in is: ', index_in
+ call error_handler(E_ERR,'get_state_meta_data',string1,source,revision,revdate)
+endif
+
+alt_index = 1 + (myindx - 1) / (NgridLon * NgridLat)
+remainder = myindx - (alt_index-1) * NgridLon * NgridLat
+lat_index = 1 + (remainder - 1) / NgridLon
+remainder2 = remainder - (lat_index - 1) * NgridLon
+lon_index = remainder2
+
+location = set_location(LON(lon_index), LAT(lat_index), ALT(alt_index), VERTISHEIGHT)
+
+if (present(var_type)) then
+ var_type = progvar(nf)%dart_kind
+endif
+
+end subroutine get_state_meta_data
+
+
+
+function get_model_time_step()
+!------------------------------------------------------------------
+!
+! Returns the the time step of the model; the smallest increment
+! in time that the model is capable of advancing the state in a given
+! implementation. This interface is required for all applications.
+
+type(time_type) :: get_model_time_step
+
+if ( .not. module_initialized ) call static_init_model
+
+get_model_time_step = model_timestep
+
+end function get_model_time_step
+
+
+
+subroutine static_init_model()
+!------------------------------------------------------------------
+!
+! Called to do one time initialization of the model.
+!
+! All the grid information comes from the initialization of
+! the dart_model_mod module.
+
+! Local variables - all the important ones have module scope
+
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
+character(len=NF90_MAX_NAME) :: varname
+character(len=paramname_length) :: kind_string
+integer :: ncid, VarID, numdims, dimlen, varsize
+integer :: iunit, io, ivar, i, index1, indexN
+integer :: ss, dd
+integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID
+logical :: shapeok
+
+if ( module_initialized ) return ! only need to do this once.
+
+! 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)
+
+! Get the model variables in a restricted scope setting.
+
+nLons = get_model_nLons()
+nLats = get_model_nLats()
+nAlts = get_model_nAlts()
+nSpecies = get_nSpecies()
+nSpeciesTotal = get_nSpeciesTotal()
+nIons = get_nIons()
+nSpeciesAll = get_nSpeciesAll()
+
+if ((debug > 0) .and. do_output() ) then
+ write(*,*)
+ write(*,*)'nLons is ',nLons
+ write(*,*)'nLats is ',nLats
+ write(*,*)'nAlts is ',nAlts
+ write(*,*)'nSpecies is ',nSpecies
+ write(*,*)'nSpeciesTotal is ',nSpeciesTotal
+ write(*,*)'nIons is ',nIons
+ write(*,*)'nSpeciesAll is ',nSpeciesAll
+endif
+
+! Read the model variable list to populate DART state vector
+! Once parsed, the values will be recorded for posterity
+call find_namelist_in_file('input.nml', 'model_vars_nml', iunit)
+read(iunit, nml = model_vars_nml, iostat = io)
+call check_namelist_read(iunit, io, 'model_vars_nml')
+
+!---------------------------------------------------------------
+! Set the time step ... causes model namelists to be read.
+! Ensures model_timestep is multiple of 'dynamics_timestep'
+
+call set_calendar_type( calendar ) ! comes from model_mod_nml
+
+model_timestep = set_model_time_step()
+
+call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
+
+write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds'
+call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate)
+
+!---------------------------------------------------------------
+! 1) get grid dimensions
+! 2) allocate space for the grids
+! 3) read them from the block restart files, could be stretched ...
+
+call get_grid_info(NgridLon, NgridLat, NgridAlt, nBlocksLon, nBlocksLat, &
+ LatStart, LatEnd, LonStart)
+
+if( debug > 0 ) then
+ write(*,*)'grid dims are ',NgridLon,NgridLat,NgridAlt
+endif
+
+allocate( LON( NgridLon ))
+allocate( LAT( NgridLat ))
+allocate( ALT( NgridAlt ))
+
+call get_grid(model_restart_dirname, nBlocksLon, nBlocksLat, &
+ nLons, nLats, nAlts, LON, LAT, ALT )
+
+!---------------------------------------------------------------
+! Compile the list of model variables to use in the creation
+! of the DART state vector. Required to determine model_size.
+!
+! Verify all variables are in the model restart file
+!
+! Compute the offsets into the state vector for the start of each
+! different variable type. Requires reading shapes from the model
+! restart file. As long as TIME is the LAST dimension, we're OK.
+!
+! Record the extent of the data type in the state vector.
+
+call verify_state_variables( model_state_variables, ncid, model_restart_dirname, &
+ nfields, variable_table )
+
+index1 = 1;
+indexN = 0;
+
+do ivar = 1, nfields
+
+ varname = trim(variable_table(ivar,1))
+ kind_string = trim(variable_table(ivar,2))
+ progvar(ivar)%varname = varname
+ progvar(ivar)%kind_string = kind_string
+ progvar(ivar)%dart_kind = get_raw_obs_kind_index( progvar(ivar)%kind_string )
+ progvar(ivar)%dimlens = 0
+
+ ! I would really like decode_model_indices to set the following (on a per-variable basis)
+ ! progvar(ivar)%storder
+ ! progvar(ivar)%numdims
+ ! progvar(ivar)%dimlens
+
+ call decode_model_indices( varname, progvar(ivar)%model_varname, progvar(ivar)%model_dim, &
+ progvar(ivar)%model_index, progvar(ivar)%long_name, &
+ progvar(ivar)%units)
+
+ varsize = NgridLon * NgridLat * NgridAlt
+
+ progvar(ivar)%storder = 'xyz3d'
+ progvar(ivar)%numdims = 3
+ progvar(ivar)%dimlens(1:progvar(ivar)%numdims) = (/ NgridLon, NgridLat, NgridAlt /)
+ progvar(ivar)%varsize = varsize
+ progvar(ivar)%index1 = index1
+ progvar(ivar)%indexN = index1 + varsize - 1
+ index1 = index1 + varsize ! sets up for next variable
+
+ if ( debug > 0 ) then
+ write(logfileunit,*)
+ write(logfileunit,*) trim(progvar(ivar)%varname),' variable number ',ivar
+ write(logfileunit,*) ' storage ',trim(progvar(ivar)%storder)
+ write(logfileunit,*) ' long_name ',trim(progvar(ivar)%long_name)
+ write(logfileunit,*) ' units ',trim(progvar(ivar)%units)
+ write(logfileunit,*) ' numdims ',progvar(ivar)%numdims
+ write(logfileunit,*) ' dimlens ',progvar(ivar)%dimlens(1:progvar(ivar)%numdims)
+ write(logfileunit,*) ' varsize ',progvar(ivar)%varsize
+ write(logfileunit,*) ' index1 ',progvar(ivar)%index1
+ write(logfileunit,*) ' indexN ',progvar(ivar)%indexN
+ write(logfileunit,*) ' dart_kind ',progvar(ivar)%dart_kind
+ write(logfileunit,*) ' kind_string ',trim(progvar(ivar)%kind_string)
+ write(logfileunit,*) ' model_varname ',trim(progvar(ivar)%model_varname)
+ write(logfileunit,*) ' model_dim ',progvar(ivar)%model_dim
+ write(logfileunit,*) ' model_index ',progvar(ivar)%model_index
+
+ write( * ,*)
+ write( * ,*) trim(progvar(ivar)%varname),' variable number ',ivar
+ write( * ,*) ' storage ',trim(progvar(ivar)%storder)
+ write( * ,*) ' long_name ',trim(progvar(ivar)%long_name)
+ write( * ,*) ' units ',trim(progvar(ivar)%units)
+ write( * ,*) ' numdims ',progvar(ivar)%numdims
+ write( * ,*) ' dimlens ',progvar(ivar)%dimlens(1:progvar(ivar)%numdims)
+ write( * ,*) ' varsize ',progvar(ivar)%varsize
+ write( * ,*) ' index1 ',progvar(ivar)%index1
+ write( * ,*) ' indexN ',progvar(ivar)%indexN
+ write( * ,*) ' dart_kind ',progvar(ivar)%dart_kind
+ write( * ,*) ' kind_string ',trim(progvar(ivar)%kind_string)
+ write( * ,*) ' model_varname ',trim(progvar(ivar)%model_varname)
+ write( * ,*) ' model_dim ',progvar(ivar)%model_dim
+ write( * ,*) ' model_index ',progvar(ivar)%model_index
+ endif
+
+enddo
+
+model_size = progvar(nfields)%indexN
+
+if ( debug > 0 ) then
+ write(logfileunit,'("grid: NgridLon, NgridLat, NgridAlt =",3(1x,i5))') NgridLon, NgridLat, NgridAlt
+ write( * ,'("grid: NgridLon, NgridLat, NgridAlt =",3(1x,i5))') NgridLon, NgridLat, NgridAlt
+ write(logfileunit, *)'model_size = ', model_size
+ write( * , *)'model_size = ', model_size
+endif
+
+allocate( ens_mean(model_size) )
+
+end subroutine static_init_model
+
+
+
+subroutine end_model()
+!------------------------------------------------------------------
+!
+! Does any shutdown and clean-up needed for model. Can be a NULL
+! INTERFACE if the model has no need to clean up storage, etc.
+
+if (allocated(LON)) deallocate(LON, LAT, ALT)
+
+end subroutine end_model
+
+
+
+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)
+
+! FIXME: put an error handler call here - we cannot initialize the model
+! this way and it would be an error if filter called it.
+
+end subroutine init_time
+
+
+
+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
+
+! FIXME: put an error handler call here - we cannot initialize the model
+! this way and it would be an error if filter called it.
+
+end subroutine init_conditions
+
+
+
+function nc_write_model_atts( ncFileID ) result (ierr)
+!------------------------------------------------------------------
+! TJH -- Writes the model-specific attributes to a netCDF file.
+! This includes coordinate variables and some metadata, but NOT
+! the model state vector.
+!
+! assim_model_mod:init_diag_output uses information from the location_mod
+! to define the location dimension and variable ID. All we need to do
+! is query, verify, and fill ...
+!
+! Typical sequence for adding new dimensions,variables,attributes:
+! NF90_OPEN ! open existing netCDF dataset
+! NF90_redef ! put into define mode
+! NF90_def_dim ! define additional dimensions (if any)
+! NF90_def_var ! define variables: from name, type, and dims
+! NF90_put_att ! assign attribute values
+! NF90_ENDDEF ! end definitions: leave define mode
+! NF90_put_var ! provide values for variable
+! NF90_CLOSE ! close: save updated netCDF dataset
+
+integer, intent(in) :: ncFileID ! netCDF file identifier
+integer :: ierr ! return value of function
+
+integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
+
+!----------------------------------------------------------------------
+! variables if we just blast out one long state vector
+!----------------------------------------------------------------------
+
+integer :: StateVarDimID ! netCDF pointer to state variable dimension (model size)
+integer :: MemberDimID ! netCDF pointer to dimension of ensemble (ens_size)
+integer :: TimeDimID ! netCDF pointer to time dimension (unlimited)
+
+integer :: StateVarVarID ! netCDF pointer to state variable coordinate array
+integer :: StateVarID ! netCDF pointer to 3D [state,copy,time] array
+
+!----------------------------------------------------------------------
+! variables if we parse the state vector into prognostic variables.
+!----------------------------------------------------------------------
+
+! for the dimensions and coordinate variables
+integer :: NLONDimID
+integer :: NLATDimID
+integer :: NALTDimID
+
+! for the prognostic variables
+integer :: ivar, VarID
+
+!----------------------------------------------------------------------
+! variables for the namelist output
+!----------------------------------------------------------------------
+
+character(len=129), allocatable, dimension(:) :: textblock
+integer :: LineLenDimID, nlinesDimID, nmlVarID
+integer :: nlines, linelen
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list