[Dart-dev] [6149] DART/branches/development/models/sqg: The model_mod: get_model_time_step() routine now returns the 'model stride',
nancy at ucar.edu
nancy at ucar.edu
Fri May 17 13:38:24 MDT 2013
Revision: 6149
Author: thoar
Date: 2013-05-17 13:38:23 -0600 (Fri, 17 May 2013)
Log Message:
-----------
The model_mod:get_model_time_step() routine now returns the 'model stride',
(as opposed to the dynamical timestep of the model). Consequently,
the model can be advanced in async=0 mode with no special modifications
to the obs_model_mod:adv_1step interface.
A new pair of namelist variables (assimilation_period_[seconds,days])
has been added to allow for control over the model stride.
HOWEVER ... I have not checked to make sure the model stride is a
multiple of the dynamical timestep of the model. As committed to the
repository, the dynamical timestep is 0.0108 * 1000000/30 = 360 seconds
and the model stride is 6 hours - so we're OK for this case.
Modified Paths:
--------------
DART/branches/development/models/sqg/fft.f90
DART/branches/development/models/sqg/model_mod.f90
DART/branches/development/models/sqg/spectral_mod.f90
DART/branches/development/models/sqg/sqg.f90
DART/branches/development/models/sqg/sqg_mod.f90
DART/branches/development/models/sqg/work/input.nml
DART/branches/development/models/sqg/write_sqg_restart.f90
-------------- next part --------------
Modified: DART/branches/development/models/sqg/fft.f90
===================================================================
--- DART/branches/development/models/sqg/fft.f90 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/fft.f90 2013-05-17 19:38:23 UTC (rev 6149)
@@ -1,11 +1,3 @@
-! ===========================================================
-! <next few lines under version control, D O N O T E D I T>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
-! ===========================================================
-
! ======================================================================
! NIST Guide to Available Math Software.
! Source for module FFT from package GO.
@@ -561,3 +553,10 @@
stop
end
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
Modified: DART/branches/development/models/sqg/model_mod.f90
===================================================================
--- DART/branches/development/models/sqg/model_mod.f90 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/model_mod.f90 2013-05-17 19:38:23 UTC (rev 6149)
@@ -6,14 +6,6 @@
module model_mod
!========================================================================
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
-!========================================================================
-
-!========================================================================
! Assimilation interface for:
! Uniform-PV two-surface QG+1 model in spectral form (Hakim 2000)
!========================================================================
@@ -21,8 +13,8 @@
!================= m o d u l e i n f o r m a t i o n ==================
use types_mod, only : r8, MISSING_R8
-use time_manager_mod, only : time_type, set_time, get_time, &
- operator(<), operator(+)
+use time_manager_mod, only : time_type, set_time, get_time, set_calendar_type, &
+ print_time, operator(<), operator(+)
use location_mod, only : location_type, get_close_maxdist_init, &
get_close_obs_init, get_close_obs, set_location, &
get_location, set_location_missing
@@ -76,12 +68,16 @@
!---------------------------------------------------------------
! Namelist with default values
-logical :: output_state_vector = .false.
-real(r8) :: channel_center = 45.0_r8
-real(r8) :: channel_width = 40.0_r8
-logical :: debug = .false.
+logical :: output_state_vector = .false.
+real(r8) :: channel_center = 45.0_r8
+real(r8) :: channel_width = 40.0_r8
+logical :: debug = .false.
+integer :: assimilation_period_days = 0
+integer :: assimilation_period_seconds = 60*60*6
+
namelist /model_nml/ output_state_vector, &
channel_center, channel_width, &
+ assimilation_period_days, assimilation_period_seconds, &
debug
real, dimension(mmax,nmax), parameter :: Rblank = 0.0
@@ -92,8 +88,8 @@
! Define the location of the state variables in module storage
type(location_type), allocatable :: state_loc(:)
-type(time_type) :: time_step
-character(len=129) :: errstring
+type(time_type) :: time_step, model_stride
+character(len=129) :: msgstring1,msgstring2
type model_static
logical :: top, bot
@@ -134,6 +130,10 @@
if (do_nml_file()) write(nmlfileunit, nml=model_nml)
if (do_nml_term()) write( * , nml=model_nml)
+call set_calendar_type('Gregorian') ! TJH addition
+
+model_stride = set_model_stride()
+
! evenly spaced longitudes, latitudes and two levels
allocate( sqg_static%lons(2*kmax) )
allocate( sqg_static%lats(2*lmax) )
@@ -281,7 +281,7 @@
!========================================================================
-subroutine adv_1step(x, current_time, target_time)
+subroutine adv_1step(x, current_time)
!-------------------------------------------------
! Does a time advance for sQG model with state vector as
! input and output.
@@ -291,12 +291,11 @@
real(r8), intent(inout) :: x(:)
type(time_type), intent(in) :: current_time
-type(time_type), intent(in), optional :: target_time
real(r8) :: cxB,cyB,cxT,cyT
integer :: j
integer :: itime, cdays, cseconds, tdays, tseconds
-type(time_type) :: ctime
+type(time_type) :: ctime, target_time
real, allocatable, dimension(:,:,:) :: thxy
real, allocatable, dimension(:,:) :: thxyB, thxyT
@@ -355,9 +354,10 @@
thspB1 = 0.0; thspB2 = 0.0
thspT1 = 0.0; thspT2 = 0.0
-! advance from current_time to target_time:
+! advance from current_time
itime = 1
ctime = current_time
+target_time = current_time + model_stride
do while ( ctime < target_time )
! save old stream-function for Ekman calculation:
@@ -408,7 +408,8 @@
if ( sqg_static%top ) call tadv(thspT, tthspT, thspT1, thspT2, sqg_static%dco, first)
! zero out k=K, l=L modes on the bottom boundary:
- thspB(kmax,:) = 0.0 ; thspB(lmax,:) = 0.0
+ thspB(kmax,:) = 0.0
+ thspB(lmax,:) = 0.0
itime = itime + 1
ctime = ctime + time_step
@@ -477,8 +478,6 @@
type(time_type), intent(out) :: time
-! for now, just set to 0
-print*,'init_time'
time = set_time(0,0)
end subroutine init_time
@@ -514,8 +513,8 @@
if ( (level < 1) .or. (level > 2) ) then
istatus = 1
obs_val = MISSING_R8
- write(errstring,*)'level ',level,' must be 1 or 2'
- call error_handler(E_ERR,'model_mod:model_interpolate', errstring, source, revision, revdate)
+ write(msgstring1,*)'level ',level,' must be 1 or 2'
+ call error_handler(E_ERR,'model_mod:model_interpolate', msgstring1, source, revision, revdate)
endif
! Get globally defined lat / lon grid specs
@@ -564,10 +563,10 @@
! Do the weighted average for interpolation
if ( debug ) write(*,*) 'fracts ', lon_fract, lat_fract
do i = 1, 2
- a(i) = lon_fract * val(2, i) + (1.0 - lon_fract) * val(1, i)
+ a(i) = lon_fract * val(2, i) + (1.0_r8 - lon_fract) * val(1, i)
end do
-obs_val = lat_fract * a(2) + (1.0 - lat_fract) * a(1)
+obs_val = lat_fract * a(2) + (1.0_r8 - lat_fract) * a(1)
end subroutine model_interpolate
@@ -589,8 +588,8 @@
!! (it is out for performance reasons, but if you get any strange values, this
!! is a good first check to re-enable.)
!if (indx < 1 .or. indx > size(x)) then
-! write(errstring,*)'index ',indx,' not between 1 and ', size(x), ' (should not be possible)'
-! call error_handler(E_ERR,'model_mod:get_val', errstring, source, revision, revdate)
+! write(msgstring1,*)'index ',indx,' not between 1 and ', size(x), ' (should not be possible)'
+! call error_handler(E_ERR,'model_mod:get_val', msgstring1, source, revision, revdate)
!endif
get_val = x(indx)
@@ -601,11 +600,12 @@
function get_model_time_step()
!-----------------------------
-! Returns the time step of the model
+! Returns the time step of a single model advance,
+! NOT the dynamical timestep of the model.
type(time_type) :: get_model_time_step
-get_model_time_step = time_step
+get_model_time_step = model_stride
end function get_model_time_step
@@ -626,8 +626,8 @@
! avoid out-of-range queries
if ( index_in > model_size ) then
- write(errstring,*)'index_in ',index_in,' must be between 1 and ', model_size
- call error_handler(E_ERR,'model_mod:get_state_meta_data', errstring, source, revision, revdate)
+ write(msgstring1,*)'index_in ',index_in,' must be between 1 and ', model_size
+ call error_handler(E_ERR,'model_mod:get_state_meta_data', msgstring1, source, revision, revdate)
endif
lev_index = (index_in-1) / ((2*lmax)*(2*kmax)) + 1
@@ -725,9 +725,9 @@
"nc_write_model_atts", "inq_dimid time," // trim(filename))
if ( TimeDimID /= unlimitedDimId ) then
- write(errstring,*)"Time Dimension ID ",TimeDimID, &
+ write(msgstring1,*)"Time Dimension ID ",TimeDimID, &
" should equal Unlimited Dimension ID",unlimitedDimID
- call error_handler(E_ERR,"nc_write_model_atts", errstring, source, revision, revdate)
+ call error_handler(E_ERR,"nc_write_model_atts", msgstring1, source, revision, revdate)
endif
!-------------------------------------------------------------------------------
@@ -1005,10 +1005,9 @@
real(r8), dimension(:), intent(out) :: statevec
if ( product(shape(theta)) /= product(shape(statevec)) ) then
- write(errstring,*)'size mismatch : product(shape(theta) /= product(shape(statevec))'
- call error_handler(E_MSG,'sqg_to_dart',errstring,source,revision,revdate)
- write(errstring,*)'size mismatch : ', product(shape(theta)), ' /= ', product(shape(statevec))
- call error_handler(E_ERR,'sqg_to_dart',errstring,source,revision,revdate)
+ write(msgstring1,*)'size mismatch : product(shape(theta) /= product(shape(statevec))'
+ write(msgstring2,*)'size mismatch : ', product(shape(theta)), ' /= ', product(shape(statevec))
+ call error_handler(E_ERR,'sqg_to_dart',msgstring1,source,revision,revdate,text2=msgstring2)
else
statevec = reshape(theta,shape(statevec))
endif
@@ -1025,10 +1024,9 @@
real, dimension(:,:,:), intent(out) :: theta
if ( product(shape(statevec)) /= product(shape(theta)) ) then
- write(errstring,*)'size mismatch : product(shape(statevec)) /= product(shape(theta))'
- call error_handler(E_MSG,'dart_to_sqg',errstring,source,revision,revdate)
- write(errstring,*)'size mismatch : ', product(shape(statevec)), ' /= ', product(shape(theta))
- call error_handler(E_ERR,'dart_to_sqg',errstring,source,revision,revdate)
+ write(msgstring1,*)'size mismatch : product(shape(statevec)) /= product(shape(theta))'
+ write(msgstring2,*)'size mismatch : ', product(shape(statevec)), ' /= ', product(shape(theta))
+ call error_handler(E_ERR,'dart_to_sqg',msgstring1,source,revision,revdate,text2=msgstring2)
else
theta = reshape(statevec,shape(theta))
endif
@@ -1037,8 +1035,6 @@
!========================================================================
-!========================================================================
-
function get_model_static_data()
!-------------------------------
! Returns static info for sQG
@@ -1051,7 +1047,43 @@
!========================================================================
+function set_model_stride()
+!------------------------------------------------------------------
+! Defines the minimum amount of time to advance the model in one 'go'.
+! This is NOT the dynamical timestep of the model. It is usually a
+! MULTIPLE of the dynamical timestep - since most models stop that way.
+!
+! If we can advance the model for 6hour chunks, for example -
+!
+! Also : All observations +/- half this timestep are assimilated.
+! In essence, this defines the minimum window used for assimilation.
+
+type(time_type) :: set_model_stride
+
+! Check the user input
+if ((assimilation_period_seconds < 0) .or. (assimilation_period_days < 0)) then
+ write(msgstring1,*)'model_nml:assimilation_period_[seconds,days] must both be positive.'
+ write(msgstring2,*)'they are : ', assimilation_period_seconds, assimilation_period_days
+ call error_handler(E_ERR,'set_model_stride',msgstring1,source,revision,revdate,text2=msgstring2)
+elseif ((assimilation_period_seconds == 0) .and. (assimilation_period_days == 0)) then
+ write(msgstring1,*)'at least one of model_nml:assimilation_period_[seconds,days] must be positive.'
+ write(msgstring2,*)'they are : ', assimilation_period_seconds, assimilation_period_days
+ call error_handler(E_ERR,'set_model_stride',msgstring1,source,revision,revdate,text2=msgstring2)
+else
+ ! FIXME ... ensure that stride is a multiple of 'time_step'
+ set_model_stride = set_time(assimilation_period_seconds, assimilation_period_days)
+endif
+
+end function set_model_stride
+
!========================================================================
! End of model_mod
!========================================================================
end module model_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
Modified: DART/branches/development/models/sqg/spectral_mod.f90
===================================================================
--- DART/branches/development/models/sqg/spectral_mod.f90 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/spectral_mod.f90 2013-05-17 19:38:23 UTC (rev 6149)
@@ -1,16 +1,7 @@
!============================================================
-! <next few lines under version control, D O N O T E D I T>
-! $Date$
-! $Author$
-! $Revision$
-! $Id$
-!============================================================
-
-!============================================================
! a module to be used with spectral surface-based models.
!============================================================
-!============================================================
MODULE spectral_mod
implicit none
@@ -77,4 +68,12 @@
type(derivative_operators) :: d_oper
END MODULE spectral_mod
+
!============================================================
+
+! <next few lines under version control, do not edit>
+! $Date$
+! $Author$
+! $Revision$
+! $Id$
+
Modified: DART/branches/development/models/sqg/sqg.f90
===================================================================
--- DART/branches/development/models/sqg/sqg.f90 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/sqg.f90 2013-05-17 19:38:23 UTC (rev 6149)
@@ -1,11 +1,3 @@
-!============================================================
-! <next few lines under version control, D O N O T E D I T>
-! $Date$
-! $Author$
-! $Revision$
-! $Id$
-!============================================================
-
!========================================================================
PROGRAM sqg
@@ -19,3 +11,10 @@
END PROGRAM sqg
!========================================================================
+
+! <next few lines under version control, do not edit>
+! $Date$
+! $Author$
+! $Revision$
+! $Id$
+
Modified: DART/branches/development/models/sqg/sqg_mod.f90
===================================================================
--- DART/branches/development/models/sqg/sqg_mod.f90 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/sqg_mod.f90 2013-05-17 19:38:23 UTC (rev 6149)
@@ -1,16 +1,7 @@
-!============================================================
-! <next few lines under version control, D O N O T E D I T>
-! $Date$
-! $Author$
-! $Revision$
-! $Id$
-!============================================================
-
!========================================================================
! Uniform-PV two-surface QG+1 model in spectral form (Hakim 2000)
!========================================================================
-!========================================================================
MODULE sqg_mod
USE netcdf
@@ -213,7 +204,7 @@
real, dimension(mmax,nmax), intent(in) :: thbyB,thbyT
real, dimension(mmax,nmax), intent(in) :: ulinB,ulinT
real, intent(in) :: lam
- logical, intent(in) :: first,bot,top
+ logical, intent(in) :: first,bot,top
real, dimension(mmax,nmax), intent(out) :: thxBr,thxTr,thyBr,thyTr
real, dimension(mmax,nmax), intent(out) :: uBr,uTr,vBr,vTr
real, dimension(mmax,nmax), intent(out) :: sblre
@@ -1293,7 +1284,7 @@
if ( it .eq. 0 ) then
! Create a new NetCDF file
- call nc_check( nf90_create(output_file, NF90_CLOBBER .or. NF90_64BIT_OFFSET, ncid), 'write_diag', 'create ' // trim(output_file) )
+ call nc_check( nf90_create(output_file, ior(NF90_CLOBBER,NF90_64BIT_OFFSET), ncid), 'write_diag', 'create ' // trim(output_file) )
! Define dimensions
call nc_check( nf90_def_dim(ncid, "nx", 2*kmax, vardim(1)), 'write_diag', 'def_dim, nx ' // trim(output_file))
@@ -1425,8 +1416,8 @@
if ( verbose .gt. 1 ) print*,'max topo = ',maxval(abs(hspR))
call xy_to_sp(cmplx(hspR,0.),hspC,2*kmax,2*lmax,kmax,lmax)
if ( verbose .gt. 1 ) print*,'max topo spectral = ',maxval(abs(hspC))
- call invert(-hspC,Cblank,Rblank,Rblank,Rblank,Rblank,Rblank, &
- hv,Rblank,hu,Cblank,Cblank,Rblank,Rblank,Rblank, &
+ call invert(-hspC,Cblank,Rblank,Rblank,Rblank,Rblank,Rblank, &
+ hv,Rblank,hu,Cblank,Cblank,Rblank,Rblank,Rblank, &
Rblank,.TRUE.,.TRUE.,.TRUE.,lam,Cblank,Cblank,Rblank)
! hu and hv have the tropo winds due to topography
if ( verbose .gt. 1 ) print*,'max tropo winds due to topography: ',maxval(abs(hu)),maxval(abs(hv))
@@ -1512,3 +1503,10 @@
END MODULE sqg_mod
!========================================================================
+
+! <next few lines under version control, do not edit>
+! $Date$
+! $Author$
+! $Revision$
+! $Id$
+
Modified: DART/branches/development/models/sqg/work/input.nml
===================================================================
--- DART/branches/development/models/sqg/work/input.nml 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/work/input.nml 2013-05-17 19:38:23 UTC (rev 6149)
@@ -1,7 +1,8 @@
&perfect_model_obs_nml
start_from_restart = .true.,
output_restart = .true.,
- async = 1,
+ async = 0,
+ adv_ens_command = "model called as a subroutine",
init_time_days = 145731,
init_time_seconds = 0,
first_obs_days = -1,
@@ -13,7 +14,6 @@
restart_out_file_name = "perfect_restart",
obs_seq_in_file_name = "obs_seq.in",
obs_seq_out_file_name = "obs_seq.out",
- adv_ens_command = "./advance_model.csh",
output_timestamps = .false.,
trace_execution = .false.,
output_forward_op_errors = .false.,
@@ -22,15 +22,17 @@
/
&model_nml
- output_state_vector = .false.,
- channel_center = 45.0,
- channel_width = 40.0,
- debug = .false.,
+ output_state_vector = .false.,
+ channel_center = 45.0,
+ channel_width = 40.0,
+ assimilation_period_days = 0,
+ assimilation_period_seconds = 21600,
+ debug = .false.
/
&filter_nml
- async = 1,
- adv_ens_command = "./advance_model.csh",
+ async = 0,
+ adv_ens_command = "model called as a subroutine",
ens_size = 25,
start_from_restart = .false.,
output_restart = .true.,
Modified: DART/branches/development/models/sqg/write_sqg_restart.f90
===================================================================
--- DART/branches/development/models/sqg/write_sqg_restart.f90 2013-05-17 17:59:34 UTC (rev 6148)
+++ DART/branches/development/models/sqg/write_sqg_restart.f90 2013-05-17 19:38:23 UTC (rev 6149)
@@ -1,11 +1,3 @@
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! <next few lines under version control, D O N O T E D I T>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
!============================================================
! read in the model state from a netCDF file and write out
! to a DART readable restart file
@@ -25,7 +17,6 @@
implicit none
integer :: model_unit, model_size, model_days, model_seconds
- integer :: j
real(r8), allocatable, dimension(:) :: model_state
real, allocatable, dimension(:,:) :: thxyB, thxyT
@@ -48,19 +39,16 @@
print*, 'Enter model time in DAYS, SECONDS [eg. 0 0]: '
read*, model_days, model_seconds
- ! open a output file to write to:
model_unit = open_restart_write('sqgRestart')
- ! set model time in the output restart file:
model_time = set_time(model_seconds, model_days)
- ! get model size:
model_size = get_model_size()
- ! get model static data:
sqg_static = get_model_static_data()
! allocate space for all variables:
allocate( model_state(model_size) )
allocate( theta(2*kmax,2*lmax,2) )
- allocate( thxyB(2*kmax,2*lmax) ) ; allocate( thxyT(2*kmax,2*lmax) )
+ allocate( thxyB(2*kmax,2*lmax) )
+ allocate( thxyT(2*kmax,2*lmax) )
! read in the theta fields from the netCDF file:
call init(infile,thxyB,thxyT)
@@ -71,7 +59,8 @@
endif
! convert theta into DART state-vector:
- theta(:,:,1) = thxyB ; theta(:,:,2) = thxyT
+ theta(:,:,1) = thxyB
+ theta(:,:,2) = thxyT
call sqg_to_dart(theta,model_state)
! write to the file:
@@ -135,3 +124,10 @@
end subroutine toggle_base_state
end program write_sqg_restart
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
More information about the Dart-dev
mailing list