!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