[Dart-dev] [7867] DART/trunk/models/tiegcm: Fixed a potentially bad error in that the variables were still being read from

nancy at ucar.edu nancy at ucar.edu
Thu Apr 9 09:14:11 MDT 2015


Revision: 7867
Author:   thoar
Date:     2015-04-09 09:14:11 -0600 (Thu, 09 Apr 2015)
Log Message:
-----------
Fixed a potentially bad error in that the variables were still being read from
the restart file even when 'secondary' was specified.
Cleaned up the logic specifying the model time.
This also allows for a more accurate 'location' of the f10_7 parameter,
which is best described as local noon. This cannot be determined 'once' by
static_init_model, for example - and the model_mod (in general) does not
have access to the current time. Turns out - every time the model is advanced
the ens_mean_for_model() is recomputed, so the model_time is re-read
from the restart file in ens_mean_for_model().

Also - added the model_mod.html - even though it still needs work.
All the bases are covered, but the overview section, references, future plans, etc
could use help.

Modified 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/work/input.nml

Added Paths:
-----------
    DART/trunk/models/tiegcm/model_mod.html

-------------- next part --------------
Modified: DART/trunk/models/tiegcm/dart_to_model.f90
===================================================================
--- DART/trunk/models/tiegcm/dart_to_model.f90	2015-04-08 23:29:37 UTC (rev 7866)
+++ DART/trunk/models/tiegcm/dart_to_model.f90	2015-04-09 15:14:11 UTC (rev 7867)
@@ -30,7 +30,7 @@
                              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
+                             dart_vector_to_tiegcm
 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, &
@@ -60,7 +60,6 @@
 ! 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                :: iunit, io, file_unit, x_size
@@ -103,8 +102,7 @@
 
 ! write fields to the binary TIEGCM restart file
 
-file_name = get_restart_file_name()
-call update_TIEGCM_restart(x_state, file_name, model_time)
+call dart_vector_to_tiegcm(x_state, model_time)
 
 if ( advance_time_present ) then
    ! update TIEGCM namelist variables used in advance_model.csh

Modified: DART/trunk/models/tiegcm/model_mod.f90
===================================================================
--- DART/trunk/models/tiegcm/model_mod.f90	2015-04-08 23:29:37 UTC (rev 7866)
+++ DART/trunk/models/tiegcm/model_mod.f90	2015-04-09 15:14:11 UTC (rev 7867)
@@ -75,9 +75,8 @@
           ens_mean_for_model
 
 !TIEGCM specific routines
-public ::   read_TIEGCM_restart, &
-          update_TIEGCM_restart, &
-          get_restart_file_name, &
+public :: tiegcm_to_dart_vector, &
+          dart_vector_to_tiegcm, &
           get_f107_value,        &
           test_interpolate
 
@@ -284,6 +283,9 @@
 ! Might as well use the Gregorian Calendar
 call set_calendar_type('Gregorian')
 
+! Convert the last year/doy/hour/minute to a dart time.
+state_time = get_state_time(tiegcm_restart_file_name)
+
 ! 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
@@ -630,8 +632,8 @@
 else
 
    write(string1,*)trim(progvar(ivar)%varname),'has unsupported rank',progvar(ivar)%rank
-   write(string2,*)trim(progvar(ivar)%varname),'or unsupported vertical system ',which_vert,&
-                   'for that rank.'
+   write(string2,*)trim(progvar(ivar)%varname),'or unsupported vertical system ', &
+                   which_vert, 'for that rank.'
    call error_handler(E_ERR,'model_interpolate', string1, &
               source, revision, revdate, text2=string2)
 
@@ -744,7 +746,8 @@
    height    = get_height(ivar, lon_index, lat_index, lev_index)
 
    if (do_output() .and. (debug > 3)) then
-      absindx = get_index( ivar, indx1=lon_index, indx2=lat_index, indx3=lev_index, knownindex=index_in)
+      absindx = get_index( ivar, indx1=lon_index, indx2=lat_index, &
+                           indx3=lev_index, knownindex=index_in)
    endif
 
    location  = set_location(lons(lon_index), lats(lat_index), height, VERTISHEIGHT)
@@ -1088,7 +1091,8 @@
 !-------------------------------------------------------------------------------
 
 call nc_check(nf90_sync(ncid), 'nc_write_model_atts', 'sync')
-if (do_output() .and. (debug > 1)) write (*,*) 'nc_write_model_atts: netCDF file ', ncid, ' is synched '
+if (do_output() .and. (debug > 1)) &
+     write (*,*) 'nc_write_model_atts: netCDF file ', ncid, ' is synched '
 
 ierr = 0 ! If we got here, things went well.
 
@@ -1178,8 +1182,8 @@
       DimCheck : do i = 1,ncNdims
 
          write(string1,'(A,i2,A)') 'inquire dimension ',i,trim(string2)
-         call nc_check(nf90_inquire_dimension(ncid,dimIDs(i),name=dimnames(i),len=dimlen), &
-               'nc_write_model_vars', trim(string1))
+         call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), name=dimnames(i), &
+                      len=dimlen), 'nc_write_model_vars', trim(string1))
 
          ! Dont care about the length of these two, setting to 1 later.
          if (dimIDs(i) == CopyDimID) cycle DimCheck
@@ -1203,10 +1207,14 @@
       where(dimIDs == TimeDimID) mycount = 1
 
       if (do_output() .and. (debug > 3)) then
