!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ASSIMCON ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE assimcon(nx,ny,nz,x,y,z,zp, & 1,3
ubar,vbar,pbar,ptbar,rhostr,qvbar, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
j1,j2,j3, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinates the input of control parameter options for a thermo-
! dynamic recovery and variational adjustment. The following is a
! summary of the available options and a list of possible combinations.
!
! recovopt Dynamic retrieval option:
!
! = 0, Do NOT retrieve p' or T'
! = 1, Retrieve p' and T' (Implied insertion of the
! interpolated velocities at the
! recovery time.)
!
! varopt Variational adjustment option:
!
! = 0, NO variational adjustment
! = 1, Perform variational adjustment
!
! insrtopt Insertion option:
!
! = 0, Do NOT insert velocities (NOTE: velocities will be inserted
! in conjunction with dynamic
! retrieval if recovopt=1)
! = 1, Insert velocities (At all available data times -
! See NOTE below.)
!
! A matrix describing these assimilation options is as follows:
!
! recovopt varopt insrtopt Comments
! ----------------------------------------------------------------------
! 0 0 0 Nothing - Exit assimilation mode
! 0 0 1 Direct insertion Only
! 0 1 0 NOT ALLOWED - Exit assimilation mode
! 0 1 1 Var. adj. + direct insertion (no recovery)
! 1 0 1 Recovery + direct insertion (no var. adj.)
! 1 1 1 Recovery, adj., & direct insertion
! 1 1 0 Recovery + var. adj. (no direct insertion)
! 1 0 0 Recovery only
!
! NOTE: The insertion option (insrtopt=1) implies that the velocities
! are to be inserted directly into the model stream at the model
! time corresponding to that of the data. This inserted velocity
! data may (varopt=1) or may not (varopt=0) be 'adjusted'.
! When a recovery is performed (recovopt=1), the 3 input data
! files (at times t1,t2,& t3) are interpolated to model time
! steps m1,m2,& m3. In effect, this is also an insertion since
! the model's prognostic variables at these times are overwritten.
! Thus, even if the insertion option is not chosen (i.e., insrtopt=0)
! there will still be an 'insertion' at the recovery time if
! the recovery option is turned on (recovopt=1).
!
! t1 t2 t3
! --x-----------|--x--|-----------x--
! m1 m2 m3
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Steven Lazarus and Alan Shapiro
!
! 2/24/1993.
!
! MODIFICATION HISTORY:
!
! 03/08/96 (Limin Zhao)
! Added an adjustment to the location of radar station so that it
! will never be on grid points.
!
! 03/19/96 (Limin Zhao)
! Fixed a bug on assimtim. They need be saved after they are read
! in first time.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space (m)
!
! OUTPUT:
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! pbar Base state pressure (Pascal)
! ptbar Base state potential temperature (K)
! rhostr Base state density (kg/m**3) * j3
! qvbar Base state water vapor specific humidity (kg/kg)
!
! u x velocity component at 3 time levels
! v y velocity component at 3 time levels
! w z velocity component at 3 time levels
!
! ptprt Perturbation potential temperature at 3 time levels
! pprt Perturbation pressure at 3 time levels
! qv Water vapor specific humidity (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rain water mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
! NOTE: Since ASSIMCON reads input data file headers only, the
! following are dummy arrays.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: nt ! The no. of t-levels of t-dependent arrays
PARAMETER (nt=3)
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: x (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z (nz) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL :: ubar (nx,ny,nz) ! Base state x-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state y-velocity (m/s)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: rhostr(nx,ny,nz) ! Base state air density (kg/m**3) times j3
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity (kg/kg)
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg)
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian defined as
! - d( zp )/d( x )
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian defined as
! - d( zp )/d( y )
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian defined as
! d( zp )/d( z )
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
REAL :: tem3 (nx,ny,nz) ! Temporary work array.
REAL :: tem4 (nx,ny,nz) ! Temporary work array.
REAL :: tem5 (nx,ny,nz) ! Temporary work array.
REAL :: tem6 (nx,ny,nz) ! Temporary work array.
REAL :: tem7 (nx,ny,nz) ! Temporary work array.
REAL :: tem8 (nx,ny,nz) ! Temporary work array.
REAL :: tem9 (nx,ny,nz) ! Temporary work array.
REAL :: tem10 (nx,ny,nz) ! Temporary work array.
REAL :: tem11 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: ichk ! Flag indicating whether the next data set is
SAVE ichk ! available (0=not available,1=available)
DATA ichk /1/
INTEGER :: iassim ! Flag indicating whether model time (curtim)
SAVE iassim ! exceeds the input data time (assimtim).
DATA iassim /1/ ! (0=yes, exit assimilation mode)
INTEGER :: icount ! Flag controling the reading of the batch file.
SAVE icount ! (0=do not read batch file,1=read batch file).
DATA icount /1/
INTEGER :: iread ! Counter indicating the input data file and data
SAVE iread ! time. Used to flag whether a variational adj.
DATA iread /1/ ! is to be performed this time step.
INTEGER :: igo ! Counter indicating the input data file and data
SAVE igo ! time. Used to flag whether a direct insertion
DATA igo /1/ ! is to be performed this time step.
INTEGER :: dtcount
SAVE dtcount
DATA dtcount /0/
INTEGER :: isrc ! Flag coordinating the data ingest
INTEGER :: istat ! Status of input data file
INTEGER :: itag ! If test flag (itag=ii,igo, or iread)
REAL :: assimtim(100) ! Time of (all) input data files
SAVE assimtim
REAL :: t1 ! Time of first data file
REAL :: t2 ! Time of second data file
REAL :: t3 ! Time of third data file
REAL :: temscl ! Grid scale used to calculate cdvdmp
REAL :: radloc
REAL :: stor1,stor2,stor3,stor4
INTEGER :: stor5
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model execution
INCLUDE 'assim.inc' ! Recovery/assimilation parameters
INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
! Routines called:
!
!-----------------------------------------------------------------------
!
! external assimbat
EXTERNAL assimrd
EXTERNAL retrint
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! If iassim is zero, then one of four things have occured:
! 1) The next available data file to be used in the assimilation
! arrived late (i.e., after the current model time exceeds that
! of the input data). Exit assimilation mode.
! 2) The wrong combination of options for recovopt,velocity,and
! insrtopt were chosen (see below). Exit assimilation mode.
! 3) The number of files (nvf) has been exceeded.
! 4) The assimilation mode has been turned off (recovopt,varopt,
! instopt= 0).
!
!-----------------------------------------------------------------------
!
WRITE(6,*)'code in assimcon, and curtim is ',curtim
IF ( iassim == 0) RETURN
!
!-----------------------------------------------------------------------
!
! Enter this IF block if this is the first entry into ASSIMCON
! (icount=1).
!
!-----------------------------------------------------------------------
!
IF ( icount == 1) THEN
!
!-----------------------------------------------------------------------
!
! Read in user-supplied recovery/assimilation parameters (ASSIMBAT.F)
! upon first entry into ASSIMCON only (icount=1). The parameter
! icount is then set to zero to avoid reading the batch file on
! subsequent visits to this routine.
!
!-----------------------------------------------------------------------
!
! CALL initassim(nx,ny,nz)
icount = 0
!
!-----------------------------------------------------------------------
!
! Check if the radar station located on grid point. If so, shift
! a small distance. Here we choose to shift (grid space/20).
!
!-----------------------------------------------------------------------
!
! radloc = xshift - dx*int(xshift/dx)
! IF (radloc.eq.0) THEN
radloc = xshift - dx*nint(xshift/dx)
IF (radloc < dx*1.0E-3) THEN
xshift = xshift + 0.05*dx
WRITE(6,*)'radar locates at x grid point'
END IF
! radloc = yshift - dy*int(yshift/dy)
! IF (radloc.eq.0) THEN
radloc = yshift - dy*nint(yshift/dy)
IF (radloc < dy*1.0E-3) THEN
yshift = yshift + 0.05*dy
WRITE(6,*)'radar locates at y grid point'
END IF
! radloc = zshift - dz*int(zshift/dz)
! IF (radloc.eq.0) THEN
radloc = zshift - dz*nint(zshift/dz)
IF (radloc < dz*1.0E-3) THEN
zshift = zshift + 0.05*dz
WRITE(6,*)'radar locates at z grid point'
END IF
WRITE(6,*)'adjusted radar locations xshift,yshift,zshift=', &
xshift,yshift,zshift
!
!-----------------------------------------------------------------------
!
! Check flags to see what combination of assimilation options
! have been chosen. If an unavailable option is chosen or the
! assimilation mode has been turned off print message and exit.
! (See table in Purpose Section above)
!
!-----------------------------------------------------------------------
!
IF ( recovopt == 0.AND.insrtopt == 0) GO TO 100
!
!-----------------------------------------------------------------------
!
! Reset nstep so that the I/O corresponds to the restart time and not
! t=0. This is a temporary fix until the option is introduced into
! ARPS proper.
!
!-----------------------------------------------------------------------
!
nstep = INT( ((curtim-tstart)+0.1*dtbig)/dtbig ) + 1
PRINT *,'in ASSIMCON, nstep =',nstep
!
!-----------------------------------------------------------------------
!
! Initialize the counter, ii, to 3. ii is declared in the common
! block of ASSIM.INC and is used by this routine and ASSIMRD.F only.
! NOTE: ii is updated in ASSIMCON only.
!
!-----------------------------------------------------------------------
!
ii = 3
!
!-----------------------------------------------------------------------
!
! Read three input data file headers. This is done so that the
! recovery can access the times of these data files. The time
! for each file is brought back via the array 'assimtim' which
! is then used to set flags to determine when to perform a
! variational adjustment, and/or direct insertion and/or the first
! recovery.
!
! NOTE: This package assumes that AT LEAST 3 input data files exist
! prior to model execution.
!
! The flag isrc coordinates the ingestion of data, where:
!
! isrc = 1, Flags ASSIMRD to read 3 file headers only
! isrc = 2, Flags ASSIMRD to read 1 file header only
! isrc = 3, Flags ASSIMRD to read 1 complete data file
! (called in ASSIMVEL.F)
! isrc = 4, Flags ASSIMRD to read 3 complete data files
!
! Note: tem11 is used to read rhobar
!
!-----------------------------------------------------------------------
!
isrc = 1 ! Flags assimrd to read 3 file headers only
itag = ii ! Counter indicating the input file name and time
PRINT *, 'ASSIMCON: CALL assimrd_1 ', itag
CALL assimrd
(nx,ny,nz,x,y,z,zp, &
isrc,itag,istat,assimtim, &
ubar,vbar,pbar,ptbar,rhostr,qvbar,tem11, &
u,v,w,pprt,ptprt,qv,qc,qr,qi,qs,qh,j1,j2,j3, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,tem10)
IF (istat /= 0) THEN
iassim = 0
WRITE(6,'(/5x,a,/5x,a,/5x,a,/5x,a)') &
'WARNING: The assimilation package assumes that 3 input', &
'data files exist prior to the start of the assimilation.', &
'Please provide at least 3 files and restart assimilation.', &
'Exit assimilation mode.'
RETURN
END IF
END IF
!
!-----------------------------------------------------------------------
!
! If ichk=0, then the data file needed for the next 'assimilation'
! does not yet exist and thus we set irecov, ivar, and insrt to zero
! during that model time step. This is done since irecov, ivar, &
! insrt cannot be flagged using an 'unknown' time (assimtim).
! If the data file 'arrives' late (i.e., after the model time
! exceeds that of the input data file), the parameter iassim is set
! to 0 and the assimilation mode is turned off for the remainder of
! the simulation. If the data arrive on time, ichk is set to 1, and
! we proceed below to determine whether a recovery/assimilation is
! to be performed during this model time step.
!
!-----------------------------------------------------------------------
!
!
! assimtim(1) = 0.0 ! for 18:36, 18:41, 18:47Z
! assimtim(2) = 350.0
! assimtim(3) = 701.0
! assimtim(1) = 30.0 ! for SSW NOV 1998 tests
! assimtim(2) = 378.0 ! comment out 12-5-98 SSW
! assimtim(3) = 732.0
WRITE(6,*)'assim time are:',assimtim(1),assimtim(2),assimtim(3)
IF ( ichk == 0) THEN
irecov = 0
ivar = 0
insrt = 0
IF ( ii > nvf.OR.igo > nvf.OR.iread > nvf) iassim = 0
ELSE
!
!-----------------------------------------------------------------------
!
! If recovopt equals 1, set flags to perform a retrieval.
!
!-----------------------------------------------------------------------
!
IF ( recovopt == 1) THEN
!
!-----------------------------------------------------------------------
!
! Determine whether a recovery/assimilation is to be performed
! during this model time step.
!
!-----------------------------------------------------------------------
!
WRITE(6,*) 'curtim = ',curtim
WRITE(6,*) 'assim time = ',assimtim(1),assimtim(2),assimtim(3)
IF (curtim < assimtim(ii).OR.ii-1 >= nvf) THEN
irecov = 0 ! Flag tells RETRPTPR.F not to perform a recovery
ELSE
ii = ii + 1
irecov = 1 ! Flag tells RETRPTPR.F to perform a recovery
END IF
END IF
!
!-----------------------------------------------------------------------
!
! If varopt equals 1, set flags to perform a variational adjustment.
!
!-----------------------------------------------------------------------
!
IF (varopt == 1) THEN
!
!-----------------------------------------------------------------------
!
! Determine whether a variational adjustment is to be performed
! during this model time step.
!
!-----------------------------------------------------------------------
!
IF (curtim < assimtim(iread) .OR. iread-1 >= nvf) THEN
ivar = 0 ! Flag tells ASSIMVEL.F not to perform a variational adjustment
ELSE
iread = iread + 1
ivar = 1 ! Flag tells ASSIMVEL.F to perform a variational adjustment
END IF
END IF
!
!-----------------------------------------------------------------------
!
! If insrtopt equals 1, set flags to perform a direct insertion.
!
!-----------------------------------------------------------------------
!
IF (insrtopt == 1) THEN
!
!-----------------------------------------------------------------------
!
! Determine whether a direct insertion is to be performed
! during this model time step.
!
!-----------------------------------------------------------------------
!
IF (curtim < assimtim(igo) .OR. igo-1 >= nvf) THEN
insrt = 0 ! Flag tells ASSIMVEL.F not to perform an insertion
ELSE
igo = igo + 1
insrt = 1 ! Flag tells ASSIMVEL.F to perform an insertion
END IF
END IF
END IF
WRITE(6,*) 'irecov,ivar,insrt: ',irecov,ivar,insrt
WRITE(6,*) 'ii,iread,igo: ',ii,iread,igo
!
!-----------------------------------------------------------------------
!
! Read the next available data file. This is done in order to set the
! flags, irecov, ivar, and insrt for the next recovery, variational
! adjustment, and/or insertion. Since it is assumed that at least 3
! data files exist prior to the assimilation period, we do not read
! the 4th file header until the model time exceeds the time of the
! 3rd file (t3).
!
! t1 t2 t3 t4
! --x--------------x--------------x--------------x--
! |
! Begin search for 4th File
!
! Similarly, we do not begin looking for the next input data file
! header [at time, t(n)] until the model time exceeds the time of the
! preceeding input data file [at time, t(n-1)].
!
! t(n-1) t(n)
! --x--------------x--
! |
! Begin search for nth File
!
!-----------------------------------------------------------------------
!
IF ( ii > 3.OR.iread > 3.OR.igo > 3) THEN
IF ( irecov == 1.OR.ivar == 1.OR.insrt == 1.OR.ichk == 0) THEN
isrc = 2 ! Flags assimrd to read 1 file header only
itag = MAX(ii,iread,igo)
PRINT *, 'ASSIMCON: CALL assimrd_2 ', &
itag, irecov, ivar, insrt
CALL assimrd(nx,ny,nz,x,y,z,zp, &
isrc,itag,istat,assimtim, &
ubar,vbar,pbar,ptbar,rhostr,qvbar,tem11, &
u,v,w,pprt,ptprt,qv,qc,qr,qi,qs,qh,j1,j2,j3, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,tem10)
IF (istat /= 0 .OR. itag > nvf) THEN
ichk = 0
ELSE
ichk = 1
!
!-----------------------------------------------------------------------
!
! Check to ensure that the input of the next data file occurs no more
! than one model time step (dtbig) past the next recovery time
! assimtim(ii). If it does, exit assimilation mode by setting
! iassim = 0.
! NOTE: This precludes any further assimilation even if subsequent
! data files become available.
!
!-----------------------------------------------------------------------
!
IF (curtim-assimtim(itag) > dtbig) THEN
iassim = 0
irecov = 0
ivar = 0
insrt = 0
WRITE(6,'(/5x,a,/5x,a,/5x,a)') &
'WARNING: The successful read of the input data file occured', &
'after the insertion/adjustment/recovery window', &
'Exit assimilation mode.'
RETURN
END IF
END IF
END IF
END IF
RETURN
100 WRITE(6,'(/5x,a,/5x,a,/10x,a,i2,/10x,a,i2,/10x,a,i2,/5x,a,/5x,a)') &
'WARNING: You have either chosen an option that is not available', &
'or have turned the assimilation mode off. The input option was', &
'recovopt=',recovopt,'varopt=',varopt,'insrtopt=',insrtopt, &
'See ASSIMCON for appropriate parameters.', &
'Exit assimilation mode.'
iassim = 0
RETURN
END SUBROUTINE assimcon