! !################################################################## !################################################################## !###### ###### !###### 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) 137,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,2 !----------------------------------------------------------------------- ! ! 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, f_walltime !----------------------------------------------------------------------- ! ! Miscellaneous local variables: ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- INCLUDE 'globcst.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ use_acct = 1 wall_init = f_walltime() cpu_init = f_cputime() acct_cpu_time = 0.0 acct_wall_time = 0.0 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) 72,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, f_walltime !----------------------------------------------------------------------- ! ! Miscellaneous local variables: ! !----------------------------------------------------------------------- REAL :: cpu0, wall0 !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- INCLUDE 'globcst.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ IF (use_acct == 0) RETURN cpu0 = 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, f_walltime !----------------------------------------------------------------------- ! ! Miscellaneous local variables: ! !----------------------------------------------------------------------- REAL :: cpu0, wall0 !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- INCLUDE 'globcst.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ IF (use_acct == 0) RETURN cpu0 = 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, f_walltime !----------------------------------------------------------------------- ! ! Miscellaneous local variables: ! !----------------------------------------------------------------------- REAL :: cpu0, wall0 !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- INCLUDE 'globcst.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ IF (use_acct == 0) RETURN cpu0 = 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, f_walltime !----------------------------------------------------------------------- ! ! Miscellaneous local variables: ! !----------------------------------------------------------------------- REAL :: cpu0, wall0 REAL :: time INTEGER :: acct !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- INCLUDE 'globcst.inc' INCLUDE 'mp.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ IF (use_acct == 0) RETURN cpu0 = 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 CALL mptotal(total_wall) CALL mptotal(total_cpu) DO acct = 1,max_acct time = acct_cpu_time(acct) CALL mptotal(time) acct_cpu_time(acct) = time time = acct_wall_time(acct) CALL mptotal(time) acct_wall_time(acct) = 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,2 !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Output the timing statistics for the ARPS model. ! !----------------------------------------------------------------------- ! ! AUTHOR: Gene Bassett ! 2000/09/14 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! OUTPUT: ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! ! Variable Declarations: ! !----------------------------------------------------------------------- IMPLICIT NONE ! Force explicit declarations !----------------------------------------------------------------------- ! ! Miscellaneous local variables: ! !----------------------------------------------------------------------- REAL :: cpu0, wall0, tc, tw !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- INCLUDE 'globcst.inc' INCLUDE 'mp.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !----------------------------------------------------------------------- ! ! Print out CPU usage. ! !----------------------------------------------------------------------- ! ! 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 == 0.0) total_cpu =1.0 IF(total_wall == 0.0) total_wall=1.0 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, '%', & '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 mpbarrier RETURN END SUBROUTINE acct_report_arps REAL FUNCTION f_walltime() 13 ! !----------------------------------------------------------------------- ! ! Function to measure wall clock time. ! !----------------------------------------------------------------------- ! INTEGER :: wall_int INTEGER :: ticspersec,max_count DOUBLE PRECISION :: walltime INTEGER :: rollover,lastval DATA rollover /0/ DATA lastval /0/ SAVE rollover,lastval CALL system_clock(wall_int,ticspersec,max_count) IF (wall_int < lastval) THEN ! assume that the clock has rolled over once WRITE (6,*) "F_WALLTIME: wall clock assumed to have rolled over once." rollover = rollover + 1 ENDIF walltime = dble(wall_int)/dble(ticspersec) & + (dble(max_count)/dble(ticspersec))*dble(rollover) lastval = wall_int f_walltime = walltime RETURN END FUNCTION f_walltime !################################################################## !################################################################## !###### ###### !###### 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