[Dart-dev] [4367] DART/trunk/models/tiegcm: This is the initial import of the TIEGCM interface. Tomoko an
nancy at ucar.edu
nancy at ucar.edu
Mon May 17 14:42:20 MDT 2010
Revision: 4367
Author: thoar
Date: 2010-05-17 14:42:20 -0600 (Mon, 17 May 2010)
Log Message:
-----------
This is the initial import of the TIEGCM interface. Tomoko and I are
working on it on bluefire (IBM). At this point, perfect_model seems to
be working correctly. We are working on the model_mod_check routine
to ensure we can put a synthetic observation wherever we want, then
we will test a single assimilation with no model advance. The model
can be compiled as a single threaded executable, so we are checking
async==2 first. After that works, we will explore async==4, the
anticipated method for production runs.
Added Paths:
-----------
DART/trunk/models/tiegcm/dart_to_model.f90
DART/trunk/models/tiegcm/model_mod.f90
DART/trunk/models/tiegcm/model_to_dart.f90
DART/trunk/models/tiegcm/shell_scripts/advance_model.csh
DART/trunk/models/tiegcm/shell_scripts/run_filter.csh
DART/trunk/models/tiegcm/shell_scripts/run_filter_async4.csh
DART/trunk/models/tiegcm/shell_scripts/run_perfect_model_obs.csh
DART/trunk/models/tiegcm/work/input.nml
DART/trunk/models/tiegcm/work/mkmf_create_fixed_network_seq
DART/trunk/models/tiegcm/work/mkmf_create_obs_sequence
DART/trunk/models/tiegcm/work/mkmf_dart_to_model
DART/trunk/models/tiegcm/work/mkmf_filter
DART/trunk/models/tiegcm/work/mkmf_model_mod_check
DART/trunk/models/tiegcm/work/mkmf_model_to_dart
DART/trunk/models/tiegcm/work/mkmf_obs_diag
DART/trunk/models/tiegcm/work/mkmf_perfect_model_obs
DART/trunk/models/tiegcm/work/mkmf_preprocess
DART/trunk/models/tiegcm/work/mkmf_wakeup_filter
DART/trunk/models/tiegcm/work/path_names_create_fixed_network_seq
DART/trunk/models/tiegcm/work/path_names_create_obs_sequence
DART/trunk/models/tiegcm/work/path_names_dart_to_model
DART/trunk/models/tiegcm/work/path_names_filter
DART/trunk/models/tiegcm/work/path_names_model_mod_check
DART/trunk/models/tiegcm/work/path_names_model_to_dart
DART/trunk/models/tiegcm/work/path_names_obs_diag
DART/trunk/models/tiegcm/work/path_names_perfect_model_obs
DART/trunk/models/tiegcm/work/path_names_preprocess
DART/trunk/models/tiegcm/work/path_names_wakeup_filter
DART/trunk/models/tiegcm/work/quickbuild.csh
-------------- next part --------------
Added: DART/trunk/models/tiegcm/dart_to_model.f90
===================================================================
--- DART/trunk/models/tiegcm/dart_to_model.f90 (rev 0)
+++ DART/trunk/models/tiegcm/dart_to_model.f90 2010-05-17 20:42:20 UTC (rev 4367)
@@ -0,0 +1,167 @@
+! DART software - Copyright \xA9 2004 - 2010 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 TIEGCM and DART
+!
+! method: Read DART state vector ("proprietary" format)
+! Reform state vector back into TIEGCM fields.
+! Replace those fields on the TIEGCM restart file with the new values,
+!
+! Replace 'mtime' variable in the TIEGCM restart file
+! with the 'valid time' of the DART state vector.
+! Write out updated namelist variables (e.g., model_time, adv_to_time)
+! in a file called 'namelist_update'
+!
+!----------------------------------------------------------------------
+
+use types_mod, only : r8
+use utilities_mod, only : get_unit, initialize_utilities, E_ERR, &
+ error_handler, timestamp, do_output
+use model_mod, only : model_type, get_model_size, init_model_instance, &
+ vector_to_prog_var, update_TIEGCM_restart, &
+ static_init_model
+use assim_model_mod, only : assim_model_type, aread_state_restart, &
+ open_restart_read, close_restart
+use time_manager_mod, only : time_type, get_time, get_date, set_calendar_type, &
+ print_time, print_date, set_date, set_time, &
+ operator(*), operator(+), operator(-), &
+ operator(>), operator(<), operator(/), &
+ operator(/=), operator(<=)
+
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+type(model_type) :: var
+type(time_type) :: model_time, adv_to_time, jan1, tbase, target_time
+real(r8), allocatable :: x_state(:)
+integer :: file_unit, x_size, ens_member, io
+character (len = 128) :: file_name = 'tiegcm_restart_p.nc', file_in = 'temp_ic'
+character (len = 128) :: file_namelist_out = 'namelist_update'
+integer :: model_doy, model_hour, model_minute ! current tiegcm mtime (doy,hour,mitute)
+integer :: adv_to_doy, adv_to_hour, adv_to_minute ! advance tiegcm mtime
+integer :: target_doy, target_hour, target_minute ! forecast time interval in mtime format
+integer :: utsec, year, month, day, sec, model_year
+
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+
+call initialize_utilities(progname='dart_to_model', output_flag=.true.)
+
+call static_init_model() ! reads input.nml, etc., sets the table
+x_size = get_model_size() ! now that we know how big state vector is ...
+allocate(x_state(x_size)) ! allocate space for the (empty) state vector
+
+! Open the DART model state ...
+! Read in the time to which TIEGCM must advance.
+! Read in the valid time for the model state
+! Read in state vector from DART
+
+file_unit = open_restart_read(file_in)
+
+call aread_state_restart(model_time, x_state, file_unit, adv_to_time)
+call close_restart(file_unit)
+
+if (do_output()) &
+ call print_time(model_time,'time for restart file '//trim(file_in))
+if (do_output()) &
+ call print_date(model_time,'date for restart file '//trim(file_in))
+
+if (do_output()) &
+ call print_time(adv_to_time,'time for restart file '//trim(file_in))
+if (do_output()) &
+ call print_date(adv_to_time,'date for restart file '//trim(file_in))
+
+
+! Parse the vector into TIEGCM fields (prognostic variables)
+call init_model_instance(var, model_time)
+call vector_to_prog_var(x_state, var)
+deallocate(x_state)
+
+!----------------------------------------------------------------------
+! This program writes out parameters to a file called 'namelist_update'
+! for TIEGCM namelist update used in advance_model.csh
+!----------------------------------------------------------------------
+
+file_unit = get_unit()
+open(unit = file_unit, file = file_namelist_out)
+
+! write fields to the binary TIEGCM restart file
+call update_TIEGCM_restart(file_name, var)
+
+
+! Get updated TIEGCM namelist variables
+!
+call get_date(model_time, model_year, month, day, model_hour, model_minute, sec)
+jan1 = set_date(model_year,1,1)
+tbase = model_time - jan1 ! total time since the start of the year
+call get_time(tbase, utsec, model_doy)
+model_doy = model_doy + 1 ! add jan1 back in
+
+call get_date(adv_to_time, year, month, day, adv_to_hour, adv_to_minute, sec)
+jan1 = set_date(year,1,1)
+tbase = adv_to_time - jan1 ! total time since the start of the year
+call get_time(tbase, utsec, adv_to_doy)
+adv_to_doy = adv_to_doy + 1 ! add jan1 back in
+
+! Calculate number of hours to advance tiegcm
+target_time = adv_to_time - model_time
+call get_time(target_time, sec, day)
+target_doy = day
+target_hour = sec/3600
+target_minute = (sec - target_hour*3600)/60
+
+!START_YEAR
+write(file_unit, *, iostat = io ) model_year
+if (io /= 0 )then
+ call error_handler(E_ERR,'dart_to_model:','cannot write model_year to STDOUT', &
+ source,revision,revdate)
+endif
+!START_DAY
+write(file_unit, *, iostat = io ) model_doy
+if (io /= 0 )then
+ call error_handler(E_ERR,'dart_to_model:','cannot write model_day to STDOUT', &
+ source,revision,revdate)
+endif
+!SOURCE_START, START, SECSTART
+write(file_unit, *, iostat = io ) model_doy,',',model_hour,',',model_minute
+if (io /= 0 )then
+ call error_handler(E_ERR,'dart_to_model:','cannot write mtime (day/hour/minute) to STDOUT', &
+ source,revision,revdate)
+endif
+!STOP, SECSTOP
+write(file_unit, *, iostat = io ) adv_to_doy,',',adv_to_hour,',',adv_to_minute
+if (io /= 0 )then
+ call error_handler(E_ERR,'dart_to_model:','cannot write adv_to_time mtime (day/hour/minute) to STDOUT', &
+ source,revision,revdate)
+endif
+!HIST, SAVE, SECHIST, SECSAVE
+write(file_unit, *, iostat = io ) target_doy,',',target_hour,',',target_minute
+if (io /= 0 )then
+ call error_handler(E_ERR,'dart_to_model:','cannot write target_time mtime (day/hour/minute) to STDOUT', &
+ source,revision,revdate)
+endif
+
+close(file_unit)
+
+!----------------------------------------------------------------------
+! When called with 'end', timestamp will also call finalize_utilities()
+!----------------------------------------------------------------------
+call timestamp(string1=source, pos='end')
+
+end program dart_to_model
Property changes on: DART/trunk/models/tiegcm/dart_to_model.f90
___________________________________________________________________
Added: mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author URL Id
Added: svn:eol-style
+ native
Added: DART/trunk/models/tiegcm/model_mod.f90
===================================================================
--- DART/trunk/models/tiegcm/model_mod.f90 (rev 0)
+++ DART/trunk/models/tiegcm/model_mod.f90 2010-05-17 20:42:20 UTC (rev 4367)
@@ -0,0 +1,1842 @@
+! DART software - Copyright \xA9 2004 - 2010 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$
+
+!------------------------------------------------------------------
+!
+! Interface for HAO-TIEGCM
+!
+!------------------------------------------------------------------
+
+! DART Modules
+use types_mod, only : r8, digits12
+use time_manager_mod, only : time_type, set_calendar_type, set_time_missing, &
+ set_time, get_time, print_time, &
+ set_date, get_date, 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, get_location, query_location, &
+ get_dist, vert_is_height
+use utilities_mod, only : file_exist, open_file, close_file, &
+ error_handler, E_ERR, E_MSG, E_WARN, nmlfileunit, &
+ do_output, find_namelist_in_file, check_namelist_read, &
+ do_nml_file, do_nml_term, nc_check, &
+ register_module
+use obs_kind_mod, only : KIND_U_WIND_COMPONENT, &! just for definition
+ KIND_V_WIND_COMPONENT, &! just for definition
+ KIND_TEMPERATURE, &! just for definition
+ KIND_ELECTRON_DENSITY, &! for Ne obs
+ KIND_DENSITY ! for neutral density observation
+ ! (for now [O]; minor species ignored)
+use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian
+
+use typesizes
+use netcdf
+
+implicit none
+private
+
+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
+!TIEGCM specific routines
+public :: model_type, &
+ init_model_instance, &
+ end_model_instance, &
+ prog_var_to_vector, &
+ vector_to_prog_var, &
+ read_TIEGCM_restart, &
+ update_TIEGCM_restart, &
+ read_TIEGCM_definition, &
+ read_TIEGCM_secondary
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+!------------------------------------------------------------------
+! define model parameters
+
+
+integer :: nilev, nlev, nlon, nlat
+real(r8),dimension(:), allocatable :: lons, lats, levs, ilevs
+real(r8),dimension(:,:,:),allocatable :: ZGtiegcm !! auxiliary variable(geopotential height[cm])
+real :: TIEGCM_missing_value !! global attribute
+
+integer :: time_step_seconds = 120 ! tiegcm time step
+integer :: time_step_days = 0 ! tiegcm time step
+character (len=19) :: restart_file = 'tiegcm_restart_p.nc'
+character (len=11) :: secondary_file = 'tiegcm_s.nc'
+
+integer, parameter :: TYPE_local_NE = 0
+integer, parameter :: TYPE_local_TN = 1
+integer, parameter :: TYPE_local_TN_NM = 2
+integer, parameter :: TYPE_local_UN = 3
+integer, parameter :: TYPE_local_UN_NM = 4
+integer, parameter :: TYPE_local_VN = 5
+integer, parameter :: TYPE_local_VN_NM = 6
+integer, parameter :: TYPE_local_O1 = 7
+integer, parameter :: TYPE_local_O1_NM = 8
+
+type(time_type) :: time_step
+integer :: model_size
+
+type model_type
+ real(r8), pointer :: vars_3d(:,:,:,:)
+ type(time_type) :: valid_time
+end type model_type
+
+integer :: state_num_3d = 9
+ ! -- interface levels --
+ ! NE
+ ! -- midpoint levels --
+ ! O1 O1_NM (O2 O2_NM NO NO_NM excluded for now)
+ ! -- midpoint levels; top slot missing --
+ ! TN TN_NM UN UN_NM VN VN_NM
+logical :: output_state_vector = .false.
+ ! .true. results in a "state-vector" netCDF file
+ ! .false. results in a "prognostic-var" netCDF file
+namelist /model_nml/ output_state_vector, state_num_3d
+
+logical :: first_pert_call = .true.
+type(random_seq_type) :: random_seq
+
+
+
+!------------------------------------------------------------------
+
+character(len = 129) :: msgstring
+logical, save :: module_initialized = .false.
+
+contains
+
+!==================================================================
+
+
+subroutine static_init_model()
+!------------------------------------------------------------------
+!
+! Called to do one time initialization of the model. As examples,
+! might define information about the model size or model timestep.
+! In models that require pre-computed static data, for instance
+! spherical harmonic weights, these would also be computed here.
+! Can be a NULL INTERFACE for the simplest models.
+
+ integer :: i
+ integer :: iunit, io
+
+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 namelist entry for model_mod from input.nml
+!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")
+
+!if (do_nml_file()) write(nmlfileunit, nml=model_nml)
+!if (do_nml_term()) write( * , nml=model_nml)
+
+! Reading in TIEGCM grid definition etc from TIEGCM restart file
+call read_TIEGCM_definition(restart_file)
+
+! Reading in Geopotential Height from TIEGCM secondary output file
+call read_TIEGCM_secondary(secondary_file)
+
+! Compute overall model size
+model_size = nlon * nlat * nlev * state_num_3d
+
+if (do_output()) write(*,*)'nlon = ',nlon
+if (do_output()) write(*,*)'nlat = ',nlat
+if (do_output()) write(*,*)'nlev = ',nlev
+if (do_output()) write(*,*)'n3D = ',state_num_3d
+
+! Might as well use the Gregorian Calendar
+call set_calendar_type('Gregorian')
+
+! The time_step in terms of a time type must also be initialized.
+time_step = set_time(time_step_seconds, time_step_days)
+
+
+
+end subroutine static_init_model
+
+
+
+subroutine init_conditions(x)
+!------------------------------------------------------------------
+! 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
+
+end subroutine init_conditions
+
+
+
+subroutine adv_1step(x, time)
+!------------------------------------------------------------------
+! 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 of filter or if the
+! program integrate_model is to be used to advance the model
+! state as a separate executable. If one of these options
+! is not going to be 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
+
+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, itype, obs_val, istatus)
+!------------------------------------------------------------------
+!
+! Given a state vector, a location, and a model state variable type,
+! interpolates the state variable field to that location and returns
+! the value in obs_val. The istatus variable should be returned as
+! 0 unless there is some problem in computing the interpolation in
+! which case an alternate value should be returned. The itype variable
+! is a model specific integer that specifies the type of field (for
+! instance temperature, zonal wind component, etc.). In low order
+! models that have no notion of types of variables, this argument can
+! be ignored. For applications in which only perfect model experiments
+! with identity observations (i.e. only the value of a particular
+! state variable is observed), this can be a NULL INTERFACE.
+
+real(r8), intent(in) :: x(:)
+type(location_type), intent(in) :: location
+integer, intent(in) :: itype
+real(r8), intent(out) :: obs_val
+integer, intent(out) :: istatus
+
+integer :: i, vstatus, which_vert
+integer :: lat_below, lat_above, lon_below, lon_above
+real(r8) :: lon_fract, temp_lon, lat_fract
+real(r8) :: lon, lat, height, lon_lat_lev(3)
+real(r8) :: bot_lon, top_lon, delta_lon, bot_lat, top_lat, delta_lat
+real(r8) :: val(2,2), a(2)
+
+if ( .not. module_initialized ) call static_init_model
+
+! Default for successful return
+istatus = 0
+vstatus = 0
+
+! Get the position, determine if it is model level or pressure in vertical
+lon_lat_lev = get_location(location)
+lon = lon_lat_lev(1) ! degree
+lat = lon_lat_lev(2) ! degree
+if(vert_is_height(location)) then
+ height = lon_lat_lev(3)
+else
+ which_vert = nint(query_location(location))
+ write(msgstring,*) 'vertical coordinate type:',which_vert,' cannot be handled'
+ call error_handler(E_ERR,'model_interpolate',msgstring,source,revision,revdate)
+endif
+
+! Get lon and lat grid specs
+bot_lon = lons(1) ! 0
+top_lon = lons(nlon) ! 355
+delta_lon = abs((lons(1)-lons(2))) ! 5
+bot_lat = lats(1) ! -87.5
+top_lat = lats(nlat) ! 87.5
+delta_lat = abs((lats(1)-lats(2))) ! 5
+
+
+! Compute bracketing lon indices: DART [0, 360] TIEGCM [0, 355]
+if(lon >= bot_lon .and. lon <= top_lon) then ! 0 <= lon <= 355
+ lon_below = int((lon - bot_lon) / delta_lon) + 1
+ lon_above = lon_below + 1
+ lon_fract = (lon - lons(lon_below)) / delta_lon
+elseif (lon < bot_lon) then ! lon < 0
+ temp_lon = lon + 360.0
+ lon_below = int((temp_lon - bot_lon) / delta_lon) + 1
+ lon_above = lon_below + 1
+ lon_fract = (temp_lon - lons(lon_below)) / delta_lon
+
+elseif (lon >= (top_lon+delta_lon)) then ! 360 <= lon
+ temp_lon = lon - 360.0
+ lon_below = int((temp_lon - bot_lon) / delta_lon) + 1
+ lon_above = lon_below + 1
+ lon_fract = (temp_lon - lons(lon_below)) / delta_lon
+else ! 355 < lon < 360 at wraparound point
+ lon_below = nlon
+ lon_above = 1
+ lon_fract = (lon - top_lon) / delta_lon
+endif
+
+! compute neighboring lat rows: TIEGCM [-87.5, 87.5] DART [-90, 90]
+! NEED TO BE VERY CAREFUL ABOUT POLES; WHAT'S BEING DONE IS NOT GREAT!
+if(lat >= bot_lat .and. lat <= top_lat) then ! -87.5 <= lat <= 87.5
+ lat_below = int((lat - bot_lat) / delta_lat) + 1
+ lat_above = lat_above + 1
+ lat_fract = (lat - lats(lat_below) ) / delta_lat
+else if(lat < bot_lat) then ! South of bottom lat
+ lat_below = 1
+ lat_above = 1
+ lat_fract = 1.
+else ! North of top lat
+ lat_below = nlat
+ lat_above = nlat
+ lat_fract = 1.
+endif
+
+! Now, need to find the values for the four corners
+ call get_val(val(1, 1), x, lon_below, lat_below, height, itype, vstatus)
+ if (vstatus /= 1) call get_val(val(1, 2), x, lon_below, lat_above, height, itype, vstatus)
+ if (vstatus /= 1) call get_val(val(2, 1), x, lon_above, lat_below, height, itype, vstatus)
+ if (vstatus /= 1) call get_val(val(2, 2), x, lon_above, lat_above, height, itype, vstatus)
+
+
+! istatus meaning return expected obs? assimilate?
+! 0 obs and model are fine; yes yes
+! 1 fatal problem; no no
+! 2 exclude valid obs yes no
+!
+istatus = vstatus
+if(istatus /= 1) then
+ do i = 1, 2
+ a(i) = lon_fract * val(2, i) + (1.0 - lon_fract) * val(1, i)
+ end do
+ obs_val = lat_fract * a(2) + (1.0 - lat_fract) * a(1)
+else
+ obs_val = 0.
+endif
+
+!print*, 'model_interpolate', lon, lat, height,obs_val
+
+
+end subroutine model_interpolate
+
+
+
+subroutine get_val(val, x, lon_index, lat_index, height, obs_kind, istatus)
+!------------------------------------------------------------------
+!
+real(r8), intent(out) :: val
+real(r8), intent(in) :: x(:)
+integer, intent(in) :: lon_index, lat_index
+real(r8), intent(in) :: height
+integer, intent(in) :: obs_kind
+integer, intent(out) :: istatus
+
+integer :: var_type
+integer :: k, lev_top, lev_bottom
+real(r8) :: zgrid, delta_z
+real(r8) :: val_top, val_bottom, frac_lev
+
+! No errors to start with
+istatus = 0
+
+! To find a layer height: what's the unit of height [cm] ?
+h_loop:do k = 1, nlev
+ zgrid = ZGtiegcm(lon_index,lat_index,k) ![cm]
+ if (height <= zgrid) then
+ lev_top = k
+ lev_bottom = lev_top -1
+ delta_z = zgrid - ZGtiegcm(lon_index,lat_index,lev_bottom)
+ frac_lev = (zgrid - height)/delta_z
+ exit h_loop
+ endif
+enddo h_loop
+
+if (obs_kind == KIND_DENSITY) then
+
+ var_type = TYPE_local_O1
+ val_top = x(get_index(lat_index, lon_index, lev_top, var_type))
+ val_bottom = x(get_index(lat_index, lon_index, lev_bottom, var_type))
+
+elseif (obs_kind == KIND_ELECTRON_DENSITY) then
+
+ var_type = TYPE_local_NE
+ val_top = x(get_index(lat_index, lon_index, lev_top, var_type))
+ val_bottom = x(get_index(lat_index, lon_index, lev_bottom, var_type))
+
+else
+
+ istatus = 1
+ val = 0.
+ return
+
+end if
+
+val = frac_lev * val_bottom + (1.0 - frac_lev) * val_top
+
+end subroutine get_val
+
+
+
+function get_index(lat_index, lon_index, lev_index, var_type)
+!------------------------------------------------------------------
+!
+integer, intent(in) :: lat_index, lon_index, lev_index, var_type
+integer :: get_index
+
+get_index = 1 + var_type + (lev_index -1)*state_num_3d &
+ + (lat_index -1)*state_num_3d*nlev &
+ + (lon_index -1)*state_num_3d*nlev*nlat
+
+end function get_index
+
+
+
+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 = time_step
+
+end function get_model_time_step
+
+
+
+subroutine get_state_meta_data(index_in, location, var_type)
+!------------------------------------------------------------------
+!
+! Given an integer index into the state vector structure, returns the
+! associated location. A second intent(out) optional argument kind
+! can be returned if the model has more than one type of field (for
+! instance temperature and zonal wind component). This interface is
+! required for all filter applications as it is required for computing
+! the distance between observations and state variables.
+
+integer, intent(in) :: index_in
+type(location_type), intent(out) :: location
+integer, optional, intent(out) :: var_type
+
+integer :: indx, num_per_col, col_num, col_elem
+integer :: lon_index, lat_index, lev_index
+real(r8) :: lon, lat, lev
+integer :: local_var_type, var_type_temp
+
+if ( .not. module_initialized ) call static_init_model
+
+! Easier to compute with a 0 to size -1 index
+indx = index_in -1
+
+! Compute number of items per column
+num_per_col = nlev * state_num_3d
+
+! What column is this index in
+col_num = indx / num_per_col
+col_elem = indx - col_num * num_per_col
+
+! what lon and lat index for this column
+lon_index = col_num /nlat
+lat_index = col_num - lon_index * nlat
+
+! Now figure out which beast in column this is
+lev_index = col_elem / state_num_3d
+
+! Get actual lon lat values from static_init_model arrays
+lon = lons(lon_index + 1)
+lat = lats(lat_index + 1)
+
+! Find which var_type this element is
+var_type_temp = mod(col_elem, state_num_3d)
+
+if (var_type_temp == 0) then !NE
+ local_var_type = KIND_ELECTRON_DENSITY
+else if (var_type_temp == 1) then !TN
+ local_var_type = KIND_TEMPERATURE
+else if (var_type_temp == 2) then !TN_NM
+ local_var_type = KIND_TEMPERATURE
+else if (var_type_temp == 3) then !UN
+ local_var_type = KIND_U_WIND_COMPONENT
+else if (var_type_temp == 4) then !UN_NM
+ local_var_type = KIND_U_WIND_COMPONENT
+else if (var_type_temp == 5) then !VN
+ local_var_type = KIND_V_WIND_COMPONENT
+else if (var_type_temp == 6) then !VN_NM
+ local_var_type = KIND_V_WIND_COMPONENT
+else if (var_type_temp == 7) then !O1
+ local_var_type = KIND_DENSITY
+else if (var_type_temp == 8) then !O1_NM
+ local_var_type = KIND_DENSITY
+else
+ write(msgstring,*)"unknown var_type for index ",index_in
+ call error_handler(E_ERR,"get_state_meta_data", msgstring, source, revision, revdate)
+endif
+
+if (var_type_temp == 0) then !NE defined at interface levels
+ lev = ilevs(lev_index + 1)
+else
+ lev = levs(lev_index + 1) !TN UN VN O1 defined at midpoints
+endif
+
+location = set_location(lon,lat,lev,2) ! 2 == pressure (hPa) 3 == height
+
+! If the type is wanted, return it
+if(present(var_type)) var_type = local_var_type
+
+end subroutine get_state_meta_data
+
+
+
+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 ( .not. module_initialized ) call static_init_model
+
+end subroutine end_model
+
+
+
+function nc_write_model_atts( ncFileID ) result (ierr)
+!------------------------------------------------------------------
+! TJH 24 Oct 2006 -- Writes the model-specific attributes to a netCDF file.
+! This includes coordinate variables and some metadata, but NOT
+! the model state vector. We do have to allocate SPACE for the model
+! state vector, but that variable gets filled as the model advances.
+!
+! As it stands, this routine will work for ANY model, with no modification.
+!
+! The simplest possible netCDF file would contain a 3D field
+! containing the state of 'all' the ensemble members. This requires
+! three coordinate variables -- one for each of the dimensions
+! [model_size, ensemble_member, time]. A little metadata is useful,
+! so we can also create some 'global' attributes.
+! This is what is implemented here.
+!
+! Once the simplest case is working, this routine (and nc_write_model_vars)
+! can be extended to create a more logical partitioning of the state vector,
+! fundamentally creating a netCDF file with variables that are easily
+! plotted. The bgrid model_mod is perhaps a good one to view, keeping
+! in mind it is complicated by the fact it has two coordinate systems.
+! There are stubs in this template, but they are only stubs.
+!
+! TJH 29 Jul 2003 -- for the moment, all errors are fatal, so the
+! return code is always '0 == normal', since the fatal errors stop execution.
+!
+! 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
+
+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
+
+integer :: TNVarID, TN_NMVarID, UNVarID, UN_NMVarID, VNVarID, VN_NMVarID
+integer :: O1VarID, O1_NMVarID
+integer :: NEVarID
+integer :: lonDimID, latDimID, levDimID, ilevDimID
+integer :: lonVarID, latVarID, levVarID, ilevVarID
+
+! we are going to need these to record the creation date in the netCDF file.
+! This is entirely optional, but nice.
+
+character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
+character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
+character(len=5) :: crzone ! needed by F90 DATE_AND_TIME intrinsic
+integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic
+character(len=NF90_MAX_NAME) :: str1
+
+integer :: i
+
+if ( .not. module_initialized ) call static_init_model
+
+!-------------------------------------------------------------------------------
+! make sure ncFileID refers to an open netCDF file,
+! and then put into define mode.
+!-------------------------------------------------------------------------------
+
+ierr = -1 ! assume things go poorly
+
+call nc_check(nf90_Inquire(ncFileID, nDimensions, nVariables, &
+ nAttributes, unlimitedDimID), 'nc_write_model_atts','inquire')
+call nc_check(nf90_Redef(ncFileID),'nc_write_model_atts','redef')
+
+!-------------------------------------------------------------------------------
+! We need the dimension ID for the number of copies/ensemble members, and
+! we might as well check to make sure that Time is the Unlimited dimension.
+! Our job is create the 'model size' dimension.
+!-------------------------------------------------------------------------------
+
+call nc_check(nf90_inq_dimid(ncid=ncFileID, name="copy", dimid=MemberDimID),&
+ 'nc_write_model_atts', 'copy dimid')
+call nc_check(nf90_inq_dimid(ncid=ncFileID, name="time", dimid= TimeDimID),&
+ 'nc_write_model_atts', 'time dimid')
+
+if ( TimeDimID /= unlimitedDimId ) then
+ write(msgstring,*)"Time Dimension ID ",TimeDimID, &
+ " should equal Unlimited Dimension ID",unlimitedDimID
+ call error_handler(E_ERR,"nc_write_model_atts", msgstring, source, revision, revdate)
+endif
+
+!-------------------------------------------------------------------------------
+! Define the model size / state variable dimension / whatever ...
+!-------------------------------------------------------------------------------
+call nc_check(nf90_def_dim(ncid=ncFileID, name="StateVariable", &
+ len=model_size, dimid = StateVarDimID),&
+ 'nc_write_model_atts', 'state def_dim')
+
+!-------------------------------------------------------------------------------
+! Write Global Attributes
+!-------------------------------------------------------------------------------
+
+call DATE_AND_TIME(crdate,crtime,crzone,values)
+write(str1,'(''YYYY MM DD HH MM SS = '',i4,5(1x,i2.2))') &
+ values(1), values(2), values(3), values(5), values(6), values(7)
+
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "creation_date" ,str1 ),&
+ 'nc_write_model_atts', 'creation put')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_source" ,source ),&
+ 'nc_write_model_atts', 'source put')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_revision",revision),&
+ 'nc_write_model_atts', 'revision put')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_revdate" ,revdate ),&
+ 'nc_write_model_atts', 'revdate put')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model","TIEGCM" ),&
+ 'nc_write_model_atts', 'model put')
+
+!-------------------------------------------------------------------------------
+! Here is the extensible part. The simplest scenario is to output the state vector,
+! parsing the state vector into model-specific parts is complicated, and you need
+! to know the geometry, the output variables (PS,U,V,T,Q,...) etc. We're skipping
+! complicated part.
+!-------------------------------------------------------------------------------
+
+if ( output_state_vector ) then
+
+ !----------------------------------------------------------------------------
+ ! Create a variable for the state vector
+ !----------------------------------------------------------------------------
+
+ ! Define the state vector coordinate variable and some attributes.
+ call nc_check(nf90_def_var(ncid=ncFileID,name="StateVariable", xtype=nf90_int, &
+ dimids=StateVarDimID, varid=StateVarVarID), &
+ 'nc_write_model_atts', 'statevariable def_var')
+ call nc_check(nf90_put_att(ncFileID, StateVarVarID, "long_name", "State Variable ID"), &
+ 'nc_write_model_atts', 'statevariable long_name')
+ call nc_check(nf90_put_att(ncFileID, StateVarVarID, "units", "indexical"), &
+ 'nc_write_model_atts', 'statevariable units')
+ call nc_check(nf90_put_att(ncFileID, StateVarVarID, "valid_range", (/ 1, model_size /)), &
+ 'nc_write_model_atts', 'statevariable valid_range')
+
+ ! Define the actual (3D) state vector, which gets filled as time goes on ...
+ call nc_check(nf90_def_var(ncid=ncFileID, name="state", xtype=nf90_real, &
+ dimids = (/ StateVarDimID, MemberDimID, unlimitedDimID /), &
+ varid=StateVarID), 'nc_write_model_atts', 'state def_var')
+ call nc_check(nf90_put_att(ncFileID, StateVarID, "long_name", "model state or fcopy"), &
+ 'nc_write_model_atts', 'state long_name')
+
+ ! Leave define mode so we can fill the coordinate variable.
+ call nc_check(nf90_enddef(ncfileID), 'nc_write_model_atts', 'state enddef')
+
+ ! Fill the state variable coordinate variable
+ call nc_check(nf90_put_var(ncFileID, StateVarVarID, (/ (i,i=1,model_size) /) ), &
+ 'nc_write_model_atts', 'state put_var')
+
+else
+
+ !----------------------------------------------------------------------------
+ ! We need to process the prognostic variables.
+ !----------------------------------------------------------------------------
+
+ ! This block is a stub for something more complicated.
+ ! Usually, the control for the execution of this block is a namelist variable.
+ ! Take a peek at the bgrid model_mod.f90 for a (rather complicated) example.
+
+ !----------------------------------------------------------------------------
+ ! Define the dimensions IDs
+ !----------------------------------------------------------------------------
+
+ call nc_check(nf90_def_dim(ncid=ncFileID, name="lon", &
+ & len = nlon, dimid = lonDimID), 'nc_write_model_atts')
+ call nc_check(nf90_def_dim(ncid=ncFileID, name="lat", &
+ & len = nlat, dimid = latDimID), 'nc_write_model_atts')
+ call nc_check(nf90_def_dim(ncid=ncFileID, name="lev", &
+ & len = nlev, dimid = levDimID), 'nc_write_model_atts')
+ call nc_check(nf90_def_dim(ncid=ncFileID, name="ilev", &
+ & len = nilev, dimid = ilevDimID), 'nc_write_model_atts')
+
+ !----------------------------------------------------------------------------
+ ! Create the (empty) Variables and the Attributes
+ !----------------------------------------------------------------------------
+
+ call nc_check(nf90_def_var(ncFileID, name="lon", &
+ & xtype=nf90_double, dimids=lonDimID, varid=lonVarID),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, lonVarID, &
+ & "long_name", "geographic longitude (-west, +east)"),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, lonVarID, "units", "degrees_east"),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, lonVarID, "valid_range", &
+ & (/ -180.0_r8, 180.0_r8 /)),'nc_write_model_atts')
+
+ call nc_check(nf90_def_var(ncFileID, name="lat", &
+ & xtype=nf90_double, dimids=latDimID, varid=latVarID),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, latVarID, &
+ & "long_name", "geographic latitude (-south +north)"),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, latVarID, "units", "degrees_north"),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, latVarID, "valid_range", &
+ & (/ -90.0_r8, 90.0_r8 /)),'nc_write_model_atts')
+
+ call nc_check(nf90_def_var(ncFileID, name="lev", &
+ & xtype=nf90_double, dimids=levDimID, varid=levVarID),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, levVarID, "long_name", "midpoint levels"),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, levVarID, "units", "hPa"),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, levVarID, "positive", "up"),&
+ 'nc_write_model_atts')
+
+ call nc_check(nf90_def_var(ncFileID, name="ilev", &
+ & xtype=nf90_double, dimids=ilevDimID, varid=ilevVarID),&
+ 'nc_write_model_atts')
+ call nc_check(nf90_put_att(ncFileID, ilevVarID, "long_name", "interface levels"),&
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list