!wdt update acct variables & f_walltime to double precision
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RANARY ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE ranary(nx,ny,nx1,nx2,ny1,ny2,iseed,amplit, rantp) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Generate a 2-D array of machine-independent random numbers
! between -amplit and +amplit with average value equal to zero.
! The input parameter, iseed, must be a negative integer number
! which needs to be defined only once in the calling subroutine
! program.
!
!-----------------------------------------------------------------------
!
! AUTHOR: V.Wong and X.Song
! 7/20/1992.
!
! Modification history:
! 7/25/1992. (MX)
! Modified to allow assignment on a portion of the array.
!
! 12/11/1992 (MX)
! Bug fix to the index bounds in loop 10.
!
! 9/1/94 (J. Levit & Y. Lu)
! Cleaned up documentation
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of random numbers in the x-direction
! ny Number of random numbers in the y-direction
! nx1,nx2 Begin and end 1st array indicies specifying a
! subdomain in which the average is set to zero
! ny1,ny2 Begin and end 2nd array indicies specifying a
! subdomain in which the average is set to zero
! iseed an arbitrary negative integer as a seed for a
! sequence of random numbers
! amplit The generated numbers stay within the range
! [-amplit,amplit]
!
! OUTPUT:
!
! rantp The array for storing random numbers
!
! WORK ARRAY:
!
! rantp Temporary working array
!
!
!-----------------------------------------------------------------------
!
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: i,j
INTEGER :: nx,ny ! Number of random numbers in each direction
INTEGER :: nx1,nx2 ! Begin and end 1st array indecies specifying
! a subdomain in which the average is set to
! zero
INTEGER :: ny1,ny2 ! Begin and end 2nd array indecies specifying
! a subdomain in which the average is set to
! zero
INTEGER :: iseed ! The seed for random number generation
REAL :: amplit ! The generated numbers stay within
! [-amplit, amplit]
REAL :: rantp (nx,ny) ! Temporary working array for storing the
! numbers
REAL :: ran3 ! The function to generate random number.
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
REAL :: wv,ranave,rantpmax
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Generate random numbers between 0 and 1, for a given iseed.
!
!-----------------------------------------------------------------------
!
ranave=0.0
DO j=ny1,ny2
DO i=nx1,nx2
rantp(i,j)=ran3(iseed)
ranave=ranave+rantp(i,j)
END DO
END DO
ranave=ranave/(FLOAT( (nx2-nx1+1)*(ny2-ny1+1) ))
!
!-----------------------------------------------------------------------
!
! Adjust the random numbers so that the new plane average value
! equals zero and the random numbers stay within the range
! [-amplit, amplit].
!
!-----------------------------------------------------------------------
!
rantpmax=0.0
DO j=ny1,ny2
DO i=nx1,nx2
rantpmax=AMAX1(rantpmax,ABS(ranave-rantp(i,j)))
END DO
END DO
wv=amplit/rantpmax
DO j=ny1,ny2
DO i=nx1,nx2
rantp(i,j)=(rantp(i,j)-ranave)*wv
END DO
END DO
RETURN
END SUBROUTINE ranary
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION RAN3 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION ran3(iseed)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Generates a random number between 0 and 1 by feeding
! a negative integer iseed.
!
! Reference: "Seminumerical Algorithms" by Donald Knuth
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! iseed an arbitrary negative integer as a seed for a
! sequence of random numbers
!
! OUTPUT:
!
! ran3 A random number between 0 and 1.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: iseed ! The seed for random number generation
REAL :: ran3 ! The function to generate random number.
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: mbig,mseed,mz,k,ii,inext,inextp,i,iff,mk,mj
REAL :: fac
PARAMETER (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
INTEGER :: ma(55)
SAVE iff, inext, inextp, ma
DATA iff /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!----------------------------------------------------------------------
!
! Initialize the sequence of random numbers between 0 and 1,
! using iseed.
!
!----------------------------------------------------------------------
!
IF (iseed < 0.OR.iff == 0) THEN
iff=1
mj=mseed-IABS(iseed)
mj=MOD(mj,mbig)
ma(55)=mj
mk=1
DO i=1,54
ii=MOD(21*i,55)
ma(ii)=mk
mk=mj-mk
IF (mk < mz) mk=mk+mbig
mj=ma(ii)
END DO
DO k=1,4
DO i=1,55
ma(i)=ma(i)-ma(1+MOD(i+30,55))
IF (ma(i) < mz) ma(i)=ma(i)+mbig
END DO
END DO
inext=0
inextp=31
iseed=1
END IF
!
!----------------------------------------------------------------------
!
! Start to generate a random number.
!
!----------------------------------------------------------------------
!
inext=inext+1
IF (inext == 56) inext=1
inextp=inextp+1
IF (inextp == 56) inextp=1
mj=ma(inext)-ma(inextp)
IF (mj < mz) mj=mj+mbig
ma(inext)=mj
ran3=mj*fac
RETURN
END FUNCTION ran3
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SWAP4BYTE ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE swap4byte(a,n) 2
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reverse order of bytes in integer*4 words
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: n
INTEGER*4 a(n)
INTEGER :: i
INTEGER*4 itema
CHARACTER (LEN=1) :: ctema(4)
CHARACTER (LEN=1) :: ctemb
EQUIVALENCE ( itema,ctema(1) )
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i = 1,n
itema = a(i)
ctemb = ctema(4)
ctema(4) = ctema(1)
ctema(1) = ctemb
ctemb = ctema(3)
ctema(3) = ctema(2)
ctema(2) = ctemb
a(i) = itema
END DO
RETURN
END SUBROUTINE swap4byte
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FINDVARTYPE ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE findvartype(iendn,itypec,lw) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Find the type of endian, type of character set and length of
! machine word. Renamed version of subroutine w3fi04 from nmcdecode.f90.
!
!-----------------------------------------------------------------------
!
!SUBROUTINE w3fi04(iendn,itypec,lw)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM: W3FI04 FIND WORD SIZE, ENDIAN, CHARACTER SET
! PRGMNR: JONES,R.E. ORG: W/NMC42 DATE: 94-10-07
!
! ABSTRACT: SUBROUTINE COMPUTES WORD SIZE, THE TYPE OF CHARACTER
! SET, ASCII OR EBCDIC, AND IF THE COMPUTER IS BIG-ENDIAN, OR
! LITTLE-ENDIAN.
!
! PROGRAM HISTORY LOG:
! 94-10-07 R.E.JONES
!
! USAGE: CALL W3FI04 (IENDN, ITYPEC, LW)
!
! OUTPUT ARGUMENT LIST:
! IENDN - INTEGER FOR BIG-ENDIAN OR LITTLE-ENDIAN
! = 0 BIG-ENDIAN
! = 1 LITTLE-ENDIAN
! = 2 CANNOT COMPUTE
! ITYPEC - INTEGER FOR TYPE OF CHARACTER SET
! = 0 ASCII CHARACTER SET
! = 1 EBCDIC CHARACTER SET
! = 2 NOT ASCII OR EBCDIC
! LW - INTEGER FOR WORDS SIZE OF COMPUTER IN BYTES
! = 4 FOR 32 BIT COMPUTERS
! = 8 FOR 64 BIT COMPUTERS
!
! ATTRIBUTES:
! LANGUAGE: SiliconGraphics 3.5 FORTRAN 77
! MACHINE: SiliconGraphics IRIS-4D/25, 35, INDIGO, INDY
!
!$$$
!
INTEGER :: itest1
INTEGER :: itest2
INTEGER :: itest3
INTEGER :: iendn
INTEGER :: itypec
INTEGER :: lw
!
CHARACTER (LEN=8) :: ctest1
CHARACTER (LEN=8) :: ctest2
CHARACTER (LEN=1) :: ctest3(8)
CHARACTER (LEN=1) :: BLANK
!
EQUIVALENCE (ctest1,itest1),(ctest2,itest2)
!
EQUIVALENCE (itest3,ctest3(1))
!
DATA ctest1/'12345678'/
DATA itest3/z'01020304'/
DATA BLANK /' '/
!
SAVE
!
! TEST FOR TYPE OF CHARACTER SET
! BLANK IS 32 (20 HEX) IN ASCII, 64 (40 HEX) IN EBCDEC
!
IF (ICHAR(BLANK) == 32) THEN
itypec = 0
ELSE IF (ICHAR(BLANK) == 64) THEN
!
! COMPUTER IS PROBABLY AN IBM360, 370, OR 390 WITH
! A 32 BIT WORD SIZE, AND BIG-ENDIAN.
!
itypec = 1
ELSE
itypec = 2
END IF
!
! TEST FOR WORD SIZE, SET LW TO 4 FOR 32 BIT COMPUTER,
! 8 FOR FOR 64 BIT COMPUTERS
!
itest2 = itest1
IF (ctest1 == ctest2) THEN
!
! COMPUTER MAY BE A CRAY, OR COULD BE DEC VAX ALPHA
! OR SGI WITH R4000, R4400, R8800 AFTER THEY CHANGE
! FORTRAN COMPILERS FOR 64 BIT INTEGER.
!
lw = 8
ELSE
lw = 4
END IF
!
! IF CHARACTER SET IS NOT ASCII OR EBCDIC SET IENDN = 2
! CAN NOT TEST FOR ENDIAN TYPE
!
IF (itypec == 2) THEN
iendn = 2
RETURN
END IF
!
! USING ITEST3 WITH Z'01020304' EQUIVALNCED TO CTEST3
! ON A 32 BIT BIG-ENDIAN COMPUTER 03 IS IN THE 3RD
! BYTE OF A 4 BYTE WORD. ON A 32 BIT LITTLE-ENDIAN
! COMPUTER IT IS IN 2ND BYTE.
! ON A 64 BIT COMPUTER Z'01020304' IS RIGHT ADJUSTED IN
! A 64 BIT WORD, 03 IS IN THE 7TH BYTE. ON A LITTLE-
! ENDIAN 64 BIT COMPUTER IT IS IN THE 2ND BYTE.
!
IF (lw == 4) THEN
IF (ICHAR(ctest3(3)) == 3) THEN
iendn = 0
ELSE IF (ICHAR(ctest3(3)) == 2) THEN
iendn = 1
ELSE
iendn = 2
END IF
ELSE IF (lw == 8) THEN
IF (ICHAR(ctest3(7)) == 3) THEN
iendn = 0
ELSE IF (ICHAR(ctest3(2)) == 3) THEN
iendn = 1
ELSE
iendn = 2
END IF
ELSE
iendn = 2
END IF
!
RETURN
END SUBROUTINE findvartype
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FIX_LAKE_ELIV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE fix_lake_eliv(h,lat,lon,nx,ny) 2
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Fill in Great Lakes. The original data set corresponds to the
! bottom of the lakes.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Richard Carpenter.
! 2000/02/03
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT/OUPUT :
!
! h Height of the grid point.
!
! INPUT :
!
! lat Latitude of the grid point.
! lon Longitude of the grid point.
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
!
!
! OUTPUT:
!
! WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx, ny ! x- y-Dimensions.
REAL :: h(nx,ny) ! Height
REAL :: lat(nx,ny) ! Latitude
REAL :: lon(nx,ny) ! Longitude
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: nlake, maxlake
PARAMETER (maxlake=12)
REAL :: lakelat1(maxlake), lakelat2(maxlake), &
lakelon1(maxlake), lakelon2(maxlake), lakeelev(maxlake)
REAL :: latpnt,lonpnt
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
PRINT *, 'ARPSTRN: Great Lakes terrain'
!
!-----------------------------------------------------------------------
!
! Define the lat/lon boxes which contain each lake
!
!-----------------------------------------------------------------------
!
k = 0
! Lake Superior
k = k + 1
lakelat1(k) = 46.4
lakelat2(k) = 49.0
lakelon1(k) = -92.5
lakelon2(k) = -84.1
lakeelev(k) = 182.9
! Lake Michigan
k = k + 1
lakelat1(k) = 41.0
lakelat2(k) = 46.3
lakelon1(k) = -88.1
lakelon2(k) = -84.7
lakeelev(k) = 176.5
! Lake Huron
k = k + 1
lakelat1(k) = 43.0
lakelat2(k) = 46.4
lakelon1(k) = -84.7
lakelon2(k) = -79.5
lakeelev(k) = 176.5
! Lake St. Clair
k = k + 1
lakelat1(k) = 42.2
lakelat2(k) = 42.8
lakelon1(k) = -83.0
lakelon2(k) = -82.2
lakeelev(k) = 174.6
! Lake Erie
k = k + 1
lakelat1(k) = 41.1
lakelat2(k) = 43.0
lakelon1(k) = -83.7
lakelon2(k) = -78.6
lakeelev(k) = 173.7
! Lake Ontario
k = k + 1
lakelat1(k) = 43.0
lakelat2(k) = 45.6
lakelon1(k) = -80.0
lakelon2(k) = -76.0
lakeelev(k) = 74.1
nlake = k
PRINT *, 'ARPSTRN: Number of Great Lakes: ', nlake
!
!-----------------------------------------------------------------------
!
! Do no allow the elevation in a box to go below the lake's elevation.
!
!-----------------------------------------------------------------------
!
DO j=1,ny
DO i=1,nx
latpnt = lat(i,j)
lonpnt = lon(i,j)
IF(lonpnt > 180) lonpnt = lonpnt - 360
DO k=1,nlake
IF (latpnt > lakelat1(k) .AND. latpnt < lakelat2(k) .AND. &
lonpnt > lakelon1(k) .AND. lonpnt < lakelon2(k) .AND. &
h(i,j) < lakeelev(k) ) THEN
! write (*,'(a,2i4,2f8.1,i4,2f7.1)')
! : "FIXELEVL",i,j,lakeelev(k),h(i,j),k,latpnt,lonpnt
h(i,j) = lakeelev(k)
END IF
END DO
! write (*,'(a,2i4,f8.1,12x,2f7.1)')
! : "FIXELEV ",i,j,h(i,j),latpnt,lonpnt
END DO
END DO
RETURN
END SUBROUTINE fix_lake_eliv
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTMNSND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE extmnsnd(nx,ny,nx_ext,ny_ext,nz_ext,lvlprof, & 2,5
iextmn,iextmx,jextmn,jextmx,kextmn,kextmx, &
xscl,yscl,x_ext,y_ext, &
hgt_ext,pln_ext,t_ext, &
rhs_ext,u_ext,v_ext, &
zsnd,plsnd,psnd,tsnd,ptsnd,rhssnd,qvsnd, &
rhosnd,usnd,vsnd)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Finds a mean sounding of all external data within the ARPS
! domain. This is used to define the model base state.
! The data are interpolated to constant height levels and
! averaged in the horizontal. The final sounding is hydrostatic.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! August, 1995 Replaces MNSOUND
!
! MODIFICATION HISTORY:
! 12/20/95 Keith Brewster Added code to insure at least one
! point goes into creating the mean sounding.
!
! 2000/04/05 (Gene Bassett)
! Added qvsnd array.
!
! 04/02/2001 (K. Brewster)
! Corrected calculations for qvsnd array that were also impacting
! rhssnd output array.
!
! Also, added parameter rhmin=0.05 (5 percent) as a minimum relative
! humidity in place of hardcoded check value of 1.0E-04.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Sizing variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny
INTEGER :: nx_ext,ny_ext,nz_ext
INTEGER :: lvlprof
!
!-----------------------------------------------------------------------
!
! Limits of calculations
!
!-----------------------------------------------------------------------
!
INTEGER :: iextmn,iextmx
INTEGER :: jextmn,jextmx
INTEGER :: kextmn,kextmx
!
!-----------------------------------------------------------------------
!
! ARPS grid variables
!
!-----------------------------------------------------------------------
!
REAL :: xscl(nx)
REAL :: yscl(ny)
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext,ny_ext)
REAL :: y_ext(nx_ext,ny_ext)
!
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height MSL
REAL :: t_ext(nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: pln_ext(nx_ext,ny_ext,nz_ext) ! ln of total pressure
REAL :: rhs_ext(nx_ext,ny_ext,nz_ext) ! RHstar=SQRT(1.-RH)
REAL :: u_ext(nx_ext,ny_ext,nz_ext) ! u wind component
REAL :: v_ext(nx_ext,ny_ext,nz_ext) ! v wind component
!
!-----------------------------------------------------------------------
!
! Input sounding variables
!
!-----------------------------------------------------------------------
!
REAL :: zsnd(lvlprof)
!
!-----------------------------------------------------------------------
!
! Output sounding variables
!
!-----------------------------------------------------------------------
!
REAL :: plsnd(lvlprof) ! ln of pressure
REAL :: psnd(lvlprof) ! Pressure (Pa)
REAL :: tsnd(lvlprof) ! Temperature (K)
REAL :: ptsnd(lvlprof) ! Potential temperature (K)
REAL :: rhssnd(lvlprof) ! RHstar
REAL :: qvsnd(lvlprof) ! Specific humidity (kg/kg)
REAL :: rhosnd(lvlprof) ! Density
REAL :: usnd(lvlprof) ! u wind component
REAL :: vsnd(lvlprof) ! v wind component
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Lapse rate and a constant for hydrostatic integration
!
!-----------------------------------------------------------------------
!
REAL :: gamma,pconst
PARAMETER (gamma = 0.0068, & ! 6.8 degrees per km
pconst = (g/rd))
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
REAL, PARAMETER :: rhmin=0.05
REAL, PARAMETER :: rhmax=1.0
INTEGER :: i,j,k,kext,imid,jmid
INTEGER :: knt,knttot,ktop,kbot
REAL :: wlow,whigh,tvbar,xmid,ymid,dist2,distmn,accept
REAL :: c1,c2,pres,qvsat,qvbot,qvtop,rh,epsilon
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
epsilon=0.01*(zsnd(2)-zsnd(1))
!
!-----------------------------------------------------------------------
!
! Zero-out all sounding arrays
! ptsnd array is used temporarily to store a count.
!
!-----------------------------------------------------------------------
!
DO k=1,lvlprof
ptsnd(k)=0.
plsnd(k)=0.
tsnd(k)=0.
rhssnd(k)=0.
qvsnd(k)=0.
usnd(k)=0.
vsnd(k)=0.
END DO
!
!-----------------------------------------------------------------------
!
! Find the external data grid location that is closest
! to the middle of the grid.
! For the case where no external grid point lies within the
! grid this grid column will be used as the mean sounding.
!
!-----------------------------------------------------------------------
!
xmid=0.5*(xscl(1)+xscl(nx-1))
ymid=0.5*(yscl(1)+yscl(ny-1))
distmn=(x_ext(1,1)-xmid)*(x_ext(1,1)-xmid) + &
(y_ext(1,1)-ymid)*(y_ext(1,1)-ymid)
imid=1
jmid=1
DO j=1,ny_ext
DO i=1,nx_ext
dist2=(x_ext(i,j)-xmid)*(x_ext(i,j)-xmid) + &
(y_ext(i,j)-ymid)*(y_ext(i,j)-ymid)
IF(dist2 < distmn) THEN
distmn=dist2
imid=i
jmid=j
END IF
END DO
END DO
iextmn=MIN(iextmn,imid)
iextmx=MAX(iextmx,imid)
jextmn=MIN(jextmn,jmid)
jextmx=MAX(jextmx,jmid)
!
!-----------------------------------------------------------------------
!
! Look at all external points possibly within domain
!
!-----------------------------------------------------------------------
!
knt=0
knttot=0
DO j=jextmn,jextmx
DO i=iextmn,iextmx
!
!-----------------------------------------------------------------------
!
! Is this point within the ARPS domain?
! Since the ARPS grid is Cartesian, need only compare
! the external grid coordinates to the x and y limits
! of the ARPS grid.
!
!-----------------------------------------------------------------------
!
knttot=knttot+1
IF((x_ext(i,j) >= xscl(1) .AND. x_ext(i,j) <= xscl(nx) .AND. &
y_ext(i,j) >= yscl(1) .AND. y_ext(i,j) <= yscl(ny)) .OR. &
(i == imid .AND. j == jmid)) THEN
knt=knt+1
!
!-----------------------------------------------------------------------
!
! Interpolate external data in vertical onto profile
! arrays.
!
!-----------------------------------------------------------------------
!
DO k=1,lvlprof
IF(zsnd(k) >= hgt_ext(i,j,kextmx)) EXIT
IF((hgt_ext(i,j,1)-zsnd(k)) < epsilon) THEN
DO kext=kextmn,kextmx-1
IF(hgt_ext(i,j,kext) >= zsnd(k)) EXIT
END DO
whigh=(zsnd(k)-hgt_ext(i,j,kext-1))/ &
(hgt_ext(i,j,kext)-hgt_ext(i,j,kext-1))
wlow=1.-whigh
ptsnd(k)=ptsnd(k)+1.
plsnd(k)=plsnd(k)+ &
whigh*pln_ext(i,j,kext)+wlow*pln_ext(i,j,kext-1)
tsnd(k)=tsnd(k)+ &
whigh*t_ext(i,j,kext) + wlow*t_ext(i,j,kext-1)
rhssnd(k)=rhssnd(k)+ &
whigh*rhs_ext(i,j,kext)+wlow*rhs_ext(i,j,kext-1)
usnd(k)=usnd(k)+ &
whigh*u_ext(i,j,kext) + wlow*u_ext(i,j,kext-1)
vsnd(k)=vsnd(k)+ &
whigh*v_ext(i,j,kext) + wlow*v_ext(i,j,kext-1)
END IF
END DO
END IF
END DO
END DO
WRITE(6,'(/a,i6,a,/a,i6,a/)') &
' extmnsnd found ',knt,' points within ARPS domain', &
' of ',knttot,' checked.'
!
!-----------------------------------------------------------------------
!
! Find lowest height with data
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a)') ' Finding range of mean sounding data ...'
accept=0.3*knt
DO k=1,lvlprof-1
WRITE(6,'(a,f10.2,a,f6.0,a,f10.0)') ' z = ',zsnd(k), &
' knt = ',ptsnd(k),' accept = ',accept
IF(ptsnd(k) > accept) EXIT
END DO
kbot=k
!
!-----------------------------------------------------------------------
!
! Find highest height with data
!
!-----------------------------------------------------------------------
!
DO k=lvlprof,2,-1
WRITE(6,'(a,f10.2,a,f6.0,a,f10.0)') ' z = ',zsnd(k), &
' knt = ',ptsnd(k),' accept = ',accept
IF(ptsnd(k) > accept) EXIT
END DO
ktop=k
!
WRITE(6,'(a,f10.2,a,f10.2,a)') &
' Height of external data for mean spans from ', &
zsnd(kbot),' to ',zsnd(ktop),' meters.'
!
!-----------------------------------------------------------------------
!
! Divide through to find average
!
!-----------------------------------------------------------------------
!
DO k=kbot,ktop
plsnd(k)=plsnd(k)/ptsnd(k)
tsnd(k)=tsnd(k)/ptsnd(k)
rhssnd(k)=rhssnd(k)/ptsnd(k)
usnd(k)=usnd(k)/ptsnd(k)
vsnd(k)=vsnd(k)/ptsnd(k)
END DO
!
!-----------------------------------------------------------------------
!
! Set variables "below-ground"
! Use a constant lapse rate, gamma.
! plsnd is a sort of first-guess log(pressure), needed for qv
! calculation from rhstar.
!
!-----------------------------------------------------------------------
!
pres = EXP(plsnd(kbot))
rh=MAX(rhmin,rhmax-(rhssnd(kbot)*rhssnd(kbot)))
qvsat = f_qvsat
( pres, tsnd(kbot) )
qvbot=rh*qvsat
c1=g/(rd*gamma)
DO k=kbot-1,1,-1
tsnd(k)=tsnd(kbot)-gamma*(zsnd(k)-zsnd(kbot))
plsnd(k)=plsnd(kbot)+ &
c1*ALOG((tsnd(kbot)-gamma*(zsnd(k)-zsnd(kbot)))/tsnd(kbot))
psnd(k) = EXP(plsnd(k))
qvsat = f_qvsat
( psnd(k), tsnd(k) )
rh=qvbot/qvsat
rhssnd(k)=SQRT(MAX(0.,(rhmax-rh)))
usnd(k)=usnd(kbot)
vsnd(k)=vsnd(kbot)
END DO
!
!-----------------------------------------------------------------------
!
! Set variables "above-top"
! Use a constant temperature.
! We're assuming stratosphere here.
!
!-----------------------------------------------------------------------
!
pres = EXP(plsnd(ktop))
rh=MAX(rhmin,rhmax-(rhssnd(ktop)*rhssnd(ktop)))
qvsat = f_qvsat
( pres, tsnd(ktop) )
qvtop=rh*qvsat
c2=g/rd
DO k=ktop+1,lvlprof
tsnd(k)=tsnd(ktop)
plsnd(k)=plsnd(ktop)-c2*(zsnd(k)-zsnd(ktop))/tsnd(ktop)
psnd(k) = EXP(plsnd(k))
qvsat = f_qvsat
( psnd(k), tsnd(k) )
rh=qvtop/qvsat
rhssnd(k)=SQRT(MAX(0.,(rhmax-rh)))
usnd(k)=usnd(ktop)
vsnd(k)=vsnd(ktop)
END DO
!
!-----------------------------------------------------------------------
!
! Calculate qv profile from RH-star.
! Temporarily use rhosnd to store virtual temperature
!
!-----------------------------------------------------------------------
!
DO k=1,lvlprof
pres = EXP(plsnd(k))
rh=MAX(rhmin,rhmax-(rhssnd(k)*rhssnd(k)))
qvsat = f_qvsat
( pres, tsnd(k) )
qvsnd(k)=rh*qvsat
rhosnd(k) = tsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
END DO
!
!-----------------------------------------------------------------------
!
! Make sure that the sounding is hydrostatic by integrating
! from kbot. rhosnd is really virtual temperature here.
!
!-----------------------------------------------------------------------
!
psnd(kbot)=EXP(plsnd(kbot))
DO k=kbot-1,1,-1
tvbar=0.5*(rhosnd(k+1)+rhosnd(k))
psnd(k)=psnd(k+1)*EXP(pconst*(zsnd(k+1)-zsnd(k))/tvbar)
END DO
!
DO k=kbot+1,lvlprof
tvbar=0.5*(rhosnd(k-1)+rhosnd(k))
psnd(k)=psnd(k-1)*EXP(pconst*(zsnd(k-1)-zsnd(k))/tvbar)
END DO
!
!-----------------------------------------------------------------------
!
! Derived variable calculations
! Compute density from virtual temperature.
! Compute potential temperature.
! Compute log of pressure.
!
!-----------------------------------------------------------------------
!
DO k=1,lvlprof
rhosnd(k)=psnd(k)/(rd*rhosnd(k))
ptsnd(k)=tsnd(k)*((p0/psnd(k))**rddcp)
! pres=exp(plsnd(k))
! print *,'psnd old, psnd new: ',pres,psnd(k),(psnd(k)-pres)
plsnd(k)=ALOG(psnd(k))
END DO
!
RETURN
END SUBROUTINE extmnsnd
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ARPSSTOP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE arpsstop(comment,outerr) 148,1
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Close down any necessary options (e.g. MPI) and stop execution of
! the program.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/04/21
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! comment Comment string.
! outerr Output error indicator (0-send to standard out,
! 1-send to standard error)
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
CHARACTER (LEN=*) :: comment ! comment
INTEGER :: outerr ! error code (0-standard out, 1-standard error)
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (outerr == 0) THEN
WRITE (6,*) trim(comment)
ELSE
WRITE (0,*) trim(comment)
END IF
IF (mp_opt > 0) THEN
CALL mpexit
(outerr)
END IF
STOP
END SUBROUTINE arpsstop
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRTCOMMENT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE wrtcomment(comment,outerr) 3
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write out a message to standard our or standard error.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/04/21
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! comment Comment string.
! outerr Output error indicator (0-send to standard out,
! 1-send to standard error)
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
CHARACTER (LEN=*) :: comment ! comment
INTEGER :: outerr ! error code (0-standard out, 1-standard error)
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (outerr == 0) THEN
IF (myproc == 0) WRITE (6,*) trim(comment)
ELSE
WRITE (0,*) trim(comment)
END IF
RETURN
END SUBROUTINE wrtcomment
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ACCT_INIT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE acct_init 1,3
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize cpu and wall clock timers.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/09/14
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
REAL :: f_cputime
DOUBLE PRECISION :: f_walltime
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
use_acct = 1
CALL init_walltime
wall_init = f_walltime
()
cpu_init = dble(f_cputime
())
acct_cpu_time = 0d0
acct_wall_time = 0d0
current_acct = 0
interrupt_acct = 0
RETURN
END SUBROUTINE acct_init
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SET_ACCT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE set_acct(account) 75,3
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Start cpu and wall timers for the specified account number after
! adding the values of the current account to its totals.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/09/14
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! account Integer pointer to the account to start timing.
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
INTEGER :: account
REAL :: f_cputime
DOUBLE PRECISION :: f_walltime
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
DOUBLE PRECISION :: cpu0, wall0
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (use_acct == 0) RETURN
cpu0 = dble(f_cputime
())
wall0 = f_walltime
()
IF (current_acct /= 0) THEN ! close off current account
acct_cpu_time(current_acct) = &
acct_cpu_time(current_acct) + cpu0 - cpu_acct
acct_wall_time(current_acct) = &
acct_wall_time(current_acct) + wall0 - wall_acct
IF (interrupt_acct /= 0) THEN ! reset the interrupt
CALL acct_interrupt
(interrupt_acct)
ENDIF
ENDIF
current_acct = account
cpu_acct = cpu0
wall_acct = wall0
RETURN
END SUBROUTINE set_acct
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ACCT_INTERRUPT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE acct_interrupt(account) 210,3
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interrupt accounting for the main timer to start tracting another
! account.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/09/14
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! account Integer pointer to the account to start timing. A value
! of zero will turn off all but to total cpu/wall timers.
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
INTEGER :: account
REAL :: f_cputime
DOUBLE PRECISION :: f_walltime
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
DOUBLE PRECISION :: cpu0, wall0
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (use_acct == 0) RETURN
cpu0 = dble(f_cputime
())
wall0 = f_walltime
()
IF (interrupt_acct /= 0) THEN ! stop an existing interrupt
CALL acct_stop_inter
ENDIF
interrupt_acct = account
cpu_inter = cpu0
wall_inter = wall0
RETURN
END SUBROUTINE acct_interrupt
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ACCT_STOP_INTER ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE acct_stop_inter 132,2
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Stop the interrupt timer and adjust the main timer accordingly.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/09/14
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
REAL :: f_cputime
DOUBLE PRECISION :: f_walltime
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
DOUBLE PRECISION :: cpu0, wall0
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (use_acct == 0) RETURN
cpu0 = dble(f_cputime
())
wall0 = f_walltime
()
IF (interrupt_acct /= 0) THEN
acct_cpu_time(current_acct) = &
acct_cpu_time(current_acct) - (cpu0 - cpu_inter)
acct_wall_time(current_acct) = &
acct_wall_time(current_acct) - (wall0 - wall_inter)
acct_cpu_time(interrupt_acct) = &
acct_cpu_time(interrupt_acct) + cpu0 - cpu_inter
acct_wall_time(interrupt_acct) = &
acct_wall_time(interrupt_acct) + wall0 - wall_inter
ENDIF
interrupt_acct = 0
RETURN
END SUBROUTINE acct_stop_inter
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ACCT_REPORT_ARPS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE acct_finish 1,8
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Stop cpu/wall timers and calculate totals.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/09/14
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
REAL :: f_cputime
DOUBLE PRECISION :: f_walltime
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
DOUBLE PRECISION :: cpu0, wall0
REAL :: time
INTEGER :: acct
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (use_acct == 0) RETURN
cpu0 = dble(f_cputime
())
wall0 = f_walltime
()
total_cpu = cpu0 - cpu_init
total_wall = wall0 - wall_init
IF (interrupt_acct /= 0) THEN ! stop an existing interrupt
CALL acct_stop_inter
END IF
IF (current_acct /= 0) THEN
CALL set_acct
(0)
END IF
!-----------------------------------------------------------------------
!
! Sum up the CPU times used by all processors if on a multi-processor
! system.
!
!-----------------------------------------------------------------------
IF (mp_opt > 0) THEN
time = real(total_wall)
CALL mptotal
(time)
total_wall = time
time = real(total_cpu)
CALL mptotal
(time)
total_cpu = time
DO acct = 1,max_acct
time = real(acct_cpu_time(acct))
CALL mptotal
(time)
acct_cpu_time(acct) = dble(time)
time = real(acct_wall_time(acct))
CALL mptotal
(time)
acct_wall_time(acct) = dble(time)
END DO
END IF
RETURN
END SUBROUTINE acct_finish
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ACCT_REPORT_ARPS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE acct_report_arps 1,4
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Output the timing statistics for the ARPS model.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/09/14
!
! MODIFICATION HISTORY:
!
! 13 March 2002 (Eric Kemp)
! Added time accounting for WRF BMJ scheme.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
REAL :: tc, tw
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Print out CPU usage.
!
!-----------------------------------------------------------------------
!
CALL flush
(6) !wdt
!wdt update
! A barrier before and after printing stats avoids spurious output
! in the body of the stats output.
CALL mpbarrier
1000 FORMAT &
(a/a,25(/1x,a,2x,1p,e13.6,a,2x,0p,f6.2,a,4x,1p,e13.6,a,2x,0p,f6.2,a))
IF (myproc == 0) THEN
IF(total_cpu == 0d0) total_cpu =1d0
IF(total_wall == 0d0) total_wall=1d0
tc = 100.0/total_cpu
tw = 100.0/total_wall
WRITE (6,'(/a/)') "ARPS CPU Summary:"
WRITE (6,1000) &
' Process CPU time WALL CLOCK time', &
' ----------------- ---------------------- ----------------------',&
'Initialization :',acct_cpu_time(init_acct), 's', &
acct_cpu_time(init_acct)*tc, '%', acct_wall_time(init_acct), 's', &
acct_wall_time(init_acct)*tw, '%', &
'Data output :',acct_cpu_time(output_acct), 's', &
acct_cpu_time(output_acct)*tc,'%', acct_wall_time(output_acct), 's', &
acct_wall_time(output_acct)*tw,'%', &
'Wind advection :',acct_cpu_time(advuvw_acct), 's', &
acct_cpu_time(advuvw_acct)*tc,'%', acct_wall_time(advuvw_acct), 's', &
acct_wall_time(advuvw_acct)*tw,'%', &
'Scalar advection:',acct_cpu_time(advs_acct), 's', &
acct_cpu_time(advs_acct)*tc, '%', acct_wall_time(advs_acct), 's', &
acct_wall_time(advs_acct)*tw, '%', &
'Coriolis force :',acct_cpu_time(coriol_acct), 's', &
acct_cpu_time(coriol_acct)*tc,'%', acct_wall_time(coriol_acct), 's', &
acct_wall_time(coriol_acct)*tw,'%', &
'Buoyancy term :',acct_cpu_time(buoy_acct), 's', &
acct_cpu_time(buoy_acct)*tc, '%', acct_wall_time(buoy_acct), 's', &
acct_wall_time(buoy_acct)*tw,'%', &
'Misc Large tstep:',acct_cpu_time(tinteg_acct), 's', &
acct_cpu_time(tinteg_acct)*tc,'%', acct_wall_time(tinteg_acct), 's', &
acct_wall_time(tinteg_acct)*tw,'%', &
'Small time steps:',acct_cpu_time(smlstp_acct), 's', &
acct_cpu_time(smlstp_acct)*tc,'%', acct_wall_time(smlstp_acct), 's', &
acct_wall_time(smlstp_acct)*tw,'%', &
'Radiation :',acct_cpu_time(rad_acct), 's', &
acct_cpu_time(rad_acct)*tc, '%', acct_wall_time(rad_acct), 's', &
acct_wall_time(rad_acct)*tw, '%', &
'Soil model :',acct_cpu_time(soil_acct), 's', &
acct_cpu_time(soil_acct)*tc, '%', acct_wall_time(soil_acct), 's', &
acct_wall_time(soil_acct)*tw,'%', &
'Surface physics :',acct_cpu_time(sfcphy_acct), 's', &
acct_cpu_time(sfcphy_acct)*tc,'%', acct_wall_time(sfcphy_acct), 's', &
acct_wall_time(sfcphy_acct)*tw,'%', &
'Turbulence :',acct_cpu_time(tmix_acct), 's', &
acct_cpu_time(tmix_acct)*tc, '%', acct_wall_time(tmix_acct), 's', &
acct_wall_time(tmix_acct)*tw,'%', &
'Comput. mixing :',acct_cpu_time(cmix_acct), 's', &
acct_cpu_time(cmix_acct)*tc, '%', acct_wall_time(cmix_acct), 's', &
acct_wall_time(cmix_acct)*tw,'%', &
'Rayleigh damping:',acct_cpu_time(raydmp_acct), 's', &
acct_cpu_time(raydmp_acct)*tc,'%', acct_wall_time(raydmp_acct), 's', &
acct_wall_time(raydmp_acct)*tw,'%', &
'TKE src terms :',acct_cpu_time(tkesrc_acct), 's', &
acct_cpu_time(tkesrc_acct)*tc,'%', acct_wall_time(tkesrc_acct), 's', &
acct_wall_time(tkesrc_acct)*tw,'%', &
'Gridscale precp.:',acct_cpu_time(qpfgrd_acct), 's', &
acct_cpu_time(qpfgrd_acct)*tc,'%', acct_wall_time(qpfgrd_acct), 's', &
acct_wall_time(qpfgrd_acct)*tw,'%', &
'Kuo cumulus :',acct_cpu_time(kuocum_acct), 's', &
acct_cpu_time(kuocum_acct)*tc,'%', acct_wall_time(kuocum_acct), 's', &
acct_wall_time(kuocum_acct)*tw,'%', &
'Kain-Fritsch :',acct_cpu_time(kfcum_acct), 's', &
acct_cpu_time(kfcum_acct)*tc, '%', acct_wall_time(kfcum_acct), 's', &
acct_wall_time(kfcum_acct)*tw, '%', &
'BMJ cumulus :',acct_cpu_time(bmjcum_acct),'s', &
acct_cpu_time(bmjcum_acct)*tc, '%', acct_wall_time(bmjcum_acct), 's', &
acct_wall_time(bmjcum_acct)*tw, '%', &
'Warmrain microph:',acct_cpu_time(warmrain_acct), 's', &
acct_cpu_time(warmrain_acct)*tc,'%',acct_wall_time(warmrain_acct),'s',&
acct_wall_time(warmrain_acct)*tw,'%', &
'Lin ice microph :',acct_cpu_time(linice_acct), 's', &
acct_cpu_time(linice_acct)*tc,'%', acct_wall_time(linice_acct), 's', &
acct_wall_time(linice_acct)*tw,'%', &
'NEM ice microph :',acct_cpu_time(nemice_acct), 's', &
acct_cpu_time(nemice_acct)*tc,'%', acct_wall_time(nemice_acct), 's', &
acct_wall_time(nemice_acct)*tw,'%', &
'Hydrometero fall:',acct_cpu_time(qfall_acct), 's', &
acct_cpu_time(qfall_acct)*tc, '%', acct_wall_time(qfall_acct), 's', &
acct_wall_time(qfall_acct)*tw, '%', &
'Bound.conditions:',acct_cpu_time(bc_acct), 's', &
acct_cpu_time(bc_acct)*tc, '%', acct_wall_time(bc_acct), 's', &
acct_wall_time(bc_acct)*tw, '%', &
'Message passing :',acct_cpu_time(mp_acct), 's', &
acct_cpu_time(mp_acct)*tc, '%', acct_wall_time(mp_acct), 's', &
acct_wall_time(mp_acct)*tw, '%', &
'Miscellaneous :',acct_cpu_time(misc_acct), 's', &
acct_cpu_time(misc_acct)*tc, '%', acct_wall_time(misc_acct), 's', &
acct_wall_time(misc_acct)*tw,'%'
WRITE(6,'(/1x,a,1p,e13.6,a,2x,0p,7x,4x,1p,e13.6,a,2x,0p,6x)') &
'Entire model :',total_cpu,'s',total_wall,'s'
WRITE(6,'(/1x,a,1p,e13.6,a,2x,0p,7x,4x,1p,e13.6,a,2x,0p,6x)') &
'Without Init/IO :', &
total_cpu-acct_cpu_time(init_acct)-acct_cpu_time(output_acct),'s', &
total_wall-acct_wall_time(init_acct)-acct_wall_time(output_acct),'s'
END IF
CALL flush
(6) !wdt
!wdt update
CALL mpbarrier
RETURN
END SUBROUTINE acct_report_arps
!wdt update
SUBROUTINE init_walltime 1
INTEGER :: lastval, rollover
COMMON /walltime/ lastval, rollover
lastval = 0
rollover = 0
RETURN
END SUBROUTINE init_walltime
DOUBLE PRECISION FUNCTION f_walltime() 5
!
!-----------------------------------------------------------------------
!
! Function to measure wall clock time.
!
!-----------------------------------------------------------------------
!
INTEGER :: lastval, rollover
INTEGER :: wall_int
INTEGER :: ticspersec,max_count
COMMON /walltime/ lastval, rollover
CALL system_clock(wall_int,ticspersec,max_count)
IF (wall_int < lastval .and. lastval > max_count*0.01) THEN
! Assume that the clock has rolled over once. The test for max_count*0.1
! is for Absoft which sometimes seems to have wall_int go backwards a
! little at times.
WRITE (6,*) "F_WALLTIME: wall clock assumed to have rolled over once."
rollover = rollover + 1
ENDIF
lastval = wall_int
f_walltime = dble(wall_int)/dble(ticspersec) &
+ (dble(max_count)/dble(ticspersec))*dble(rollover)
RETURN
END FUNCTION f_walltime
!wdt update - subroutine moved from another file to here:
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FIX_SOIL_NSTYP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE fix_soil_nstyp(nx,ny,nstypin,nstyp, & 7
tsfc,tsoil,wetsfc,wetdp,wetcanp)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Adjust soil variables if nstyp differs between what is defined in
! ARPS and what is in the data file.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/08/11
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx
! ny
! nstypin Number of soil types in the data to be imported
! nstyp Number of soil types to use in the arps
!
! INPUT/OUTPUT:
!
! tsfc Temperature at surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx
INTEGER :: ny
INTEGER :: nstypin,nstyp
REAL :: tsfc (nx,ny,0:nstyp) ! Temperature at surface (K)
REAL :: tsoil (nx,ny,0:nstyp) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyp) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyp) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyp) ! Canopy water amount
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
INTEGER :: is
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Only one soil type, use average values.
!
!-----------------------------------------------------------------------
IF ( nstypin == 1 ) THEN
tsfc (1:nx-1,1:ny-1,1) = tsfc (1:nx-1,1:ny-1,0)
tsoil (1:nx-1,1:ny-1,1) = tsoil (1:nx-1,1:ny-1,0)
wetsfc (1:nx-1,1:ny-1,1) = wetsfc (1:nx-1,1:ny-1,0)
wetdp (1:nx-1,1:ny-1,1) = wetdp (1:nx-1,1:ny-1,0)
wetcanp(1:nx-1,1:ny-1,1) = wetcanp(1:nx-1,1:ny-1,0)
END IF
!-----------------------------------------------------------------------
!
! Use average values for soil types beyond that provided in the data.
!
!-----------------------------------------------------------------------
IF ( nstypin < nstyp ) THEN
DO is=nstypin+1,nstyp
tsfc (1:nx-1,1:ny-1,is) = tsfc (1:nx-1,1:ny-1,0)
tsoil (1:nx-1,1:ny-1,is) = tsoil (1:nx-1,1:ny-1,0)
wetsfc (1:nx-1,1:ny-1,is) = wetsfc (1:nx-1,1:ny-1,0)
wetdp (1:nx-1,1:ny-1,is) = wetdp (1:nx-1,1:ny-1,0)
wetcanp(1:nx-1,1:ny-1,is) = wetcanp(1:nx-1,1:ny-1,0)
ENDDO
END IF
RETURN
END SUBROUTINE fix_soil_nstyp
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FIX_STYPFRCT_NSTYP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE fix_stypfrct_nstyp(nx,ny,nstypin,nstyp,stypfrct) 6
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Adjust stypfrct if the number of soil types read in differs from
! the number specified for arps to use, nstyp.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/08/11
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx
! ny
! nstypin Number of soil types in the data file
! nstyp Number of soil types to use in the arps
!
! INPUT/OUTPUT:
!
! stypfrct Soil type fraction
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx
INTEGER :: ny
INTEGER :: nstypin,nstyp
REAL :: stypfrct(nx,ny,nstyp)
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
INTEGER :: i,j,is
REAL :: frctot
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Readjust stypfrct if only part of the soil type were read in.
!
!-----------------------------------------------------------------------
IF (nstypin > nstyp) THEN
DO j=1,ny
DO i=1,nx
frctot = 0.0
DO is=1,nstyp
frctot = frctot + stypfrct(i,j,is)
END DO
IF (frctot == 0) THEN
frctot = 1.
stypfrct(i,j,1) = 1.
ENDIF
DO is=1,nstyp
stypfrct(i,j,is) = stypfrct(i,j,is) / frctot
END DO
END DO
END DO
ENDIF
!-----------------------------------------------------------------------
!
! Fill in any unset values.
!
!-----------------------------------------------------------------------
IF ( nstypin < nstyp ) THEN
stypfrct(1:nx-1,1:ny-1,nstypin+1:nstyp) = 0.
END IF
!-----------------------------------------------------------------------
!
! Make sure the the fraction totals up to 1 (or 0)
!
!-----------------------------------------------------------------------
DO j=1,ny
DO i=1,nx
frctot = 0.
DO is = 1,nstyp
frctot = frctot + stypfrct(i,j,is)
END DO
IF (frctot /= 0) THEN
DO is = 1,nstyp
stypfrct(i,j,is) = stypfrct(i,j,is) / frctot
END DO
END IF
END DO
END DO
RETURN
END SUBROUTINE fix_stypfrct_nstyp
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE REMAP_SOIL_VARS ######
!###### ######
!###### Developed by ######
!###### Weather Decision Technologies ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE remap_soil_vars(nx,ny,nstypin,nstyp, & 1
tsfc_in,tsoil_in,wetsfc_in,wetdp_in,wetcanp_in,soiltyp_in, &
tsfcin,tsoilin,wsfcin,wdpin,wcanpin, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,soiltyp)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Rearrange soil variables to be consistent with any difference between
! the soil types to be used and the ones in the data set.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/10/27
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx
! ny
! nstypin Number of soil types in the data file
! nstyp Number of soil types to use in the arps
! tsfcin
! tsoilin
! wsfcin
! wdpin
! wcanpin
! tsfc_in Temperature at surface in data set (K)
! tsoil_in Deep soil temperature in data set (K)
! wetsfc_in Surface soil moisture in data set
! wetdp_in Deep soil moisture in data set
! wetcanp_in Canopy water amount in data set
! soiltyp_in Soil type in data set
!
! OUTPUT:
!
! tsfc Temperature at surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
! soiltyp Soil type
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx
INTEGER :: ny
INTEGER :: nstypin,nstyp
INTEGER :: tsfcin,tsoilin,wsfcin,wdpin,wcanpin
REAL :: tsfc_in (nx,ny,0:nstypin) ! Temperature at surface (K)
REAL :: tsoil_in (nx,ny,0:nstypin) ! Deep soil temperature (K)
REAL :: wetsfc_in (nx,ny,0:nstypin) ! Surface soil moisture
REAL :: wetdp_in (nx,ny,0:nstypin) ! Deep soil moisture
REAL :: wetcanp_in(nx,ny,0:nstypin) ! Canopy water amount
INTEGER :: soiltyp_in(nx,ny,nstypin) ! Soil type in model domain
REAL :: tsfc (nx,ny,0:nstyp) ! Temperature at surface (K)
REAL :: tsoil (nx,ny,0:nstyp) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyp) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyp) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyp) ! Canopy water amount
INTEGER :: soiltyp(nx,ny,nstyp) ! Soil type in model domain
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
INTEGER :: i,j,k,is,isin
INTEGER :: ifound_styp
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Match the soil type in the old data set. If not found then set to
! average soil properties (except water).
!
!-----------------------------------------------------------------------
DO j=1,ny
DO i=1,nx
DO is=1,nstyp
ifound_styp = 0
DO isin=1,nstypin
IF (soiltyp(i,j,is) == soiltyp_in(i,j,isin)) THEN
IF (tsfcin == 1) tsfc(i,j,is) = tsfc_in(i,j,isin)
IF (tsoilin == 1) tsoil(i,j,is) = tsoil_in(i,j,isin)
IF (wsfcin == 1) wetsfc(i,j,is) = wetsfc_in(i,j,isin)
IF (wdpin == 1) wetdp(i,j,is) = wetdp_in(i,j,isin)
IF (wcanpin == 1) wetcanp(i,j,is) = wetcanp_in(i,j,isin)
ifound_styp = 1
EXIT
END IF
END DO
IF (ifound_styp == 0) THEN
IF (tsfcin == 1) tsfc(i,j,is) = tsfc_in(i,j,0)
IF (tsoilin == 1) tsoil(i,j,is) = tsoil_in(i,j,0)
IF (wsfcin == 1) wetsfc(i,j,is) = wetsfc_in(i,j,0)
IF (wdpin == 1) wetdp(i,j,is) = wetdp_in(i,j,0)
IF (wcanpin == 1) wetcanp(i,j,is) = wetcanp_in(i,j,0)
IF (tsfcin == 1 .and. soiltyp(i,j,is) == 13) THEN
tsfc(i,j,is) = tsoil_in(i,j,0)
END IF
END IF
END DO
END DO
END DO
END SUBROUTINE remap_soil_vars
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INTERP_SOIL ######
!###### ######
!###### Developed by ######
!###### Weather Decision Technologies ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE intrp_soil(nx,ny,nx1,ny1,nstyp,nstyp1,wx,wy,ix,jy, & 2
tsfc,tsoil,wetsfc,wetdp,wetcanp, &
soiltyp,stypfrct,vegtyp, &
tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1, &
soiltyp1,stypfrct1,vegtyp1)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interpolate soil properties onto another grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/11/01
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx
! ny
! nx1
! ny1
! xw Weight factor in x-direction
! yw Weight factor in y-direction
! ix Old grid index (lower left) for new grid point
! jy Old grid index (lower left) for new grid point
! nstyp Number of soil types on original grid
! nstyp1 Number of soil types to use on new grid
! tsfc Temperature at surface in data set (K)
! tsoil Deep soil temperature in data set (K)
! wetsfc Surface soil moisture in data set
! wetdp Deep soil moisture in data set
! wetcanp Canopy water amount in data set
! soiltyp Soil type in data set
! stypfrct Soil type fraction
! vegtyp Vegetation type in data set
!
! OUTPUT:
!
! tsfc1 Temperature at surface in data set (K)
! tsoil1 Deep soil temperature in data set (K)
! wetsfc1 Surface soil moisture in data set
! wetdp1 Deep soil moisture in data set
! wetcanp1 Canopy water amount in data set
! soiltyp1 Soil type in data set
! stypfrct1 Soil type fraction
! vegtyp1 Vegetation type in data set
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nx1,ny1
INTEGER :: nstyp,nstyp1
REAL :: wx(nx1,ny1) ! Weight factor in x-direction
REAL :: wy(nx1,ny1) ! Weight factor in y-direction
INTEGER :: ix(nx1,ny1) ! Old grid index (lower left) for new grid point
INTEGER :: jy(nx1,ny1) ! Old grid index (lower left) for new grid point
REAL :: tsfc (nx,ny,0:nstyp) ! Temperature at surface (K)
REAL :: tsoil (nx,ny,0:nstyp) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyp) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyp) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyp) ! Canopy water amount
INTEGER :: soiltyp(nx,ny,nstyp) ! Soil type in model domain
REAL :: stypfrct(nx,ny,nstyp)
INTEGER :: vegtyp(nx,ny)
REAL :: tsfc1 (nx1,ny1,0:nstyp1) ! Temperature at surface (K)
REAL :: tsoil1 (nx1,ny1,0:nstyp1) ! Deep soil temperature (K)
REAL :: wetsfc1 (nx1,ny1,0:nstyp1) ! Surface soil moisture
REAL :: wetdp1 (nx1,ny1,0:nstyp1) ! Deep soil moisture
REAL :: wetcanp1(nx1,ny1,0:nstyp1) ! Canopy water amount
INTEGER :: soiltyp1(nx1,ny1,nstyp1) ! Soil type in model domain
REAL :: stypfrct1(nx1,ny1,nstyp1)
INTEGER :: vegtyp1(nx1,ny1)
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
INCLUDE 'arpssfc.inc'
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
! Arrays start at zero in case no soil type defined (i.e. soiltyp=0)
REAL :: tsfc1sum (0:nsoiltyp) ! Ground sfc. temperature (K)
REAL :: tsoil1sum (0:nsoiltyp) ! Deep soil temperature (K)
REAL :: wetsfc1sum (0:nsoiltyp) ! Surface soil moisture
REAL :: wetdp1sum (0:nsoiltyp) ! Deep soil moisture
REAL :: wetcanp1sum (0:nsoiltyp) ! Canopy water amount
REAL :: stypfrct1sum(0:nsoiltyp) ! Frction of soil type
INTEGER :: i,j,i1,j1,is,ii
REAL :: weight,maxweight,frctot
REAL :: totweight(0:nsoiltyp)
INTEGER :: maxtype
INTEGER :: soiltype
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
tsfc1 = 0
tsoil1 = 0
wetsfc1 = 0
wetdp1 = 0
wetcanp1 = 0
soiltyp1 = 0
stypfrct1 = 0
vegtyp1 = 0
IF (nstyp == 1) THEN
stypfrct = 1 ! make sure stypfrct is set for nstyp=1
! copy level 0 to level 1 if level 1 is undefined:
IF (tsfc(1,1,1) <= 0) tsfc(:,:,1) = tsfc(:,:,0)
IF (tsoil(1,1,1) <= 0) tsoil(:,:,1) = tsoil(:,:,0)
IF (wetsfc(1,1,1) < 0) wetsfc(:,:,1) = wetsfc(:,:,0)
IF (wetdp(1,1,1) < 0) wetdp(:,:,1) = wetdp(:,:,0)
IF (wetcanp(1,1,1) < 0) wetcanp(:,:,1) = wetcanp(:,:,0)
ENDIF
DO j1=1,ny1-1
DO i1=1,nx1-1
tsfc1sum = 0.
tsoil1sum = 0.
wetsfc1sum = 0.
wetdp1sum = 0.
wetcanp1sum = 0.
stypfrct1sum = 0.
!vegtyp1sum = 0.
maxweight = 0.
totweight = 0.
vegtyp1(i1,j1) = 0
DO j=jy(i1,j1),jy(i1,j1)+1
DO i=ix(i1,j1),ix(i1,j1)+1
weight = wx(i1,j1)*(1.+ix(i1,j1)-i)*wy(i1,j1)*(1+jy(i1,j1)-j) &
+ (1.-wx(i1,j1))*(i-ix(i1,j1))*wy(i1,j1)*(1+jy(i1,j1)-j) &
+ wx(i1,j1)*(1.+ix(i1,j1)-i)*(1.-wy(i1,j1))*(j-jy(i1,j1)) &
+ (1.-wx(i1,j1))*(i-ix(i1,j1))*(1.-wy(i1,j1))*(j-jy(i1,j1))
DO is=1,nstyp
IF (stypfrct(i,j,is) > 0) THEN
soiltype = soiltyp(i,j,is)
IF (soiltype < 1 .or. soiltype > nsoiltyp) soiltype = 0
tsfc1sum(soiltyp(i,j,is)) = tsfc1sum(soiltyp(i,j,is)) &
+ weight*stypfrct(i,j,is)*tsfc(i,j,is)
tsoil1sum(soiltyp(i,j,is)) = tsoil1sum(soiltyp(i,j,is)) &
+ weight*stypfrct(i,j,is)*tsoil(i,j,is)
wetsfc1sum(soiltyp(i,j,is)) = wetsfc1sum(soiltyp(i,j,is)) &
+ weight*stypfrct(i,j,is)*wetsfc(i,j,is)
wetdp1sum(soiltyp(i,j,is)) = wetdp1sum(soiltyp(i,j,is)) &
+ weight*stypfrct(i,j,is)*wetdp(i,j,is)
wetcanp1sum(soiltyp(i,j,is)) = wetcanp1sum(soiltyp(i,j,is)) &
+ weight*stypfrct(i,j,is)*wetcanp(i,j,is)
stypfrct1sum(soiltyp(i,j,is)) = stypfrct1sum(soiltyp(i,j,is)) &
+ weight*stypfrct(i,j,is)
totweight(soiltyp(i,j,is)) = totweight(soiltyp(i,j,is)) + weight
END IF
END DO
! use the vegtyp of the closest point
IF (weight > maxweight) THEN
vegtyp1(i1,j1) = vegtyp(i,j)
maxweight = weight
END IF
!!use most popular vegtyp weighted by distance
!IF (vegtyp(i,j) > 0 .and. vegtyp(i,j) <= nvegtyp) THEN
! vegtyp1sum(vegtyp(i,j)) = vegtyp1sum(vegtyp(i,j)) + weight
!END IF
END DO
END DO
!!use most popular vegtyp weighted by distance
!maxweight = 0
!vegtyp1(i1,j1) = 0
!DO ii = 1,nvegtyp
! IF (vegtyp1sum(ii) > maxweight) THEN
! vegtyp1(i1,j1) = ii
! maxweight = vegtyp1sum(ii)
! END IF
!END DO
! get the nstyp1 largest weighted soil types
DO is = 1,nstyp1
maxtype = -1
maxweight = 0.
DO ii = 0,nsoiltyp
IF (stypfrct1sum(ii) > maxweight) THEN
maxweight = stypfrct1sum(ii)
maxtype = ii
END IF
END DO
IF (maxtype /= -1) THEN
soiltyp1(i1,j1,is) = maxtype
tsfc1(i1,j1,is) = tsfc1sum(maxtype)/stypfrct1sum(maxtype)
tsoil1(i1,j1,is) = tsoil1sum(maxtype)/stypfrct1sum(maxtype)
wetsfc1(i1,j1,is) = wetsfc1sum(maxtype)/stypfrct1sum(maxtype)
wetdp1(i1,j1,is) = wetdp1sum(maxtype)/stypfrct1sum(maxtype)
wetcanp1(i1,j1,is) = wetcanp1sum(maxtype)/stypfrct1sum(maxtype)
stypfrct1(i1,j1,is) = stypfrct1sum(maxtype)/totweight(maxtype)
stypfrct1sum(maxtype) = 0.
ELSE
IF (is == 1) THEN
WRITE (6,*) &
"INTRP_SOIL: WARNING, no soil type found, variables not assigned!"
ELSE
stypfrct1(i1,j1,is) = 0.
ENDIF
ENDIF
END DO
! renormalize
frctot = 0.
DO is = 1,nstyp1
frctot = frctot + stypfrct1(i1,j1,is)
END DO
IF (frctot /= 0) THEN
DO is = 1,nstyp1
stypfrct1(i1,j1,is) = stypfrct1(i1,j1,is) / frctot
END DO
ELSE
stypfrct1(i1,j1,1) = 1.
IF (nstyp1 .gt. 1) stypfrct1(i1,j1,2:nstyp1) = 0.
END IF
tsfc1(i1,j1,0) = 0.
tsoil1(i1,j1,0) = 0.
wetsfc1(i1,j1,0) = 0.
wetdp1(i1,j1,0) = 0.
wetcanp1(i1,j1,0) = 0.
DO is = 1,nstyp1
tsfc1(i1,j1,0) = tsfc1(i1,j1,0) &
+ stypfrct1(i1,j1,is)*tsfc1(i1,j1,is)
tsoil1(i1,j1,0) = tsoil1(i1,j1,0) &
+ stypfrct1(i1,j1,is)*tsoil1(i1,j1,is)
wetsfc1(i1,j1,0) = wetsfc1(i1,j1,0) &
+ stypfrct1(i1,j1,is)*wetsfc1(i1,j1,is)
wetdp1(i1,j1,0) = wetdp1(i1,j1,0) &
+ stypfrct1(i1,j1,is)*wetdp1(i1,j1,is)
wetcanp1(i1,j1,0) = wetcanp1(i1,j1,0) &
+ stypfrct1(i1,j1,is)*wetcanp1(i1,j1,is)
END DO
!FIXME: remove? GMB
! tsfc1(i1,j1,0) = tsfc1(i1,j1,0)
! tsoil1(i1,j1,0) = tsoil1(i1,j1,0)
! wetsfc1(i1,j1,0) = wetsfc1(i1,j1,0)
! wetdp1(i1,j1,0) = wetdp1(i1,j1,0)
! wetcanp1(i1,j1,0) = wetcanp1(i1,j1,0)
END DO
END DO
END SUBROUTINE intrp_soil
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE LD2DGRID ######
!###### ######
!###### Developed by ######
!###### Weather Decision Technologies ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE ld2dgrid(dxin,dyin,ctrlatin,ctrlonin, &
mprojin,trlat1in,trlat2in,trlonin,sclfctin)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Load 2-D grid parameters into the ARPS common blocks (useful for
! C programs using ARPS subroutines).
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2001/08/03
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
REAL :: dxin,dyin,ctrlatin,ctrlonin,trlat1in,trlat2in,trlonin,sclfctin
INTEGER :: mprojin
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
dx = dxin
dy = dyin
ctrlat = ctrlatin
ctrlon = ctrlonin
mapproj = mprojin
trulat1 = trlat1in
!wdt update
trulat2 = trlat2in
trulon = trlonin
sclfct = sclfctin
RETURN
END SUBROUTINE ld2dgrid
!wdt testing:
SUBROUTINE test_dump(data,filename,nx,ny,nz,stagdim,incfake),1
! write out a unformatted dump to debug MP version
implicit none
include 'mp.inc'
integer nx,ny,nz
integer stagdim,incfake
real data(nx,ny,nz)
character*(*) filename
character*80 tempname
integer lenfl
integer count
common /save_count/ count
!return
count = count + 1
write(tempname, '(a,a,i4.4)') trim(filename),'_c',count
if (mp_opt > 0) then
write(tempname, '(a,a,2i2.2)') &
trim(adjustl(tempname)),'_',loc_x,loc_y
endif
open (11,file=tempname,form='unformatted')
write (11) nx,ny,nz,stagdim,incfake
write (11) data
close (11)
CALL test_check_in
(tempname)
RETURN
END
SUBROUTINE test_dump2(data,filename,nx,ny,nz,ib,ie,jb,je,kb,ke)
! write out a unformatted dump to debug MP version
implicit none
include 'mp.inc'
integer nx,ny,nz
integer ib,ie,jb,je,kb,ke ! compare for i=ib,ie ...
real data(nx,ny,nz)
character*(*) filename
character*80 tempname
integer lenfl
integer count
common /save_count/ count
!return
count = count + 1
write(tempname, '(a,a,i4.4)') trim(filename),'_c',count
if (mp_opt > 0) then
write(tempname, '(a,a,2i2.2)') &
trim(adjustl(tempname)),'_',loc_x,loc_y
endif
open (11,file=tempname,form='unformatted')
write (11) nx,ny,nz,-1,-1
write (11) ib,ie,jb,je,kb,ke
write (11) data
close (11)
RETURN
END
subroutine test_check_in(comment) 1,4
implicit none
include 'mp.inc'
character*(*) comment
!return
CALL flush
(6)
CALL flush
(0)
write (*,*) "XXX checking in:",myproc," ",trim(comment)
CALL flush
(6)
CALL mpbarrier
return
end