[Dart-dev] [7775] DART/trunk: This represents a big change to the TIEGCM model_mod.
nancy at ucar.edu
nancy at ucar.edu
Fri Mar 27 16:44:01 MDT 2015
Revision: 7775
Author: thoar
Date: 2015-03-27 16:44:00 -0600 (Fri, 27 Mar 2015)
Log Message:
-----------
This represents a big change to the TIEGCM model_mod.
The composition of the DART state vector can now be specified at run time
through the model_nml:variables construct. This has been implemented for more
than a year on the tiegcm branch and is being merged onto the trunk.
This change requires that static_init_model() have access to a tiegcm_restart_p.nc
and tiegcm_s.nc file as well as a tiegcm.nml.
Since the tiegcm branch cannot be directly merged onto the trunk,
these changes are duplicates of those made to the branch. The revision
history only exists on the tiegcm branch unfortunately.
I am backing out a change made to location_mod.f90, that was not necessary.
The print statement has enough precision to handle up to 999 km, TIEGCM
never gets quite that high.
Modified Paths:
--------------
DART/trunk/location/threed_sphere/location_mod.f90
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_perfect_model_obs.csh
DART/trunk/models/tiegcm/work/input.nml
DART/trunk/obs_def/obs_def_upper_atm_mod.f90
DART/trunk/obs_kind/DEFAULT_obs_kind_mod.F90
DART/trunk/observations/SSUSI/convert_f16_edr_dsk.f90
DART/trunk/utilities/model_mod_check.f90
Added Paths:
-----------
DART/trunk/models/tiegcm/model_mod.nml
-------------- next part --------------
Modified: DART/trunk/location/threed_sphere/location_mod.f90
===================================================================
--- DART/trunk/location/threed_sphere/location_mod.f90 2015-03-27 21:08:43 UTC (rev 7774)
+++ DART/trunk/location/threed_sphere/location_mod.f90 2015-03-27 22:44:00 UTC (rev 7775)
@@ -803,7 +803,7 @@
case (VERTISPRESSURE)
write(charstring, '(A,1X,F12.7,A)') trim(string1), loc%vloc / 100.0_r8, ' hPa'
case (VERTISHEIGHT)
- write(charstring, '(A,1X,F20.7,A)') trim(string1), loc%vloc / 1000.0_r8, ' km'
+ write(charstring, '(A,1X,F12.7,A)') trim(string1), loc%vloc / 1000.0_r8, ' km'
case (VERTISSCALEHEIGHT)
write(charstring, '(A,1X,F12.7,A)') trim(string1), loc%vloc, ' scale ht'
case default
Modified: DART/trunk/models/tiegcm/dart_to_model.f90
===================================================================
--- DART/trunk/models/tiegcm/dart_to_model.f90 2015-03-27 21:08:43 UTC (rev 7774)
+++ DART/trunk/models/tiegcm/dart_to_model.f90 2015-03-27 22:44:00 UTC (rev 7775)
@@ -4,37 +4,40 @@
!
! $Id$
+!----------------------------------------------------------------------
+
+!> purpose: interface between TIEGCM and DART
+!>
+!> method: Read DART state vector.
+!> 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'
+!>
+!> 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'.
+
program dart_to_model
-!----------------------------------------------------------------------
-! 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, E_MSG, &
- error_handler, finalize_utilities, 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
+ error_handler, finalize_utilities, do_output, &
+ find_namelist_in_file, check_namelist_read, &
+ logfileunit
+use model_mod, only : get_model_size, static_init_model, get_f107_value, &
+ get_restart_file_name, update_TIEGCM_restart
+use assim_model_mod, only : 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, &
+ 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
@@ -43,77 +46,113 @@
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate = "$Date$"
-type(model_type) :: var
+!-----------------------------------------------------------------------
+! namelist parameters with default values.
+!-----------------------------------------------------------------------
+
+character(len=256) :: file_in = 'dart_restart'
+character(len=256) :: file_namelist_out = 'namelist_update'
+logical :: advance_time_present = .false.
+
+namelist /dart_to_model_nml/ file_in, file_namelist_out, advance_time_present
+
+!-----------------------------------------------------------------------
+! global storage
+!-----------------------------------------------------------------------
+
+character(len=256) :: file_name
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
+integer :: iunit, io, file_unit, x_size
-!----------------------------------------------------------------------
-!----------------------------------------------------------------------
+!======================================================================
call initialize_utilities(progname='dart_to_model', output_flag=.true.)
-call static_init_model() ! reads input.nml, etc., sets the table
+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")
+
+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.
+! Open the DART model state ...
+! (possibly) 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)
+if ( advance_time_present ) then
+ call aread_state_restart(model_time, x_state, file_unit, adv_to_time)
+else
+ call aread_state_restart(model_time, x_state, file_unit)
+endif
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()) then
+ call print_time(model_time,'time of restart file '//trim(file_in))
+ call print_date(model_time,'date of restart file '//trim(file_in))
+endif
-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))
+if (do_output() .and. advance_time_present) then
+ call print_time(adv_to_time,'advance-to-time for '//trim(file_in))
+ call print_date(adv_to_time,'advance-to-date for '//trim(file_in))
+endif
+! write fields to the binary TIEGCM restart file
-! Parse the vector into TIEGCM fields (prognostic variables)
-call init_model_instance(var, model_time)
-call vector_to_prog_var(x_state, var)
+file_name = get_restart_file_name()
+call update_TIEGCM_restart(x_state, file_name, model_time)
+
+if ( advance_time_present ) then
+ ! update TIEGCM namelist variables used in advance_model.csh
+ call write_tiegcm_time_control(file_namelist_out, model_time, adv_to_time)
+endif
+
deallocate(x_state)
-!----------------------------------------------------------------------
-! This program writes out parameters to a file called 'namelist_update'
-! for TIEGCM namelist update used in advance_model.csh
-!----------------------------------------------------------------------
+write( * ,*)''
+write(logfileunit,*)''
+call error_handler(E_MSG,'dart_to_model','Finished successfully.',source,revision,revdate)
+call finalize_utilities()
-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)
+!======================================================================
+contains
+!======================================================================
-! Get updated TIEGCM namelist variables
-!
+subroutine write_tiegcm_time_control(filename, model_time, adv_to_time)
+character(len=*), intent(in) :: filename
+type(time_type), intent(in) :: model_time
+type(time_type), intent(in) :: adv_to_time
+
+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
+
+real(r8) :: f10_7
+
+! Write updated TIEGCM namelist variables to a text file.
+! It is up to advance_model.csh to update the TIEGCM namelist.
+
+file_unit = get_unit()
+open(unit = file_unit, file = trim(filename))
+
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
+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
+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
+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
+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
@@ -122,57 +161,66 @@
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
-!F107
-if (size(var%vars_1d) > 0) then
- write(file_unit, *, iostat = io ) var%vars_1d(1)
- if (io /= 0 )then
- call error_handler(E_ERR,'dart_to_model:','cannot write f107 (var%vars_1d) to STDOUT', &
- source,revision,revdate)
- endif
-else
- write(file_unit, *, iostat = io ) 'NA'
- if (io /= 0 )then
- call error_handler(E_ERR,'dart_to_model:','cannot write f107 (var%vars_1d) to STDOUT', &
- source,revision,revdate)
- endif
-endif
+! Write the information to a text file so we can grep the desired strings and
+! then update the right parts of the tiegcm.nml - without having to write the
+! whole tiegcm.nml.
+! NOTE: The free-format write implicitly puts a blank at the start of the 'line'.
+! This is REQUIRED for the logic in advance_model.csh to correctly parse
+! the namelist_update file.
+write(file_unit, *, iostat=io) 'START_YEAR = ',model_year
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write START_YEAR to '//trim(filename),source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'START_DAY = ',model_doy
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write START_DAY to '//trim(filename),source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'SOURCE_START = ',model_doy,',',model_hour,',',model_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write SOURCE_START to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'START = ',model_doy,',',model_hour,',',model_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write START to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'SECSTART = ',model_doy,',',model_hour,',',model_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write SECSTART to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'STOP = ',adv_to_doy,',',adv_to_hour,',',adv_to_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write STOP to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'SECSTOP = ',adv_to_doy,',',adv_to_hour,',',adv_to_minute
+if (io /= 0 ) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write SECSTOP to '//trim(filename),source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'HIST = ', target_doy,',',target_hour,',',target_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write HIST to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'SAVE = ', target_doy,',',target_hour,',',target_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write SAVE to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'SECHIST = ', target_doy,',',target_hour,',',target_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write SECHIST to '//trim(filename), source,revision,revdate)
+
+write(file_unit, *, iostat=io) 'SECSAVE = ', target_doy,',',target_hour,',',target_minute
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write SECSAVE to '//trim(filename), source,revision,revdate)
+
+f10_7 = get_f107_value(x_state)
+write(file_unit, '('' F107 = '',f14.5)', iostat=io) f10_7
+if (io /= 0) call error_handler(E_ERR,'dart_to_model:', &
+ 'cannot write F107 to '//trim(filename),source,revision,revdate)
+
close(file_unit)
-call error_handler(E_MSG,'dart_to_model','Finished successfully.',source,revision,revdate)
-call finalize_utilities()
+end subroutine write_tiegcm_time_control
-
end program dart_to_model
! <next few lines under version control, do not edit>
Modified: DART/trunk/models/tiegcm/model_mod.f90
===================================================================
--- DART/trunk/models/tiegcm/model_mod.f90 2015-03-27 21:08:43 UTC (rev 7774)
+++ DART/trunk/models/tiegcm/model_mod.f90 2015-03-27 22:44:00 UTC (rev 7775)
@@ -6,47 +6,57 @@
module model_mod
-!------------------------------------------------------------------
+!-------------------------------------------------------------------------------
!
-! Interface for HAO-TIEGCM
+! Interface for HAO-TIEGCM
!
-!------------------------------------------------------------------
-! DART Modules
-use types_mod, only : r8, digits12, missing_r8, i4, PI
+!-------------------------------------------------------------------------------
+
+use types_mod, only : r4, r8, MISSING_R8, MISSING_R4, PI, &
+ earth_radius, gravity, obstypelength
+
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, &
+ 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, loc_get_close_obs => get_close_obs,&
+ get_close_obs_init, loc_get_close_obs => get_close_obs,&
set_location, get_location, query_location, &
- get_dist, vert_is_height, horiz_dist_only, &
- get_close_type, &
- VERTISPRESSURE, VERTISHEIGHT, vert_is_pressure
-use utilities_mod, only : file_exist, open_file, close_file, &
- error_handler, E_ERR, E_MSG, E_WARN, nmlfileunit, &
+ get_dist, vert_is_height, horiz_dist_only, &
+ get_close_type, vert_is_undef, VERTISUNDEF, &
+ VERTISPRESSURE, VERTISHEIGHT, vert_is_level
+
+use utilities_mod, only : file_exist, open_file, close_file, logfileunit, &
+ 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, &! neutral density obs
- KIND_PRESSURE, &! neutral density obs
- KIND_ELECTRON_DENSITY, &! Ne obs
- KIND_ATOMIC_OXYGEN_MIXING_RATIO, &! neutral density obs
- KIND_MOLEC_OXYGEN_MIXING_RATIO, &! neutral density obs
- KIND_1D_PARAMETER, &! just for definition
- KIND_GEOPOTENTIAL_HEIGHT
+ do_nml_file, do_nml_term, nc_check, register_module, &
+ file_to_text, find_textfile_dims, to_upper
+
+use obs_kind_mod, only : KIND_U_WIND_COMPONENT, &
+ KIND_V_WIND_COMPONENT, &
+ KIND_TEMPERATURE, &! neutral temperature obs
+ KIND_PRESSURE, &! neutral pressure obs
+ KIND_MOLEC_OXYGEN_MIXING_RATIO, &! neutral composition obs
+ KIND_1D_PARAMETER, &
+ KIND_GEOPOTENTIAL_HEIGHT, &
+ KIND_GEOMETRIC_HEIGHT, &
+ KIND_VERTICAL_TEC, &! total electron content
+ get_raw_obs_kind_index, get_raw_obs_kind_name
+
use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian
-use mpi_utilities_mod,only : my_task_id
+
+use mpi_utilities_mod,only : my_task_id
+
use typesizes
use netcdf
implicit none
private
+!DART mandatory public interfaces
public :: get_model_size, &
adv_1step, &
get_state_meta_data, &
@@ -63,111 +73,153 @@
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, &
- read_TIEGCM_namelist
+public :: read_TIEGCM_restart, &
+ update_TIEGCM_restart, &
+ get_restart_file_name, &
+ get_f107_value, &
+ test_interpolate
! version controlled file description for error handling, do not edit
character(len=256), parameter :: source = &
- "$URL$"
-character(len=32 ), parameter :: revision = "$Revision$"
-character(len=128), parameter :: revdate = "$Date$"
+ '$URL$'
+character(len=32 ), parameter :: revision = '$Revision$'
+character(len=128), parameter :: revdate = '$Date$'
-!------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! namelist with default values
+! output_state_vector = .true. results in a "state-vector" netCDF file
+! output_state_vector = .false. results in a "prognostic-var" netCDF file
+
+! IMPORTANT: Change output file names in tiegcm.nml to match these names
+! i.e. OUTPUT='tiegcm_restart_p.nc'
+! SECOUT='tiegcm_s.nc'
+character(len=256) :: tiegcm_restart_file_name = 'tiegcm_restart_p.nc'
+character(len=256) :: tiegcm_secondary_file_name = 'tiegcm_s.nc'
+character(len=256) :: tiegcm_namelist_file_name = 'tiegcm.nml'
+logical :: output_state_vector = .false.
+integer :: debug = 0
+logical :: estimate_f10_7 = .false.
+integer :: assimilation_period_seconds = 3600
+
+integer, parameter :: max_num_variables = 30
+integer, parameter :: max_num_columns = 6
+character(len=NF90_MAX_NAME) :: variables(max_num_variables * max_num_columns) = ' '
+
+namelist /model_nml/ output_state_vector, tiegcm_restart_file_name, &
+ tiegcm_secondary_file_name, tiegcm_namelist_file_name, &
+ variables, debug, estimate_f10_7, assimilation_period_seconds
+
+!-------------------------------------------------------------------------------
+! Everything needed to describe a variable
+
+type progvartype
+ private
+ character(len=NF90_MAX_NAME) :: varname
+ character(len=NF90_MAX_NAME) :: long_name
+ character(len=NF90_MAX_NAME) :: units
+ character(len=obstypelength), dimension(NF90_MAX_VAR_DIMS) :: dimnames
+ integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens ! ntime, [nlev,] nlat, nlon
+ integer :: rank
+ integer :: varsize ! prod(dimlens(1:rank))
+ integer :: index1 ! location in dart state vector of first occurrence
+ integer :: indexN ! location in dart state vector of last occurrence
+ character(len=obstypelength) :: verticalvar
+ character(len=obstypelength) :: kind_string
+ integer :: dart_kind
+ integer :: xtype
+ integer :: rangeRestricted
+ real(r8) :: minvalue
+ real(r8) :: maxvalue
+ real(r8) :: missingR8
+ real(r4) :: missingR4
+ logical :: has_missing_value
+ logical :: update ! is this a state variable (updateable)
+ character(len=NF90_MAX_NAME) :: origin
+end type progvartype
+
+type(progvartype), dimension(max_num_variables) :: progvar
+integer :: nfields ! number of tiegcm variables in DART state
+
+!-------------------------------------------------------------------------------
! define model parameters
integer :: nilev, nlev, nlon, nlat
real(r8),dimension(:), allocatable :: lons, lats, levs, ilevs, plevs, pilevs
real(r8) :: TIEGCM_missing_value !! global attribute
real(r8) :: TIEGCM_reference_pressure
-integer :: time_step_seconds
-integer :: time_step_days
+integer :: time_step_seconds
+integer :: time_step_days
type(time_type) :: time_step
-
- ! IMPORTANT: Change output file names in
- ! tiegcm.nml to match with the names below
- ! i.e. OUTPUT='tiegcm_restart_p.nc'
- ! SECOUT='tiegcm_s.nc'
-character (len=19) :: restart_file = 'tiegcm_restart_p.nc'
-character (len=11) :: secondary_file = 'tiegcm_s.nc'
-character (len=10) :: namelist_file = 'tiegcm.nml'
- ! 3d TIEGCM variables are packed into
- ! DART state vector in the following order
-integer, parameter :: TYPE_local_ZG = 0
-integer, parameter :: TYPE_local_TN = 1
-integer, parameter :: TYPE_local_TN_NM = 2
-integer, parameter :: TYPE_local_O1 = 3
-integer, parameter :: TYPE_local_O1_NM = 4
-integer, parameter :: TYPE_local_O2 = 5
-integer, parameter :: TYPE_local_O2_NM = 6
-integer, parameter :: TYPE_local_UN = 7
-integer, parameter :: TYPE_local_UN_NM = 8
-integer, parameter :: TYPE_local_VN = 9
-integer, parameter :: TYPE_local_VN_NM = 10
-integer, parameter :: TYPE_local_NE = 11
+! Codes for interpreting the columns of the variable_table
+integer, parameter :: VT_VARNAMEINDX = 1 ! ... variable name
+integer, parameter :: VT_KINDINDX = 2 ! ... DART kind
+integer, parameter :: VT_MINVALINDX = 3 ! ... minimum value if any
+integer, parameter :: VT_MAXVALINDX = 4 ! ... maximum value if any
+integer, parameter :: VT_ORIGININDX = 5 ! ... file of origin
+integer, parameter :: VT_STATEINDX = 6 ! ... update (state) or not
-type model_type
- real(r8), pointer :: vars_3d(:,:,:,:)
- real(r8), pointer :: vars_1d(:)
- type(time_type) :: valid_time
-end type model_type
+character(len=obstypelength) :: variable_table(max_num_variables, max_num_columns)
-logical :: only_neutral_density = .true.
- ! .true. excludes UN VN NE (state_num_3d = 7)
- ! .false. includes UN VN NE (state_num_3d = 12)
-integer :: state_num_3d = 7
- ! -- interface levels --
- ! NE ZG
- ! -- midpoint levels --
- ! O1 O1_NM O2 O2_NM
- ! -- midpoint levels; top slot missing --
- ! TN TN_NM UN UN_NM VN VN_NM
+! include_vTEC = .true. vTEC must be calculated from other vars
+! include_vTEC = .false. just ignore vTEC altogether
-integer :: state_num_1d = 0
-logical :: estimate_parameter = .false.
- ! IMPORTANT: 1 D model parameters (e.g., F107) are read in from "tiegcm.nml"
- ! When "estimate_parameter = .true.", "state_num_1d" should be greater than or
- ! equal to 1 so that 1 D model parameters will be included in the state vector
- ! (note "estimate_parameter" option is still under
- ! development by Tomoko Matsuo as of June 24, 2011)
-
-integer :: model_size
-real(r8), allocatable :: ens_mean(:)
-
- ! FOR NOW OBS LOCATIONS ARE EXPECTED GIVEN IN HEIGHT [m],
- ! AND SO VERTICAL LOCALIZATION COORDINATE IS *always* HEIGHT
- ! (note that gravity adjusted geopotential height (ZG)
- ! read in from "tiegcm_s.nc")
-!integer :: vert_localization_coord = VERTISHEIGHT
+logical :: include_vTEC = .true.
+logical :: include_vTEC_in_state = .false.
-logical :: output_state_vector = .false.
- ! .true. results in a "state-vector" netCDF file
- ! .false. results in a "prognostic-var" netCDF file
-logical :: first_pert_call = .true.
-type(random_seq_type) :: random_seq
-!------------------------------------------------------------------
+! IMPORTANT: 1 D model parameters (e.g., F107) are read in from "tiegcm.nml"
+! (note "estimate_f10_7" option is still under
+! development by Tomoko Matsuo as of June 24, 2011)
-character(len = 129) :: msgstring, msgstring2, msgstring3
-logical, save :: module_initialized = .false.
+real(r8) :: f10_7
+type(time_type) :: state_time ! module-storage declaration of current model time
-namelist /model_nml/ output_state_vector, state_num_3d, state_num_1d
+integer :: model_size
+real(r8), allocatable :: ens_mean(:)
+! FOR NOW OBS LOCATIONS ARE EXPECTED GIVEN IN HEIGHT [m],
+! AND SO VERTICAL LOCALIZATION COORDINATE IS *always* HEIGHT
+! (note that gravity adjusted geopotential height (ZG)
+! read in from "tiegcm_s.nc" *WARNING* ZG is 'cm', DART is mks)
+
+real(r8), allocatable, dimension(:,:,:) :: ZG
+integer :: ivarZG
+
+logical :: first_pert_call = .true.
+type(random_seq_type) :: random_seq
+character(len=512) :: string1, string2, string3
+logical, save :: module_initialized = .false.
+
+! Codes for restricting the range of a variable
+integer, parameter :: BOUNDED_NONE = 0 ! ... unlimited range
+integer, parameter :: BOUNDED_BELOW = 1 ! ... minimum, but no maximum
+integer, parameter :: BOUNDED_ABOVE = 2 ! ... maximum, but no minimum
+integer, parameter :: BOUNDED_BOTH = 3 ! ... minimum and maximum
+
+interface prog_var_to_vector
+ module procedure var1d_to_vector
+ module procedure var2d_to_vector
+ module procedure var3d_to_vector
+ module procedure var4d_to_vector
+end interface
+
+interface apply_attributes
+ module procedure apply_attributes_1D
+ module procedure apply_attributes_2D
+ module procedure apply_attributes_3D
+ module procedure apply_attributes_4D
+end interface
+
+!===============================================================================
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.
@@ -175,11 +227,11 @@
! spherical harmonic weights, these would also be computed here.
! Can be a NULL INTERFACE for the simplest models.
- integer :: i
- integer :: iunit, io
+integer :: iunit, io
+real(r8) :: total_steps
if (module_initialized) return ! only need to do this once
-
+
! Print module information to log file and stdout.
call register_module(source, revision, revdate)
@@ -187,89 +239,149 @@
! 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")
+! 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)
+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)
+if (do_output()) then
+ write( * ,*)'static_init_model: debug level is ',debug
+ write(logfileunit,*)'static_init_model: debug level is ',debug
+endif
-! Reading in TIEGCM namelist input file (just for definition)
-call read_TIEGCM_namelist(namelist_file)
+! Read in TIEGCM namelist input file (just for definition)
+! Read in TIEGCM grid definition etc from TIEGCM restart file
+! Read in TIEGCM auxiliary variables from TIEGCM 'secondary' file
-! Compute overall model size
-model_size = nlon * nlat * nlev * state_num_3d + state_num_1d
+call read_TIEGCM_namelist(tiegcm_namelist_file_name)
+call read_TIEGCM_definition(tiegcm_restart_file_name)
+call read_TIEGCM_secondary(tiegcm_secondary_file_name)
-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
-if (do_output()) write(*,*) 'n1D = ', state_num_1d
-if (do_output()) write(*,*) 'model_size = ', model_size
+! error-check and convert namelist input to variable_table
+! and the 'progvar' database in the scope of the module
+call verify_variables(variables, nfields)
+! Compute overall model size
+
+model_size = progvar(nfields)%indexN
+
+if (do_output() .and. (debug > 0)) then
+ write(*,*) 'nlon = ', nlon
+ write(*,*) 'nlat = ', nlat
+ write(*,*) 'nlev = ', nlev
+ write(*,*) 'nilev = ', nilev
+ write(*,*) 'model_size = ', model_size
+ if (estimate_f10_7) &
+ write(*,*) 'estimating f10.7 ... init value ',f10_7
+endif
+
allocate (ens_mean(model_size))
+ens_mean = MISSING_R8
! 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)
+! Ensure assimilation_period is a multiple of the dynamical timestep
+! The time is communicated to TIEGCM through their "STOP" variable,
+! which is an array of length 3 corresponding to day-of-year, hour, minute
+! SO - there is some combination of 'STEP' and assimilation_period_seconds
+! that must be an integer number of minutes.
+time_step_seconds = time_step_seconds + time_step_days*86400
+
+if (assimilation_period_seconds < time_step_seconds) then
+ write(string1,*)'assimilation_period_seconds must be >= STEP'
+ write(string2,*)' input.nml: assimilation_period_seconds ',assimilation_period_seconds
+ write(string3,*)'tiegcm.nml: STEP ',time_step_seconds
+ call error_handler(E_ERR,'static_init_model',string1, &
+ source, revision, revdate, text2=string2,text3=string3)
+endif
+
+total_steps = real(assimilation_period_seconds,r8)/real(time_step_seconds,r8)
+
+if ( time_step_seconds*nint(total_steps) /= assimilation_period_seconds) then
+ write(string1,*)'assimilation_period_seconds must be an integer number of tiegcm "STEP"s'
+ write(string2,*)' input.nml: assimilation_period_seconds ',assimilation_period_seconds
+ write(string3,*)'tiegcm.nml: STEP ',time_step_seconds
+ call error_handler(E_ERR,'static_init_model',string1, &
+ source, revision, revdate, text2=string2,text3=string3)
+endif
+
+if ( mod(assimilation_period_seconds,60) /= 0 ) then
+ write(string1,*)'assimilation_period_seconds must be an integer number of tiegcm "STEP"s'
+ write(string2,*)'assimilation_period_seconds=',assimilation_period_seconds, &
+ ' STEP=',time_step_seconds
+ write(string3,*)'AND must be an integer number of minutes because of tiegcm "STOP"'
+ call error_handler(E_ERR,'static_init_model',string1, &
+ source, revision, revdate, text2=string2,text3=string3)
+endif
+
+time_step = set_time(assimilation_period_seconds, 0)
+
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
+! 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.
+! If this option is not to be used this can be a NULL INTERFACE.
real(r8), intent(out) :: x(:)
if ( .not. module_initialized ) call static_init_model
+write(string1,*) 'no good way to specify initial conditions'
+call error_handler(E_ERR,'init_conditions',string1,source,revision,revdate)
+
+x(:) = MISSING_R8 ! just to silence compiler messages
+
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
+! 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
+! 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
+! a separate model-specific executable), this can be a
! NULL INTERFACE.
real(r8), intent(inout) :: x(:)
type(time_type), intent(in) :: time
+write(string1,*) 'TIEGCM cannot be advanced as a subroutine from within DART.'
+write(string2,*) 'check your input.nml setting of "async" and try again.'
+call error_handler(E_ERR,'adv_1step',string1,source,revision,revdate)
+
+! just to silence compiler messages
+
+x(:) = MISSING_R8
+call print_time(time)
+
end subroutine adv_1step
+!-------------------------------------------------------------------------------
+
function get_model_size()
-!------------------------------------------------------------------
-!
! Returns the size of the model as an integer. Required for all
! applications.
@@ -282,99 +394,137 @@
end function get_model_size
+!-------------------------------------------------------------------------------
+
subroutine init_time(time)
-!------------------------------------------------------------------
-!
-! Companion interface to init_conditions. Returns a time that is somehow
+! Companion interface to init_conditions. Returns a time that is somehow
! appropriate for starting up a long integration of the model.
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list