SUBROUTINE SFLX(ICE,SATFLG,DT,Z,NSOIL, NROOT, SLDPTH, i F,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2,TBOT,CH,CHKFF, i VEGTYP,SOITYP,SHDFAC, o ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX, 2 T1,CMC,SMC,STC,SNODEP,SOILW,SOILM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE AND ACCUMULATE EVAPOTRANSPIRATION VALUES CC ======= BASED UPON A SOIL HYDROLOGY MODEL FROM OREGON STATE CC UNIVERSITY. THE ROUTINE ALSO CALCULATES SOIL TEMPS CC AND MOISTURE VALUES AT TWO SOIL DEPTHS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Argument list in the CALL SFLX C Prepared by Fei Chen, Ken Mitchell in 07/97 C 1. Calling Statement C CALL SFLX(ICE,SATFLG,DT,Z,NSOIL,NROOT, SLDPTH, C & F,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2,TBOT,CH,CHKFF, C & VEGTYP,SOITYP,SHDFRC, C & ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX, C & T1,CMC,SMC,STC,SNODEP,SOILW,SOILM) C 2. Input C *** General Parameters *** C ICE: SEA-ICE FLAG (=1: sea-ice, =0: land) C SATFLG: ATMOSPHERIC SATURATION FLAG (=0.0: Saturation (ETP=0.0 in SFLX) C =1.0: otherwise) C DT: TIMESTEP (SEC) C Z: HT OF FIRST MODEL LEVEL ABOVE GROUND (M) C NSOIL: NUMBER OF SOIL LAYERS C NROOT: NUMBER OF ROOT-ZONE LAYERS C SLDPTH: THE THICKNESS OF EACH SOIL LAYER (M) C *** Atmospheric variables *** C F: NET DOWNWARD RADIATION (W M-2; positive, if downward) C F=SOLDN-SOLUP+LONGDN=(1-ALBEDO)*SOLDN+LONGDN C SOLDN: SOLAR DOWNWARD RADIATION (W M-2; positive, if downward) C SFCPRS: PRESSURE AT 1ST MDL LVL ABV SKIN (PASCALS) C PRCP: PRECIP RATE (KG M-2 S-1) C SFCTMP: AIR TEMPERATURE AT 1ST MDL LVL ABV SKIN (K) C TH2: AIR POTENTIAL TEMPERATURE EXTRAPOLATED TO THE SURFACE (K) C Q2: MIXING RATIO AT 1ST MDL LVL ABV SKIN (KG KG-1) C Q2SAT: SAT MIXING RATIO AT 1ST MDL LVL ABV SKIN (KG KG-1) C DQSDT2: SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=T2 (KG KG-1 K-1) C TBOT: BOTTOM SOIL TEMPERATURE (local yearly-mean air temperature) C CH: DRAG COEF FOR HEAT AND MOISTURE (M S-1) C CHKFF: DRAG COEF FOR HEAT AND MOISTURE MULTIPLIED BY RHO*Cp C *** Canopy/Soil Characteristics *** C VEGTYP: Vegetation type C SOITYP: Soil type C SHDFRC: PLANT SHADING FACTOR (or GREENESS FRACTION) (UNITLESS) C *** State Variables *** C T1: SKIN (GRND SFC) TEMPERATURE (K) C CMC: CANOPY MOISTURE CONTENT (M) C SMC(NSOIL): SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION) C STC(NSOIL): SOIL TEMP (K) C SNODEP: WATER EQUIVALENT SNOW DEPTH (M) C *** Note: The above state variables are also output *** C 3. OUTPUT C ETP: POTENTIAL EVAPORATION (W M-2) C ETA: ACTUAL HEAT FLUX (W M-2: negative, if upward from surface) C H: SENSIBLE HEAT FLUX (W M-2: negative, if upward from surface) C S: SOIL HEAT FLUX (W M-2: negative, if downward from surface) C RUNOFF1: GRND SFC RUNOFF (M S-1) C RUNOFF2: UNDERGROUND RUNOFF (M S-1) C Q1: EFFECTIVE MIXING RATIO AT GRND SFC ( KG KG-1) C SNMAX: SNOW MELT (M) (WATER EQUIVALENT) C SOILW: AVAILABLE WTAER CONTENT (M) C SOILM: MAX WATER CONTENT (M) C #include "sp.h" PARAMETER (NSOLD=4) INTEGER ICE INTEGER NSOIL,VEGTYP,SOITYP,NROOT LOGICAL SNOWNG LOGICAL FRZGRA REAL B REAL CFACTR REAL CH REAL CM REAL CMC REAL CMCMAX REAL DF1 REAL DKSAT REAL DT REAL DWSAT REAL EPSCA REAL ESD REAL ETA REAL ETP REAL F REAL F1 REAL PC REAL PRCP REAL Q1 REAL Q2 REAL Q2SAT REAL RCH REAL RR REAL S REAL SFCPRS REAL SFCTMP REAL SHDFAC REAL SMC(NSOLD),STC(NSOLD),SLDPTH(NSOLD) REAL SMCDRY REAL SMCMAX REAL SMCREF REAL SMCWLT REAL SNODEP REAL T1 REAL T1MIN REAL T1V REAL T24 REAL T2V REAL TBOT REAL TH2 REAL TH2V REAL Z REAL Z0 REAL ZSOIL ( NSOLD ) C Some universal veg/soil parameters PARAMETER(CFACTR=0.5,CMCMAX=0.5E-3,TOPT=298.0,RSMAX=5000.0) !!$omp threadlocal /rite/ COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOF, & DEW,RIB C C Initialization RUNOFF1 = 0.0 RUNOFF2 = 0.0 RUNOFF3 = 0.0 SNMAX = 0.0 C C Specify the VEG/SOIL parameters C CALL REDPRM(VEGTYP,SOITYP, o CMCMAX,TOPT,RSMAX,ALBEDO,Z0,SHDFAC,NROOT,RCMIN,RGL,HS, o B,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF,SMCDRY,F1) IF(ICE.EQ.1) THEN DO KZ = 1, NSOIL ZSOIL(KZ)=-3.*FLOAT(KZ)/FLOAT(NSOIL) ENDDO ELSE ZSOIL(1)=-SLDPTH(1) DO KZ = 2, NSOIL ZSOIL(KZ)=-SLDPTH(KZ)+ZSOIL(KZ-1) ENDDO ENDIF C C print*,' 1st DT=', DT C print*,' *********** Input Parameters *************' C print*,' ' C print*,' SOILTP=', SOILTP,' VEGTYP=', VEGTYP, ' NSOIL=', NSOIL, C & ' LAND=', LAND,' Z0=', Z0,' B=', B,' F1=', F1,' CFACTR=', CFACTR C print*,' PC=', PC,' SMCMAX=', SMCMAX,' SMCREF=', SMCREF, C & ' SMCWLT=', SMCWLT, ' SMCDRY=', SMCDRY,'DKSAT=',DKSAT, C & ' DWSAT=', DWSAT,'DT=',DT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C EXECUTABLE CODE BEGINS HERE...INITIALIZE MISC VARIABLES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC S = 0.0 SNOWNG = .FALSE. FRZGRA = .FALSE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT THE INPUT SNOW DEPTH (METERS) FROM A SNOWBOARD TYPE C MEASUREMENT TO A LIQUID EQUIVALENT SNOW DEPTH (ESD, IN METERS) C USING A 10 TO 1 RATIO. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ESD = SNODEP * 0.1 ESD = SNODEP IF(ICE.EQ.1) THEN ESD = 0.01 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC VIRTUAL TEMPS AND POTENTIAL TEMPS AT GRND (SUB 1) AND AT C THE 1ST MDL LVL ABV THE GRND (SUB 2). EXPON IS CP DIVD BY R. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC T1V = T1 * (1.0 + 0.61 * Q2 ) T2V = SFCTMP * ( 1.0 + 0.61 * Q2 ) C TH2 = SFCTMP + ( 0.0097545 * Z ) TH2V = TH2 * (1.0 + 0.61 * Q2 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DETERMINE IF IT'S PRECIPING AND WHAT KIND OF PRECIP IT IS. C IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING! C IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND C TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( PRCP .GT. 0.0 ) THEN IF ( SFCTMP .LE. 273.16 ) THEN SNOWNG = .TRUE. ELSE IF ( SFCTMP .LE. 273.16 ) FRZGRA = .TRUE. ENDIF ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF EITHER PRCP FLAG IS SET, CONVERT THE PRCP RATE FROM C KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH (METERS) AND ACCUM IT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( ( SNOWNG ) .OR. ( FRZGRA ) ) THEN ESD =ESD + PRCP * DT * 0.001 PRCP1 = 0.0 ELSE PRCP1 = PRCP ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Check: if ESD < 1.E-6 m. The snow cover will be ignored. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(ESD.LT.1.E-6) THEN SNMAX=ESD ESD = 0.0 SNOWNG = .FALSE. ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC GRND HEAT FLUX FOR SNOW PACK CASE (WILL EFFECTIVELY C CALC THE HEAT FLUX FROM THE SNOW PACK). REMEMBER, IF C T1 > 273.16 K, THERE WILL BE SNOW MELT AND THE SKIN TEMP WILL C BE 273.16 K. CONST (0.35) IS THERMAL DIFFUSIVITY FOR SNOW. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( ESD .GT. 0.0 ) THEN T1MIN = AMIN1 ( T1, 273.16 ) S = 0.35 * ( T1MIN - STC(1) ) / ( 10.0 * ESD ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SET THE MINIMUM FLUX WITHIN SNOW TO -200 W M-2 C FOR NUMERICAL REASONS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC S = MIN ( MAX ( S,-100.0 ), 100.0 ) ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE GRND HEAT FLUX FOR CASE WITH NO SNOW PACK CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL TDFCND ( DF1, SMC(1), B, F1 ) S = DF1 * ( STC(1) - T1 ) / ( 0.5 * ZSOIL(1) ) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION C AND OTHER PARTIAL PRODUCTS AND SUMS FOR LATER CALCULATIONS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL PENMAN ( SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,F,T24,S,Q2, & Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,DQSDT2) IF(SATFLG.LE.0.5) ETP = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALL CANRES to calculate the canopy resistance and convert it C into PC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(SHDFAC.GT.1.E-6) THEN CALL CANRES(SOLDN,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, & SMCWLT,SMCREF,RCMIN,RC,PC,NROOT,Q2SAT,DQSDT2, & TOPT,RSMAX,RGL,HS) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF POINT IS TO BE FULLY PROCESSED.... C BASED ON WHETHER OR NOT THERE IS A SNOW PACK, CALL SNOPAC OR C NOPAC TO CALCULATE AN ESTIMATE OF ACTUAL EVAPOTRANSPIRATION, C UPDATE SOIL MOISTURE, AND UPDATE SOIL TEMPERATURE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( ESD .EQ. 0.0 ) THEN CALL NOPAC ( ETP, ETA, PRCP, SMC, SMCMAX, SMCWLT, & SMCREF,SMCDRY, CMC, CMCMAX, NSOIL, DT, SHDFAC, & Q1,Q2,T1,SFCTMP,T24,TH2,F,F1,S,STC, & EPSCA, B, PC, RCH, RR, CFACTR, & ZSOIL, DKSAT, DWSAT, TBOT,RUNOFF1,RUNOFF2, & RUNOFF3, EDIR1, EC1, ETT1,NROOT,ICE) ELSE CALL SNOPAC ( ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, & SMCREF, SMCDRY, CMC, CMCMAX, NSOIL, DT, Q1, & Q2,T1,SFCTMP,T24,TH2,F,F1,S,STC,EPSCA,SFCPRS, & B, PC, RCH, RR, CFACTR, ESD, & ZSOIL, DWSAT, DKSAT, TBOT, SHDFAC,RUNOFF1, & RUNOFF2,RUNOFF3,EDIR1,EC1,ETT1,NROOT,SNMAX,ICE) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RECONVERT THE LIQUID EQUIVALENT SNOW DEPTH (METERS) TO AN C OUTPUT SNOW DEPTH (METERS) USING A 10 TO 1 RATIO. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SNODEP = ESD * 10.0 SNODEP = ESD CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC SUB-PRODUCTS THEN USE TO CALC SENSIBLE HEAT FLUX. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC HEAT = -CHKFF * ( TH2 - T1 ) T14 = T1 * T1 * T1 * T1 RR = T14 * 6.48E-8 / ( SFCPRS * CH ) + 1.0 H = HEAT C C Convert ETA from KG M-2 S-1 to W M-2 C ETA = ETA*2.5E+6 ETP = ETP*2.5E+6 C Convert the sign of soil heat flux so that: C S>0: warm the surface (night time) C S<0: cool the surface (day time) S=-1.0*S C C Convert RUNOFF3 from M to M S-1 C RUNOFF3 = RUNOFF3/DT RUNOFF2 = RUNOFF2+RUNOFF3 SOILM=-1.0*SMC(1)*ZSOIL(1) DO K = 2, NSOIL SOILM=SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K)) ENDDO SOILWM=-1.0*(SMCMAX-SMCWLT)*ZSOIL(1) SOILWW=-1.0*(SMC(1)-SMCWLT)*ZSOIL(1) DO K = 2, NROOT SOILWM=SOILWM+(SMCMAX-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K)) SOILWW=SOILWW+(SMC(K)-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K)) ENDDO SOILW=SOILWW/SOILWM C RETURN END C C SUBROUTINE CANRES(SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, & SMCWLT,SMCREF,RCMIN,RC,PC,NROOT,Q2SAT,DQSDT2, & TOPT,RSMAX,RGL,HS) C ********************************************************************** C SUBROUTINE CANRES C ----------------- C THIS ROUTINE CALCULATES THE CANOPY RESISTANCE WHICH DEPENDS C ON INCOMING SOLAR RADIATION, AIR TEMPERATURE, ATMOSPHERIC C WATER VAPOR PRESSURE DEFICIT AT THE LOWEST MODEL LEVEL, ANS SOIL C MOISTURE. C ---------------------------------------------------------------------- C SOURCE: JARVIS (1976), Jacquemin and Noilhan (1990 BLM) C ---------------------------------------------------------------------- C ---------------------------------------------------------------------- C INPUT: SOLAR: incoming solar radiation C CH: surface exchange coefficient for heat and moisture C SFCTMP: air temperature at 1st level C Q2: air humidity at 1st level C SFCPRS: surface pressure C SMC: volumetric soil moisture C ZSOIL: soil depth C NSOIL: number of soil layers C SMCWLT: wilting point C SMCREF: reference soil moisture (field capacity) C C OUTPUT: PC: plant coefficient C RC: canopy resistance C ---------------------------------------------------------------------- C ********************************************************************** PARAMETER ( NSOLD = 4 ) DIMENSION ZSOIL(NSOIL), SMC(NSOIL), PART(NSOLD) INTEGER DATE, DDTIME PARAMETER (SIGMA=5.67E-8, RD=287.04, CP=1004.6, SLV=2.5E6) C ---------------------------------------------------------------------- C FOR LOW CROPS USE THESE VALUES FROM JACQUEMIN AND NOILHAN (90 BLM) C ---------------------------------------------------------------------- RCS = 0.0 RCT = 0.0 RCQ = 0.0 RCSOIL = 0.0 RC = 0.0 C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO INCOMING SOLAR RADIATION. C ---------------------------------------------------------------------- FF = 0.55*2.0*SOLAR/RGL RCS = (FF + RCMIN/RSMAX) / (1.0 + FF) RCS = AMAX1(RCS,0.0001) C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO THE TEMPERATURE AT THE FIRST MODEL LEVEL. C ---------------------------------------------------------------------- RCT = 1.0 - 0.0016*((TOPT-SFCTMP)**2.0) RCT = AMAX1(RCT,0.0001) C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT THE FIRST MODEL LEVEL. C ---------------------------------------------------------------------- P = SFCPRS QS = Q2SAT C RCQ expression from Ssib RCQ = 1.0/(1.0+HS*(QS-Q2)) RCQ = AMAX1(RCQ,0.01) C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY. C DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP. C ---------------------------------------------------------------------- GX = (SMC(1) - SMCWLT) / (SMCREF - SMCWLT) IF (GX .GT. 1.) GX = 1. IF (GX .LT. 0.) GX = 0. PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX DO 10 K = 2, NROOT GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT) IF (GX .GT. 1.) GX = 1. IF (GX .LT. 0.) GX = 0. PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX 10 CONTINUE DO 20 K = 1, NROOT RCSOIL = RCSOIL+PART(K) 20 CONTINUE RCSOIL = AMAX1(RCSOIL,0.0001) C ---------------------------------------------------------------------- C DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. C CONVERT CANOPY RESISTANCE (RC) TO PLANT COEFFICIENT (PC). C ---------------------------------------------------------------------- RC = RCMIN/(RCS*RCT*RCQ*RCSOIL) TAIR4 = SFCTMP**4 ST1 = (4.*SIGMA*RD)/CP SLVCP = SLV/CP RR = ST1*TAIR4/(SFCPRS*CH) + 1.0 DELTA = SLVCP*DQSDT2 PC = (RR+DELTA)/(RR*(1.+RC*CH)+DELTA) RETURN END C FUNCTION DEVAP ( ETP1, SMC, ZSOIL, SHDFAC, SMCMAX, B, & DKSAT, DWSAT, SMCDRY, SMCREF, SMCWLT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC NAME: DIRECT EVAPORATION (DEVAP) FUNCTION VERSION: N/A CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL B REAL DEVAP REAL DKSAT REAL DWSAT REAL ETP1 REAL FX REAL SHDFAC REAL SMC REAL SMCDRY REAL SMCMAX REAL WCND REAL WDF REAL ZSOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GET SOIL H2O DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY GIVEN SMC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL WDFCND ( WDF, WCND, SMC, SMCMAX, B, DKSAT, DWSAT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FX REPRESENTS SOIL MOISTURE FLUX / ATMOSPHERIC DEMAND CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FX = ( -2. * WDF * ( SMC - SMCDRY ) / ZSOIL - WCND ) / ETP1 C*** test of BETA Methode, F. Chen 7/11/96 *** FX = (SMC - SMCWLT) / (SMCREF - SMCWLT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FX > 1 REPRESENTS DEMAND CONTROL C FX < 1 REPRESENTS FLUX CONTROL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FX = MAX ( MIN ( FX, 1. ) ,0. ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DEVAP = FX * ( 1.0 - SHDFAC ) * ETP1 RETURN END SUBROUTINE HRT ( RHSTS, STC, SMC, SMCMAX, NSOIL, ZSOIL, YY, ZZ1, & TBOT, B, F1, DF1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY CC ======= TERM OF THE SOIL THERMAL DIFFUSION EQUATION. ALSO TO CC COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE CC TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD = 4 ) INTEGER NSOIL REAL AI ( NSOLD ) REAL BI ( NSOLD ) REAL CAIR REAL CH2O REAL CI ( NSOLD ) REAL CSOIL REAL DDZ REAL DDZ2 REAL DENOM REAL DFMAX1 REAL DF1 REAL DFMAXN REAL DTSDZ REAL DTSDZ2 REAL F1 REAL HCPCT INTEGER K REAL MXSMC REAL RHSTS ( NSOIL ) REAL S REAL SMC ( NSOIL ) REAL SMCMAX REAL STC ( NSOIL ) REAL TBOT REAL YY REAL ZBOT REAL ZSOIL ( NSOIL ) REAL ZZ1 !!$omp threadlocal /abci/ COMMON /ABCI/ AI, BI, CI PARAMETER(CAIR=1004.6,CH2O=4.2 E6,CSOIL=1.26E6,ZBOT=-3.0) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE THERMAL DIFFUSIVITY FOR THE WETTER OF THE TOP SOIL C LAYER OR THE NEXT LOWER SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MXSMC = MAX ( SMC(1), SMC(2) ) CALL TDFCND ( DFMAX1, MXSMC, B, F1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC HCPCT = SMC(1)*CH2O + (1.0-SMCMAX)*CSOIL + (SMCMAX-SMC(1))*CAIR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ = 1.0 / ( -0.5 * ZSOIL(2) ) AI(1) = 0.0 CI(1) = ( DFMAX1 * DDZ ) / ( ZSOIL(1) * HCPCT ) BI(1) = -CI(1) + DF1 / ( 0.5 * ZSOIL(1) * ZSOIL(1)*HCPCT*ZZ1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL C LAYERS. RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT C AND FLUX TO CALC RHSTS FOR THE TOP SOIL LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) ) S = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 ) RHSTS(1) = ( DFMAX1 * DTSDZ - S ) / ( ZSOIL(1) * HCPCT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE DDZ2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 40 K = 2, NSOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THIS SOIL LAYER'S HEAT CAPACITY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC HCPCT = SMC(K)*CH2O +(1.0-SMCMAX)*CSOIL +(SMCMAX-SMC(K))*CAIR IF ( K .NE. NSOIL ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE THERMAL DIFFUSIVITY FOR THE WETTER OF THE C CURRENT OR THE NEXT LOWER SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MXSMC = MAX ( SMC(K), SMC(K+1) ) CALL TDFCND ( DFMAXN, MXSMC, B, F1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) ) DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1)) CI(K) = -DFMAXN * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT) ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE THERMAL DIFFUSIVITY FOR THE LOWEST SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL TDFCND ( DFMAXN, SMC(K), B, F1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTSDZ2 = (STC(K)-TBOT) / (.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SET MATRIX COEF, CI TO ZERO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CI(K) = 0. END IF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT RHSTS(K) = ( DFMAXN * DTSDZ2 - DFMAX1 * DTSDZ ) / DENOM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC AI(K) = - DFMAX1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT) BI(K) = -(AI(K) + CI(K)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RESET VALUES OF DFMAX1, DTSDZ, AND DDZ FOR LOOP TO NEXT SOIL LYR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DFMAX1 = DFMAXN DTSDZ = DTSDZ2 DDZ = DDZ2 40 CONTINUE RETURN END SUBROUTINE HRTICE (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, & TBOT, B, F1, DF1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY CC ======= TERM OF THE SOIL THERMAL DIFFUSION EQUATION. ALSO TO CC COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE CC TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD = 4 ) INTEGER NSOIL REAL AI ( NSOLD ) REAL BI ( NSOLD ) REAL CAIR REAL CH2O REAL CI ( NSOLD ) REAL CSOIL REAL DDZ REAL DDZ2 REAL DENOM REAL DFMAX1 REAL DF1 REAL DFMAXN REAL DTSDZ REAL DTSDZ2 REAL F1 REAL HCPCT INTEGER K REAL MXSMC REAL RHSTS ( NSOIL ) REAL S REAL SMC ( NSOIL ) REAL SMCMAX REAL STC ( NSOIL ) REAL TBOT REAL YY REAL ZBOT REAL ZSOIL ( NSOIL ) REAL ZZ1 !!$omp threadlocal /abci/ COMMON /ABCI/ AI, BI, CI PARAMETER(CAIR=1004.6,CH2O=4.2 E6,CSOIL=1.26E6) ZBOT=ZSOIL(NSOIL) DFMAX1=2.2 HCPCT=1880.0*917.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ = 1.0 / ( -0.5 * ZSOIL(2) ) AI(1) = 0.0 CI(1) = ( DFMAX1 * DDZ ) / ( ZSOIL(1) * HCPCT ) BI(1) = -CI(1) + DF1/( 0.5 * ZSOIL(1) * ZSOIL(1) * HCPCT * ZZ1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL C LAYERS. RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT C AND FLUX TO CALC RHSTS FOR THE TOP SOIL LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) ) S = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 ) RHSTS(1) = ( DFMAX1 * DTSDZ - S ) / ( ZSOIL(1) * HCPCT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE DDZ2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 40 K = 2, NSOIL IF ( K .NE. NSOIL ) THEN DFMAXN=2.2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) ) DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1)) CI(K) = -DFMAXN * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT) ELSE DFMAXN=2.2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTSDZ2 = (STC(K)-TBOT)/(.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SET MATRIX COEF, CI TO ZERO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CI(K) = 0. END IF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT RHSTS(K) = ( DFMAXN * DTSDZ2 - DFMAX1 * DTSDZ ) / DENOM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC AI(K) = - DFMAX1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT) BI(K) = -(AI(K) + CI(K)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RESET VALUES OF DFMAX1, DTSDZ, AND DDZ FOR LOOP TO NEXT SOIL LYR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DFMAX1 = DFMAXN DTSDZ = DTSDZ2 DDZ = DDZ2 40 CONTINUE RETURN END SUBROUTINE HSTEP ( STC, RHSTS, DT, NSOIL ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD=4 ) INTEGER K INTEGER NSOIL REAL AI ( NSOLD ) REAL BI ( NSOLD ) REAL CI ( NSOLD ) REAL DT REAL RHSTS ( NSOIL ) REAL STC ( NSOIL ) !!$omp threadlocal /abci/ COMMON /ABCI/ AI, BI, CI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 10 K = 1 , NSOIL RHSTS(K) = RHSTS(K) * DT AI(K) = AI(K) * DT BI(K) = 1. + BI(K) * DT CI(K) = CI(K) * DT 10 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SOLVE THE TRI-DIAGONAL MATRIX EQUATION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ROSR12 ( CI,AI,BI,CI,RHSTS,RHSTS,NSOIL ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 20 K = 1 , NSOIL STC(K) = STC(K) + CI(K) 20 CONTINUE RETURN END SUBROUTINE NOPAC ( ETP, ETA, PRCP, SMC, SMCMAX, SMCWLT, & SMCREF,SMCDRY,CMC,CMCMAX, NSOIL, DT, SHDFAC, & Q1, Q2, T1, SFCTMP, T24, TH2, F, F1, S, STC, & EPSCA, B, PC, RCH, RR, CFACTR, ZSOIL, DKSAT, & DWSAT, TBOT, RUNOFF1, RUNOFF2,RUNOFF3, EDIR1, & EC1, ETT1, NROOT, ICE) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE CC ======= SOIL MOISTURE CONTENT AND SOIL HEAT CONTENT VALUES FOR CC THE CASE WHEN NO SNOW PACK IS PRESENT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOIL REAL B REAL BETA REAL CFACTR REAL CMC REAL CMCMAX REAL CP REAL DEW REAL DF1 REAL DKSAT REAL DRIP REAL DT REAL DWSAT REAL EC REAL EDIR REAL EPSCA REAL ETA REAL ETA1 REAL ETP REAL ETP1 REAL ETT REAL F REAL F1 REAL FLX1 REAL FLX2 REAL FLX3 REAL PC REAL PRCP REAL PRCP1 REAL Q2 REAL RCH REAL RHO REAL RIB REAL RR REAL RUNOFF REAL S REAL SFCTMP REAL SHDFAC REAL SIGMA REAL SMC ( NSOIL ) REAL SMCDRY REAL SMCMAX REAL SMCREF REAL SMCWLT REAL STC ( NSOIL ) REAL T1 REAL T24 REAL TBOT REAL TH2 REAL YY REAL YYNUM REAL ZSOIL ( NSOIL ) REAL ZZ1 C !!$omp threadlocal /rite/ COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF, & DEW,RIB PARAMETER(CP=1004.6, SIGMA=5.67E-8) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C EXECUTABLE CODE BEGINS HERE..... C CONVERT ETP FROM KG M-2 S-1 TO MS-1 AND INITIALIZE DEW. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PRCP1 = PRCP * 0.001 ETP1 = ETP * 0.001 DEW = 0.0 IF ( ETP .GT. 0.0 ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT PRCP FROM KG M-2 S-1 TO M S-1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC, & CMCMAX,SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, & EDIR1, EC1, ETT1, SFCTMP,Q2,NROOT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ETA = ETA1 * 1000.0 ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW C AND REINITIALIZE ETP1 TO ZERO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DEW = -ETP1 ETP1 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT PRCP FROM KG M-2 S-1 TO M S-1 AND ADD DEW AMT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PRCP1 = PRCP1 + DEW C CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC, & CMCMAX,SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, & EDIR1, EC1, ETT1, SFCTMP, Q2, NROOT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ETA = ETA1 * 1000.0 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C BASED ON ETP AND E VALUES, DETERMINE BETA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( ETP .LE. 0.0 ) THEN BETA = 0.0 IF ( ETP .LT. 0.0 ) THEN BETA = 1.0 ETA = ETP ENDIF ELSE BETA = ETA / ETP ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC MIXING RATIO AT GRND LEVEL (SKIN) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Q1 = Q2 + ETA * CP / RCH CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR, C CALC ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN C CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL TDFCND ( DF1, SMC(1), B, F1 ) YYNUM = F - SIGMA * T24 YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1 CALL SHFLX ( S,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,TBOT, & B,F1,DF1, ICE ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SET FLX1, AND FLX3 TO ZERO SINCE THEY ARE NOT USED. FLX2 C WAS SIMILARLY INITIALIZED IN THE PENMAN ROUTINE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FLX1 = 0.0 FLX3 = 0.0 C RETURN END SUBROUTINE PENMAN(SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,F,T24,S,Q2, & Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,DQSDT2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. CC ======= VARIOUS PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND CC PASSED BACK TO THE CALLING ROUTINE FOR LATER USE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC LOGICAL SNOWNG LOGICAL FRZGRA REAL A REAL BETA REAL CH REAL CP REAL CPH2O REAL CPICE REAL DELTA REAL DEW REAL DRIP REAL EC REAL EDIR REAL ELCP REAL EPSCA REAL ETP REAL ETT REAL F REAL FLX1 REAL FLX2 REAL FLX3 REAL FNET REAL LSUBC REAL LSUBF REAL PRCP REAL Q2 REAL Q2SAT REAL R REAL RAD REAL RCH REAL RHO REAL RIB REAL RR REAL RUNOFF REAL S REAL SFCPRS REAL SFCTMP REAL SIGMA REAL T24 REAL T2V REAL TH2 !!$omp threadlocal /rite/ COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF, & DEW,RIB PARAMETER(CP=1004.6,CPH2O=4.218E+3,CPICE=2.106E+3,R=287.04, & ELCP=2.4888E+3,LSUBF=3.335E+5,LSUBC=2.5E+6,SIGMA=5.67E-8) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C EXECUTABLE CODE BEGINS HERE...DEFINE FUNCTION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC TT(T) = T * T FLX2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DELTA = ELCP * DQSDT2 T24 = TT( SFCTMP ) * TT( SFCTMP ) RR = T24 * 6.48E-8 / ( SFCPRS * CH ) + 1.0 RHO = SFCPRS / ( R * T2V ) RCH = RHO * CP * CH CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT C EFFECTS CAUSED BY FALLING PRECIPITATION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( .NOT. SNOWNG ) THEN IF ( PRCP .GT. 0.0 ) RR = RR + CPH2O * PRCP / RCH ELSE RR = RR + CPICE * PRCP / RCH ENDIF FNET = F - SIGMA * T24 - S CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO C ICE ON IMPACT IN THE CALCULATION OF FLX2 AND FNET. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( FRZGRA ) THEN FLX2 = -LSUBF * PRCP FNET = FNET - FLX2 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FINISH PENMAN EQUATION CALCULATIONS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC RAD = FNET / RCH + TH2 - SFCTMP A = ELCP * ( Q2SAT - Q2 ) EPSCA = ( A * RR + RAD * DELTA ) / ( DELTA + RR ) ETP = EPSCA * RCH / LSUBC RETURN END SUBROUTINE ROSR12 ( P, A, B, C, D, DELTA, NSOIL ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN CC ======= BELOW: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER K INTEGER KK INTEGER NSOIL REAL P(*) REAL A(*) REAL B(*) REAL C(*) REAL D(*) REAL DELTA(*) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C(NSOIL) = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SOLVE THE COEFS FOR THE 1ST SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC P(1) = -C(1) / B(1) DELTA(1) = D(1) / B(1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 10 K = 2 , NSOIL P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) ) DELTA(K) = (D(K)-A(K)*DELTA(K-1))*(1.0/(B(K)+A(K)*P(K-1))) 10 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SET P TO DELTA FOR LOWEST SOIL LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC P(NSOIL) = DELTA(NSOIL) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ADJUST P FOR SOIL LAYERS 2 THRU NSOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 20 K = 2 , NSOIL KK = NSOIL - K + 1 P(KK) = P(KK) * P(KK+1) + DELTA(KK) 20 CONTINUE RETURN END SUBROUTINE SHFLX ( S,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,TBOT, & B,F1,DF1,ICE ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SOIL HEAT FLUX. SOIL TEMPERATURE CONTENT CC ======= (AN ARRAY CONTAINING SOIL TEMPERATURES) IS ALSO UPDATED. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD = 4 ) INTEGER NSOIL REAL B REAL DF1 REAL DT REAL F1 REAL RHSTS ( NSOLD ) REAL S REAL SMC ( NSOIL ) REAL SMCMAX REAL STC ( NSOIL ) REAL T1 REAL TBOT REAL YY REAL ZSOIL ( NSOIL ) REAL ZZ1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(ICE.EQ.1) THEN CALL HRTICE(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,B, * F1,DF1) ELSE CALL HRT(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,B,F1,DF1) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C HSTEP ROUTINE CALCS/UPDATES SOIL TEMPS BASED ON RHSTS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL HSTEP ( STC,RHSTS,DT,NSOIL ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE GRND (SKIN) TEMPERATURE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE SFC SOIL HEAT FLUX CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC S = DF1 * (STC(1) - T1) / (0.5 * ZSOIL(1)) RETURN END SUBROUTINE SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC,CMCMAX, & SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, EDIR1, EC1, & ETT1, SFCTMP,Q2,NROOT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CC ======= CONTENT (SMC - A PER UNIT VOLUME MEASUREMENT) IS A CC DEPENDENT VARIABLE THAT IS UPDATED WITH PROGNOSTIC EQNS. CC THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD = 4 ) INTEGER K INTEGER NSOIL REAL B REAL BETA REAL CFACTR REAL CMC REAL CMCMAX REAL DEW REAL DKSAT REAL DRIP REAL DT REAL DWSAT REAL EC REAL EDIR REAL ET ( NSOLD ) REAL ETA1 REAL ETP1 REAL ETT REAL EXCESS REAL FLX1 REAL FLX2 REAL FLX3 REAL PC REAL PCPDRP REAL PRCP1 REAL RHO REAL RHSCT REAL RHSTT ( NSOLD ) REAL RIB REAL RUNOF REAL RUNOFF REAL SHDFAC REAL SMC ( NSOIL ) REAL SMCA ( NSOLD ) REAL SMCDRY REAL SMCFG ( NSOLD ) REAL SMCMAX REAL SMCREF REAL SMCWLT REAL TRHSCT REAL ZSOIL ( NSOIL ) !!$omp threadlocal /rite/ COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOF, & DEW,RIB CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C EXECUTABLE CODE BEGINS HERE....IF THE POTENTIAL EVAPOTRANS- C PIRATION IS GREATER THAN ZERO... CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DUMMY=0. EDIR = 0. EC = 0. ETT = 0. DO 10 K = 1, NSOIL ET ( K ) = 0. 10 CONTINUE C IF ( ETP1 .GT. 0.0 ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EDIR = DEVAP ( ETP1, SMC(1), ZSOIL(1), SHDFAC, SMCMAX, & B, DKSAT, DWSAT, SMCDRY,SMCREF, SMCWLT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT C TRANSPIRATION, AND ACCUMULATE IT FOR ALL SOIL LAYERS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ETT = 0. IF(SHDFAC.GT.0.0) THEN CALL TRANSP ( ET, NSOIL, ETP1, SMC, CMC, ZSOIL, SHDFAC, & SMCWLT, CMCMAX, PC, CFACTR, SMCREF, SFCTMP,Q2,NROOT) DO 20 K = 1 , NSOIL ETT = ETT + ET ( K ) 20 CONTINUE ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE CANOPY EVAPORATION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1 C******** EC should be limited by the total amount of available C water on the canopy. Modified by F.CHEN on 10/18/94 C******** CMC2MS=CMC/DT EC=MIN(CMC2MS,EC) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EDIR1=EDIR EC1=EC ETT1=ETT ETA1 = EDIR + ETT + EC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC RHSCT = SHDFAC * PRCP1 - EC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT RHSCT (A RATE) TO TRHSCT (AN AMT) AND ADD IT TO EXISTING C CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP C AND WILL FALL TO THE GRND. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DRIP = 0. TRHSCT = DT * RHSCT EXCESS = CMC + TRHSCT IF ( EXCESS .GT. CMCMAX ) DRIP = EXCESS - CMCMAX CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT C GOES INTO THE SOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PCPDRP = (1. - SHDFAC) * PRCP1 + DRIP / DT CC print*,' **************** SMLX ******************' CC print*,' PCPDRP=', PCPDRP, ' EDIR=', EDIR,' ET=', ET, CC * 'SMC(1)=', SMC(1), 'SMC(2)=', SMC(2), ' PRCP1=', PRCP1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE C TENDENCY EQUATIONS. ITERATE TWICE THRU THE CALLS. THIS ELIMINATES C 2-DELTA-T OSCILLATIONS IN THE SMC VALUES EVEN WITH TIMESTEPS OF C UP TO 1 HOUR. THE 1ST SET OF CALLS YIELDS 1ST GUESS SMC VALUES, C THE 2ND SET OF CALLS YIELDS FINAL SMC VALUES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( PCPDRP .GT. 0.0 ) THEN CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMC,NSOIL,PCPDRP,ZSOIL, & DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT) CALL SSTEP ( SMCFG,SMC,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, & CMCMAX, RUNOFF3, ZSOIL) DO 30 K = 1, NSOIL SMCA(K) = ( SMC(K) + SMCFG(K) ) * 0.5 30 CONTINUE CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMCA,NSOIL,PCPDRP,ZSOIL, & DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT) CALL SSTEP ( SMC,SMC,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, & CMCMAX, RUNOFF3, ZSOIL) ELSE CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMC,NSOIL,PCPDRP,ZSOIL, & DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT) CALL SSTEP ( SMC,SMC,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, & CMCMAX, RUNOFF3, ZSOIL) ENDIF RUNOF = RUNOFF RETURN END SUBROUTINE SNOPAC ( ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, & SMCREF, SMCDRY, CMC, CMCMAX, NSOIL, DT, Q1, & Q2,T1,SFCTMP,T24,TH2,F,F1,S,STC,EPSCA,SFCPRS, & B, PC, RCH, RR, CFACTR, ESD, & ZSOIL, DWSAT, DKSAT, TBOT, SHDFAC, RUNOFF1, & RUNOFF2,RUNOFF3,EDIR1,EC1,ETT1,NROOT,SNMAX,ICE) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE CC ======= SOIL MOISTURE CONTENT AND SOIL HEAT CONTENT VALUES FOR CC THE CASE WHEN A SNOW PACK IS PRESENT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOIL LOGICAL SNOWNG REAL B REAL BETA REAL CFACTR REAL CMC REAL CMCMAX REAL CP REAL CPH2O REAL CPICE REAL DENOM REAL DEW REAL DF1 REAL DKSAT REAL DRIP REAL DT REAL DWSAT REAL EC REAL EDIR REAL EPSCA REAL ESD REAL ETA REAL ETA1 REAL ETP REAL ETP1 REAL ETP2 REAL ETT REAL EX REAL F REAL FLX1 REAL FLX2 REAL FLX3 REAL F1 REAL LSUBF REAL LSUBC REAL LSUBS REAL PC REAL PRCP REAL PRCP1 REAL Q1 REAL Q2 REAL RCH REAL RHO REAL RIB REAL RR REAL RUNOFF REAL S REAL S1 REAL SFCTMP REAL SHDFAC REAL SIGMA REAL SMC ( NSOIL ) REAL SMCDRY REAL SMCMAX REAL SMCREF REAL SMCWLT REAL SNMAX REAL STC ( NSOIL ) REAL T1 REAL T11 REAL T12 REAL T12A REAL T12B REAL T24 REAL TBOT REAL TH2 REAL YY REAL ZSOIL( NSOIL ) REAL ZZ1 !!$omp threadlocal /rite/ COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF, & DEW,RIB PARAMETER(CP=1004.6,CPH2O=4.218E+3,CPICE=2.106E+3, & LSUBF=3.335E+5,LSUBC=2.5E+6,LSUBS=2.83E+6,SIGMA=5.67E-8) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C EXECUTABLE CODE BEGINS HERE... C CONVERT ETP FM KG M-2 S-1 TO M S-1 AND THEN TO AN AMOUNT C (M) AND CALL IT AN EFFECTIVE SNOWPACK REDUCTION AMT, ETP2 (M). C THIS IS THE AMT THE SNOWPACK WOULD BE REDUCED DUE TO EVAPORATION C FROM THE SNOW SFC DURING THE TIMESTEP. EVAPORATION WILL C PROCEED AT THE POTNTL RATE UNLESS THE SNOW DEPTH IS LESS THAN C THE EXPECTED SNOWPACK REDUCTION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PRCP1=PRCP1*0.001 ETP2 = ETP * 0.001 * DT BETA = 1.0 IF(ICE.NE.1) THEN IF ( ESD .LT. ETP2 ) THEN BETA = ESD / ETP2 ENDIF ENDIF DEW = 0.0 IF ( ETP .LT. 0.0 ) THEN DEW = -ETP * 0.001 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC AN 'EFFECTIVE SNOW-GRND SFC TEMP' BASED ON HEAT FLUXES C BTWN THE SNOW PACK AND THE SOIL AND ON NET RADIATION. THE C CONSTANT 0.35 (W M-1 K-1) IS THE THERMAL DIFFUSIVITY OF SNOW, CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC T12A = ((F - SIGMA * T24) / RCH+TH2-SFCTMP-BETA*EPSCA) / RR T12B = 0.35 * STC(1) / ( 10.0 * ESD * RR * RCH ) DENOM = 1.0 + 0.35 / ( 10.0 * ESD * RR * RCH ) T12 = (SFCTMP + T12A + T12B ) / DENOM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RECALCULATE THE HEAT FLUX INTO OR OUT OF THE SOIL, C AND SET THE EFFECTIVE POTENTIAL EVAPOTRANSP TO ZERO. C THE CONSTANT 0.35 (W M-1 K-1) IS THE THERMAL DIFFUSIVITY OF SNOW. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC S = 0.35 * ( T12 - STC(1) ) / ( 10.0 * ESD ) S = MIN ( MAX ( S,-100.0 ), 100.0 ) IF ( T12 .LE. 273.16 ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BLO FREEZING, C NO SNOW MELT WILL OCCUR. SET THE SKIN TEMP TO THIS C EFFECTIVE TEMP AND SET THE EFFECTIVE PRECIP TO ZERO. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ESD = AMAX1 ( 0.0 , ESD - ETP2 ) T1 = T12 FLX3 = 0.0 EX = 0.0 SNMAX = 0.0 ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW C MELT WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNMAX. C REVISE THE EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE C IT WOULD HAVE CHGD DUE TO THE LATENT HEAT RELEASED BY THE C MELTING. CALC THE LATENT HEAT RELEASED, FLX3. SET THE EFFECTIVE C PRECIP, PRCP1 TO THE SNOW MELT RATE, EX FOR USE IN SMFLX. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C This is the old version of snow melt C EX = 0.001 * RCH * RR * ( T12 - 273.16 ) / LSUBF C SNMAX = EX * DT C EX = AMIN1 ( ESD / DT, EX ) C ESD = AMAX1 ( ESD - SNMAX, 0.0 ) C T1 = AMAX1 ( 273.16, (T12-(LSUBF * EX * 1000.0 / (RCH*RR)))) C FLX3 = LSUBF * EX * 1000.0 C PRCP1 = PRCP1 + EX C Below is the correction of evaporation considering snow melting C F. Chen (96/3/18) T1 = 273.16 QSAT = (0.622*6.11E2)/(SFCPRS-0.378*6.11E2) ETP = RCH*(QSAT-Q2)/CP ETP2 = ETP*0.001*DT BETA = 1.0 IF ( ESD .LE. ETP2 ) THEN BETA = ESD / ETP2 ESD = 0.0 SNMAX = 0.0 EX = 0.0 ELSE ESD = MAX ( 0.0 , ESD - ETP2 ) ETP3 = ETP*LSUBC S = 0.35 * ( T1 - STC(1) ) / ( 10.0 * ESD ) S = MIN ( MAX ( S,-100.0 ), 100.0 ) SEH = RCH*(T1-TH2) T14 = T1*T1 T14 = T14*T14 FLX3 = F - SIGMA*T14 - S - SEH - ETP3 IF(FLX3.LE.0.0) FLX3=0.0 EX = FLX3*0.001/LSUBF SNMAX = EX * DT ENDIF IF(SNMAX.LT.ESD) THEN ESD = ESD - SNMAX ELSE EX = ESD/DT SNMAX = ESD ESD = 0.0 FLX3 = EX*1000.0*LSUBF C CALL TDFCND ( DF1, SMC(1), B, F1 ) C YYNUM = F - SIGMA * T24 C YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-(FLX3+ETP3)) / RR C ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0 C T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1 C S = DF1 * ( STC(1) - T1 ) / ( 0.5 * ZSOIL(1) ) ENDIF PRCP1 = PRCP1 + EX ENDIF C SET THE EFFECTIVE POTENTIAL EVAPOTRANSP TO ZERO ETP1 = 0.0 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SMFLX RETURNS SOIL MOISTURE VALUES AND PRELIMINARY VALUES OF C EVAPOTRANSPIRATION. IN THIS, THE SNOW PACK CASE, THE PRELIM- C VALUES (ETA1) ARE NOT USED IN SUBSEQUENT CALCULATION OF E. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(ICE.NE.1) THEN CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT, & SMCREF,SHDFAC,CMCMAX,SMCDRY,CFACTR,RUNOFF1,RUNOFF2, & RUNOFF3, EDIR1, EC1, ETT1,SFCTMP,Q2,NROOT) ENDIF ETA = BETA * ETP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC/UPDATE THE GRND SFC (SKIN) MIXING RATIO. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Q1 = Q2 + ETA * CP / RCH CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE SOIL TEMP DIFFUSIVITY/CONDUCTIVITY. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL TDFCND ( DF1, SMC(1), B, F1 ) IF(ICE.EQ.1) DF1=2.2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C THE 'ADJUSTED TOP SOIL LYR TEMP' ( YY ) AND THE 'ADJUSTED SOIL C HEAT FLUX' ( ZZ1 ) ARE SET TO THE TOP SOIL LYR TEMP, AND 1, C RESPECTIVELY. THESE ARE CLOSE-ENOUGH APPROXIMATIONS BECAUSE C THE SFC HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE C THE FLUX AT THE SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT C SINCE WE WILL NOT USE ITS VALUE AS REVISED BY SHFLX. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ZZ1 = 1.0 YY = STC(1)-0.5*S*ZSOIL(1)*ZZ1/DF1 T11 = T1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE THE SOIL HEAT FLUX C ( S1 ) AND THE SKIN TEMP ( T11 ) OUTPUT FROM THIS CALL ARE NOT C USED IN ANY SUBSEQUENT CALCULATIONS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL SHFLX ( S1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL,TBOT, & B,F1,DF1,ICE ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF PRECIP IS FALLING, CALC HEAT FLUX FROM SNOW SFC TO NEWLY C ACCUMULATING PRECIP. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FLX1 = 0.0 IF ( SNOWNG ) THEN FLX1 = CPICE * PRCP * ( T1 - SFCTMP ) ELSE IF ( PRCP .GT. 0.0 ) FLX1 = CPH2O * PRCP * ( T1 - SFCTMP ) ENDIF RETURN END SUBROUTINE SRT (RHSTT,RUNOFF,EDIR,ET,SMC,SMCA,NSOIL,PCPDRP, & ZSOIL,DWSAT,DKSAT,SMCMAX,B, RUNOFF1, RUNOFF2,DT, SMCWLT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY CC ======= TERM OF THE SOIL WATER DIFFUSION EQUATION. ALSO TO CC COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE CC TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD = 4 ) INTEGER K INTEGER NSOIL REAL AI ( NSOLD ) REAL B REAL BI ( NSOLD ) REAL CI ( NSOLD ) REAL DMAX ( NSOLD ) REAL DDZ REAL DDZ2 REAL DENOM REAL DENOM2 REAL DKSAT REAL DSMDZ REAL DSMDZ2 REAL DWSAT REAL EDIR REAL ET ( NSOIL ) REAL INFMAX REAL MXSMC REAL MXSMC2 REAL NUMER REAL PCPDRP REAL PDDUM REAL RHSTT ( NSOIL ) REAL RUNOFF REAL SMC ( NSOIL ) REAL SMCA ( NSOIL ) REAL SMCMAX REAL WCND REAL WCND2 REAL WDF REAL WDF2 REAL ZSOIL ( NSOIL ) !!$omp threadlocal /abci/ COMMON /ABCI/ AI, BI, CI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ************* Include the infiltration formule from Schaake and Koren C Model C CC MODIFIED BY Q DUAN CC REPLACE THE MAXIMUM INFILTRATION EQUATION IN OSU MODEL CC BY SCHAAKE/KOREN EXPRESSION CC CC INITIALIZE PARAMETERS AND VARIABLES CC IOHINF=1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PDDUM = PCPDRP RUNOFF1 = 0.0 IF ( PCPDRP .NE. 0.0 ) THEN CC++ Modified by Q. Duan, 5/16/94 C IF (IOHINF .EQ. 1) THEN C REFKDT = 3.0 C REFDK = 3.4341E-6 CFEI Modify the value of REFKDT and REFDK according to recent C PILPS2d and GSWP results REFDK=2.E-6 is the sat. DK. value C for the soil type 2 (Plast PILPS2Fei GSWP runs unssed REFKDT=0.75, REFDK=2.E-6). REFKDT=3.0 REFDK = 2.0E-6 DT1 = DT/86400. SMCAV = SMCMAX - SMCWLT DMAX(1)=-ZSOIL(1)*SMCAV DMAX(1)=DMAX(1)*(1.0 - (SMCA(1)-SMCWLT)/SMCAV) DD=DMAX(1) DO KS=2,NSOIL DMAX(KS)=(ZSOIL(KS-1)-ZSOIL(KS))*SMCAV DMAX(KS)=DMAX(KS)*(1.0 - (SMCA(KS)-SMCWLT)/SMCAV) DD=DD+DMAX(KS) ENDDO KDT = REFKDT * DKSAT/REFDK CC VAL = (1.-EXP(-KDT*SQRT(DT1))) C remove the SQRT by F. Chen, 07/08/96 VAL = (1.-EXP(-KDT*DT1)) DDT = DD*VAL PX = PCPDRP*DT IF(PX.LT.0.0) PX = 0.0 INFMAX = (PX*(DDT/(PX+DDT)))/DT INFMAX = MIN(INFMAX, DKSAT) C ELSE C INFMAX = -DWSAT*(SMCMAX-SMCA(1))/(0.5*ZSOIL(1))+DKSAT C IF (SMCA(1).GE.SMCMAX) INFMAX = DKSAT C C END IF C C*************** RUNOFF1: SURFACE RUNOFF ******* C IF ( PCPDRP .GT. INFMAX ) THEN RUNOFF1 = PCPDRP - INFMAX PDDUM = INFMAX END IF END IF C C********* Below is the original OSU model C C PDDUM = PCPDRP C print*,' ************** in SRT *******************' C print*,' PCPDRP=', PCPDRP C RUNOFF1 = 0.0 C IF ( PCPDRP .NE. 0.0 ) THEN C INFMAX = - DWSAT * (SMCMAX - SMCA(1))/(0.5 * ZSOIL(1)) + DKSAT C IF ( SMCA(1) .GE. SMCMAX ) INFMAX = DKSAT C IF ( PCPDRP .GT. INFMAX ) THEN C RUNOFF1 = PCPDRP - INFMAX C C*************** RUNOFF1: SURFACE RUNOFF ******* C C PDDUM = INFMAX C END IF C END IF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DETERMINE SOIL MOISTURE DIFFUSIVITY (WDF) & CONDUCTIVITY (WCND) C COEFS FOR TOP SOIL LYR USING THE HIGHER OF SMC FOR THE TOP C SOIL LAYER OR SMC FOR THE NEXT LOWER SOIL LAYER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MXSMC = MAX( SMCA(1), SMCA(2) ) CALL WDFCND ( WDF,WCND,MXSMC,SMCMAX,B,DKSAT,DWSAT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ = 1. / ( -.5 * ZSOIL(2) ) AI(1) = 0.0 BI(1) = WDF * DDZ / ( -ZSOIL(1) ) CI(1) = -BI(1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL C MOISTURE GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DSMDZ = ( SMC(1) - SMC(2) ) / ( -.5 * ZSOIL(2) ) RHSTT(1) = (WDF * DSMDZ + WCND - PDDUM + EDIR + ET(1))/ZSOIL(1) SSTT = WDF*DSMDZ + WCND + EDIR + ET(1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE DDZ2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 10 K = 2 , NSOIL DENOM2 = ( ZSOIL(K-1) - ZSOIL(K) ) IF ( K .NE. NSOIL ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE SOIL WATER DIFFUSIVITY FOR THE WETTER OF THE C CURRENT OR THE NEXT LOWER SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C MXSMC2 = MAX ( SMCA(K), SMCA(K+1) ) C ************ PILPS TEST ************* MXSMC2 = SMCA(K) CALL WDFCND ( WDF2,WCND2,MXSMC2,SMCMAX,B,DKSAT,DWSAT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DENOM = ( ZSOIL(K-1) - ZSOIL(K+1) ) DSMDZ2 = ( SMC(K) - SMC(K+1) ) / ( DENOM * 0.5 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 2.0 / DENOM CI(K) = -WDF2 * DDZ2 / DENOM2 ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC C CONDUCTIVITY FOR THIS LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL WDFCND( WDF2,WCND2,SMCA(NSOIL),SMCMAX,B,DKSAT,DWSAT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DSMDZ2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SET MATRIX COEF CI TO ZERO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CI(K) = 0.0 END IF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NUMER = (WDF2 * DSMDZ2) + WCND2 - (WDF * DSMDZ) - WCND + ET(K) RHSTT(K) = NUMER / (-DENOM2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC AI(K) = -WDF * DDZ / DENOM2 BI(K) = -( AI(K) + CI(K) ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(K.EQ.NSOIL) THEN C*************** RUNOFF2: ground water runoff *********** RUNOFF2 = WCND2 ENDIF IF ( K .NE. NSOIL ) THEN WDF = WDF2 WCND = WCND2 DSMDZ = DSMDZ2 DDZ = DDZ2 END IF 10 CONTINUE RETURN END SUBROUTINE SSTEP ( SMCOUT, SMCIN, CMC, RHSTT, RHSCT, DT, & NSOIL, SMCMAX, CMCMAX, RUNOFF3, ZSOIL ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE/UPDATE THE SOIL MOISTURE CONTENT VALUES CC ======= AND THE CANOPY MOISTURE CONTENT VALUES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER ( NSOLD = 4 ) INTEGER K INTEGER NSOIL REAL AI ( NSOLD ) REAL BI ( NSOLD ) REAL CI ( NSOLD ) REAL CMC REAL CMCMAX REAL DT REAL RHSCT REAL RHSTT ( NSOIL ) REAL SMCIN ( NSOIL ) REAL SMCMAX REAL SMCOUT ( NSOIL ) REAL ZSOIL(NSOIL) !!$omp threadlocal /abci/ COMMON /ABCI/ AI, BI, CI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE C TRI-DIAGONAL MATRIX ROUTINE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 10 K = 1 , NSOIL RHSTT(K) = RHSTT(K) * DT AI(K) = AI(K) * DT BI(K) = 1. + BI(K) * DT CI(K) = CI(K) * DT 10 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ROSR12 ( CI, AI, BI, CI, RHSTT, RHSTT, NSOIL ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A C NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ****************** RUNOFF3: Runoff within soil layers ******* RUNOFF3 = 0. DO 20 K = 1 , NSOIL SMCOUT(K) = SMCIN(K) + CI(K) IF(SMCOUT(K).GT.SMCMAX) THEN IF(K.EQ.1) THEN DDZ = -ZSOIL(1) ELSE KK11 = K -1 DDZ = -ZSOIL(K) + ZSOIL(KK11) ENDIF RUNOFF3 = RUNOFF3 - ((SMCMAX-SMCOUT(K))*DDZ) ENDIF SMCOUT(K) = MAX ( MIN ( SMCOUT(K), SMCMAX ), 0.02 ) 20 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT RHSCT TO AN 'AMOUNT' VALUE AND ADD TO PREVIOUS C CMC VALUE TO GET NEW CMC VALUE. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CMC = CMC + DT * RHSCT c CMC = MIN ( MAX ( CMC, 0.), CMCMAX ) c c code changed 8-26-98 by Manikin per suggestions of Black/Mitchell c NOTE CALCULATION OF RHSCT (and DRIP) PREVIOUSLY IN SUBROUTINE SMFLX c IF (CMC.LT.1.E-20) CMC=0.0 CMC = MIN(CMC,CMCMAX) RETURN END SUBROUTINE TDFCND ( DF, SMC, B, F1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF CC ======= THE SOIL FOR A GIVEN POINT AND TIME. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL B REAL DF REAL F1 REAL PF REAL SMC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT C OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ccc SMC=0.27 IF ( SMC .GT. 0. ) THEN PF = F1 - B * ALOG10 ( SMC ) ELSE PF = 5.2 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF THE PARTIAL SUM/PRODUCT IS LOWER THAN OR EQUAL TO 5.1, COMPUTE C A THERMAL CONDUCTIVITY/DIFFUSIVITY. CAP ITS VALUE AT 1.9. C OTHERWISE SET THE THERMAL CONDUCTIVITY/DIFFUSIVITY TO A CONSTANT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( PF .LE. 5.1 ) THEN DF = EXP ( -( 2.7 + PF ) ) * 420.0 DF = AMIN1( DF, 1.9 ) ELSE DF = 0.1744 END IF C *** Snesitivity test *** C DF = 1.2*DF RETURN END SUBROUTINE TRANSP ( ET, NSOIL, ETP1, SMC, CMC, ZSOIL, & SHDFAC,SMCWLT,CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE TRANSPIRATION FROM THE VEGTYP FOR THIS PT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER K INTEGER NSOIL REAL CFACTR REAL CMC REAL CMCMAX REAL ET ( NSOIL ) REAL ETP1 REAL ETP1A REAL GX REAL PC REAL SHDFAC REAL SMC ( NSOIL ) REAL SMCREF REAL SMCWLT REAL ZSOIL ( NSOIL ) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 20 K = 1, NSOIL ET(K) = 0. 20 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC AN 'ADJUSTED' POTNTL TRANSPIRATION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ETP1A = SHDFAC * PC * ETP1 * (1.0 - (CMC /CMCMAX) ** CFACTR) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE PLANT TRANSP. FIRST, CALC GX AND ET FOR THE 1ST C (TOP) SOIL LAYER. GX IS THE SIMPLE TRANSPIRATION FUNCTION, A C RATIO OF THE ACTUAL SMC SURPLUS TO THE REFERENCE SMC SURPLUS. C USE THE ABSOLUTE DEPTH OF THE 1ST LAYER (ZSOIL(1)) TO GET C THE PROPORTION OF THE POTENTIAL IN THE FINAL CALCULATION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GX = ( SMC(1) - SMCWLT ) / ( SMCREF - SMCWLT ) C GX = MAX ( MIN ( GX, 1. ), 0. ) C Test canopy resistance C GX = 1.0 C ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE, C BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE C ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 10 K = 2, NROOT C GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT ) C GX = MAX ( MIN ( GX, 1. ), 0. ) C test canopy resistance C GX = 1.0 C ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A 10 CONTINUE RETURN END SUBROUTINE WDFCND ( WDF,WCND,SMC,SMCMAX,B,DKSAT,DWSAT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SOIL WATER DIFFUSIVITY AND SOIL CC ======= HYDRAULIC CONDUCTIVITY. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL B REAL DKSAT REAL DWSAT REAL EXPON REAL FACTR2 REAL SMC REAL SMCMAX REAL WCND REAL WDF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SMC = SMC SMCMAX = SMCMAX FACTR2 = SMC / SMCMAX CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EXPON = B + 2.0 WDF = DWSAT * FACTR2 ** EXPON CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EXPON = ( 2.0 * B ) + 3.0 WCND = DKSAT * FACTR2 ** EXPON RETURN END C SUBROUTINE REDPRM(VEGTYP,SOITYP, o CMCMAX,TOPT,RSMAX,ALBEDO,Z0,SHDFAC,NROOT,RCMIN,RGL,HS, o B,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF,SMCDRY,F1) C C This subroutine reads the soil and vegetation parameters C required for the execusion of the OSU land-surface scheme. C by F. Chen 3/15/95 C INTEGER NROOT,VEGTYP,SOITYP C Set-up soil parameters CALL PRMSOI(SOITYP,B,SMCDRY,F1,SMCMAX,SMCREF, & PSISAT,DKSAT,DWSAT,SMCWLT) C C Set-up vegetation parameters CALL PRMVEG(VEGTYP,CMCMAX,TOPT,RSMAX,ALBEDO,Z0, & SHDFAC,NROOT,RCMIN,RGL,HS) C RETURN END C SUBROUTINE PRMVEG(VEGTYP,CMCMAX,TOPT,RSMAX,ALBEDO,Z0, & SHDFAC,NROOT,RCMIN,RGL,HS) C C Set-up vegetation parameters for a given vegetaion type C C Input: C VEGTYP: Vegetation type C Ouput: C Vegetation parameters: C ALBEDO: SFC albedo C CMXTBL: MAX CNPY Capacity C Z0: Roughness length C SHDFAC: Plant shade factor C NROOT: Rooting depth C RCMIN: Mimimum stomatal resistance C RSMAX and RGL: Parameters used in radiation stress function C HS: Parameter used in vapor pressure deficit function C TOPT: Parameter used in temperature stress function C CMCMAX: Maximum canopy water capacity C C ************************************************************************* C C SSiB Vegetation Types (Dorman and Sellers, 1989; JAM) C C 1: Broadleaf-evergreen trees (tropical forest) C 2: Broadleaf-deciduous tress C 3: Broadleaf and needleleaf tress (mixed forest) C 4: Needleleaf-evergreen trees C 5: Needleleaf-deciduous tress (larch) C 6: Broadleaf tress with groundcover (savanna) C 7: Groundcover only (perennial) C 8: Broadleaf shrubs with perennial groundcover C 9: Broadleaf shrubs with bare soil C 10: Dwarf trees and shrubs with groundcover (tundra) C 11: Bare soil C 12: Cultivations (The same parameters for the Type 7) C 13: Glacial C C*************************************************************************** C INTEGER VEGTYP,NROTBL(13) REAL ALBTBL(13),Z0TBL(13),SHDTBL(13),RSMTBL(13), * RGLTBL(13),HSTBL(13) C DATA ALBTBL/0.11, 0.19, 0.16, 0.13, 0.19, 0.19, 0.19, * 0.29, 0.29, 0.14, 0.15, 0.19, 0.15/ DATA Z0TBL/2.653,0.826, 0.563, 1.089, 0.854, 0.856, 0.075, * 0.238,0.065, 0.076, 0.011, 0.075, 0.011/ DATA SHDTBL/0.90, 0.80, 0.80, 0.85, 0.40, 0.50, 0.80, * 0.20, 0.15, 0.30, 0.0, 0.80, 0.0/ DATA NROTBL/13*3/ DATA RSMTBL/150.0,100.0, 125.0, 150.0, 100.0, 70.0, 40.0, * 300.0,400.0, 150.0, 999.0, 40.0, 999.0/ DATA RGLTBL/30.0, 30.0, 30.0, 30.0, 30.0, 65.0, 100.0, * 100.0,100.0, 100.0, 999.0, 100.0, 999.0/ DATA HSTBL/41.69, 54.53, 51.93, 47.35, 47.35, 54.53, 36.35, * 42.00, 42.00, 42.00, 999.0, 36.35, 999.0/ SAVE C ALBEDO = ALBTBL(VEGTYP) C Z0 = Z0TBL(VEGTYP) C NROOT = NROTBL(VEGTYP) C SHDFAC = SHDTBL(VEGTYP) RCMIN = RSMTBL(VEGTYP) RGL = RGLTBL(VEGTYP) HS = HSTBL(VEGTYP) IF(VEGTYP.EQ.11) SHDFAC=0.0 C RETURN END C C SUBROUTINE PRMSOI(SOLTYP,SBB,SDRYSMC,SF11,SMAXSMC,SREFSMC, & SSATPSI,SSATDK,SSATDW,SWLTSMC) C C Set-up soil Parameters for given soil type C C Input: C SOLTYP: Soil type C Ouput: C C Soil parameters: C SMAXSMC: MAX soil moisture content C SREFSMC: Reference soil moisture C SWLTSMC: Wilting PT soil moisture contents C SDRYSMC: Air dry soil moist content limits C SSATPSI: SAT soil potential coefs. C SSATDK: SAT soil diffusivity/conductivity coefs. C SBB: Soil diffusivity/conductivity coef. C SSATDW: SAT soil diffusivity/conductivity coefs. C SF11: Soil diffusivity/conductivity coef. C C ************************************************************************* C C SOIL TYPES Zobler (1986) Cosby et al (1984) C 1 COARSE LOAMY SAND C 2 MEDIUM SILTY CLAY LOAM C 3 FINE LIGHT CLAY C 4 COARSE-MEDIUM SANDY LOAM C 5 COARSE-FINE SANDY CLAY C 6 MEDIUM-FINE CLAY LOAM C 7 COARSE-MED-FINE SANDY CLAY LOAM C 8 ORGANIC LOAM C 9 LAND ICE LOAMY SAND C C*************************************************************************** C INTEGER SOLTYP REAL BB(9),DRYSMC(9),F11(9),MAXSMC(9),REFSMC(9),SATPSI(9), & SATDK(9),SATDW(9),WLTSMC(9) C DATA MAXSMC/0.421, 0.464, 0.468, 0.434, 0.406, 0.465, 0.404, & 0.439, 0.421/ DATA DRYSMC/0.07, 0.14, 0.22, 0.08, 0.18, 0.16, 0.120, & 0.10, 0.07/ DATA SATPSI/0.04, 0.62, 0.47, 0.14, 0.10, 0.26, 0.14, & 0.36, 0.04/ DATA SATDK /1.41E-5, 0.20E-5, 0.10E-5, 0.52E-5, 0.72E-5, & 0.25E-5, 0.45E-5, 0.34E-5, 1.41E-5/ DATA BB /4.26, 8.72, 11.55, 4.74, 10.73, 8.17, 6.77, & 5.25, 4.26/ DATA REFSMC/0.283, 0.387, 0.412, 0.312, 0.338, 0.382, 0.315, & 0.329, 0.283/ DATA WLTSMC/0.029, 0.119, 0.139, 0.047, 0.020, 0.103, 0.069, & 0.066, 0.029/ DATA SATDW /5.71E-6, 2.33E-5, 1.16E-5, 7.95E-6, 1.90E-5, & 1.14E-5, 1.06E-5, 1.46E-5, 5.71E-6/ DATA F11 /-0.999, -1.116, -2.137, -0.572, -3.201, -1.302, & -1.519, -0.329, -0.999/ C NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC) C F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0 C REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 m/s= 0.5 mm/day C REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1) C WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB) (Wetzel and Chang, 1987) C WLTSMC=WLTSMC1-0.5*WLTSMC1 SAVE C SOIL PARAMETERS SBB = BB(SOLTYP) SDRYSMC = DRYSMC(SOLTYP) SF11 = F11(SOLTYP) SMAXSMC = MAXSMC(SOLTYP) SREFSMC = REFSMC(SOLTYP) SSATPSI = SATPSI(SOLTYP) SSATDK = SATDK(SOLTYP) SSATDW = SATDW(SOLTYP) SWLTSMC = WLTSMC(SOLTYP) C RETURN END SUBROUTINE DCOEF ( Z, Z0, T1V, TH2V, SFCSPD, CM, CH ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC NAME: DETERMINE COEFFICIENTS (DCOEF) VERSION: N/A CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL A REAL AH REAL AM REAL BETA REAL B1 REAL B2 REAL CH REAL CM REAL CTP REAL CUP REAL CUS REAL DEW REAL DRIP REAL EC REAL EDIR REAL ETT REAL EXMCH REAL FLX1 REAL FLX2 REAL FLX3 REAL G REAL PR REAL RHO REAL RIB REAL RUNOFF REAL SFCSPD REAL TH2V REAL T1V REAL VK REAL Z REAL Z0 REAL Z0H !!$omp threadlocal /rite/ COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RHO,RUNOFF, & DEW,RIB DATA B1 / 9.4 / DATA B2 / 15.0 / DATA CUS / 7.4 / DATA EXMCH / -1. / DATA G / 9.806 / C C DATA PR / .74 / C Set PR=1 the same as in Ek's version for PILPS 1994 C DATA PR / 1.0 / DATA VK / .4 / C SFCSPD=MAX(SFCSPD, 0.01) C Set Z0H as function of Z0 as in Beljaars and Betts (1992?) C Z0H = Z0/10. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC A FRICTION VELOCITY FOR USE IN CALCULATING THE DRAG C COEF FOR MOMENTUM AND ONE FOR USE IN CALCULATING THE DRAG C COEF FOR HEAT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC A = VK / ALOG( Z / Z0 ) AM = A * A AH = A * VK / ALOG( Z / Z0H ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC A BULK RICHARDSON NUMBER. CONSTRAIN ITS VALUE IN THE C STABLE CASE TO A MAXIMUM OF 1.0 TO AVOID CREATING EXCHANGE C COEFFICIENTS THAT APPROACH ZERO. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC RIB = G * Z * ( TH2V - T1V ) / ( TH2V * SFCSPD * SFCSPD ) RIB = MIN( RIB, 1.0 ) IF ( RIB .GE. 0. ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF THE RICHARDSON NUMBER IS .GE. ZERO, THE AIR IS STABLY C STRATIFIED. CALC THE DRAG COEFFICIENTS USING A METHOD C DEVELOPED BY MAHRT (MONTHLY WEATHER REVIEW, 1987). THE C SMALLEST ALLOWABLE CM VALUE IS 1.0E-6. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CM = AM * SFCSPD * EXP( EXMCH * RIB ) IF ( CM .LT. 1.0E-6) THEN CM = 1.0E-6 CH = 1.0E-6 ELSE CH = ( CM * AH / AM ) / PR ENDIF ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C OTHERWISE, THE AIR IS UNSTABLY STRATIFIED AND THE DRAG C COEFFICIENTS WILL BE CALCULATED AS FOLLOWS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CUP = CUS * AM * B1 * SQRT( -RIB * Z / Z0 ) CTP = CUS * AH * B1 * SQRT( -RIB * Z / Z0H ) CM = ( 1.0 - (B1 * RIB)/(1.0 + CUP) ) * AM * SFCSPD CH = ( 1.0 - (B2 * RIB)/(1.0 + CTP) ) * AH * SFCSPD / PR END IF RETURN END