-         write(string1,*)'nc_write_model_vars',trim(progvar(ivar)%varname),' file id ',ncid
-         write(string2,*)'nc_write_model_vars '//trim(varname)//' start is ',mystart(1:ncNdims)
-         write(string3,*)'nc_write_model_vars '//trim(varname)//' count is ',mycount(1:ncNdims)
-         call error_handler(E_MSG,'nc_write_model_vars',string1,text2=string2,text3=string3)
+         write(string1,*)'nc_write_model_vars',trim(progvar(ivar)%varname), &
+                         ' file id ',ncid
+         write(string2,*)'nc_write_model_vars '//trim(varname)//' start is ', &
+                         mystart(1:ncNdims)
+         write(string3,*)'nc_write_model_vars '//trim(varname)//' count is ', &
+                         mycount(1:ncNdims)
+         call error_handler(E_MSG, 'nc_write_model_vars', string1, &
+                    text2=string2, text3=string3)
       endif
 
       ! Now that we have the hyperslab indices ... it is easy.
@@ -1381,6 +1389,14 @@
 
 ens_mean = filter_ens_mean
 
+! The most current model time must be read from 
+! tiegcm_restart_p.nc after each forecast step. 
+! state_time is used by get_state_meta_data to determine
+! the location of the f10_7 parameter (local noon)
+! The ens_mean_for_model() is called after each forecast step.
+
+state_time = get_state_time(tiegcm_restart_file_name)
+
 end subroutine ens_mean_for_model
 
 
@@ -1389,105 +1405,62 @@
 !===============================================================================
 
 
-subroutine read_TIEGCM_restart(file_name, statevec, model_time)
+subroutine tiegcm_to_dart_vector(statevec, model_time)
 !
 ! Read TIEGCM restart file fields and pack it into a DART vector
 !
-character(len=*),       intent(in)  :: file_name
 real(r8), dimension(:), intent(out) :: statevec
 type(time_type),        intent(out) :: model_time
 
-integer :: ncid, VarID
+integer :: ncid, ncid1, ncid2, VarID, i, ivar
 integer :: time_dimlen, dimlen, ncNdims
 integer :: LonDimID, LatDimID, LevDimID, IlevDimID, TimeDimID
 integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, mystart, mycount
+character(len=512) :: filename
 
 real(r8), allocatable, dimension(:)       :: temp1D
 real(r8), allocatable, dimension(:,:)     :: temp2D
 real(r8), allocatable, dimension(:,:,:)   :: temp3D
 real(r8), allocatable, dimension(:,:,:,:) :: temp4D
 
-integer :: i, ivar
-
 if ( .not. module_initialized ) call static_init_model
 
-if( .not. file_exist(file_name)) then
-   write(string1,*)trim(file_name)//' does not exist.'
-   call error_handler(E_ERR,'read_TIEGCM_restart',string1,source,revision,revdate)
-endif
+! All variables come from either the restart or secondary file.
+! Just open them both. 
+! Since the restart and secondary files have already had their metadata
+! checked and compared favorably, we do not need to check sizes again.
 
