subroutine rnevp(irun,tl,ql,rain,pl,plebot,pletop,clfrac,dt, 1 rcon) C*********************************************************************** C Subroutine rnevp C ---------------- C Purpose: Calculate the amount of re-evaporation of rain in a C single level as a function of the rain rate and C moisture deficit amount. [Sud and Molod, MWR 1988] C C Input: C irun The number of grid points to calculate C tl (irun) Temperature in deg K C ql (irun) Specific Humidity in kg/kg C rain (irun) Input rain rate in mm per time step C pl (irun) Mid-level pressure in mb C plebot (irun) Bottom edge-level pressure in mb C pletop (irun) Top edge-level pressure in mb C clfrac (irun) Fraction of layer into which re-evaporation occurs C dt Time step over which input rain occured in seconds C C Output: C rcon (irun) Resulting rain amount after re-evaporation in mm per C time step C*********************************************************************** implicit none integer irun real tl(irun),ql(irun),rain(irun),pl(irun),plebot(irun) real pletop(irun),clfrac(irun),rcon(irun) real dt real qeq(irun) real actevap,arearat,deltaq,mass,massinv,potevap real teq,qsteq,dqdt,exparg,expapp real dum,getcon,elocp,gamfac,gravcon,relax integer I,N real expcoef(4) data expcoef/1., -2.507213e-1, 2.92732e-2, -3.8278e-3/ CFPP$ EXPAND (QSAT) gravcon = 100./getcon('GRAVITY') elocp = getcon('LATENT HEAT COND')/getcon('CP') gamfac = getcon('LATENT HEAT COND') * getcon('EPS') * elocp . / getcon('RGAS') C Iteration loop to find equillibrium saturation specific humidity do I = 1,irun teq = tl(i) qeq(i) = ql(i) do N=1,2 if(N.eq.1)relax=0.5 if(N.gt.1)relax=1. call qsat ( teq,pl(i),qsteq,dqdt,.true. ) deltaq=(qsteq-qeq(i))/(1.+elocp*dqdt) qeq(i)=qeq(i)+deltaq*relax teq=teq-deltaq*elocp*relax enddo enddo do I=1,irun mass = (plebot(i)-pletop(i)) * gravcon massinv = 1./mass potevap=(qeq(i)-ql(i))*mass rcon(I)=rain(I) if ((rcon(I).GT.0.).AND.(potevap.GT.0.)) then CCC exparg= -1.E-4*dt*sqrt(rcon(i)*3600.*(sqrt(pl(I)*0.001))/dt) exparg= -1.04E-4*dt*((rcon(i)*3600.*(sqrt(pl(I)*0.001))/dt) 1 **0.578) CCC arearat = 1.-(exp(exparg)) expapp = ((expcoef(4) * exparg + expcoef(3)) * exparg 1 + expcoef(2)) * exparg + expcoef(1) arearat = 1.-(1./(expapp*expapp*expapp*expapp)) actevap = potevap*arearat*min(1.,clfrac(i)) IF(actevap.GE.rcon(i)) actevap = rcon(i) rcon(I)=rcon(I)-actevap ql(I) = ql(I)+actevap*massinv tl(I) = tl(I)-actevap*massinv*elocp endif enddo return end SUBROUTINE QSAT (TT,P,Q,DQDT,LDQDT) IMPLICIT NONE REAL TT, P, Q, DQDT LOGICAL LDQDT REAL AIRMW, H2OMW, LICE, RUNIV, RGAS PARAMETER ( AIRMW = 28.97 ) PARAMETER ( H2OMW = 18.01 ) PARAMETER ( LICE = 2.834e6 ) PARAMETER ( RUNIV = 8314.3 ) PARAMETER ( RGAS = RUNIV/AIRMW) REAL TMIN, TMAX, TICE, TUSEICE PARAMETER ( TMIN = 223.15 ) PARAMETER ( TMAX = 323.15 ) PARAMETER ( TICE = 273.16 ) PARAMETER ( TUSEICE = -40.00 ) REAL ESFAC, ERFAC PARAMETER ( ESFAC = H2OMW/AIRMW ) PARAMETER ( ERFAC = (1.0-ESFAC)/ESFAC ) REAL B0, B1, B2, B3, B4, B5, B6 PARAMETER (B6= 6.136820929E-11*ESFAC) PARAMETER (B5= 2.034080948E-8 *ESFAC) PARAMETER (B4= 3.031240396E-6 *ESFAC) PARAMETER (B3= 2.650648471E-4 *ESFAC) PARAMETER (B2= 1.428945805E-2 *ESFAC) PARAMETER (B1= 4.436518521E-1 *ESFAC) PARAMETER (B0= 6.107799961E+0 *ESFAC) REAL BI0, BI1, BI2, BI3, BI4, BI5, BI6 PARAMETER (BI6= 1.838826904E-10*ESFAC) PARAMETER (BI5= 4.838803174E-8 *ESFAC) PARAMETER (BI4= 5.824720280E-6 *ESFAC) PARAMETER (BI3= 4.176223716E-4 *ESFAC) PARAMETER (BI2= 1.886013408E-2 *ESFAC) PARAMETER (BI1= 5.034698970E-1 *ESFAC) PARAMETER (BI0= 6.109177956E+0 *ESFAC) REAL C1, C2, C3, C4, C5, C6 PARAMETER (C1= B1 ) PARAMETER (C2= B2*2.) PARAMETER (C3= B3*3.) PARAMETER (C4= B4*4.) PARAMETER (C5= B5*5.) PARAMETER (C6= B6*6.) REAL CI1, CI2, CI3, CI4, CI5, CI6 PARAMETER (CI1= BI1 ) PARAMETER (CI2= BI2*2.) PARAMETER (CI3= BI3*3.) PARAMETER (CI4= BI4*4.) PARAMETER (CI5= BI5*5.) PARAMETER (CI6= BI6*6.) REAL T, D, W, QX, DQX T = AMIN1(TT,TMAX) - TICE DQX = 0. QX = 0. IF(T.GT.0.) THEN QX = (T*(T*(T*(T*(T*(T*B6+B5)+B4)+B3)+B2)+B1)+B0) IF (LDQDT) THEN DQX = (T*(T*(T*(T*(T*C6+C5)+C4)+C3)+C2)+C1) ENDIF ELSEIF(T.LT.TUSEICE) THEN qx = esfac * 6.107*exp( esfac*lice/rgas*(1./tice - 1./tt) ) c QX = (T*(T*(T*(T*(T*(T*BI6+BI5)+BI4)+BI3)+BI2)+BI1)+BI0) IF (LDQDT) THEN dqx = qx * esfac*lice/(rgas*tt*tt) c DQX = (T*(T*(T*(T*(T*CI6+CI5)+CI4)+CI3)+CI2)+CI1) ENDIF ELSE W = (TUSEICE - T)/TUSEICE QX = W *(T*(T*(T*(T*(T*(T*B6+B5)+B4)+B3)+B2)+B1)+B0) & + (1.-W)*(T*(T*(T*(T*(T*(T*BI6+BI5)+BI4)+BI3)+BI2)+BI1)+BI0) IF (LDQDT) THEN DQX = W *(T*(T*(T*(T*(T*C6+C5)+C4)+C3)+C2)+C1) & + (1.-W)*(T*(T*(T*(T*(T*CI6+CI5)+CI4)+CI3)+CI2)+CI1) ENDIF ENDIF D = (P-ERFAC*QX) IF(D.LT.0.) THEN Q = 1.0 IF (LDQDT) DQDT = 0. ELSE D = 1.0 / D Q = AMIN1(QX * D,1.0) IF (LDQDT) DQDT = (1.0 + ERFAC*Q) * D * DQX ENDIF RETURN END FUNCTION GETCON (NAME) C*********************************************************************** C C GENERIC FUNCTION GETCON IS A REPOSITORY OF GLOBAL VARIABLES, I.E. C A MEMORY FOR SCALAR VALUES NEEDED THROUGHOUT A LARGE PROGRAM. C C THIS SPECIFIC FUNCTION, GETCON, REMEMBERS FLOATING POINT VALUES. C THE FUNCTION IS CALLED WITH A CHARACTER NAME TO INTERROGATE A VALUE. C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** C PARAMETER (MAXCON=36) CHARACTER*(*) NAME CHARACTER*16 ANAME(MAXCON) DIMENSION ACON (MAXCON) C COMPUTATIONAL CONSTANTS C ----------------------- PARAMETER ( VECMAX = 65535.5 ) PARAMETER ( CALTOJ = 4184. ) PARAMETER ( UNDEF = 1.E15 ) C ASTRONOMICAL CONSTANTS C ---------------------- PARAMETER ( OB = 23.45 ) PARAMETER ( PER = 102. ) PARAMETER ( ECC = .0167 ) PARAMETER ( AE = 6371E3 ) PARAMETER ( EQNX = 80.5 ) PARAMETER ( SOLS = 176.5 ) PARAMETER ( S0 = 1365.0 ) C TERRESTRIAL CONSTANTS C --------------------- PARAMETER ( GRAV = 9.81 ) PARAMETER ( SRFPRS = 984.7 ) PARAMETER ( PIMEAN = 984.7 ) PARAMETER ( PSTD = 1000.0 ) PARAMETER ( TSTD = 280.0 ) PARAMETER ( SDAY = 86400.0 ) PARAMETER ( SSALB = 0.99 ) PARAMETER ( CO2 = 330.0 ) C THERMODYNAMIC CONSTANTS C ----------------------- PARAMETER ( CPD = 1004.16 ) PARAMETER ( CPV = 1869.46 ) PARAMETER ( ALHL = 2.499E6 ) PARAMETER ( ALHS = 2.845E6 ) PARAMETER ( STFBOL = 5.67E-8 ) PARAMETER ( AIRMW = 28.97 ) PARAMETER ( H2OMW = 18.01 ) PARAMETER ( RUNIV = 8314.3 ) PARAMETER ( RGAS = RUNIV/AIRMW) PARAMETER ( RKAP = RGAS/CPD ) PARAMETER ( HEATW = 597.2 ) PARAMETER ( HEATI = 680.0 ) C TURBULENCE CONSTANTS C -------------------- PARAMETER ( VKRM = 0.41 ) C MOISTURE CONSTANTS C ------------------ PARAMETER ( TICE = 273.16 ) PARAMETER ( EPS = 0.622 ) PARAMETER ( VIRTCON= 0.609 ) PARAMETER ( EPSFAC = EPS*HEATW/RGAS*CALTOJ ) DATA ANAME(1 ),ACON(1 ) / 'CP ', CPD / DATA ANAME(2 ),ACON(2 ) / 'RGAS ', RGAS / DATA ANAME(3 ),ACON(3 ) / 'KAPPA ', RKAP / DATA ANAME(4 ),ACON(4 ) / 'LATENT HEAT COND', ALHL / DATA ANAME(5 ),ACON(5 ) / 'GRAVITY ', GRAV / DATA ANAME(6 ),ACON(6 ) / 'STEFAN-BOLTZMAN ', STFBOL / DATA ANAME(7 ),ACON(7 ) / 'VON KARMAN ', VKRM / DATA ANAME(8 ),ACON(8 ) / 'EARTH RADIUS ', AE / DATA ANAME(9 ),ACON(9 ) / 'OBLIQUITY ', OB / DATA ANAME(10),ACON(10) / 'ECCENTRICITY ', ECC / DATA ANAME(11),ACON(11) / 'PERIHELION ', PER / DATA ANAME(12),ACON(12) / 'VERNAL EQUINOX ', EQNX / DATA ANAME(13),ACON(13) / 'SUMMER SOLSTICE ', SOLS / DATA ANAME(14),ACON(14) / 'MAX VECT LENGTH ', VECMAX / DATA ANAME(15),ACON(15) / 'MOL WT H2O ', H2OMW / DATA ANAME(16),ACON(16) / 'MOL WT AIR ', AIRMW / DATA ANAME(17),ACON(17) / 'CPV ', CPV / DATA ANAME(18),ACON(18) / 'CPD ', CPD / DATA ANAME(19),ACON(19) / 'UNIV GAS CONST ', RUNIV / DATA ANAME(20),ACON(20) / 'LATENT HEAT SBLM', ALHS / DATA ANAME(21),ACON(21) / 'FREEZING-POINT ', TICE / DATA ANAME(23),ACON(23) / 'CALTOJ ', CALTOJ / DATA ANAME(24),ACON(24) / 'EPS ', EPS / DATA ANAME(25),ACON(25) / 'HEATW ', HEATW / DATA ANAME(26),ACON(26) / 'EPSFAC ', EPSFAC / DATA ANAME(27),ACON(27) / 'VIRTCON ', VIRTCON/ DATA ANAME(28),ACON(28) / 'PIMEAN ', PIMEAN / DATA ANAME(29),ACON(29) / 'SDAY ', SDAY / DATA ANAME(30),ACON(30) / 'HEATI ', HEATI / DATA ANAME(31),ACON(31) / 'S0 ', S0 / DATA ANAME(32),ACON(32) / 'PSTD ', PSTD / DATA ANAME(33),ACON(33) / 'TSTD ', TSTD / DATA ANAME(34),ACON(34) / 'SSALB ', SSALB / DATA ANAME(35),ACON(35) / 'UNDEF ', UNDEF / DATA ANAME(36),ACON(36) / 'CO2 ', CO2 / DO 10 I=1,MAXCON IF(NAME.EQ.ANAME(I)) THEN GETCON = ACON(I) RETURN ENDIF 10 CONTINUE 900 PRINT *,' CANNOT FIND FLOATING POINT CONSTANT - ',NAME PRINT *,' GETCON - CANNOT FIND CONSTANT REQUESTED' STOP END