PROGRAM arpsensbc,44
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Generate one set of perturbation external BOUNDARY files from two
! sets of ARPS external boundary files and write it out for the use in
! ENSEMBLE forecast.
! The idea is based on the SLAF (Scaled Lagged Average Forecast).
! The procedure is as follows:
! 1, read file a;
! 2, read in file b;
! 3, find the difference bwtween a and b, store it in a. a=a-b
! 4, read in the base file c or use b as the control file (iread=0).
! 5, generate the perturbation c=c+a/n (b is used as c in code)
! 6, n is input as the variable iorder; it can be ... -2,-1,+1,+2 ...
!
! AUTHOR: Dingchen Hou
! History: Apr. 30, 1998: developed initially for SAMEX98.
! Sep. 02, 1999: modified to change the input file to namelist
! format.
! Feb-Apr, 2002: (F. KONG) Major modifications to:
! - add random perturbation option (with irandom = 1,
! and iseed provided in arpsensbc.input)
! - read in domain config directly from data
! This version allows to produce BC for both SLAF and
! BMG
!-----------------------------------------------------------------------
!
! DATA ARRAYS READ IN:
! (and another set with b repalcing a)
!
! au x component of velocity (m/s)
! av y component of velocity (m/s)
! aww vertical component of velocity in Cartesian
! coordinates (m/s).
!
! apt potential temperature (K)
! ap pressure (Pascal)
!
! aqv water vapor mixing ratio (kg/kg)
! aqc Cloud water mixing ratio (kg/kg)
! aqr Rainwater mixing ratio (kg/kg)
! aqi Cloud ice mixing ratio (kg/kg)
! aqs Snow mixing ratio (kg/kg)
! aqh Hail mixing ratio (kg/kg)
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'exbc.inc'
INCLUDE 'grid.inc'
!-----------------------------------------------------------------------
!
! Dimension nx, ny, nz
!
!-----------------------------------------------------------------------
INTEGER :: nch,idummy,hdmpfm
INTEGER :: nx, ny, nz
!
!-----------------------------------------------------------------------
REAL, ALLOCATABLE :: au (:,:,:) ! u-velocity (m/s)
REAL, ALLOCATABLE :: av (:,:,:) ! v-velocity (m/s)
REAL, ALLOCATABLE :: aww (:,:,:) ! w-velocity (m/s)
REAL, ALLOCATABLE :: apt (:,:,:) ! potential temperature (K)
REAL, ALLOCATABLE :: ap (:,:,:) ! pressure (Pascal)
REAL, ALLOCATABLE :: aqv (:,:,:) ! water vapor specific humidity
REAL, ALLOCATABLE :: aqc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: aqr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: aqi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: aqs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: aqh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: bu (:,:,:) ! u-velocity (m/s)
REAL, ALLOCATABLE :: bv (:,:,:) ! v-velocity (m/s)
REAL, ALLOCATABLE :: bww (:,:,:) ! w-velocity (m/s)
REAL, ALLOCATABLE :: bpt (:,:,:) ! potential temperature (K)
REAL, ALLOCATABLE :: bp (:,:,:) ! pressure (Pascal)
REAL, ALLOCATABLE :: bqv (:,:,:) ! water vapor specific humidity
REAL, ALLOCATABLE :: bqc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: bqr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: bqi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: bqs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: bqh (:,:,:) ! Hail mixing ratio (kg/kg)
!-----------------------------------------------------------------------
!
! Work Arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: filename1(10),filename2(10), &
filename3(10),filnmout(10)
INTEGER :: lenfil1(10),lenfil2(10),lenfil3(10),lenfilo(10)
INTEGER :: iorder,irandom,iseed,iread,ierr
INTEGER :: ireturn
INTEGER :: i,j,k
CHARACTER (LEN=15) :: ctime
INTEGER :: is,isets, istat
INTEGER :: sd_id,clipxy, clipz
!
!------------------------------------------------------------------------
!
! NAMELIST blocks for input parameters
!
!------------------------------------------------------------------------
! NAMELIST /grid_dims/ nx, ny, nz
! NAMELIST /ptbcpara/ iorder,iread,dx,dy,dz,ctrlat,ctrlon
NAMELIST /ptbcpara/ iorder,irandom,iseed,iread,ngbrz,exbcdmp,exbcfmt
NAMELIST /input_data/ isets,filename1,filename2,filename3
NAMELIST /outpt_data/ filnmout
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
READ (5,ptbcpara,END=999)
WRITE(6,'(/a)') 'Namelist block ptbcpara sucessfully read in.'
WRITE(6,ptbcpara)
READ (5,input_data,END=999)
WRITE(6,'(/a)') 'Namelist block input_data sucessfully read in.'
! WRITE(6,input_data)
READ (5,outpt_data,END=999)
WRITE(6,'(/a)') 'Namelist block output_data sucessfully read in.'
! WRITE(6,outpt_data)
DO is=1,isets
WRITE(6,'(/a,i5)')' The data set No.', is
lenfil1(is) = LEN_trim(filename1(is))
lenfil2(is) = LEN_trim(filename2(is))
lenfil3(is) = LEN_trim(filename3(is))
lenfilo(is) = LEN_trim(filnmout(is))
WRITE(6,'(/a,1x,a)') ' data file 1 ', filename1(is)(1:lenfil1(is))
WRITE(6,'(/a,1x,a)') ' data file 2 ', filename2(is)(1:lenfil2(is))
WRITE(6,'(/a,1x,a)') ' data file 3 ', filename3(is)(1:lenfil3(is))
WRITE(6,'(/a,1x,a)') ' output file ', filnmout(is)(1:lenfilo(is))
END DO
GO TO 1001
999 CONTINUE
PRINT *,'ERROR in reading input file, Program Terminated!'
STOP
1001 CONTINUE
!
!-----------------------------------------------------------------------
!
! Read in the dimensions etc.
!
!-----------------------------------------------------------------------
IF (exbcfmt == 0) THEN
! lenfil1(1)=len_trim(filename1(1))
CALL getunit
( nch )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(filename1(1)), '-F f77 -N ieee', ierr)
OPEN(nch,FILE=trim(filename1(1)),STATUS='old', &
FORM='unformatted',IOSTAT=istat)
IF ( istat /= 0 ) THEN
PRINT *, 'Problem to open BC file, STOP'
STOP
END IF
READ(nch,ERR=99) nx,ny,nz,dx,dy,dz,ctrlat,ctrlon
READ(nch,ERR=99) ctime
READ(nch,ERR=99) ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd, &
qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd, &
qhbcrd,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
GO TO 900
99 CONTINUE
WRITE (6,*) "READEXBC: ERROR reading data ", &
"from file ",trim(filename1(1)(1:lenfil1(1))),"STOP"
STOP
900 CONTINUE
CLOSE (nch)
CALL retunit( nch )
ELSE ! HDF4 format
CALL hdfopen(trim(filename1(1)), 1, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "READEXBC: ERROR opening ", &
trim(filename1(1))," for reading."
ierr = 1
! GO TO 900
STOP
END IF
CALL hdfrdc(sd_id,15,"ctime",ctime,istat)
CALL hdfrdi(sd_id,"nx",nx,istat)
CALL hdfrdi(sd_id,"ny",ny,istat)
CALL hdfrdi(sd_id,"nz",nz,istat)
CALL hdfrdr(sd_id,"dx",dx,istat)
CALL hdfrdr(sd_id,"dy",dy,istat)
CALL hdfrdr(sd_id,"dz",dz,istat)
CALL hdfrdr(sd_id,"dzmin",dzmin,istat)
CALL hdfrdi(sd_id,"strhopt",strhopt,istat)
CALL hdfrdr(sd_id,"zrefsfc",zrefsfc,istat)
CALL hdfrdr(sd_id,"dlayer1",dlayer1,istat)
CALL hdfrdr(sd_id,"dlayer2",dlayer2,istat)
CALL hdfrdr(sd_id,"zflat",zflat,istat)
CALL hdfrdr(sd_id,"strhtune",strhtune,istat)
CALL hdfrdi(sd_id,"mapproj",mapproj,istat)
CALL hdfrdr(sd_id,"trulat1",trulat1,istat)
CALL hdfrdr(sd_id,"trulat2",trulat2,istat)
CALL hdfrdr(sd_id,"trulon",trulon,istat)
CALL hdfrdr(sd_id,"sclfct",sclfct,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlat,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlon,istat)
CALL hdfrdi(sd_id,"ubcflg",ubcrd,istat)
CALL hdfrdi(sd_id,"vbcflg",vbcrd,istat)
CALL hdfrdi(sd_id,"wbcflg",wbcrd,istat)
CALL hdfrdi(sd_id,"ptbcflg",ptbcrd,istat)
CALL hdfrdi(sd_id,"prbcflg",prbcrd,istat)
CALL hdfrdi(sd_id,"qvbcflg",qvbcrd,istat)
CALL hdfrdi(sd_id,"qcbcflg",qcbcrd,istat)
CALL hdfrdi(sd_id,"qrbcflg",qrbcrd,istat)
CALL hdfrdi(sd_id,"qibcflg",qibcrd,istat)
CALL hdfrdi(sd_id,"qsbcflg",qsbcrd,istat)
CALL hdfrdi(sd_id,"qhbcflg",qhbcrd,istat)
CALL hdfrdi(sd_id, 'clipxy', clipxy,istat)
IF (istat == 0 .AND. clipxy < ngbrz) THEN
WRITE (6,*) "READEXBC: ERROR, clipxy (ngbrz) in exbc file too small"
ierr = 1
! GO TO 900
STOP
END IF
CALL hdfrdi(sd_id, 'clipz', clipz,istat)
IF (istat == 0 .AND. clipz > rayklow) THEN
WRITE (6,*) "READEXBC: ERROR, clipz (rayklow) in exbc file too large"
ierr = 1
! GO TO 900
STOP
END IF
!wdt end block
! alternate dump format ...
END IF
PRINT *, 'nx,ny,nz:',nx,ny,nz
PRINT *, 'dx,dy,dz:',dx,dy,dz
PRINT *, 'ctrlat,ctrlon:',ctrlat,ctrlon
PRINT *, 'ctime:',ctime
!-----------------------------------------------------------------------
!
! Allocate the variables and initialize the them to zero
!
!-----------------------------------------------------------------------
ALLOCATE(au(nx,ny,nz),STAT=istat)
au=0
ALLOCATE(av(nx,ny,nz),STAT=istat)
av=0
ALLOCATE(aww(nx,ny,nz),STAT=istat)
aww=0
ALLOCATE(apt(nx,ny,nz),STAT=istat)
apt=0
ALLOCATE(ap(nx,ny,nz),STAT=istat)
ap=0
ALLOCATE(aqv(nx,ny,nz),STAT=istat)
aqv=0
ALLOCATE(aqc(nx,ny,nz),STAT=istat)
aqc=0
ALLOCATE(aqr(nx,ny,nz),STAT=istat)
aqr=0
ALLOCATE(aqi(nx,ny,nz),STAT=istat)
aqi=0
ALLOCATE(aqs(nx,ny,nz),STAT=istat)
aqs=0
ALLOCATE(aqh(nx,ny,nz),STAT=istat)
aqh=0
ALLOCATE(bu(nx,ny,nz),STAT=istat)
bu=0
ALLOCATE(bv(nx,ny,nz),STAT=istat)
bv=0
ALLOCATE(bww(nx,ny,nz),STAT=istat)
bww=0
ALLOCATE(bpt(nx,ny,nz),STAT=istat)
bpt=0
ALLOCATE(bp(nx,ny,nz),STAT=istat)
bp=0
ALLOCATE(bqv(nx,ny,nz),STAT=istat)
bqv=0
ALLOCATE(bqc(nx,ny,nz),STAT=istat)
bqc=0
ALLOCATE(bqr(nx,ny,nz),STAT=istat)
bqr=0
ALLOCATE(bqi(nx,ny,nz),STAT=istat)
bqi=0
ALLOCATE(bqs(nx,ny,nz),STAT=istat)
bqs=0
ALLOCATE(bqh(nx,ny,nz),STAT=istat)
bqh=0
ALLOCATE(tem1(nx,ny,nz), STAT=istat)
tem1=0
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
is=0
101 CONTINUE
is=is+1
ierr=0
PRINT *,filename1(is)(1:lenfil1(is)),lenfil1(is)
! CALL readexbc(nx,ny,nz,dx,dy,dz,ctrlat,ctrlon, &
! filename1(is)(1:lenfil1(is)),lenfil1(is), ctime, &
! ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd, &
! qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd, &
! au,av,aww,apt,ap,aqv,aqc,aqr,aqi,aqs,aqh, ierr)
CALL readexbc(nx,ny,nz, &
filename1(is)(1:lenfil1(is)),lenfil1(is), ctime, &
au,av,aww,apt,ap,aqv,aqc,aqr,aqi,aqs,aqh, ierr)
IF (ierr /= 0) THEN
PRINT *,'ierr1=', ierr,' file',filename1(is)(1:lenfil1(is)), &
' not read, STOP!'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! Read in the verification data.
!
!-----------------------------------------------------------------------
!
ierr=0
! CALL readexbc(nx,ny,nz,dx,dy,dz,ctrlat,ctrlon, &
! filename2(is)(1:lenfil2(is)),lenfil2(is), ctime, &
! ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd, &
! qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd, &
! bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
CALL readexbc(nx,ny,nz, &
filename2(is)(1:lenfil2(is)),lenfil2(is), ctime, &
bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
IF (ierr /= 0) THEN
PRINT *,'ierr2=', ierr,' file',filename2(is)(1:lenfil2(is)), &
' not read, STOP!'
STOP
END IF
IF (irandom == 0) THEN
!
!-----------------------------------------------------------------------
!
! Find difference = first data - 2nd data set and store in 1st set
! a=a-b
! To reduce memory requirements, the difference fields are
! written to the same arrays as the first data set.
!
!-----------------------------------------------------------------------
!
CALL difebc(nx,ny,nz,0,tem1, &
au, av, aww, apt, ap, &
aqv, aqc, aqr, aqi, aqs, aqh, &
bu, bv, bww, bpt, bp, &
bqv, bqc, bqr, bqi, bqs, bqh, &
au, av, aww, apt, ap, &
aqv, aqc, aqr, aqi, aqs, aqh, &
ireturn )
ELSE IF (irandom == 1) THEN ! add random perturbation
do k=1,nz-1
do i=1,nx
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
au(i,j,k) = 2.*0.10*10.0*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny
iseed=mod(iseed*7141+54773,259200)
av(i,j,k) = 2.*0.10*10.0*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
apt(i,j,k) = 2.*0.10*11.0*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
ap(i,j,k) = 2.*0.10*1000.0*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
aqv(i,j,k) = 2.*0.10*2e-3*(iseed/259199.-.5)
enddo
enddo
enddo
aww = 0.0
aqc = 0.0
aqr = 0.0
aqi = 0.0
aqs = 0.0
aqh = 0.0
END IF
!
!-----------------------------------------------------------------------
!
! Read in the base data array if (iread=/=0)
!
!-----------------------------------------------------------------------
!
IF (iread /= 0) THEN
ierr=0
! CALL readexbc(nx,ny,nz,dx,dy,dz,ctrlat,ctrlon, &
! filename3(is)(1:lenfil3(is)),lenfil3(is), ctime, &
! ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd, &
! qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd, &
! bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
CALL readexbc(nx,ny,nz, &
filename3(is)(1:lenfil3(is)),lenfil3(is), ctime, &
bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
IF (ierr /= 0) THEN
PRINT *,'ierr3=', ierr,' file',filename3(is)(1:lenfil3(is)), &
' not read, STOP!'
STOP
END IF
END IF
!-----------------------------------------------------------------------
! Find the ptb field = 3rd file + difference/n
! b=b+a/n
! To reduce memory requirements, the perturbation fields are
! written to the same arrays as the base state fields.
!
!-----------------------------------------------------------------------
!
IF (iorder /= 0) THEN
CALL difebc(nx,ny,nz,iorder,tem1, &
bu, bv, bww, bpt, bp, &
bqv, bqc, bqr, bqi, bqs, bqh, &
au, av, aww, apt, ap, &
aqv, aqc, aqr, aqi, aqs, aqh, &
bu, bv, bww, bpt, bp, &
bqv, bqc, bqr, bqi, bqs, bqh, &
ireturn )
DO k=1,nz
DO j=1,ny
DO i=1,nx
IF (k /= 1) THEN
IF (bp(i,j,k) >= bp(i,j,k-1)) bp(i,j,k)=bp(i,j,k-1)-1.0
END IF
IF (bqv(i,j,k) <= 1.0E-8) THEN
bqv(i,j,k)=1.0E-8
END IF
END DO
END DO
END DO
END IF
!
print *, filnmout(is)(1:lenfilo(is))
print *, ctime
print *, 'ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd:'
print *, ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd
CALL writexbc(nx,ny,nz,filnmout(is)(1:lenfilo(is)),lenfilo(is), &
ctime, &
ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd, &
qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd, &
bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh)
IF (is < isets) GO TO 101
STOP
END PROGRAM arpsensbc
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DIFEBC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE difebc(nx,ny,nz,nscale,tem1, & 2,11
uprt, vprt, wprt, ptprt, pprt, &
qvprt, qc, qr, qi, qs, qh, &
auprt, avprt, awprt, aptprt, apprt, &
aqvprt, aqc, aqr, aqi, aqs, aqh, &
duprt, dvprt, dwprt, dptprt, dpprt, &
dqvprt, dqc, dqr, dqi, dqs, dqh, &
ireturn)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE: du=u+au/nscale (nscale=/=0) or du=u-au (nscale==0)
! (this is done by the subroutine SUBTRBC)
!
! Subtract the forecast fields from the interpolated verification
! fields (names beginning with "a") and output to the difference
! fields (names beginning with "d"). The input and difference
! fields may share the same storage location. For this subroutine
! it is assumed the forecast and corresponding verification
! data are at the same physical location, however, the physical
! location may differ between variables. That is uprt and auprt
! are at the same location, but that may differ from pprt and apprt.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster Ou School of Meteorology. April 1992
!
! MODIFICATION HISTORY:
! 14 May 1992 (KB) changed from arps2.5 to arps3.0
! 03 Aug 1992 (KB) updated to account for changes in arps3.0
!
! 09/07/1995 (KB)
! Added differencing of surface (soil) fields.
!
!-----------------------------------------------------------------------
!
! INPUT :
! nx,ny,nz Array dimensions for forecast field.
!
! FORECAST FIELDS:
!
! uprt perturbation x component of velocity (m/s)
! vprt perturbation y component of velocity (m/s)
! wprt perturbation vertical component of velocity in Cartesian
! coordinates (m/s).
!
! ptprt perturbation potential temperature (K)
! pprt perturbation pressure (Pascal)
!
! qvprt perturbation water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! INTERPOLATED VERIFICATION FIELDS:
!
! auprt perturbation x component of velocity (m/s)
! avprt perturbation y component of velocity (m/s)
! awprt perturbation vertical component of velocity in Cartesian
! coordinates (m/s).
!
! aptprt perturbation potential temperature (K)
! apprt perturbation pressure (Pascal)
!
! aqvprt perturbation water vapor mixing ratio (kg/kg)
! aqc Cloud water mixing ratio (kg/kg)
! aqr Rainwater mixing ratio (kg/kg)
! aqi Cloud ice mixing ratio (kg/kg)
! aqs Snow mixing ratio (kg/kg)
! aqh Hail mixing ratio (kg/kg)
!
! OUTPUT :
!
! DIFFERENCE FIELDS (may share storage with forecast fields
! or interpolated fields in calling program):
!
! duprt perturbation x component of velocity (m/s)
! dvprt perturbation y component of velocity (m/s)
! dwprt perturbation vertical component of velocity in Cartesian
! coordinates (m/s).
!
! dptprt perturbation potential temperature (K)
! dpprt perturbation pressure (Pascal)
!
! dqvprt perturbation water vapor mixing ratio (kg/kg)
! dqc Cloud water mixing ratio (kg/kg)
! dqr Rainwater mixing ratio (kg/kg)
! dqi Cloud ice mixing ratio (kg/kg)
! dqs Snow mixing ratio (kg/kg)
! dqh Hail mixing ratio (kg/kg)
!
!
! tem1 Work array
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! 3 dimensions of array
INTEGER :: nscale
!
!-----------------------------------------------------------------------
!
! Model Arrays
!
!-----------------------------------------------------------------------
!
REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s)
REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s)
REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor specific humidity
REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
!-----------------------------------------------------------------------
!
! Verification data interpolated to model grid
!
!-----------------------------------------------------------------------
!
REAL :: auprt (nx,ny,nz) ! Perturbation u-velocity (m/s)
REAL :: avprt (nx,ny,nz) ! Perturbation v-velocity (m/s)
REAL :: awprt (nx,ny,nz) ! Perturbation w-velocity (m/s)
REAL :: aptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: apprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: aqvprt (nx,ny,nz) ! Perturbation water vapor specific humidity
REAL :: aqc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: aqr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: aqi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: aqs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: aqh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
!-----------------------------------------------------------------------
!
! Difference arrays
!
!-----------------------------------------------------------------------
!
REAL :: duprt (nx,ny,nz) ! perturbation x component of velocity (m/s)
REAL :: dvprt (nx,ny,nz) ! perturbation y component of velocity (m/s)
REAL :: dwprt (nx,ny,nz) ! perturbation vertical component of
! velocity in Cartesian coordinates (m/s)
REAL :: dptprt (nx,ny,nz) ! perturbation potential temperature (K)
REAL :: dpprt (nx,ny,nz) ! perturbation pressure (Pascal)
REAL :: dqvprt (nx,ny,nz) ! perturbation water vapor mixing ratio (kg/kg)
REAL :: dqc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: dqr (nx,ny,nz) ! Rainwater mixing ratio (kg/kg)
REAL :: dqi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: dqs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: dqh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: tem1 (nx,ny,nz) ! A work array
INTEGER :: ireturn, i,j,k
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: is,js,ks,ie,je,ke
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
is=1
js=1
ks=1
ie=nx-1
je=ny-1
ke=nz-1
!-----------------------------------------------------------------------
!
! Scalars
!
!-----------------------------------------------------------------------
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k)=0.0
END DO
END DO
END DO
PRINT *, ' ptprt: '
CALL subtrbc(nx,ny,nz, ptprt,tem1,aptprt,tem1,dptprt,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' pprt: '
CALL subtrbc(nx,ny,nz, pprt, tem1, apprt, tem1, dpprt,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qvprt: '
CALL subtrbc(nx,ny,nz, qvprt,tem1,aqvprt,tem1,dqvprt,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qc: '
CALL subtrbc(nx,ny,nz, qc, tem1, aqc, tem1, dqc,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qr: '
CALL subtrbc(nx,ny,nz, qr, tem1, aqr, tem1, dqr,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qi: '
CALL subtrbc(nx,ny,nz, qi, tem1, aqi, tem1, dqi,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qs: '
CALL subtrbc(nx,ny,nz, qs, tem1, aqs, tem1, dqs,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qh: '
CALL subtrbc(nx,ny,nz, qh, tem1, aqh, tem1, dqh,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! u wind components
!
!-----------------------------------------------------------------------
ie=nx
PRINT *, ' uprt: '
CALL subtrbc(nx,ny,nz,uprt,tem1,auprt,tem1,duprt,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! v wind components
!
!-----------------------------------------------------------------------
ie=nx-1
je=ny
PRINT *, ' vprt: '
CALL subtrbc(nx,ny,nz,vprt,tem1,avprt,tem1,dvprt,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! w wind components
!
!-----------------------------------------------------------------------
je=ny-1
ke=nz
PRINT *, ' wprt: '
CALL subtrbc(nx,ny,nz,wprt,tem1,awprt,tem1,dwprt,nscale, &
is,js,ks,ie,je,ke)
!
RETURN
END SUBROUTINE difebc
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SUBTRBC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE subtrbc(nx,ny,nz, a,abar,b,bbar,c,nscale, & 11
istr,jstr,kstr,iend,jend,kend)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Subtracts 2 three-dimensional arrays, represented by
! means plus perturbations.
!
! AUTHOR: Keith Brewster OU School of Meteorology. Feb 1992
!
! MODIFICATION HISTORY:
! 11 Aug 1992 (KB) changed from arps2.5 to arps3.0
!
!
!-----------------------------------------------------------------------
!
! INPUT:
! a perturbation data array
! abar mean data array
! b perturbation data array to subtract from a
! bbar mean data array to subtract from a
!
! OUTPUT:
! c difference array a-b
! (may share storage in calling program with array a or b)
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! 3 dimensions of array
INTEGER :: nscale
REAL :: a (nx,ny,nz) ! data array
REAL :: abar(nx,ny,nz) ! base state of data array a
REAL :: b (nx,ny,nz) ! data array to subtract from a
REAL :: bbar(nx,ny,nz) ! base state of data arrya b
REAL :: c (nx,ny,nz) ! difference array a-b
INTEGER :: istr,jstr,kstr
INTEGER :: iend,jend,kend
INTEGER :: i,j,k,imid,jmid,kmid
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
imid=nint(0.5*FLOAT(istr+iend))
jmid=nint(0.5*FLOAT(jstr+jend))
kmid=nint(0.5*FLOAT(kstr+kend))
!
!-----------------------------------------------------------------------
!
! Tell us about a sample input point
!
!-----------------------------------------------------------------------
!
IF (nscale == 0) THEN
PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), &
' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid))
ELSE
PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), &
' b= ',b(imid,jmid,kmid)
END IF
!
!-----------------------------------------------------------------------
!
! Subtraction
!
!-----------------------------------------------------------------------
!
DO k=kstr,kend
DO j=jstr,jend
DO i=istr,iend
IF (nscale == 0) THEN
c(i,j,k)=a(i,j,k)+abar(i,j,k)-(b(i,j,k)+bbar(i,j,k))
ELSE
c(i,j,k)=a(i,j,k)+b(i,j,k)/FLOAT(nscale)
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Tell us about a sample output point
!
!-----------------------------------------------------------------------
!
IF (nscale == 0) THEN
PRINT *, ' c= ',c(imid,jmid,kmid)
ELSE
PRINT *, ' c= ',c(imid,jmid,kmid)+abar(imid,jmid,kmid)
END IF
RETURN
END SUBROUTINE subtrbc