-call error_handler(E_MSG,'read_TIEGCM_restart:','reading restart ['//trim(file_name)//']')
+call nc_check(nf90_open(tiegcm_restart_file_name, NF90_NOWRITE, ncid1), &
+              'tiegcm_to_dart_vector','open '//trim(tiegcm_restart_file_name))
 
-! Open the netCDF file and then read all the static information.
+call nc_check(nf90_open(tiegcm_secondary_file_name, NF90_NOWRITE, ncid2), &
+              'tiegcm_to_dart_vector','open '//trim(tiegcm_secondary_file_name))
 
-call nc_check( nf90_open( file_name, NF90_NOWRITE, ncid ), &
-                              'read_TIEGCM_restart', 'open')
+! state_time is set by static_init_model()
+model_time = state_time
 
-!... check for matching dimensions
-call nc_check( nf90_inq_dimid(ncid, 'lon', LonDimID), &
-         'read_TIEGCM_restart', 'inq_dimid lon')
-call nc_check( nf90_inquire_dimension(ncid, LonDimID, len=dimlen), &
-         'read_TIEGCM_restart', 'inquire_dimension lon')
-if (dimlen .ne. nlon) then
-  write(string1, *) trim(file_name), ' dim_lon = ',dimlen, ' DART expects ',nlon
-  call error_handler(E_ERR,'read_TIEGCM_restart',string1,source,revision,revdate)
-endif
-
-call nc_check( nf90_inq_dimid(ncid, 'lat', LatDimID), &
-         'read_TIEGCM_restart', 'inq_dimid lat')
-call nc_check( nf90_inquire_dimension(ncid, LatDimID, len=dimlen), &
-         'read_TIEGCM_restart', 'inquire_dimension lat')
-if (dimlen .ne. nlat) then
-  write(string1, *) trim(file_name), ' dim_lat = ',dimlen, ' DART expects ',nlat
-  call error_handler(E_ERR,'read_TIEGCM_restart',string1,source,revision,revdate)
-endif
-
-call nc_check( nf90_inq_dimid(ncid, 'lev', LevDimID), &
-         'read_TIEGCM_restart', 'inq_dimid lev')
-call nc_check( nf90_inquire_dimension(ncid, LevDimID, len=dimlen), &
-         'read_TIEGCM_restart', 'inquire_dimension lev')
-if (dimlen .ne. (nlev+1)) then
-  write(string1, *) trim(file_name), ' dim_lev = ',dimlen, ' DART expects ',nlev+1
-  call error_handler(E_ERR,'read_TIEGCM_restart',string1,source,revision,revdate)
-endif
-
-call nc_check( nf90_inq_dimid(ncid, 'ilev', IlevDimID), &
-         'read_TIEGCM_restart', 'inq_dimid ilev')
-call nc_check( nf90_inquire_dimension(ncid, IlevDimID, len=dimlen), &
-         'read_TIEGCM_restart', 'inquire_dimension ilev')
-if (dimlen .ne. nilev) then
-  write(string1, *) trim(file_name), ' dim_ilev = ',dimlen, ' DART expects ',nilev
-  call error_handler(E_ERR,'read_TIEGCM_restart',string1,source,revision,revdate)
-endif
-
-call nc_check( nf90_inq_dimid(ncid, 'time', TimeDimID), &
-         'read_TIEGCM_restart', 'inq_dimid time')
-call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=time_dimlen ), &
-         'read_TIEGCM_restart', 'inquire_dimension time')
-
 ! Loop over all variables in DART state.
-! Make sure the shape of the variable matches what we expect.
-! Get the hyperslab with the most current (last) time.
-
 ReadVariable: do ivar = 1,nfields
 
-   string2 = trim(file_name)//' '//trim(progvar(ivar)%varname)
+   ! Handle the special and CALCULATE cases first.
 
    if (trim(progvar(ivar)%varname) == 'f10_7') then
       statevec( progvar(ivar)%index1 ) = f10_7
       cycle ReadVariable
    endif
 
-   ! ZG is already a module variable that is required.
+   ! ZG is already a module variable that is required and comes from
+   ! the tiegcm_s.nc (aka "secondary" file).
    if (trim(progvar(ivar)%varname) == 'ZG') then
       call prog_var_to_vector(ivar, ZG, statevec)
       ! deallocate(ZG)
-      ! FIXME After the create_vtec() rewrite, ZG should be deallocated.
+      ! FIXME After the create_vtec() rewrite, ZG could be deallocated.
       ! It can be accessed directly from the DART state
       ! maybe should also be removed from the read_tiegcm_seconday() routine
       cycle ReadVariable
    endif
 
+   ! VTEC is a 'CALCULATE'
    if (trim(progvar(ivar)%varname) == 'VTEC') then
       allocate(temp2D(nlon,nlat))
       call create_vtec(ncid, time_dimlen, temp2D)
@@ -1496,16 +1469,42 @@
       cycle ReadVariable
    endif
 
+   ! Check to see which file contains the variable.
+   if     (variable_table(ivar,VT_ORIGININDX) == 'RESTART') then
+      ncid     = ncid1
+      filename = tiegcm_restart_file_name
+
+      call get_diminfo(filename, ncid, LonDimID=LonDimID, LatDimID=LatDimID, &
+                  LevDimID=LevDimID, iLevDimID=iLevDimID, &
+                 TimeDimID=TimeDimID, ntimes=time_dimlen)
+
+   elseif (variable_table(ivar,VT_ORIGININDX) == 'SECONDARY') then
+      ncid     = ncid2
+      filename = tiegcm_secondary_file_name
+
+      call get_diminfo(filename, ncid, LonDimID=LonDimID, LatDimID=LatDimID, &
+                  LevDimID=LevDimID, iLevDimID=iLevDimID, &
+                 TimeDimID=TimeDimID, ntimes=time_dimlen)
+
+   else ! MUST be a CALCULATE we do not know how to calculate
+      string1 = 'unknown operation to calculate variable '//trim(progvar(ivar)%varname)
+      string2 = 'valid options are "RESTART", "SECONDARY", or "CALCULATE"'
+      call error_handler(E_ERR, 'tiegcm_to_dart_vector', string1, &
+                 source, revision, revdate, text2=string2)
+   endif
+
+   string2 = trim(filename)//' '//trim(progvar(ivar)%varname)
+
    call nc_check(nf90_inq_varid(ncid, trim(progvar(ivar)%varname), VarID), &
-          'read_TIEGCM_restart', 'inq_varid '//trim(string2))
+          'tiegcm_to_dart_vector', 'inq_varid '//trim(string2))
 
    call nc_check(nf90_inquire_variable(ncid,VarID,dimids=dimIDs,ndims=ncNdims), &
-            'read_TIEGCM_restart', 'inquire '//trim(string2))
+            'tiegcm_to_dart_vector', 'inquire '//trim(string2))
 
    if (ncNdims /= progvar(ivar)%rank) then
       write(string1,*)trim(string2),'has ',ncNDims,'dimensions.'
       write(string3,*)'same thing in DART has ',progvar(ivar)%rank,' dimensions.'
-      call error_handler(E_ERR, 'read_TIEGCM_restart', string1, &
+      call error_handler(E_ERR, 'tiegcm_to_dart_vector', string1, &
                       source, revision, revdate, text2=string3)
    endif
 
@@ -1515,7 +1514,7 @@
 
       write(string1,'(a,i2,A)') 'inquire dimension ',i,trim(string2)
       call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
-               'read_TIEGCM_restart', trim(string1))
+               'tiegcm_to_dart_vector', trim(string1))
 
       ! Check the shape of the variable
       if ( dimIDs(i) == LevDimID ) dimlen = dimlen - 1
@@ -1525,7 +1524,7 @@
       elseif ( dimlen /= progvar(ivar)%dimlens(i) ) then
          write(string1,*) trim(string2),' has dim/dimlen ',i,dimlen
          write(string3,*)' expected dimlen of ', progvar(ivar)%dimlens(i)
-         call error_handler(E_ERR, 'read_TIEGCM_restart', string1, &
+         call error_handler(E_ERR, 'tiegcm_to_dart_vector', string1, &
                          source, revision, revdate, text2=string3)
       endif
 
@@ -1533,6 +1532,8 @@
 
    enddo DimCheck
 
+   ! Get the hyperslab with the most current (last) time.
+
    where(dimIDs == TimeDimID) mystart = time_dimlen
    where(dimIDs == TimeDimID) mycount = 1
    where(dimIDs ==  LonDimID) mycount = nlon
@@ -1541,9 +1542,9 @@
    where(dimIDs == iLevDimID) mycount = nilev
 
    if (do_output() .and. (debug > 3)) then
-      write(*,*)'read_TIEGCM_restart:',trim(string2)
-      write(*,*)'read_TIEGCM_restart: start ',mystart(1:ncNdims)
-      write(*,*)'read_TIEGCM_restart: count ',mycount(1:ncNdims)
+      write(*,*)'tiegcm_to_dart_vector:',trim(string2)
+      write(*,*)'tiegcm_to_dart_vector: start ',mystart(1:ncNdims)
+      write(*,*)'tiegcm_to_dart_vector: count ',mycount(1:ncNdims)
       write(*,*)
    endif
 
@@ -1551,7 +1552,7 @@
       allocate(temp1D(mycount(1)))
       call nc_check(nf90_get_var(ncid, VarID, values=temp1D, &
               start = mystart(1:ncNdims), count = mycount(1:ncNdims)), &
-              'read_TIEGCM_restart', 'get_var '//trim(string2))
+              'tiegcm_to_dart_vector', 'get_var '//trim(string2))
       call apply_attributes(ivar, ncid, VarID, temp1D)
       call prog_var_to_vector(ivar, temp1D, statevec)
       deallocate(temp1D)
@@ -1560,7 +1561,7 @@
       allocate(temp2D(mycount(1), mycount(2)))
       call nc_check(nf90_get_var(ncid, VarID, values=temp2D, &
               start = mystart(1:ncNdims), count = mycount(1:ncNdims)), &
-              'read_TIEGCM_restart', 'get_var '//trim(string2))
+              'tiegcm_to_dart_vector', 'get_var '//trim(string2))
       call apply_attributes(ivar, ncid, VarID, temp2D)
       call prog_var_to_vector(ivar, temp2D, statevec)
       deallocate(temp2D)
@@ -1569,7 +1570,7 @@
       allocate(temp3D(mycount(1), mycount(2), mycount(3)))
       call nc_check(nf90_get_var(ncid, VarID, values=temp3D, &
               start = mystart(1:ncNdims), count = mycount(1:ncNdims)), &
-              'read_TIEGCM_restart', 'get_var '//trim(string2))
+              'tiegcm_to_dart_vector', 'get_var '//trim(string2))
       call apply_attributes(ivar, ncid, VarID, temp3D)
       call prog_var_to_vector(ivar, temp3D, statevec)
       deallocate(temp3D)
@@ -1578,38 +1579,34 @@
       allocate(temp4D(mycount(1), mycount(2), mycount(3), mycount(4)))
       call nc_check(nf90_get_var(ncid, VarID, values=temp4D, &
               start = mystart(1:ncNdims), count = mycount(1:ncNdims)), &
-              'read_TIEGCM_restart', 'get_var '//trim(string2))
+              'tiegcm_to_dart_vector', 'get_var '//trim(string2))
       call apply_attributes(ivar, ncid, VarID, temp4D)
       call prog_var_to_vector(ivar, temp4D, statevec)
       deallocate(temp4D)
    else
       write(string1,*) trim(string2),' has unsupported number of dims (', &
                                      progvar(ivar)%rank,')'
-      call error_handler(E_ERR, 'read_TIEGCM_restart', string1, &
+      call error_handler(E_ERR, 'tiegcm_to_dart_vector', string1, &
                          source, revision, revdate)
    endif
 
 enddo ReadVariable
 
-! Convert the last year/doy/hour/minute to a dart time.
-model_time = get_state_time(trim(file_name), ncid)
-state_time = model_time   ! state_time is scoped for entire module
+call nc_check( nf90_close(ncid1), 'tiegcm_to_dart_vector', 'close restart')
+call nc_check( nf90_close(ncid2), 'tiegcm_to_dart_vector', 'close secondary')
 
-call nc_check( nf90_close(ncid), 'read_TIEGCM_restart', 'close')
+end subroutine tiegcm_to_dart_vector
 
-end subroutine read_TIEGCM_restart
 
-
 !-------------------------------------------------------------------------------
 
 
-subroutine update_TIEGCM_restart(statevec, filename, dart_time)
+subroutine dart_vector_to_tiegcm(statevec, dart_time)
 ! Updates TIEGCM restart file fields with the reshaped variables
 ! in the DART state vector. The variable MISSING attribute has to be
 ! reinstated ... etc.
 
 real(r8), dimension(:), intent(in) :: statevec
-character(len=*),       intent(in) :: filename
 type(time_type),        intent(in) :: dart_time
 
 integer :: i, ivar
@@ -1620,47 +1617,33 @@
 character(len=NF90_MAX_NAME) :: varname
 character(len=NF90_MAX_NAME), dimension(NF90_MAX_VAR_DIMS) :: dimnames
 
-type(time_type) :: file_time
-
 if ( .not. module_initialized ) call static_init_model
 
-if( .not. file_exist(filename)) then
-   write(string1,*) trim(filename),' not available.'
-   call error_handler(E_ERR,'update_TIEGCM_restart',string1,source,revision,revdate)
-endif
+call nc_check(nf90_open(tiegcm_restart_file_name,NF90_WRITE,ncid),'dart_vector_to_tiegcm','open')
 
-if (do_output()) &
-   call error_handler(E_MSG,'update_TIEGCM_restart','opening',text2=filename)
+call get_diminfo(tiegcm_restart_file_name, ncid, LevDimID=LevDimID, &
+                 TimeDimID=TimeDimID, ntimes=timeindex)
 
-call nc_check(nf90_open(trim(filename),NF90_WRITE,ncid),'update_TIEGCM_restart','open')
-
-call SanityCheck(trim(filename), ncid, TimeDimID=TimeDimID, ntimes=timeindex)
-
-call nc_check(nf90_inq_dimid(ncid, 'lev', LevDimID), &
-                   'update_TIEGCM_restart', 'inq_dimid lev'//trim(filename))
-
-file_time = get_state_time(trim(filename), ncid, timeindex)
-
-if ( file_time /= dart_time ) then
-   call print_time(dart_time,'DART   current time',logfileunit)
-   call print_time(file_time,'TIEGCM current time',logfileunit)
-   call print_time(dart_time,'DART   current time')
-   call print_time(file_time,'TIEGCM current time')
-   write(string1,*)trim(filename),' current time /= model time. FATAL error.'
-   call error_handler(E_ERR,'update_TIEGCM_restart',string1,source,revision,revdate)
+if ( state_time /= dart_time ) then
+   call print_time( dart_time,'DART   current time',logfileunit)
+   call print_time(state_time,'TIEGCM current time',logfileunit)
+   call print_time( dart_time,'DART   current time')
+   call print_time(state_time,'TIEGCM current time')
+   write(string1,*)trim(tiegcm_restart_file_name),' DART time /= TIEGCM time. FATAL error.'
+   call error_handler(E_ERR,'dart_vector_to_tiegcm',string1,source,revision,revdate)
 endif
 
 if (do_output() .and. (debug > 0)) then
-   call print_date(file_time,'update_TIEGCM_restart: date in restart file '//trim(filename))
-   call print_time(file_time,'update_TIEGCM_restart: time in restart file '//trim(filename))
-   call print_date(file_time,'update_TIEGCM_restart: date in restart file '//trim(filename),logfileunit)
-   call print_time(file_time,'update_TIEGCM_restart: time in restart file '//trim(filename),logfileunit)
+   call print_date(state_time,'dart_vector_to_tiegcm: date in restart file '//trim(tiegcm_restart_file_name))
+   call print_time(state_time,'dart_vector_to_tiegcm: time in restart file '//trim(tiegcm_restart_file_name))
+   call print_date(state_time,'dart_vector_to_tiegcm: date in restart file '//trim(tiegcm_restart_file_name),logfileunit)
+   call print_time(state_time,'dart_vector_to_tiegcm: time in restart file '//trim(tiegcm_restart_file_name),logfileunit)
 endif
 
 UPDATE : do ivar = 1,nfields
 
    varname = trim(progvar(ivar)%varname)
-   string2 = trim(filename)//' '//trim(varname)
+   string2 = trim(tiegcm_restart_file_name)//' '//trim(varname)
 
    ! Skip the variables that are not supposed to be updated.
    if ( .not. progvar(ivar)%update ) cycle UPDATE
@@ -1669,10 +1652,10 @@
    ! At the same time, create the mystart(:) and mycount(:) arrays.
 
    call nc_check(nf90_inq_varid(ncid, varname, VarID), &
-                 'update_TIEGCM_restart', 'inq_varid '//trim(string2))
+                 'dart_vector_to_tiegcm', 'inq_varid '//trim(string2))
 
    call nc_check(nf90_inquire_variable(ncid,VarID,dimids=dimIDs,ndims=ncNdims), &
-                 'update_TIEGCM_restart', 'inquire '//trim(string2))
+                 'dart_vector_to_tiegcm', 'inquire '//trim(string2))
 
    mystart(:) = 1
    mycount(:) = 1
@@ -1680,7 +1663,7 @@
 
       write(string1,'(''inquire dimension'',i2,A)') i,trim(string2)
       call nc_check(nf90_inquire_dimension(ncid,dimIDs(i),name=dimnames(i),len=dimlen), &
-                    'update_TIEGCM_restart', string1)
+                    'dart_vector_to_tiegcm', string1)
 
       ! Dont care about the length of Time Dimension, setting to 1 later.
       if (dimIDs(i) == TimeDimID) cycle DimCheck
@@ -1689,7 +1672,7 @@
       if ( dimlen /= progvar(ivar)%dimlens(i) ) then
          write(string1,*) trim(string2),' dim/dimlen ',i,dimlen
          write(string2,*)'expected it to be ',progvar(ivar)%dimlens(i)
-         call error_handler(E_ERR, 'update_TIEGCM_restart', string1, &
+         call error_handler(E_ERR, 'dart_vector_to_tiegcm', string1, &
                       source, revision, revdate, text2=string2)
       endif
 
@@ -1701,7 +1684,7 @@
       write(string1,*) trim(string2),&
          ' required to have "Time" as the last/unlimited dimension'
       write(string2,*)' last dimension is ',trim(dimnames(ncNdims))
-      call error_handler(E_ERR, 'update_TIEGCM_restart', string1, &
+      call error_handler(E_ERR, 'dart_vector_to_tiegcm', string1, &
                    source, revision, revdate, text2=string2)
    endif
 
@@ -1709,10 +1692,10 @@
    where(dimIDs == TimeDimID) mycount = 1
 
    if (do_output() .and. (debug > 3)) then
-      write(*,*)'update_TIEGCM_restart '//trim(varname)//' start is ',mystart(1:ncNdims)
-      write(*,*)'update_TIEGCM_restart '//trim(varname)//' count is ',mycount(1:ncNdims)
+      write(*,*)'dart_vector_to_tiegcm '//trim(varname)//' start is ',mystart(1:ncNdims)
+      write(*,*)'dart_vector_to_tiegcm '//trim(varname)//' count is ',mycount(1:ncNdims)
       do i = 1,NcNdims
-         write(*,*)'update_TIEGCM_restart '//trim(varname)//' dim ',&
+         write(*,*)'dart_vector_to_tiegcm '//trim(varname)//' dim ',&
                    i,' is ',trim(dimnames(i))
       enddo
    endif
@@ -1725,27 +1708,15 @@
 
 enddo UPDATE
 
-call nc_check(nf90_sync( ncid), 'update_TIEGCM_restart', 'sync')
-call nc_check(nf90_close(ncid), 'update_TIEGCM_restart', 'close')
+call nc_check(nf90_sync( ncid), 'dart_vector_to_tiegcm', 'sync')
+call nc_check(nf90_close(ncid), 'dart_vector_to_tiegcm', 'close')
 
-end subroutine update_TIEGCM_restart
+end subroutine dart_vector_to_tiegcm
 
 
 !-------------------------------------------------------------------------------
 
 
-function get_restart_file_name()
-character(len=256) :: get_restart_file_name
-
-if ( .not. module_initialized ) call static_init_model
-get_restart_file_name = adjustl(tiegcm_restart_file_name)
-
-end function get_restart_file_name
-
-
-!-------------------------------------------------------------------------------
-
-
 function get_f107_value(x)
 real(r8), dimension(:), intent(in) :: x
 real(r8)                           :: get_f107_value
@@ -2151,7 +2122,7 @@
 
 call nc_check(nf90_open(file_name, NF90_NOWRITE, ncid), 'read_TIEGCM_secondary', 'open')
 
-call SanityCheck(file_name, ncid, TimeDimID=TimeDimID, ntimes=time_dimlen)
+call get_diminfo(file_name, ncid, TimeDimID=TimeDimID, ntimes=time_dimlen)
 
 allocate(ZG(nlon,nlat,nilev)) ! comes from module storage
 
@@ -2335,7 +2306,7 @@
 
 !-------------------------------------------------------------------------------
 ! Now that we know how many variables, etc., fill metadata structure 'progvar'
-! read_TIEGCM_restart() uses this structure to figure out what to put in the
+! tiegcm_to_dart_vector() uses this structure to figure out what to put in the
 ! DART state vector.
 !-------------------------------------------------------------------------------
 
@@ -2365,7 +2336,8 @@
       write(string1,'(''unknown option ['',a,''] for variable '',a)') &
                              trim(variable_table(ivar,VT_ORIGININDX)), trim(varname)
       string2 = 'valid options are "RESTART", "SECONDARY", or "CALCULATE"'
-      call error_handler(E_ERR,'verify_variables',string1,source,revision,revdate,text2=string2)
+      call error_handler(E_ERR, 'verify_variables', string1, &
+                 source, revision, revdate, text2=string2)
    endif
 
    progvar(ivar)%varname           = varname
@@ -2501,7 +2473,7 @@
       write(string1,'(''inquire dimension'',i2,A)') i,trim(string2)
       call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), name=dimname, len=dimlen), &
                                           'verify_variables', string1)
-      ! read_TIEGCM_restart only reads the latest time step, so no matter how many
+      ! tiegcm_to_dart_vector() only reads the latest time step, so no matter how many
       ! time steps are defined, the time dimension is only 1
       if (trim(dimname) == 'time') dimlen = 1
 
@@ -2737,7 +2709,8 @@
 !-------------------------------------------------------------------------------
 
 
-subroutine vert_interp(x, lon_index, lat_index, height, ikind, vertstagger, ivar, val, istatus)
+subroutine vert_interp(x, lon_index, lat_index, height, ikind, vertstagger, &
+                       ivar, val, istatus)
 ! returns the value at an arbitrary height on an existing horizontal grid location.
 ! istatus == 0 is a 'good' return code.
 
@@ -2961,7 +2934,7 @@
 integer, optional, intent(in) :: knownindex
 integer                       :: get_index
 
-integer :: ndim1, ndim2, ndim3, ndim4
+integer :: ndim1, ndim2, ndim3
 integer :: ix1, ix2, ix3, ix4
 
 real(r8) :: height
@@ -3642,7 +3615,7 @@
 !-------------------------------------------------------------------------------
 
 
-subroutine SanityCheck(filename, ncid, LonDimID, LatDimID, LevDimID, iLevDimId, &
+subroutine get_diminfo(filename, ncid, LonDimID, LatDimID, LevDimID, iLevDimID, &
                        TimeDimID, ntimes)
 ! Just make sure the dimensions of the netCDF file match those
 ! that were used for static_init_model()
@@ -3658,77 +3631,77 @@
 
 integer :: dimlen, DimID, DimID2
 
-call nc_check(nf90_inq_dimid(ncid, 'lon', DimID), 'SanityCheck', 'inq_dimid lon')
+call nc_check(nf90_inq_dimid(ncid, 'lon', DimID), 'get_diminfo', 'inq_dimid lon')
 call nc_check(nf90_inquire_dimension(ncid, DimID, len=dimlen), &
-   'SanityCheck', 'inquire_dimension lon')
-if (dimlen .ne. nlon) then
+   'get_diminfo', 'inquire_dimension lon')
+
+if (dimlen /= nlon) then
    write(string1, *) trim(filename), ' dim_lon = ',dimlen, ' DART expects ',nlon
-   call error_handler(E_ERR,'SanityCheck',string1,source,revision,revdate)
+   call error_handler(E_ERR,'get_diminfo',string1,source,revision,revdate)
 endif
-if (present(LonDimID)) LonDimID = dimlen
+if (present(LonDimID)) LonDimID = DimID
 
-
-call nc_check(nf90_inq_dimid(ncid, 'lat', DimID), 'SanityCheck', 'inq_dimid lat')
+call nc_check(nf90_inq_dimid(ncid, 'lat', DimID), 'get_diminfo', 'inq_dimid lat')
 call nc_check(nf90_inquire_dimension(ncid, DimID, len=dimlen), &
-   'SanityCheck', 'inquire_dimension lat')
-if (dimlen .ne. nlat) then
+   'get_diminfo', 'inquire_dimension lat')
+if (dimlen /= nlat) then
    write(string1, *) trim(filename), ' dim_lat = ',dimlen, ' DART expects ',nlat
-   call error_handler(E_ERR,'SanityCheck',string1,source,revision,revdate)
+   call error_handler(E_ERR,'get_diminfo',string1,source,revision,revdate)
 endif
-if (present(LatDimID)) LatDimID = dimlen
+if (present(LatDimID)) LatDimID = DimID
 
 ! TIEGCM has a useless top level for everything defined on 'lev'
 ! 'nlev' had a local value (i.e. in DART) that is 1 less than 'nilev'
 
-call nc_check(nf90_inq_dimid(ncid, 'lev', DimID), 'SanityCheck', 'inq_dimid lev')
+call nc_check(nf90_inq_dimid(ncid, 'lev', DimID), 'get_diminfo', 'inq_dimid lev')
 call nc_check(nf90_inquire_dimension(ncid, DimID, len=dimlen), &
-   'SanityCheck', 'inquire_dimension lev')
-if (dimlen .ne. (nlev+1)) then
+   'get_diminfo', 'inquire_dimension lev')
+if (dimlen /= (nlev+1)) then
    write(string1, *) trim(filename), ' dim_lev = ',dimlen, ' DART expects ',nlev
-   call error_handler(E_ERR,'SanityCheck',string1,source,revision,revdate)
+   call error_handler(E_ERR,'get_diminfo',string1,source,revision,revdate)
 endif
-if (present(LevDimID)) LevDimID = dimlen
+if (present(LevDimID)) LevDimID = DimID
 
 
-call nc_check(nf90_inq_dimid(ncid, 'ilev', DimID), 'SanityCheck', 'inq_dimid ilev')
+call nc_check(nf90_inq_dimid(ncid, 'ilev', DimID), 'get_diminfo', 'inq_dimid ilev')
 call nc_check(nf90_inquire_dimension(ncid, DimID, len=dimlen), &
-   'SanityCheck', 'inquire_dimension ilev')
+   'get_diminfo', 'inquire_dimension ilev')
 if (dimlen .ne. nilev) then
    write(string1, *) trim(filename), ' dim_ilev = ',dimlen, ' DART expects ',nilev
-   call error_handler(E_ERR,'SanityCheck',string1,source,revision,revdate)
+   call error_handler(E_ERR,'get_diminfo',string1,source,revision,revdate)
 endif
-if (present(iLevDimID)) iLevDimID = dimlen
+if (present(iLevDimID)) iLevDimID = DimID
 
 
-call nc_check(nf90_inq_dimid(ncid, 'time', DimID), 'SanityCheck', 'inq_dimid time')
+call nc_check(nf90_inq_dimid(ncid, 'time', DimID), 'get_diminfo', 'inq_dimid time')
 call nc_check(nf90_inquire(ncid, unlimitedDimId = DimID2), &
-   'SanityCheck', 'inquire id of unlimited dimension')
-if (DimID .ne. DimID2) then
+   'get_diminfo', 'inquire id of unlimited dimension')
+if (DimID /= DimID2) then
    write(string1, *) trim(filename), 'has an unlimited dimension ID of ',DimID2
    write(string2, *) 'and a time dimension ID of ',DimID, &
                      '. DART requires them to be the same.'
-   call error_handler(E_ERR, 'SanityCheck', string1, &
+   call error_handler(E_ERR, 'get_diminfo', string1, &
               source, revision, revdate, text2=string2)
 endif
 if (present(TimeDimID)) TimeDimID = DimID
 
 
 call nc_check(nf90_inquire_dimension(ncid, DimID, len=dimlen), &
-   'SanityCheck', 'inquire_dimension time')
+   'get_diminfo', 'inquire_dimension time')
 if (present(ntimes)) ntimes = dimlen
 
-end subroutine SanityCheck
+end subroutine get_diminfo
 
 
 
-function get_state_time(filename, ncid, lasttime)
+function get_state_time(filename, ncfileid, lasttime)
 ! Gets the latest time in the netCDF file.
 character(len=*),  intent(in) :: filename
-integer,           intent(in) :: ncid
+integer, optional, intent(in) :: ncfileid
 integer, optional, intent(in) :: lasttime
 type(time_type) :: get_state_time
 
-integer :: TimeDimID, time_dimlen, DimID, dimlen, VarID
+integer :: ncid, TimeDimID, time_dimlen, DimID, dimlen, VarID
 
 integer,  parameter         :: nmtime = 3
 integer,  dimension(nmtime) :: mtime  ! day, hour, minute
@@ -3736,6 +3709,13 @@
 integer, allocatable, dimension(:,:) :: mtimetmp
 integer, allocatable, dimension(:)   :: yeartmp
 
+if ( present(ncfileid) ) then
+   ncid = ncfileid
+else
+   call nc_check(nf90_open(filename, NF90_NOWRITE, ncid), &
+              'get_state_time','open '//trim(filename))
+endif
+
 call nc_check(nf90_inq_dimid(ncid, 'time', TimeDimID), &
         'get_state_time', 'inquire id of time')
 call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=time_dimlen ), &
@@ -3770,7 +3750,7 @@
 call nc_check(nf90_inq_varid(ncid, 'year', VarID), 'get_state_time', 'inq_varid year')
 call nc_check(nf90_get_var(ncid, VarID, values=yeartmp), 'get_state_time', 'get_var year')
 
-! pick off the latest
+! pick off the latest/last
 mtime = mtimetmp(:,time_dimlen)
 year  = yeartmp(   time_dimlen)
 
@@ -3787,6 +3767,10 @@
    call print_time(get_state_time, str=trim(filename)//':get_state_time: time ')
 endif
 
+if ( .not. present(ncfileid) ) then
+   call nc_check(nf90_close(ncid), 'get_state_time', 'close '//trim(filename))
+endif
+
 end function get_state_time
 
 
@@ -3861,7 +3845,8 @@
 
 if ( progvar(ivar)%rangeRestricted == BOUNDED_BOTH ) then
    if (maxvalue < minvalue) then
-      write(string1,*)'&model_nml state_variable input error for ',trim(progvar(ivar)%varname)
+      write(string1,*)'&model_nml state_variable input error for ', &
+                      trim(progvar(ivar)%varname)
       write(string2,*)'minimum value (',minvalue,') must be less than '
       write(string3,*)'maximum value (',maxvalue,')'
       call error_handler(E_ERR,'SetVariableLimits',string1, &

Added: DART/trunk/models/tiegcm/model_mod.html
===================================================================
--- DART/trunk/models/tiegcm/model_mod.html	                        (rev 0)
+++ DART/trunk/models/tiegcm/model_mod.html	2015-04-09 15:14:11 UTC (rev 7867)
@@ -0,0 +1,1215 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN"
+          "http://www.w3.org/TR/html4/strict.dtd">
+<HTML>
+<HEAD>
+<TITLE>module model_mod (TIEGCM)</TITLE>
+<link rel="stylesheet" type="text/css" href="../../doc/html/doc.css" />
+<link href="../../doc/images/dart.ico" rel="shortcut icon" />
+</HEAD>
+<BODY>
+<A NAME="TOP"></A>
+
+<H1>MODULE model_mod (TIEGCM)</H1>
+
+<table border=0 summary="" cellpadding=5>
+<tr>
+    <td valign=middle>
+    <img src="../../doc/images/Dartboard7.png" alt="DART project logo" height=70 />
+    </td>
+    <td>
+       <P>Jump to <a href="../../index.html">DART Documentation Main Index</a><br />
+          <small><small>version information for this file: <br />
+          <!-- version tag follows, do not edit -->
+          $Id$</small></small>
+       </P></td>
+</tr>
+</table>
+
+<A HREF="#Namelist">NAMELIST</A> /
+<A HREF="#Interface">INTERFACES</A> /
+<A HREF="#FilesUsed">FILES</A> /
+<A HREF="#References">REFERENCES</A> /
+<A HREF="#Errors">ERRORS</A> /
+<A HREF="#FuturePlans">PLANS</A> /
+<A HREF="#PrivateComponents">PRIVATE COMPONENTS</A> /
+<A HREF="#Legalese">TERMS OF USE</A>
+
+<H2>Overview</H2>
+
+<P>
+This is the DART interface to the Thermosphere Ionosphere Electrodynamic 
+General Circulation Model 
+(<a href="http://www.hao.ucar.edu/modeling/tgcm/tie.php">TIEGCM</a>).
+It is <strong>strongly</strong>
+recommended that you become familiar with running TIEGCM
+<strong>before</strong> you try to run DART/TIEGCM. The TIEGCM source code
+and restart files are <strong>not</strong> included in DART, you must get them
+from the NCAR High Altitude Observatory 
+<a href="http://www.hao.ucar.edu/modeling/tgcm/download.php">download website</a>.
+DART is designed so that the TIEGCM source code can be used with no 
+modifications, as DART runs TIEGCM as a completely separate executable.
+<br /><br />
+Certain assumptions are made about the mannner in which TIEGCM is run.
+There can only be 1 each of the [restart,secondary] files. 
+The last timestep in the restart file is the only timestep which is 
+converted to a DART state, and only the last timestep in the TIEGCM 
+restart file is ever modified by DART. In the course of an assimilation
+experiment, it is usually necessary to make a short forecast with TIEGCM.
+DART writes out an ancillary file with the information necessary to 
+advance TIEGCM to the required time. The DART script 
+<em class=program>advance_model.csh</em> reads this information and
+modifies the <em class=file>tiegcm.nml</em> such that TIEGCM stops
+at the requested time. 
+<br /><br />
+The TIEGCM variables to be converted to a DART state, and possibly updated
+by the assimilation, are specified in a Fortran namelist. It is required to
+associate the variable name with a 'generic' DART counterpart -- 
+<em class=code>NE</em> is <em class=code>KIND_ELECTRON_DENSITY</em>, 
+for example. The composition of the DART input and what gets updated in the
+TIEGCM restart file are under complete user control.
+<br /><br />
+The run scripts <em class=program>run_filter.csh</em> and 
+<em class=program>run_perfect_model_obs.csh</em> are configured to run under
+the LSF queueing system. The scripting examples exploit an 'embarassingly-simple'
+parallel paradigm in that each TIEGCM instance is a single-threaded executable 
+and all ensemble members may run simultaneously. As such, there is an advantage to
+matching the ensemble size to the number of tasks. Requesting more tasks than the
+number of ensemble members may speed up the DART portion of an 
+assimilation (ie <em class=program>filter</em>) but will not make the 
+model advances faster.  <em class=program>filter</em> may be compiled with 
+MPI and can exploit all available tasks. 
+</P>
+
+<H2>Quickstart guide to running</H2>
+
+<P>It is important to understand basic DART nomenclature and mechanisms. 
+Please take the time to read and run the 
+<a href="../../tutorial/index.pdf">DART tutorial</a>.
+<br /><br />
+Both  <em class=program>run_filter.csh</em> and 
+<em class=program>run_perfect_model_obs.csh</em> are heavily internally
+commented. Please read and understand the scripts. The overall process is to
+</P>
+<OL>
+<li>Specify resources (wall-clock time, number of nodes, tasks that sort of thing).
+<li>Set shell variables to identify the location of the DART exectuables, the TIEGCM
+executables, initial ensemble, etc.
+<li>Establish a temporary working directory for the experiment.
+<li>Populate that directory with the initial ensemble and required namlists.
+<li>Convert each TIEGCM ensemble member to a DART initial conditions file. 
+<li>Run either <em class=program>filter</em> or 
+    <em class=program>run_perfect_model_obs.csh</em>.
+</OL>
+<OL><em class=program>perfect_model_obs</em> will 
+       <li>Check for any desired observations at the current time of the model 
+           state and create the synthetic observations for all observation 
+           times in the specified assimilation window. If the model needs to
+           be advanced, it then
+       <li>creates a unique run-time directory for the model advance,
+       <li>copies the required information into that directory,
+       <li>conveys the desired forecast stopping time to TIEGCM via the
+           <em class=file>tiegcm.nml</em> and
+       <li>runs a single executable of TIEGCM.
+       <li>Steps 1-5 are repeated until the input DART observation sequence file
+           has been exhausted.
+</OL>
+<OL><em class=program>filter</em> will 
+       <li>Check for any desired observations at the current time of the model 
+           state and assimilates all the observations in the specified 
+           assimilation window. If the model needs to be advanced, it then
+       <li>creates a set of run-time directories, one for each task. A single task
+           may be responsible for advancing more than one TIEGCM instance. If so,
+           each instance is done serially, one after another.
+           See the documentation for 
+           <a href="../../doc/html/filter_async_modes.html"
+          >options for running DART in parallel</a>.
+       <li>Copy the required information into that directory.

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list