[Dart-dev] [6653] DART/trunk: The portion of the GITM incorporation that required
nancy at ucar.edu
nancy at ucar.edu
Fri Dec 6 14:32:29 MST 2013
Revision: 6653
Author: thoar
Date: 2013-12-06 14:32:28 -0700 (Fri, 06 Dec 2013)
Log Message:
-----------
The portion of the GITM incorporation that required
changes to existing routines, along with some mundane
whitespace and comment changes found along the way.
Modified Paths:
--------------
DART/trunk/location/column/location_mod.f90
DART/trunk/models/gitm/GITM2/src/ModSize.f90
DART/trunk/models/gitm/dart_to_gitm.f90
DART/trunk/models/gitm/gitm_to_dart.f90
DART/trunk/models/gitm/model_mod.f90
DART/trunk/models/gitm/model_mod.html
DART/trunk/models/gitm/shell_scripts/advance_model.csh
DART/trunk/models/gitm/work/input.nml
DART/trunk/models/gitm/work/mkmf_create_fixed_network_seq
DART/trunk/models/gitm/work/mkmf_create_obs_sequence
DART/trunk/models/gitm/work/mkmf_dart_to_gitm
DART/trunk/models/gitm/work/mkmf_filter
DART/trunk/models/gitm/work/mkmf_gitm_to_dart
DART/trunk/models/gitm/work/mkmf_model_mod_check
DART/trunk/models/gitm/work/mkmf_obs_diag
DART/trunk/models/gitm/work/mkmf_obs_seq_to_netcdf
DART/trunk/models/gitm/work/mkmf_obs_sequence_tool
DART/trunk/models/gitm/work/mkmf_perfect_model_obs
DART/trunk/models/gitm/work/mkmf_preprocess
DART/trunk/models/gitm/work/mkmf_restart_file_tool
DART/trunk/models/gitm/work/mkmf_wakeup_filter
DART/trunk/models/gitm/work/path_names_create_fixed_network_seq
DART/trunk/models/gitm/work/path_names_create_obs_sequence
DART/trunk/models/gitm/work/path_names_dart_to_gitm
DART/trunk/models/gitm/work/path_names_filter
DART/trunk/models/gitm/work/path_names_gitm_to_dart
DART/trunk/models/gitm/work/path_names_model_mod_check
DART/trunk/models/gitm/work/path_names_obs_diag
DART/trunk/models/gitm/work/path_names_obs_seq_to_netcdf
DART/trunk/models/gitm/work/path_names_obs_sequence_tool
DART/trunk/models/gitm/work/path_names_perfect_model_obs
DART/trunk/models/gitm/work/path_names_restart_file_tool
DART/trunk/models/gitm/work/quickbuild.csh
DART/trunk/models/tiegcm/model_mod.f90
DART/trunk/obs_def/obs_def_radar_mod.f90
DART/trunk/obs_def/obs_def_reanalysis_bufr_mod.f90
DART/trunk/obs_def/obs_def_upper_atm_mod.f90
DART/trunk/obs_sequence/obs_sequence_tool.nml
-------------- next part --------------
Modified: DART/trunk/location/column/location_mod.f90
===================================================================
--- DART/trunk/location/column/location_mod.f90 2013-12-06 20:14:49 UTC (rev 6652)
+++ DART/trunk/location/column/location_mod.f90 2013-12-06 21:32:28 UTC (rev 6653)
@@ -399,7 +399,7 @@
read(*, *) location%vloc
location%vloc = 100.0_r8 * location%vloc
else if(location%which_vert == VERTISHEIGHT ) then
- write(*, *) 'Vertical coordinate height (in gpm)'
+ write(*, *) 'Vertical coordinate height (in meters)'
read(*, *) location%vloc
else if(location%which_vert == VERTISSURFACE ) then
write(*, *) 'Vertical coordinate surface height'
Modified: DART/trunk/models/gitm/GITM2/src/ModSize.f90
===================================================================
--- DART/trunk/models/gitm/GITM2/src/ModSize.f90 2013-12-06 20:14:49 UTC (rev 6652)
+++ DART/trunk/models/gitm/GITM2/src/ModSize.f90 2013-12-06 21:32:28 UTC (rev 6653)
@@ -1,7 +1,3 @@
-! This code may (or may not) be part of the GITM distribution,
-! So it is not protected by the DART copyright agreement.
-!
-! DART $Id$
module ModSizeGitm
@@ -14,9 +10,3 @@
integer :: nBlocks
end module ModSizeGitm
-
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
Modified: DART/trunk/models/gitm/dart_to_gitm.f90
===================================================================
--- DART/trunk/models/gitm/dart_to_gitm.f90 2013-12-06 20:14:49 UTC (rev 6652)
+++ DART/trunk/models/gitm/dart_to_gitm.f90 2013-12-06 21:32:28 UTC (rev 6653)
@@ -7,30 +7,32 @@
program dart_to_gitm
!----------------------------------------------------------------------
-! purpose: interface between DART and the gitm model
+! purpose: interface between DART and the GITM model
!
! method: Read DART state vector and overwrite values in a gitm restart file.
! If the DART state vector has an 'advance_to_time' present, a
-! file called gitm_in.DART is created with a time_manager_nml namelist
-! appropriate to advance gitm to the requested time.
+! file called DART_GITM_time_control.txt is created with
+! information appropriate to advance gitm to the requested time.
!
-! The dart_to_gitm_nml namelist setting for advance_time_present
+! The dart_to_gitm_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, finalize_utilities, &
find_namelist_in_file, check_namelist_read, &
logfileunit, open_file, close_file
+
+use model_mod, only : static_init_model, get_model_size, &
+ get_gitm_restart_dirname, statevector_to_restart_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, statevector_to_restart_file, &
- get_model_size, get_base_time, get_gitm_restart_dirname
implicit none
@@ -40,49 +42,52 @@
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate = "$Date$"
-!------------------------------------------------------------------
-! The namelist variables
-!------------------------------------------------------------------
+!-----------------------------------------------------------------------
+! namelist parameters with default values.
+!-----------------------------------------------------------------------
character (len = 128) :: dart_to_gitm_input_file = 'dart_restart'
logical :: advance_time_present = .false.
-character(len=256) :: gitm_restart_dirname = 'gitm_restartdir'
namelist /dart_to_gitm_nml/ dart_to_gitm_input_file, &
- advance_time_present, &
- gitm_restart_dirname
+ advance_time_present
!----------------------------------------------------------------------
+! global storage
+!----------------------------------------------------------------------
-integer :: iunit, io, x_size, diff1, diff2
-type(time_type) :: model_time, adv_to_time, base_time
+integer :: iunit, io, x_size
+type(time_type) :: model_time, adv_to_time
real(r8), allocatable :: statevector(:)
-logical :: verbose = .TRUE.
+character(len=256) :: gitm_restart_dirname = 'none'
-!----------------------------------------------------------------------
+!======================================================================
call initialize_utilities(progname='dart_to_gitm')
!----------------------------------------------------------------------
+! Read the namelist.
+!----------------------------------------------------------------------
+
+call find_namelist_in_file("input.nml", "dart_to_gitm_nml", iunit)
+read(iunit, nml = dart_to_gitm_nml, iostat = io)
+call check_namelist_read(iunit, io, "dart_to_gitm_nml")
+
+!----------------------------------------------------------------------
! Call model_mod:static_init_model() which reads the gitm namelists
! to set grid sizes, etc.
!----------------------------------------------------------------------
call static_init_model()
+call get_gitm_restart_dirname(gitm_restart_dirname)
+write(*,*)
+write(*,*) 'dart_to_gitm: converting DART file ', "'"//trim(dart_to_gitm_input_file)//"'"
+write(*,*) 'to gitm restart files in directory ', "'"//trim(gitm_restart_dirname)//"'"
+
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_gitm_nml", iunit)
-read(iunit, nml = dart_to_gitm_nml, iostat = io)
-call check_namelist_read(iunit, io, "dart_to_gitm_nml")
-
-write(*,*)
-write(*,*) 'dart_to_gitm: converting DART file ', "'"//trim(dart_to_gitm_input_file)//"'"
-write(*,*) 'to gitm restart files in directory ', "'"//trim(gitm_restart_dirname)//"'"
-
!----------------------------------------------------------------------
! Reads the valid time, the state, and the target time.
!----------------------------------------------------------------------
@@ -96,14 +101,12 @@
endif
call close_restart(iunit)
-print *, 'read state vector'
!----------------------------------------------------------------------
! update the current gitm 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 statevector_to_restart_file(statevector, gitm_restart_dirname, model_time)
if ( advance_time_present ) then
@@ -162,7 +165,7 @@
! the end time comes first.
call get_date(adv_to_time,iyear,imonth,iday,ihour,imin,isec)
-write(iunit,'(''#TIMEEND'')')
+write(iunit,'(''#TIMEEND'')')
write(iunit,'(i4.4,10x,''year'' )')iyear
write(iunit,'(i2.2,12x,''month'' )')imonth
write(iunit,'(i2.2,12x,''day'' )')iday
@@ -172,7 +175,7 @@
write(iunit,*)
call get_date(model_time,iyear,imonth,iday,ihour,imin,isec)
-write(iunit,'(''#TIMESTART'')')
+write(iunit,'(''#TIMESTART'')')
write(iunit,'(i4.4,10x,''year'' )')iyear
write(iunit,'(i2.2,12x,''month'' )')imonth
write(iunit,'(i2.2,12x,''day'' )')iday
Modified: DART/trunk/models/gitm/gitm_to_dart.f90
===================================================================
--- DART/trunk/models/gitm/gitm_to_dart.f90 2013-12-06 20:14:49 UTC (rev 6652)
+++ DART/trunk/models/gitm/gitm_to_dart.f90 2013-12-06 21:32:28 UTC (rev 6653)
@@ -7,26 +7,28 @@
program gitm_to_dart
!----------------------------------------------------------------------
-! purpose: interface between gitm and DART
+! purpose: interface between the GITM model and DART
!
! method: Read gitm "restart" files of model state
! Reform fields into a DART state vector (control vector).
! Write out state vector in "proprietary" format for DART.
! The output is a "DART restart file" format.
-!
+!
! USAGE: The gitm dirname is read from the gitm_in namelist
! <edit gitm_to_dart_output_file in input.nml:gitm_to_dart_nml>
! gitm_to_dart
-!
-! author: Tim Hoar 6/24/09
!----------------------------------------------------------------------
use types_mod, only : r8
+
use utilities_mod, only : initialize_utilities, finalize_utilities, &
find_namelist_in_file, check_namelist_read
-use model_mod, only : get_model_size, restart_file_to_statevector, &
- get_gitm_restart_dirname
-use assim_model_mod, only : awrite_state_restart, open_restart_write, close_restart
+
+use model_mod, only : static_init_model, get_model_size, &
+ get_gitm_restart_dirname, restart_file_to_statevector
+
+use assim_model_mod, only : open_restart_write, awrite_state_restart, close_restart
+
use time_manager_mod, only : time_type, print_time, print_date
implicit none
@@ -42,50 +44,48 @@
!-----------------------------------------------------------------------
character(len=128) :: gitm_to_dart_output_file = 'dart_ics'
-character(len=256) :: gitm_restart_dirname = 'gitm_restartdir'
-namelist /gitm_to_dart_nml/ &
- gitm_to_dart_output_file, &
- gitm_restart_dirname
+namelist /gitm_to_dart_nml/ gitm_to_dart_output_file
!----------------------------------------------------------------------
! global storage
!----------------------------------------------------------------------
-logical :: verbose = .FALSE.
-integer :: io, iunit, x_size
+integer :: iunit, io, x_size
type(time_type) :: model_time
real(r8), allocatable :: statevector(:)
+character(len=256) :: gitm_restart_dirname = 'none'
!======================================================================
-!call initialize_utilities(progname='gitm_to_dart', output_flag=verbose)
call initialize_utilities(progname='gitm_to_dart')
!----------------------------------------------------------------------
-! Read the namelist to get the output dirname.
+! Read the namelist.
!----------------------------------------------------------------------
call find_namelist_in_file("input.nml", "gitm_to_dart_nml", iunit)
read(iunit, nml = gitm_to_dart_nml, iostat = io)
call check_namelist_read(iunit, io, "gitm_to_dart_nml") ! closes, too.
+!----------------------------------------------------------------------
+! Call model_mod:static_init_model() which reads the gitm namelists
+! to set grid sizes, etc.
+!----------------------------------------------------------------------
+
+call static_init_model()
+call get_gitm_restart_dirname(gitm_restart_dirname)
+
write(*,*)
write(*,*) 'gitm_to_dart: converting gitm restart files in directory ', &
- "'"//trim(gitm_restart_dirname)//"'"
+ "'"//trim(gitm_restart_dirname)//"'"
write(*,*) ' to DART file ', "'"//trim(gitm_to_dart_output_file)//"'"
-!----------------------------------------------------------------------
-! get to work
-!----------------------------------------------------------------------
-
x_size = get_model_size()
allocate(statevector(x_size))
-call get_gitm_restart_dirname( gitm_restart_dirname )
+call restart_file_to_statevector(gitm_restart_dirname, statevector, model_time)
-call restart_file_to_statevector(gitm_restart_dirname, statevector, model_time)
-
iunit = open_restart_write(gitm_to_dart_output_file)
call awrite_state_restart(model_time, statevector, iunit)
Modified: DART/trunk/models/gitm/model_mod.f90
===================================================================
--- DART/trunk/models/gitm/model_mod.f90 2013-12-06 20:14:49 UTC (rev 6652)
+++ DART/trunk/models/gitm/model_mod.f90 2013-12-06 21:32:28 UTC (rev 6653)
@@ -37,7 +37,8 @@
use obs_kind_mod, only : paramname_length, &
get_raw_obs_kind_index, &
- get_raw_obs_kind_name
+ get_raw_obs_kind_name, &
+ KIND_GEOPOTENTIAL_HEIGHT
use mpi_utilities_mod, only: my_task_id
@@ -75,12 +76,9 @@
! generally useful routines for various support purposes.
! the interfaces here can be changed as appropriate.
-public :: get_gridsize, &
- get_grid_val, &
- restart_file_to_statevector, &
+public :: restart_file_to_statevector, &
statevector_to_restart_file, &
get_gitm_restart_dirname, &
- get_base_time, &
get_state_time
! version controlled file description for error handling, do not edit
@@ -96,26 +94,8 @@
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) :: gitm_restart_dirname = 'gitm_restartdir'
-
-namelist /model_nml/ &
- gitm_restart_dirname, &
- output_state_vector, &
- assimilation_period_days, & ! for now, this is the timestep
- assimilation_period_seconds, &
- model_perturbation_amplitude,&
- calendar, &
- debug
-
!------------------------------------------------------------------
+! things which can/should be in the model_nml
!
! The DART state vector may consist of things like:
!
@@ -134,7 +114,7 @@
! QH long_name = "GRAUPEL MIXING RATIO" float QH(TIME, ALT, LAT, LON)
!
! The variables in the gitm restart file that are used to create the
-! DART state vector are specified in the input.nml:gitm_vars_nml namelist.
+! DART state vector are specified in input.nml:model_nml: gitm_state_variables
!
!------------------------------------------------------------------
@@ -143,8 +123,24 @@
character(len=NF90_MAX_NAME) :: gitm_state_variables(max_state_variables * num_state_table_columns ) = ' '
character(len=NF90_MAX_NAME) :: variable_table(max_state_variables, num_state_table_columns )
-namelist /gitm_vars_nml/ gitm_state_variables
+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) :: gitm_restart_dirname = 'gitm_restartdir'
+namelist /model_nml/ &
+ gitm_restart_dirname, &
+ output_state_vector, &
+ assimilation_period_days, & ! for now, this is the timestep
+ assimilation_period_seconds, &
+ model_perturbation_amplitude,&
+ calendar, &
+ debug, &
+ gitm_state_variables
+
integer :: nfields
! Everything needed to describe a variable
@@ -219,11 +215,6 @@
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
@@ -231,6 +222,7 @@
contains
+
!==================================================================
! All the REQUIRED interfaces come first - just by convention.
!==================================================================
@@ -238,7 +230,7 @@
function get_model_size()
!------------------------------------------------------------------
-! Done - TJH.
+
! Returns the size of the model as an integer.
! Required for all applications.
@@ -251,10 +243,11 @@
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
@@ -280,13 +273,16 @@
write(string1,*) 'cannot run this model with async=0'
call error_handler(E_ERR,'adv_1step',string1,source,revision,revdate)
+x = MISSING_R8 ! just to silence compiler warnings
+
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_
@@ -333,10 +329,197 @@
end subroutine get_state_meta_data
+!==================================================================
-function get_model_time_step()
+
+subroutine model_interpolate(x, location, obs_type, interp_val, istatus)
!------------------------------------------------------------------
+! PURPOSE:
!
+! For a given lat, lon, and height, interpolate the correct state value
+! to that location for the filter from the gitm state vectors
+!
+! Variables needed to be stored in the MODEL_MODULE data structure
+!
+! LON = 1D array storing the local grid center coords (degrees)
+! LAT = 1D array storing the local grid center coords (degrees)
+! ALT = 1D array storing the local grid center coords (meters)
+!
+! ERROR codes:
+!
+! ISTATUS = 99: general error in case something terrible goes wrong...
+! ISTATUS = 15: dont know what to do with vertical coord supplied
+! ISTATUS = 16: longitude illegal
+! ISTATUS = 17: latitude illegal
+! ISTATUS = 18: altitude illegal
+! ISTATUS = 20: asking to interpolate an unknown obs kind
+
+! Passed variables
+
+real(r8), intent(in) :: x(:)
+type(location_type), intent(in) :: location
+integer, intent(in) :: obs_type
+real(r8), intent(out) :: interp_val
+integer, intent(out) :: istatus
+
+! Local storage
+
+real(r8) :: loc_array(3), llon, llat, lvert, lon_fract, lat_fract, alt_fract
+integer :: base_offset, end_offset, blon(2), blat(2), balt(2), i, j, k, ier, nhgt
+real(r8) :: cube(2, 2, 2), square(2, 2), line(2)
+
+if ( .not. module_initialized ) call static_init_model
+
+! Assume failure. Set return val to missing, then the code can
+! just set istatus to something indicating why it failed, and return.
+! If the interpolation is good, the interp_val will be set to the
+! good value, and the last line here sets istatus to 0.
+! make any error codes set here be in the 10s
+
+interp_val = MISSING_R8 ! the DART bad value flag
+istatus = 99 ! unknown error
+
+! Get the individual locations values
+
+loc_array = get_location(location)
+llon = loc_array(1)
+llat = loc_array(2)
+lvert = loc_array(3)
+
+! IF (debug > 2) print *, 'requesting interpolation at ', llon, llat, lvert
+print *, 'requesting interpolation at ', llon, llat, lvert
+
+! Only height and level for vertical location type is supported at this point
+if (.not. vert_is_height(location) .and. .not. vert_is_level(location)) THEN
+ istatus = 15
+ return
+endif
+
+! Find the start and end offsets for this field in the state vector x(:)
+
+! FIXME: this fails if you ask for a kind that doesn't exist in the state vector.
+! in most cases you just want to set a bad istatus and return without stopping
+! the entire assimilation. if the error calls are commented out of the two
+! index range subroutines, the offsets will come back as 0. since the
+! return values have already been set, just give it a more specific error
+! code and return here.
+
+if (obs_type == KIND_GEOPOTENTIAL_HEIGHT ) then
+ ! ok to continue. offsets unused in this case, but
+ ! set them to something > 0 to indicate its ok.
+ base_offset = 1
+ end_offset = 1
+else
+ call get_index_range(obs_type, base_offset, end_offset)
+endif
+
+if (debug > 2) print *, 'base offset now ', base_offset
+
+! fail if this kind isn't in the state vector.
+if (base_offset <= 0) then
+ istatus = 20
+ return
+endif
+
+! Need to find bounding center indices and fractional offsets for lat, lon, alt
+call find_lon_bounds(llon, blon(1), blon(2), lon_fract, ier)
+if(ier /= 0) then
+ istatus = 16
+ return
+endif
+
+call find_lat_or_alt_bounds(llat, NgridLat, LAT, blat(1), blat(2), lat_fract, ier)
+if(ier /= 0) then
+ istatus = 17
+ return
+endif
+
+if (vert_is_height(location)) then ! call the height interpolation routine.
+
+ call find_lat_or_alt_bounds(lvert, NgridAlt, ALT, balt(1), balt(2), alt_fract, ier)
+
+else if (vert_is_level(location)) then ! set the levels and fraction.
+
+ nhgt = int(lvert)
+ if (nhgt < 1 .or. nhgt > NgridAlt) then
+ istatus = 18
+ return
+ endif
+
+ ! if we are below the top level, set the lower bound to the integer part of the
+ ! level and set the fraction between it and the next level to the fractional part.
+ ! if we are asking for the top level, set the upper bound to the requested level
+ ! and set the fraction to 1.
+ if (nhgt < NGridAlt) then
+ balt(1) = nhgt
+ balt(2) = nhgt + 1
+ alt_fract = lvert - real(nhgt,r8)
+ else
+ balt(1) = nhgt - 1
+ balt(2) = nhgt
+ alt_fract = 1.0_r8
+ endif
+ ier = 0
+else
+ ! shouldn't happen
+ istatus = 99
+ return
+endif
+if(ier /= 0) then
+ istatus = 18
+ return
+endif
+
+! if we're asking about height, we have the alt arrays directly.
+if (obs_type == KIND_GEOPOTENTIAL_HEIGHT) then
+
+ ! Interpolate to the given altitude - lat/lon doesn't matter here.
+ interp_val = (1 - alt_fract) * ALT(balt(1)) + alt_fract * ALT(balt(2))
+
+ istatus = 0
+ return
+endif
+
+! for the rest of the state vector contents - find the offset to the start
+! of the requested kind, and do a tri-linear interpolation inside the enclosing
+! cube from the grid.
+
+! Get the grid values for the first
+do i = 1, 2
+ do j = 1, 2
+ do k = 1, 2
+ cube(i, j, k) = get_grid_value(base_offset, blon(i), blat(j), balt(k), x)
+ end do
+ end do
+end do
+
+! Interpolate to the given altitude
+do i = 1, 2
+ do j = 1, 2
+ square(i, j) = (1 - alt_fract) * cube(i, j, 1) + alt_fract * cube(i, j, 2)
+ end do
+end do
+
+! Interpolate to the given latitude
+do i = 1, 2
+ line(i) = (1 - lat_fract) * square(i, 1) + lat_fract * square(i, 2)
+end do
+
+! Interpolate to the given longitude
+interp_val = (1 - lon_fract) * line(1) + lon_fract * line(2)
+
+! All good.
+istatus = 0
+
+return
+end subroutine model_interpolate
+
+
+!==================================================================
+
+
+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.
@@ -350,10 +533,11 @@
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
@@ -361,14 +545,11 @@
! 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
+character(len=NF90_MAX_NAME) :: varname
+character(len=paramname_length) :: kind_string
+integer :: varsize
+integer :: iunit, io, ivar, index1, indexN
integer :: ss, dd
-integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID
-logical :: shapeok
if ( module_initialized ) return ! only need to do this once.
@@ -410,12 +591,6 @@
write(*,*)'nSpeciesAll is ',nSpeciesAll
endif
-! Read the gitm variable list to populate DART state vector
-! Once parsed, the values will be recorded for posterity
-call find_namelist_in_file('input.nml', 'gitm_vars_nml', iunit)
-read(iunit, nml = gitm_vars_nml, iostat = io)
-call check_namelist_read(iunit, io, 'gitm_vars_nml')
-
!---------------------------------------------------------------
! Set the time step ... causes gitm namelists to be read.
! Ensures model_timestep is multiple of 'dynamics_timestep'
@@ -460,8 +635,7 @@
!
! Record the extent of the data type in the state vector.
-call verify_state_variables( gitm_state_variables, ncid, gitm_restart_dirname, &
- nfields, variable_table )
+call verify_state_variables( gitm_state_variables, nfields, variable_table)
index1 = 1;
indexN = 0;
@@ -480,8 +654,12 @@
! progvar(ivar)%numdims
! progvar(ivar)%dimlens
- call decode_gitm_indices( varname, progvar(ivar)%gitm_varname, progvar(ivar)%gitm_dim, &
- progvar(ivar)%gitm_index, progvar(ivar)%long_name, &
+ ! This routine also checks to make sure user specified accurate GITM variables
+ call decode_gitm_indices( varname, &
+ progvar(ivar)%gitm_varname, &
+ progvar(ivar)%gitm_dim, &
+ progvar(ivar)%gitm_index, &
+ progvar(ivar)%long_name, &
progvar(ivar)%units)
if (progvar(ivar)%varname == 'f107') then ! if we are dealing with f107
@@ -553,10 +731,11 @@
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.
@@ -565,10 +744,11 @@
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
@@ -590,10 +770,11 @@
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
@@ -614,7 +795,9 @@
end subroutine init_conditions
+!==================================================================
+
function nc_write_model_atts( ncFileID ) result (ierr)
!------------------------------------------------------------------
! TJH -- Writes the model-specific attributes to a netCDF file.
@@ -656,10 +839,10 @@
!----------------------------------------------------------------------
! for the dimensions and coordinate variables
-integer :: NLONDimID
-integer :: NLATDimID
-integer :: NALTDimID
-integer :: NWLDimID !alex: number of WaveLengths (for EUV spectrum, 10.7cm is 1 wl)
+integer :: LON_dim_id
+integer :: LAT_dim_id
+integer :: ALT_dim_id
+integer :: WL_dim_id !alex: number of WaveLengths (for EUV spectrum, 10.7cm is 1 wl)
! for the prognostic variables
integer :: ivar, VarID
@@ -833,16 +1016,16 @@
!----------------------------------------------------------------------------
call nc_check(nf90_def_dim(ncid=ncFileID, name='LON', &
- len = NgridLon, dimid = NLONDimID),'nc_write_model_atts', 'LON def_dim '//trim(filename))
+ len = NgridLon, dimid = LON_dim_id),'nc_write_model_atts', 'LON def_dim '//trim(filename))
call nc_check(nf90_def_dim(ncid=ncFileID, name='LAT', &
- len = NgridLat, dimid = NLATDimID),'nc_write_model_atts', 'LAT def_dim '//trim(filename))
+ len = NgridLat, dimid = LAT_dim_id),'nc_write_model_atts', 'LAT def_dim '//trim(filename))
call nc_check(nf90_def_dim(ncid=ncFileID, name='ALT', &
- len = NgridAlt, dimid = NALTDimID),'nc_write_model_atts', 'ALT def_dim '//trim(filename))
+ len = NgridAlt, dimid = ALT_dim_id),'nc_write_model_atts', 'ALT def_dim '//trim(filename))
call nc_check(nf90_def_dim(ncid=ncFileID, name='WL', &
- len = 1, dimid = NWLDimID),'nc_write_model_atts', 'WL def_dim '//trim(filename)) !Alex
+ len = 1, dimid = WL_dim_id),'nc_write_model_atts', 'WL def_dim '//trim(filename)) !Alex
!----------------------------------------------------------------------------
! Create the (empty) Coordinate Variables and the Attributes
@@ -850,7 +1033,7 @@
! Grid Longitudes
call nc_check(nf90_def_var(ncFileID,name='LON', xtype=nf90_real, &
- dimids=NLONDimID, varid=VarID),&
+ dimids=LON_dim_id, varid=VarID),&
'nc_write_model_atts', 'LON def_var '//trim(filename))
call nc_check(nf90_put_att(ncFileID, VarID, 'type', 'x1d'), &
'nc_write_model_atts', 'LON type '//trim(filename))
@@ -865,7 +1048,7 @@
! Grid Latitudes
call nc_check(nf90_def_var(ncFileID,name='LAT', xtype=nf90_real, &
- dimids=NLATDimID, varid=VarID),&
+ dimids=LAT_dim_id, varid=VarID),&
'nc_write_model_atts', 'LAT def_var '//trim(filename))
call nc_check(nf90_put_att(ncFileID, VarID, 'type', 'y1d'), &
'nc_write_model_atts', 'LAT type '//trim(filename))
@@ -880,7 +1063,7 @@
! Grid Altitudes
call nc_check(nf90_def_var(ncFileID,name='ALT',xtype=nf90_real, &
- dimids=NALTDimID,varid=VarID), &
+ dimids=ALT_dim_id,varid=VarID), &
'nc_write_model_atts', 'ALT def_var '//trim(filename))
call nc_check(nf90_put_att(ncFileID, VarID, 'type', 'z1d'), &
'nc_write_model_atts', 'ALT type '//trim(filename))
@@ -895,7 +1078,7 @@
! Grid wavelengths
call nc_check(nf90_def_var(ncFileID,name='WL', xtype=nf90_real, &
- dimids=NWLDimID, varid=VarID),&
+ dimids=WL_dim_id, varid=VarID),&
'nc_write_model_atts', 'WL def_var '//trim(filename))
call nc_check(nf90_put_att(ncFileID, VarID, 'type', 'x1d'), &
'nc_write_model_atts', 'WL type '//trim(filename))
@@ -919,8 +1102,8 @@
! match shape of the variable to the dimension IDs
- call define_var_dims(progvar(ivar), myndims, mydimids, MemberDimID, unlimitedDimID, &
- NLONDimID, NLATDimID, NALTDimID, NgridLon, NgridLat, NgridAlt, NWLDimID) !Alex
+ call define_var_dims(progvar(ivar), myndims, mydimids, MemberDimID, unlimitedDimID,&
+ LON_dim_id, LAT_dim_id, ALT_dim_id, WL_dim_id)
! define the variable and set the attributes
@@ -989,10 +1172,11 @@
end function nc_write_model_atts
+!==================================================================
+
function nc_write_model_vars( ncFileID, state_vec, copyindex, timeindex ) result (ierr)
!------------------------------------------------------------------
-!
! TJH 29 Aug 2011 -- all errors are fatal, so the
! return code is always '0 == normal', since the fatal errors stop execution.
!
@@ -1204,10 +1388,11 @@
end function nc_write_model_vars
+!==================================================================
+
subroutine pert_model_state(state, pert_state, interf_provided)
!------------------------------------------------------------------
-!
! Perturbs a model state for generating initial ensembles.
! The perturbed state is returned in pert_state.
! A model may choose to provide a NULL INTERFACE by returning
@@ -1244,15 +1429,16 @@
end subroutine pert_model_state
+!==================================================================
+
subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, &
obs_loc, obs_kind, num_close, close_ind, dist)
!------------------------------------------------------------------
-
! Given a DART location (referred to as "base") and a set of candidate
! locations & kinds (obs, obs_kind), returns the subset close to the
! "base", their indices, and their distances to the "base" ...
-
+!
! For vertical distance computations, general philosophy is to convert all
! vertical coordinates to a common coordinate. This coordinate type is defined
! in the namelist with the variable "vert_localization_coord".
@@ -1273,17 +1459,16 @@
! Initialize variables to missing status
-num_close = 0
-close_ind = -99
-dist = 1.0e9 !something big and positive (far away)
-istatus1 = 0
-istatus2 = 0
-is_in_obs_kind = 0
+num_close = 0
+close_ind = -99
+dist = 1.0e9_r8 !something big and positive (far away)
+istatus1 = 0
+istatus2 = 0
+is_in_obs_kind = 0
is_in_close_ind = 0
-f107_ind = -37 !a bad index, hopefully out of bounds of obs_kind
+f107_ind = -37 !a bad index, hopefully out of bounds of obs_kind
-
! Convert base_obs vertical coordinate to requested vertical coordinate if necessary
base_array = get_location(base_obs_loc)
@@ -1302,8 +1487,8 @@
! Loop over potentially close subset of obs priors or state variables
! This way, we are decreasing the number of distance computations that will follow.
- ! This is a horizontal-distance operation and we don't need to have the relevant vertical
- ! coordinate information yet (for obs_loc).
+ ! This is a horizontal-distance operation and we don't need to have the relevant
+ ! vertical coordinate information yet (for obs_loc).
call loc_get_close_obs(gc, base_obs_loc, base_obs_kind, obs_loc, obs_kind, &
num_close, close_ind)
@@ -1337,9 +1522,10 @@
local_obs_loc = obs_loc(t_ind)
local_obs_which = nint(query_location(local_obs_loc))
- ! Convert local_obs vertical coordinate to requested vertical coordinate if necessary.
- ! This should only be necessary for obs priors, as state location information already
- ! contains the correct vertical coordinate (filter_assim's call to get_state_meta_data).
+ ! Convert vertical coordinate to requested vertical coordinate if necessary.
+ ! Should only be necessary for obs priors, as state location information already
+ ! contains the correct vertical coordinate
+ ! (filter_assim's call to get_state_meta_data).
if (.not. horiz_dist_only) then
!fixme if (local_obs_which /= wrf%dom(1)%vert_coord) then
!fixme call vert_interpolate(ens_mean, local_obs_loc, obs_kind(t_ind), istatus2)
@@ -1348,15 +1534,17 @@
!fixme endif
endif
- ! Compute distance - set distance to a very large value if vert coordinate is missing
- ! or vert_interpolate returned error (istatus2=1)
+ ! Compute distance - set distance to a very large value if vert coordinate is
+ ! missing or vert_interpolate returned error (istatus2=1)
local_obs_array = get_location(local_obs_loc)
- if (( (.not. horiz_dist_only) .and. &
+ ! TJH FIXME ... the pbs_file script actually modifies the value of dist in
+ ! TJH FIXME ... this file and recompiles. NOT APPROPRIATE.
+ if (( (.not. horiz_dist_only) .and. &
(local_obs_array(3) == MISSING_R8)) .or. &
(istatus2 == 1) ) then
- dist(k) = 1.0e9
+ dist(k) = 1.0e9_r8
else
if (close_ind(k) .eq. f107_ind) then !check if we came across the parameter
dist(k) = 0 !changed by pbs_file script
@@ -1370,7 +1558,9 @@
end subroutine get_close_obs
+!==================================================================
+
subroutine ens_mean_for_model(filter_ens_mean)
!------------------------------------------------------------------
! If needed by the model interface, this is the current mean
@@ -1385,84 +1575,50 @@
end subroutine ens_mean_for_model
-
!==================================================================
! The remaining PUBLIC interfaces come next
!==================================================================
-
-subroutine get_gridsize(num_LON, num_LAT, num_ALT )
- integer, intent(out) :: num_LON, num_LAT, num_ALT
+subroutine restart_file_to_statevector(dirname, state_vector, model_time)
!------------------------------------------------------------------
-! public utility routine.
-
-if ( .not. module_initialized ) call static_init_model
-
- num_LON = NgridLon
- num_LAT = NgridLat
- num_ALT = NgridAlt
-
-end subroutine get_gridsize
-
-subroutine get_grid_val( lon_a, lat_a, alt_a )
- real(r8), dimension(:), intent(out) :: lon_a, lat_a, alt_a
-!------------------------------------------------------------------
-! public utility routine.
-
-if ( .not. module_initialized ) call static_init_model
-
- lon_a = LON(:)
- lat_a = LAT(:)
- alt_a = ALT(:)
-
-end subroutine get_grid_val
-
-
-
+! Reads the current time and state variables from a gitm restart
+! file and packs them into a dart state vector.
+!
! FIXME:
-! this routine needs:
-! 1. a base dirname for the restart files.
-! they will have the format 'dirname/bNNNN.rst' where NNNN has
-! leading 0s and is the block number. blocks start in the
-! southwest corner of the lat/lon grid and go west first, then
-! north and end in the northeast corner. (assuming var 'dirname')
-! the other info is in 'dirname/header.rst'
+! this routine needs:
+! 1. a base dirname for the restart files.
+! they will have the format 'dirname/bNNNN.rst' where NNNN has
+! leading 0s and is the block number. blocks start in the
+! southwest corner of the lat/lon grid and go west first, then
+! north and end in the northeast corner. (assuming var 'dirname')
+! the other info is in 'dirname/header.rst'
!
-! 2. the overall grid size, lon/lat/alt when you've read in all
-! the blocks. (nGridLon, nGridLat, nGridAlt, will compute totalVarSize)
+! 2. the overall grid size, lon/lat/alt when you've read in all
+! the blocks. (nGridLon, nGridLat, nGridAlt, will compute totalVarSize)
!
-! 3. the number of blocks in Lon and Lat (nBlocksLon, nBlocksLat)
+! 3. the number of blocks in Lon and Lat (nBlocksLon, nBlocksLat)
!
-! 4. the number of lon/lats in a single grid block (nLons, nLats, nAlts)
+! 4. the number of lon/lats in a single grid block (nLons, nLats, nAlts)
!
-! 5. the number of neutral species (and probably a mapping between
-! the species number and the variable name) (nSpeciesTotal, nSpecies)
+! 5. the number of neutral species (and probably a mapping between
+! the species number and the variable name) (nSpeciesTotal, nSpecies)
!
-! 6. the number of ion species (ditto - numbers <-> names) (nIons)
+! 6. the number of ion species (ditto - numbers <-> names) (nIons)
!
! we assume that the 'UseTopography' flag is false - that all columns
! have the same altitude arrays. this is true on earth but not on
! other planets.
!
-! in addition to reading in the state data, it fills Longitude,
-! Latitude, and Altitude arrays with the grid spacing. this grid
-! is orthogonal and rectangular but can have irregular spacing along
-! any or all of the three dimensions.
-!
+! in addition to reading in the state data, it fills Longitude,
+! Latitude, and Altitude arrays with the grid spacing. this grid
+! is orthogonal and rectangular but can have irregular spacing along
+! any or all of the three dimensions.
-subroutine restart_file_to_statevector(dirname, state_vector, model_time)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list