!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ASSIMRD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE assimrd(nx,ny,nz,x,y,z,zp, & 4,21
isrc,itag,ireturn,assimtim, &
ubar,vbar,pbar,ptbar,rhostr,qvbar,rhobar, &
u,v,w,pprt,ptprt,qv,qc,qr,qi,qs,qh,j1,j2,j3, &
tem1,tem5,tem6,tem7,tem8,tem9,tem10,tem11, &
tem12,tem13)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinates the ingestion of history data dump files in various
! formats and the ingestion of 'pre-processed' Lincoln Lab radar data.
! Other data types can be used by selecting the 'other' option for
! the 'dtyp' parameter in ASSIM.INPUT. The velocity fields can be
! used for an insertion, adjustment and/or recovery (see ASSIMCON.F).
!
! NOTE: In order to use the 'other' option, corresponding changes
! are needed in this routine. In this version, only two type
! of data are processed.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Alan Shapiro and Steven Lazarus
!
! 2/25/1993
!
! MODIFICATION HISTORY:
!
! 02/22/96 (Limin Zhao)
! Changes are made to match the new version of assimilation package.
!
! 02/26/96 (Steve Lazarus & Limin Zhao)
! Added the simple moisture parameterization (Kessler formula) for
! qr and qv.
!
! 02/27/96 (Limin Zhao)
! A critical value of qr is used to set flags at the area outside of
! rainwater regions for OSSE type of experiments.
!
! 02/29/96 (Limin Zhao)
! Added two temporary arrays to avoid overwritten of pbar and ptbar,
! which will be used in thermodynamic retrieval.
!
! 05/08/96 (Limin Zhao)
! A bug was fixed when control input parameters are set
! varopt=0,insrtopt=0 and recovopt =1.
! A correct filename is given now.
!
! 05/14/96 (Limin Zhao and Alan Shapiro)
! Modified the stratage to get qv by using background temperature
! and pressure information.
!
! 05/15/96 (Limin Zhao)
! A control parameter is hardwired to control the moisture retrieval
! switch on or off.
!
! 09/24/96 (Limin Zhao)
! Fix a bug in qr reading with isrc=4. qr should be read in for
! thermodynamic retrieval.
!
!-----------------------------------------------------------------------
!
! 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)
!
! isrc Flag, indicating the calling routine
! itag Counter indicating file name and time
!
! OUTPUT:
!
! ireturn Flag, indicating read status of input file
! = 0 Successful read
! = 1 Unsuccessful read
!
! assimtim Times of all input files
!
! 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) times j3
! qvbar Base state water vapor specific humidity (kg/kg)
!
! u x component of velocity at a given time level (m/s)
! v y component of velocity at a given time level (m/s)
! w Vertical component of Cartesian velocity at a given
! time level (m/s)
! pprt Perturbation pressure at times tpast and tpresent (Pascal)
! ptprt Perturbation potential temperature at times tpast and
! tpresent (K)
!
! 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
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem5 Temporary work array.
! tem6 Temporary work array.
! tem7 Temporary work array.
! tem8 Temporary work array.
! tem9 Temporary work array.
! tem10 Temporary work array.
! tem11 Temporary work array.
! tem12 Temporary work array.
! tem13 Temporary work array.
!
! rhobar Base state density (kg/m**3)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nt ! Number of time levels of time-dependent arrays.
INTEGER :: tim ! Index of time level.
INTEGER :: tpast ! Index of time level for the past time.
INTEGER :: tpresent ! Index of time level for the present time.
INTEGER :: tfuture ! Index of time level for the future time.
PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=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.
INTEGER :: isrc ! Flag indicating source of calling routine
INTEGER :: itag ! Counter indicating file name and time
INTEGER :: nstyps
REAL :: assimtim(100) ! Time of data input files
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 :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3)
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 :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
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 -d(zp)/d(x)
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian -d(zp)/d(y)
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian d(zp)/d(z)
REAL :: tem1 (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
REAL :: tem12 (nx,ny,nz) ! Temporary work array
REAL :: tem13 (nx,ny,nz) ! Temporary work array
INTEGER :: ireturn
REAL :: frsvnth, coef
REAL :: rad,xs,ys,zs
REAL :: xmove,ymove
REAL :: plcl
INTEGER :: in,jn
INTEGER :: klcl
INTEGER :: itimesv
!
! real dummy(153,99,43) !Note: temporary fix for extra arrays
! integer idummy(153,99,43)
!
! real dummy(133,133,37) !Note: temporary fix for extra arrays
! integer idummy(133,133,37)
!
!-----------------------------------------------------------------------
!
! Routines called:
!
!-----------------------------------------------------------------------
!
EXTERNAL dtahead
EXTERNAL dtaread
EXTERNAL radread
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: is ! Counter, used for reading input data
SAVE is
DATA is/1/
INTEGER :: i,j,k,n ! Loop index
INTEGER :: nchanl ! FORTRAN I/O channel number for history data output.
INTEGER :: lengbf,lendtf,lbasdmpf
REAL :: xor ! Coordinate (xor,yor,zor) of the model grid
REAL :: yor ! origin relative to the radar.
REAL :: zor !
REAL :: storstop ! Temporarily stores the model stop time
CHARACTER (LEN=128 ) :: saverunm ! Temporarily stores the name of this run
CHARACTER (LEN=128 ) :: filein ! Input file name for recovery
REAL :: temp,pres,qvsat
INTEGER :: imoist
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'assim.inc'
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! If isrc=1, then read the times (assimtim) from 3 data file headers.
! This is done in order to determine when either a recovery,
! adjustment, and/or an insertion should be performed. It is
! assumed that 3 data files exist at the beginning of the
! assimilation period. For this reason the counter, itag, is
! initially set to 3 as well as updated in ASSIMCON.F.
!
! The grid/base state file (inigbf) and input file names (assimdat)
! are to be provided by the user in ASSIM.INPUT. The grid/base
! state file is used when model data is ingested only.
!
! NOTE: None of the time-dependent variables are read in during
! this step.
!
!-----------------------------------------------------------------------
!
WRITE(6,*) 'code in assimrd: is,isrc=',is,isrc
IF( isrc == 1) THEN ! Called by ASSIMCON
DO i=1,itag
tim = i
lendtf=LEN(assimdat(i))
CALL strlnth
( assimdat(i), lendtf)
WRITE(6,*) 'The filename assimdat(',i,')=', &
assimdat(i)(1:lendtf)
!c
!clz
!c The option for dtyp=0 is not well considered for SOP, It will
!c be modified in future
!c
IF( dtyp == 0) THEN
saverunm = runname
CALL dtahead(nx,ny,nz, &
inifmt,inigbf,lengbf,assimdat(i), &
lendtf,assimtim(i),x,y,z,zp,tem11,tem12,tem13, &
ptprt(1,1,1,tim),pprt(1,1,1,tim),tem7, &
qc(1,1,1,tim),qr(1,1,1,tim),qi(1,1,1,tim), &
qs(1,1,1,tim),qh(1,1,1,tim),tem9, &
ubar,vbar,tem8,ptbar,pbar,rhobar,qvbar, &
ireturn,tem1,tem9,tem10)
runname = saverunm
IF(ireturn /= 0) RETURN
ELSE IF(dtyp == 1) THEN
CALL radread(assimdat(i),lendtf, &
isrc,assimtim(i),ireturn, &
nx,ny,nz,dx,dy,dz, &
xor,yor,zor, &
tem1,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11,tem12,tem13)
IF(ireturn /= 0) RETURN
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Read in a single header to recover the time. For recovopt=1
!
! NOTE: None of the time-dependent variables are read in during
! this step.
!
!-----------------------------------------------------------------------
!
ELSE IF(isrc == 2) THEN ! Called by ASSIMCON
tim = tpresent
lendtf=LEN(assimdat(itag))
CALL strlnth( assimdat(itag), lendtf)
WRITE(6,*) 'The filename assimdat(',itag,')=', &
assimdat(itag)(1:lendtf)
IF( dtyp == 0) THEN
saverunm = runname
CALL dtahead(nx,ny,nz, &
inifmt,inigbf,lengbf,assimdat(itag), &
lendtf,assimtim(itag),x,y,z,zp,tem11,tem12,tem13, &
ptprt(1,1,1,tim),pprt(1,1,1,tim),tem7, &
qc(1,1,1,tim),qr(1,1,1,tim),qi(1,1,1,tim), &
qs(1,1,1,tim),qh(1,1,1,tim),tem9, &
ubar,vbar,tem8,ptbar,pbar,rhobar,qvbar, &
ireturn,tem1,tem9,tem10)
runname = saverunm
IF( ireturn /= 0) RETURN
ELSE IF(dtyp == 1) THEN
CALL radread(assimdat(itag),lendtf, &
isrc,assimtim(itag),ireturn, &
nx,ny,nz,dx,dy,dz, &
xor,yor,zor, &
tem1,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11,tem12,tem13)
IF(ireturn /= 0) RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Read in a single data file. The radial velocities are used
! to perform a variational adjustment and/or insertion.
!
! NOTE: If model data is ingested (dtyp=0) all current model
! prognostic variables are overwritten with the exeption of
! u,v,w and qv. uprt,vprt,wprt, and qvprt are returned in
! tem4,tem5,tem6,and tem7 respectively. The total fields
! can then reconstructed from these and their respective base
! state quantities. qv is reconstructed here.
!
!-----------------------------------------------------------------------
!
ELSE IF(isrc == 3) THEN ! Called by ASSIMVEL
tim = tpresent
IF( dtyp == 0) THEN
lengbf=LEN(inigbf)
CALL strlnth( inigbf, lengbf)
WRITE(6,'(/a,a)') &
'The grid/base state file name is ',inigbf(1:lengbf)
lendtf=LEN(assimdat(is))
CALL strlnth( assimdat(is), lendtf)
WRITE(6,'(/a,a)')'The file name is ',assimdat(is)(1:lendtf)
saverunm = runname
!
!-----------------------------------------------------------------------
!
! Set the tem8, tem5 and tem6 arrays to zero. This is done to avoid
! overwriting wbar, pbar and ptbar on calls to dtaread.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem8(i,j,k) = 0.0
tem5(i,j,k) = 0.0
tem6(i,j,k) = 0.0
END DO
END DO
END DO
storstop = tstop ! Temporary. For ingestion of model data
! generated prior to 4.0
CALL ctim2abss( year, month, day, hour, minute, second, itimesv )
CALL dtaread(nx,ny,nz,nstyps, &
inifmt,nchanl,inigbf,lengbf,assimdat(is), &
lendtf,assimtim(is),x,y,z,zp,tem11,tem12,tem13, &
ptprt(1,1,1,tim),pprt(1,1,1,tim),tem7, &
qc(1,1,1,tim),qr(1,1,1,tim),qi(1,1,1,tim), &
qs(1,1,1,tim),qh(1,1,1,tim),tem9,tem9,tem9, &
ubar,vbar,tem8,ptbar,pbar,rhobar,qvbar, &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
ireturn,tem6,tem9,tem10)
CALL abss2ctim( itimesv, year, month, day, hour, minute, second )
tstop = storstop
runname = saverunm
IF(ireturn /= 0) RETURN
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
rhostr(i,j,k) = rhobar(i,j,k)*j3(i,j,k)
qv(i,j,k,tim) = tem7(i,j,k) + qvbar(i,j,k)
END DO
END DO
END DO
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
tem11(i,j,k) = tem11(i,j,k) + ubar(i,j,k)
tem12(i,j,k) = tem12(i,j,k) + vbar(i,j,k)
tem13(i,j,k) = tem13(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
!
! Interpolate input 'model' velocity components to scalar grid point.
! Where:
! tem7 = scalar u component
! tem8 = scalar v component
! tem9 = scalar w component
!
!-----------------------------------------------------------------------
!
CALL avgx(tem11, 0, nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1,tem7)
CALL avgy(tem12, 0, nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1,tem8)
CALL avgz(tem13, 0, nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1,tem9)
!
!-----------------------------------------------------------------------
!
! Calculate the radial velocity component of the input model data.
! This is for OSSE type experiments. The input model data may or
! may not match that of the current model stream.
!
! WARNING: The insertion below assumes that, the user has supplied
! the radar coordinates with respect to the lower left hand
! corner of the grid (See ASSIM.INPUT). The grid origin is
! assumed to be at (0,0,0).
!
! * (0,0,0) grid origin
! . .
! . . ymove
! . .
! radar (- ..........
! xmove
!
!
!-----------------------------------------------------------------------
!
umove=0.0
vmove=0.0
xmove= xshift - umove*(curtim-assimtim(1))
ymove= yshift - vmove*(curtim-assimtim(1))
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
xs = 0.5*(x(i)+x(i+1)) - xmove
ys = 0.5*(y(j)+y(j+1)) - ymove
zs = 0.5*(z(k)+z(k+1)) - zshift
rad = SQRT(xs**2+ys**2+zs**2)
!
!
! tem1(i,j,k) = (tem7(i,j,k)*xs + tem9(i,j,k)*ys
! : + tem9(i,j,k)*zs)/rad
tem1(i,j,k) = (tem7(i,j,k)*xs + tem8(i,j,k)*ys &
+ tem9(i,j,k)*zs)/rad
END DO
END DO
END DO
!
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
tem11(i,j,k) = tem11(i,j,k) - ubar(i,j,k)
tem12(i,j,k) = tem12(i,j,k) - vbar(i,j,k)
tem13(i,j,k) = tem13(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set flags outside of rainwater area.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
IF(qr(i,j,k,tim) < 0.0005) THEN
tem1(i,j,k) = spval
tem11(i,j,k) = spval
tem12(i,j,k) = spval
tem13(i,j,k) = spval
END IF
END DO
END DO
END DO
!
!clz for real data ingestion
!
ELSE IF(dtyp == 1) THEN
!cc
!cc
!ccNote: no overwritten for model variables.
!cc
!cc
lendtf=LEN(assimdat(is))
CALL strlnth( assimdat(is), lendtf)
WRITE(6,'(/a,a)')'The filename', assimdat(is)(1:lendtf)
CALL radread(assimdat(is),lendtf, &
isrc,assimtim(is),ireturn, &
nx,ny,nz,dx,dy,dz, &
xor,yor,zor, &
tem1,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11,tem12,tem13)
IF(ireturn /= 0) RETURN
!
!-----------------------------------------------------------------------
!
! Compute the rainwater from the radar reflectivity factor. See
! Kessler (1969). Reflectivity is assumed to be in dBz, i.e.,
! Z (mm**6/m**3), and dBz = 10log10(Z).
!
! Note: When reflectivity is too small or too large, qr might be
! reliable from the qr-Z relation, a upper limit and lower
! limit are set up here.
!
!
! Set qr max at 50 dbz for AMS99 testing (SSW 12-7-98)
!
!
!-----------------------------------------------------------------------
!
imoist=1
IF(imoist == 0) GO TO 911
WRITE(6,*) 'This run is with qr and qv'
coef = LOG(10.)/10.
frsvnth = 4./7.
DO k = 1,nz-1
DO j = 1,ny-1
DO i = 1,nx-1
IF (tem8(i,j,k) > 50.0.AND.tem8(i,j,k) /= spval) THEN
! IF (tem8(i,j,k).gt.75.0.and.tem8(i,j,k).ne.spval) THEN
! IF (tem8(i,j,k).gt.75.0.and.tem8(i,j,k).ne.spval) THEN
! tem8(i,j,k) = 75.0
tem8(i,j,k) = 50.0
END IF
! IF (tem8(i,j,k).ne.spval.and.tem8(i,j,k).gt.20.0) THEN
! IF (tem8(i,j,k).ne.spval.and.tem8(i,j,k).gt.0.0) THEN
! IF (tem8(i,j,k).ne.spval.and.tem8(i,j,k).gt.5.0) THEN
IF (tem8(i,j,k) /= spval) THEN
qr(i,j,k,tim) = (.001/rhobar(i,j,k))* &
( (1./17300.)*EXP(coef*tem8(i,j,k)) )**frsvnth
ELSE
qr(i,j,k,tim) = 0.0
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! For real data, dtyp=1, assume qv=qvs in rainwater regions with
! updraft. In downdraft area, model background mean is used.
!
! Calculate the saturation water vapor specific humidity,
! qvs, from pressure and potential temperature using Teten's formula.
!
!
! Comment out Limin's old qvsat for AMS99 testing (SSW 12-7-98)
!
!
!-----------------------------------------------------------------------
!
! CALL satmrpt(nx,ny,nz,pbar,ptbar,tem9) ! calculate qvs
! DO 195 k=1,nz-1
! DO 195 j=1,ny-1
! DO 195 i=1,nx-1
! pres = pbar(i,j,k)+pprt(i,j,k,tim)
! temp = (ptbar(i,j,k)+ptprt(i,j,k,tim))*((pres/p0)**rddcp)
!ccxxx IF(temp.ge.273.16) THEN
! qvsat=(380. / pres) * !Teten's formula
! : exp( 17.27 * (temp-273.15)/(temp-35.86) )
!ccxxx ELSE
!ccxxx qvsat=(380. / pres) *
!ccxxx : exp( 21.875* (temp-273.15)/(temp-7.5) )
!ccxxx ENDIF
!ccxxx
! IF(tem8(i,j,k).gt.20.0.and.tem8(i,j,k).ne.spval
! : .and.w(i,j,k,tim).ge.0.0) THEN
! qv(i,j,k,tim) = qvsat
! ELSE IF(tem8(i,j,k).gt.20.0.and.tem8(i,j,k).
! : ne.spval.and.w(i,j,k,tim).lt.0.0) THEN
! qv(i,j,k,tim) = 0.8*qvsat
! ENDIF
!195 CONTINUE
911 CONTINUE
END IF
is = is + 1
!
!-----------------------------------------------------------------------
!
! Read in the winds, which may (or may not) be variationally adjusted,
! at three time steps, i.e.
!
! t1 t2 t3
! ---|----------|----------|---
!
! The time interval between these files corresponds to the frequency
! of the radar data. The time between each radar observation does not
! have to be the same (see RETRINT.F).
! NOTE:
! The adjusted winds are dumped via the dump3d routine (see ASSIMVEL.F).
! Hence, for the recovery option (recovopt=1), the adjusted winds and
! associated model data at that time step are always written off in
!
! a model format as specified by the user in ARPSINIT.INPUT
!
!-----------------------------------------------------------------------
!
ELSE IF(isrc == 4) THEN ! Called by RETPTPR
tim = 1
lengbf=LEN(inigbf)
CALL strlnth( inigbf, lengbf)
WRITE(6,'(/a,a)') &
'The grid/base state file name is ',inigbf(1:lengbf)
!
!
!ccxxx
!ccxxx You have to hardwire the filenames if you choose to do
!ccxxx thermodynamic retrieval only
!ccxxx
!ccxxx adjdat(1) = 'assim4.bin000032'
!ccxxx adjdat(2) = 'assim4.bin000384'
!ccxxx adjdat(3) = 'assim4.bin000736'
!
! adjdat(1) = 'ADA_15.bin000000.9408171830'
! adjdat(2) = 'ADA_15.bin000000.9408171836'
! adjdat(3) = 'ADA_15.bin000000.9408171841'
!
! adjdat(1) = 'WANG15.bin000000.940817183548'
! adjdat(2) = 'WANG15.bin000000.940817183600'
! adjdat(3) = 'WANG15.bin000000.940817183612'
!
! adjdat(1) = 'KI1881.bin000000'
! adjdat(2) = 'KI1881.bin000351'
! adjdat(3) = 'KI1881.bin000702'
DO n = itag-3,itag-1
WRITE(6,*)'itag,n',itag,n
! IF (varopt.eq.0.and.insrtopt.eq.0) THEN
! lendtf=len(assimdat(n))
! CALL strlnth( assimdat(n), lendtf)
! filein = assimdat(n)(1:lendtf)
! write(6,'(/a,a)')
! : 'The input file name is ',filein
! ELSE
lendtf=LEN(adjdat(n))
CALL strlnth( adjdat(n), lendtf)
filein = adjdat(n)(1:lendtf)
WRITE(6,'(/a,a)') &
'The input file name is ',filein
! ENDIF
saverunm = runname
!
!-----------------------------------------------------------------------
!
! Set the tem8 array to zero. This is done to avoid overwriting
! wbar on calls to dtaread.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem8(i,j,k) = 0.0
tem7(i,j,k) = 0.0
tem6(i,j,k) = 0.0
END DO
END DO
END DO
storstop = tstop ! Temporary. For ingestion of model data
! generated prior to 4.0
CALL ctim2abss( year, month, day, hour, minute, second, itimesv )
CALL dtaread(nx,ny,nz,nstyps, &
hdmpfmt,nchanl,inigbf,lengbf,filein, &
lendtf,assimtim(n),x,y,z,zp,tem11,tem12,tem13, &
ptprt(1,1,1,tim),pprt(1,1,1,tim),tem7, &
qc(1,1,1,tim),tem6,qi(1,1,1,tim), &
qs(1,1,1,tim),qh(1,1,1,tim),tem1,tem1,tem1, &
ubar,vbar,tem8,ptbar,pbar,rhobar,qvbar, &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
tem1(1,1,1),tem1(1,1,1),tem1(1,1,1),tem1(1,1,1), &
ireturn,tem5,tem9,tem10)
CALL abss2ctim( itimesv, year, month, day, hour, minute, second )
!
!
! if( n.eq.1 ) assimtim(n) = 0.0 ! 18:30z
! if( n.eq.2 ) assimtim(n) = 360.0 ! 18:36z
! if( n.eq.3 ) assimtim(n) = 660.0 ! 18:41z
!
! if( n.eq.1 ) assimtim(n) = 0.0 ! 18:35:48z
! if( n.eq.2 ) assimtim(n) = 12.0 ! 18:36:00z
! if( n.eq.3 ) assimtim(n) = 24.0 ! 18:36:12z
!
! if( n.eq.1 ) assimtim(n) = 0.0 ! 18:35:36z for SDVR data
! if( n.eq.2 ) assimtim(n) = 351.0 ! 18:41:26z
! if( n.eq.3 ) assimtim(n) = 702.0 ! 18:47:17z
!
tstop = storstop
runname = saverunm
IF( ireturn /= 0) RETURN
!
!-----------------------------------------------------------------------
!
! The arrays tem11,tem12 and tem13 contain uprt, vprt, and wprt
! respectively. The recovery requires the total velocities,i.e.,
! u = uprt + ubar.
!
! Reconstruct qv from qvprt which is returned from the call to
! DTAREAD in the tem7 array.
!
! Construct rhostr from rhobar and j3.
!
!-----------------------------------------------------------------------
!
DO i = 1,nx
DO j = 1,ny-1
DO k = 1,nz-1
u(i,j,k,tim) = tem11(i,j,k) + ubar(i,j,k)
END DO
END DO
END DO
DO i = 1,nx-1
DO j = 1,ny
DO k = 1,nz-1
v(i,j,k,tim) = tem12(i,j,k) + vbar(i,j,k)
END DO
END DO
END DO
DO i = 1,nx-1
DO j = 1,ny-1
DO k = 1,nz
w(i,j,k,tim) = tem13(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qv(i,j,k,tim) = tem7(i,j,k) + qvbar(i,j,k)
qr(i,j,k,tim) = tem6(i,j,k)
END DO
END DO
END DO
tim = tim + 1
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
rhostr(i,j,k) = rhobar(i,j,k)*j3(i,j,k)
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE assimrd
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PALCL ######
!###### ######
!###### Steven Lazarus ######
!###### University of Oklahoma ######
!###### School of Meteorology ######
!##################################################################
!##################################################################
SUBROUTINE palcl(nx,ny,nz,in,jn,p,pienv,tinit,qinit,klcl,plcl)
!
!-----------------------------------------------------------------------
!
! Purpose:
!
! Subroutine calculates the lifting condensation level
!
!-----------------------------------------------------------------------
!
! Author: Steve Lazarus
! 6/10/93
!
! Modification History:
!
!-----------------------------------------------------------------------
!
! Input :
!
! nx,ny,nz Model grid dimensions
! in,jn A given point in the horizontally homogenous domain
! p Sounding pressure (Pa)
! pienv Sounding pressure (non-dimensional)
! tinit Sounding temperature (K)
! qinit Sounding mixing ratio (kg/kg)
!
! Cp Specific heat for dry air (kg/(m s**2))
! Rd Gas constant for dry air (kg/(m s**2))
! p0 Surface pressure (Pa)
!
! Output:
!
! klcl Index just above LCL
! plcl Pressure at LCL (Pa)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: pie, qvs
REAL :: qtest
REAL :: plcl
REAL :: plcl1
REAL :: pilcl
REAL :: pilcl1
REAL :: pavg
REAL :: piavg
INTEGER :: in,jn
INTEGER :: nx,ny,nz
REAL :: p(nx,ny,nz),pienv(nx,ny,nz),tinit(nx,ny,nz),qinit(nx,ny,nz)
INTEGER :: k,klcl
INTEGER :: iter
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Calculate non-dimensional pressure, pienv.
!
!-----------------------------------------------------------------------
!
PRINT *,'Entering subroutine palcl'
DO k=1,nz
pienv(in,jn,k) = pie( END DO
k=1
qtest = 999.
10 CONTINUE
IF(qtest > qinit(in,jn,2)) THEN
k= k + 1
qtest = qvs( p(in,jn,k),pienv(in,jn,k),tinit(in,jn,2))
GO TO 10
END IF
klcl = k
IF(klcl == 1) klcl=2
!
!-----------------------------------------------------------------------
!
! Interpolate to determine the 'exact' pressure of the LCL
!
!-----------------------------------------------------------------------
!
plcl = p(in,jn,klcl)
plcl1 = p(in,jn,klcl-1)
pilcl = pienv(in,jn,klcl)
pilcl1 = pienv(in,jn,klcl-1)
qtest = 999.
iter = 0
20 CONTINUE
IF( ABS(qtest - qinit(in,jn,2)) > .00001) THEN
iter = iter + 1
pavg = (plcl + plcl1 )*.5
piavg = (pilcl + pilcl1)*.5
qtest = qvs(pavg,piavg,tinit(in,jn,2))
IF(qtest > qinit(in,jn,1)) THEN
plcl1 = pavg
pilcl1 = piavg
ELSE
plcl = pavg
pilcl = piavg
END IF
GO TO 20
END IF
plcl = pavg
pilcl= piavg
RETURN
END SUBROUTINE palcl
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION PIE ######
!###### ######
!###### Steven Lazarus ######
!###### University of Oklahoma ######
!###### School of Meteorology ######
!##################################################################
!##################################################################
FUNCTION pie(pdim,p0,rd,cp)
!
!-----------------------------------------------------------------------
!
! Purpose:
!
! Calculate the non-dimensional pressure given the dimensional
! pressure.
!
!-----------------------------------------------------------------------
!
! Author: Steve Lazarus
! 6/10/93
!
! Modification History:
!
!-----------------------------------------------------------------------
!
! Input :
!
! pdim Sounding pressure (Pa)
! p0 Surface pressure (Pa)
! Rd Gas constant for dry air (kg/(m s**2))
! Cp Specific heat for dry air (kg/(m s**2))
!
! Output:
!
! pie Non-dimensional pressure
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable declarations.
!
!-----------------------------------------------------------------------
!
REAL :: pdim,rd,cp,p0
!
!
pie = (pdim/p0)**(rd/cp)
!
! end of function pie
END FUNCTION pie
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION QVS ######
!###### ######
!###### Steven Lazarus ######
!###### University of Oklahoma ######
!###### School of Meteorology ######
!##################################################################
!##################################################################
FUNCTION qvs(pdim,pnon,t)
!
!-----------------------------------------------------------------------
!
! Purpose:
!
! Calculate the saturation vapor pressure
!
!-----------------------------------------------------------------------
!
! Author: Steve Lazarus
! 6/10/93
!
! Modification History:
!
!-----------------------------------------------------------------------
!
! Input :
!
! pdim Sounding pressure (Pa)
! pnon Sounding pressure (nondimensional)
! t Sounding Potential temperature at level k (K)
!
! Output:
!
! qvs Saturation vapor pressure
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable declarations.
!
!-----------------------------------------------------------------------
!
REAL :: pdim,pnon,t
!
a = 7.5*LOG(10.)
b = 3.8/pdim
!
es = 610.78*EXP(a*(pnon*t - 273.16)/(pnon*t - 35.86))
qvs = 0.622*es/(pdim-es)
!
! end of function qvs
END FUNCTION qvs