PROGRAM MAIN IMPLICIT NONE CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NCEP/OH/OSU UNCOUPLED 1-D LAND-SURFACE MODEL: VERSION 1.1 07 MAR 99 CC CC THIS MAIN PROGRAM AND ITS FAMILY OF SUBROUTINES COMPRISE VERSION 1.1 CC OF THE "PUBLIC" RELEASE OF THE UNCOUPLED 1-D COLUMN VERSION OF THE CC NCEP LAND-SURFACE MODEL (LSM). THE NCEP LSM IS A DESCENDANT OF AN CC EARLY GENERATION OF THE OREGON STATE UNIVERSITY LSM, BUT INCLUDES CC SUBSTANTIAL PHYSICS EXTENSIONS AND RECODING ACCOMPLISHED ALONG THE CC WAY BY NCEP, OH, AFGWC, AND AFGL. CC CC FOR DOCUMENTATION OF THIS CODE AND INSTRUCTIONS ON ITS EXECUTION AND CC INPUT/OUTPUT FILES, SEE DOCUMENTATION FILE README_1.1 IN SAME PUBLIC CC SERVER DIRECTORY AS THIS CODE. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC PROGRAM HISTORY LOG CC CC VERSION 1.1 -- 07 MARCH 99 : PUBLIC RELEASE VERSION CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC COMING SOON: CC CC 1 - ARBITRARY NUMBER AND THICKNESS OF SOIL LAYERS CC CC 2 - IMPROVED SUBSURFACE HEAT FLUX UNDER SHALLOW SNOWPACK CC CC 3 - NAMELIST DIRECTED OPTIONAL INPUT OF NON-DEFAULT PARAMETER CC VALUES IN SUB REDPRM CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CHARACTER*72 CNTRFL, FILENAME INTEGER NSOLD PARAMETER (NSOLD=4) LOGICAL LSATURATED INTEGER ICE INTEGER IDAY INTEGER IIDAY INTEGER IMONTH INTEGER IREC1 INTEGER IREC3 INTEGER IREC5 INTEGER IRECD INTEGER NDAILY INTEGER NOUT1 INTEGER NOUT3 INTEGER NOUT5 INTEGER NREAD1 INTEGER NRECL INTEGER NROOT INTEGER NRUN INTEGER NSOIL INTEGER ISLOPTYP INTEGER ISOILTP INTEGER IVEGTYP INTEGER ITIME INTEGER jday INTEGER ITIME_ctl INTEGER jday_ctl INTEGER jday0 INTEGER IJ INTEGER INDI REAL R REAL CP PARAMETER (R = 287.04, CP = 1004.5) REAL BETA REAL DRIP REAL EC REAL EDIR REAL ETT REAL F REAL FLX1 REAL FLX2 REAL FLX3 REAL RUNOF REAL DEW REAL RIB REAL RUNOFF3 REAL SMSCDIF_O REAL SMSCDIF C LW REAL SIGMA REAL LATITUDE REAL LONGITUDE REAL RUNOFF2 REAL Q1 REAL AET REAL ALB REAL ALBEDO REAL CH REAL CM REAL CMC REAL DQSDT REAL CZIL REAL REFKDT REAL DQSDT2 REAL DT C LW REAL EMISS REAL ESAT REAL ETA REAL ETA_OBS REAL ETP REAL FUP REAL GHF REAL H REAL H_OBS REAL LE REAL LWDN REAL LW_in REAL PRCP REAL PTU REAL Par_in REAL Par_out REAL Q2 REAL Q2SAT REAL RES REAL RESTOT REAL RH REAL RNET REAL RUNOFF1 REAL Rg REAL SFCSPD REAL SFCPRS REAL SFCTMP REAL SHDFAC REAL ALBEDOM (13) REAL SHDFACM (13) REAL MLAI (13) REAL XLAI REAL SKN_IRT REAL SMC(NSOLD) REAL SMCMAX REAL SNMAX REAL SNOALB REAL SOILW REAL SOILM REAL SOLDN REAL STC(NSOLD) REAL S REAL S_OBS REAL T1 REAL T1_OBS REAL T14 REAL T1V REAL T2V C LW REAL TAK REAL TBOT REAL TH2 REAL TH2V REAL T_16cm REAL T_2cm REAL T_32cm REAL T_4cm REAL T_64cm REAL T_8cm REAL Z REAL Z0 REAL eddyuw REAL sm_20cm REAL sm_5cm REAL sm_60cm REAL u_bar REAL uprim2 REAL vprim2 REAL w_dir REAL wet REAL wprim2 C REAL PCPDAY REAL PCPSUM REAL ETADAY REAL ETSUM REAL RUNOFFSUM REAL RUNOFFDAY REAL SMCMM REAL SMCDIF REAL SMCNW REAL SMC_OBS(NSOLD) REAL ETADAY_O REAL ETSUM_O REAL SMCDIF_O REAL SMCMM_O REAL SMCNW_O C REAL SH2O(NSOLD) REAL SLDPTH(NSOLD) REAL SNOWH REAL SNEQV C************************************************************* CMIC$ TASKCOMMON RITE COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOF, & DEW,RIB,RUNOFF3 C------------------------------------------------------------- C THERE ARE 10 STEPS IN THIS MAIN PROGRAM DRIVER C C DRIVER STEP 1 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 1. READ CONTROL FILE <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C >>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<< CNTRFL = 'controlfile' CALL READCNTL (CNTRFL,NRUN,DT,NSOIL,NROOT,Z,SLDPTH, & ISOILTP,IVEGTYP,ISLOPTYP,ALBEDOM,SHDFACM,MLAI,SNOALB,ICE, & TBOT,Z0, CZIL,REFKDT,T1, & STC,SMC,SH2O,CMC,SNOWH,SNEQV,FILENAME, & LATITUDE, & LONGITUDE, & jday_ctl, & ITIME_ctl) C THIS IS FOR THE DAILY WATER BALANCE ------------------ C Initialize daily accum vars after reading control file PCPDAY = 0.0 PCPSUM = 0.0 ETADAY = 0.0 ETSUM = 0.0 ETADAY_O = 0.0 ETSUM_O = 0.0 RUNOFFSUM = 0.0 RUNOFFDAY = 0.0 SMCMM = 0.0 C ----- DO NOT INITIALIZE THIS INSIDE TIME LOOP ! ------------- C initialize variables for water balance (daily sm differences) C initial SMCMM from control file !!!!!!!!!!!!! C but only 3 layers !!!!!!!!!!!!! do ij = 1,3 SMCMM = SLDPTH(ij)*1000.*SMC(ij) + SMCMM end do C initial OBS SMCMM from control file C but only 3 layers !!!!!!!!!!!!! SMCMM_O = SMCMM SMSCDIF_O = 0.0 SMCNW_O = 0.0 SMSCDIF = 0.0 SMCNW = 0.0 C THIS PART ABOVE IS FOR THE DAILY WATER BALANCE ------- C ----- DO NOT INITIALIZE THE ABOVE INSIDE TIME LOOP ! -------- C DRIVER STEP 2 OPEN INPUT/OUTPUT FILES <<<<<<<<<<<<<<<<<<< C 2. OPEN INPUT/OUTPUT FILES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C >> OPEN FILE CONTAINING ATM. STATION DATA <<<<< NREAD1=25 OPEN(NREAD1,FILE=FILENAME,STATUS='OLD') C C ------------------------------------------------------------ C c (NOUT's must be positive for binary output - GrADS) c (NOUT's must be negative for ascii output intead) NOUT1 = -43 NOUT3 = -45 NOUT5 = -47 NDAILY= -49 C DATA WORD LENGTH (RECORD LENGTH, DEPENDS ON THE BINARY SIZE OF C THE REAL VARIABLES TO BE WRITEN FOR GrADS) NRECL=1 C HYDROLOGY IF (NOUT1 .GT. 0) THEN IREC1 = 1 OPEN(UNIT=NOUT1,FILE= & 'HYDRO.GRS', & STATUS='NEW' & ,ACCESS='DIRECT' & ,FORM='UNFORMATTED' & ,RECL=NRECL &) ELSE OPEN(UNIT=-NOUT1,FILE= & 'HYDRO.TXT', & STATUS='NEW') END IF C THERMODYNAMICS IF (NOUT3 .GT. 0) THEN IREC3 = 1 OPEN(UNIT=NOUT3,FILE= & 'THERMO.GRS', & STATUS='NEW' &,ACCESS='DIRECT' &,FORM='UNFORMATTED' &,RECL=NRECL &) ELSE OPEN(UNIT=-NOUT3,FILE= & 'THERMO.TXT', & STATUS='NEW') END IF C OBSERVATIONAL DATA IF (NOUT5 .GT. 0) THEN IREC5 = 1 OPEN(UNIT=NOUT5, & FILE='OBS_DATA.GRS', & STATUS='NEW' &, ACCESS='DIRECT' &, FORM='UNFORMATTED' &,RECL=NRECL &) ELSE OPEN(UNIT=-NOUT5, & FILE='OBS_DATA.TXT', & STATUS='NEW') END IF C DAILY ACCUMULATION VALUES IF (NDAILY .GT. 0) THEN IRECD = 1 OPEN(UNIT=NDAILY, & FILE='DAILY.GRS', & STATUS='NEW' &, ACCESS='DIRECT' &, FORM='UNFORMATTED' &,RECL=NRECL &) ELSE OPEN(UNIT=-NDAILY, & FILE='DAILY.TXT', & STATUS='NEW') END IF C C >>>>>>>> OPENING OUTPUT FILES FINISHED <<<<<<<<<<<<<<<<<<<<<<<<<<< C ------------------------------------------------------------ C Initialize Day number, summed sfc energy bal residual, C Phota Thermal Unit (PTU) C jday0 = 0 - 1 IIDAY = 0 RESTOT = 0.0 PTU = 0.0 C C READ FIRST RECORD OF FORCING TO CHECK THE START DATE <<<<<<<<<<<< C C >>>>>>>>> READPILPS-->READBND - Read Tilden Meyers BND data <<<<< CALL READBND(jday, ITIME, SFCTMP, RH, SFCPRS, Rg, * Par_in, Par_out, RNET, LW_in, GHF, PRCP, wet, SKN_IRT, * T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm, sm_5cm, * sm_20cm, sm_60cm, w_dir, u_bar, eddyuw, * uprim2, vprim2, wprim2, H, LE, * IMONTH,IDAY,NREAD1) rewind NREAD1 C CHECK DATE and time for simulation on first step IF ((ITIME_ctl .NE. ITIME) .or. (jday_ctl .NE. jday)) THEN PRINT*,'date/time specified in control file does not match & initial read from forcing data file, the program stops' PRINT*,' ' PRINT*,'Julian day in control file= ',jday_ctl PRINT*,'Julian day in forcing data file= ',jday PRINT*,' ' PRINT*,'Initial time in control file= ',ITIME_ctl PRINT*,'Initial time in forcing data file= ',ITIME stop 999 ENDIF C DRIVER STEP 3 (time step loop) <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C ---------------------------------------------------------------------| C ---------------------------------------------------------------------| C | C START TIME LOOP | C |C | DO INDI = 1,NRUN C DRIVER STEP 4 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 4. READ FORCING DATA <<<<<<<<<<<<<<<<<<<<<<<<< C C READ THE REQUIRED ATMOSPHERIC FORCING DATA, AS WELL AS NONREQUIRED C (BUT EXTREMELY USEFUL) COMPANION SIMULTANEOUS VALIDATING FLUX-SITE DATA C FROM THE EAST-CENTRAL ILLINOIS OBSERVING FLUX SITE (JUST SOUTH OF C CHAMPAIGN IL, NEAR BONDVILLE IL) OPERATED AND MAINTAINED (WITH SOME C GEWEX/GCIP SUPPORT) BY TILDEN MEYERS OF NOAA/ARL C C THE SEVEN REQUIRED SURFACE ATMOSPHERIC FORCING VARIABLES REQUIRED BY C THIS LSM AND MOST MODERN-ERA LSMs ARE: C 1) AIR TEMPERATURE, 2) AIR HUMIDITY, 3) WIND SPEED, 4) SFC PRESSURE, C 5) DOWNWARD SOLAR RAD, 6) DOWNWARD LONGWAVE RAD, AND MOST IMPORTANTLY, C 7) PRECIPITATION C C (NOTE: FOR THE UNITS REQUIRED OF THE FORCING DATA BY THIS LSM, C SEE COMMENT BLOCK AT BEGINNING OF ROUTINE "SFLX". ROUTINE C "READBND" HERE DOES THE REQUIRED UNITS CONVERSION FOR CALL SFLX) C C >>>>>>>>> READPILPS-->READBND - Read Tilden Meyers BND data <<<<< CALL READBND(jday, ITIME, SFCTMP, RH, SFCPRS, Rg, * Par_in, Par_out, RNET, LW_in, GHF, PRCP, wet, SKN_IRT, * T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm, sm_5cm, * sm_20cm, sm_60cm, w_dir, u_bar, eddyuw, * uprim2, vprim2, wprim2, H, LE, * IMONTH,IDAY,NREAD1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SAVE OBSERVATIONS FOR LATER VALIDATION OUTPUT C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Skin Temperature in Kelvin, (T1 IN SFLX) T1_OBS = SKN_IRT + 273.16 C C >>>>>>>> SIGN CONVENTIONS FOR READ FLUXES <<<<<<<<<<<<<<<<<< C C S_OBS: SOIL HEAT FLUX READ FROM DATA (S IN SFLX) C (W M-2: NEGATIVE, IF DOWNWARD FROM SURFACE) S_OBS = GHF C C H_OBS: SENSIBLE HEAT FLUX READ FROM DATA (H IN SFLX) C (W M-2: NEGATIVE, IF UPWARD FROM SURFACE) H_OBS = H C C ETA_OBS: LATENT HEAT FLUX READ FROM DATA (ETA IN SFLX) C (W M-2: NEGATIVE, IF UPWARD FROM SURFACE) ETA_OBS = LE C C SMC_OBS: Soil moisture content (volumetric) FROM DATA C (just for water balance purposes) - don't have 4th layer SMC_OBS(1) = sm_5cm SMC_OBS(2) = sm_20cm SMC_OBS(3) = sm_60cm C SMC_OBS(4) = 0.0 C C Keep the day number of the simulation IF (jday0 .NE. jday) THEN IIDAY = IIDAY + 1 jday0 = jday C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IN A FUTURE VERSION, CONSIDER IMPLEMENTING A PHYSICAL C TIME-DEPENDENT PHENOLOGY MODEL, THAT MIGHT CALCULATE A GROWING SEASON C GROWTH-STAGE INDEX, SUCH AS PTU = PHOTO THERMAL UNIT, DEPENDENT ON C ACCUMULATIVE GROWING SEASON AIR TEMPERATURE AND RADIATION. C FOR NOW, INCREMENT PTU IN A DUMMY WAY AS A STUB PLACE HOLDER C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PTU = PTU + 0.10 C C DRIVER STEP 5 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 5. INTERPOLATE DAILY LAND SURFACE CHARACTERISTICS FROM MONTHLY C VALUES: C C........................... monthly, daily ............ CALL MONTH_D(JDAY,ALBEDOM,ALB) CALL MONTH_D(JDAY,SHDFACM,SHDFAC) CALL MONTH_D(JDAY,MLAI ,XLAI) END IF C....................................................... C C >>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C *********************************************************** C C RADIATION C The following step (OBTENTION OF LWDN) has been C commented out, given that C the current forcing data file available provides the C downward long wave radiative forcing as measured by NOAA's C SURFRAD site located within 5 Km of the flux/meteorological C measurements system site. C C OBTENTION OF LWDN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C C COMPUTATION OF LWDN (INCOMING LW RADIATION) FROM TAIR AND Q: C C...............LWDN = EMISS*SIGMA*(TAK)**4. C C WHERE: TAK = AIR TEMP IN KELVIN C EMISS = 0.7 OR (IDSO AND JACKSON, 1969): C C EMISS = (1 - 0.261 EXP(-7.77*10^(-4)X(273-TAK)^2) C C NEED STEFAN-BOLTZMANN CONSTANT, SIGMA C SIGMA = 5.672 * 10^-8 W M^-2 T^-4 C C SIGMA = 5.672E-8 C TAK = SFCTMP C EMISS = 1 - 0.261*EXP((-7.77E-4)*(273-TAK)**2.) C C LWDN = EMISS*SIGMA*TAK**4. C *********************************************************** C DRIVER STEP 6 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 6. ASSIGN DOWNWARD SOLAR AND LONGWAVE RADIATION <<<<<<<<<<<<< C Using LW_in READ FROM DATA instead of LWDN calculated LWDN = LW_in SOLDN = Rg C DRIVER STEP 7 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 7. CALCULATE A SATURATION MIX RATIO (Q2SAT) <<<<<<<<<<<<<<<< C C....SET SATURATION FLAG (SEE COMMENT IN SFLX) LSATURATED = .FALSE. C NEED Q2 (FROM REL.HUMID.) USE SUBROUTINE QDATAP CALL QDATAP(SFCTMP,SFCPRS,RH,Q2,Q2SAT, ESAT) IF (Q2 .LT. 0.1E-5) Q2=0.1E-5 IF (Q2 .GE. Q2SAT) THEN Q2 = Q2SAT*0.99 LSATURATED = .TRUE. ENDIF C CALCULATE SLOPE OF SAT SPECIFIC HUMIDITY CURVE FOR PENMAN: DQSDT2 DQSDT2 = DQSDT (SFCTMP, SFCPRS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC TH2 = SFCTMP + ( 0.0098 * Z ) T2V = SFCTMP * (1.0 + 0.61 * Q2 ) T1V = T1 * (1.0 + 0.61 * Q2 ) TH2V = TH2 * (1.0 + 0.61 * Q2 ) C C DRIVER STEP 8 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 8. CALCULATE CH (EXCHANGE COEFFICIENT) <<<<<<<<<<<<<<<<<<<<< C SFCSPD = sqrt(u*u+v*v), if individual (u,v) components provided, C but Tilden Meyers data provides wind speed as u_bar. SFCSPD = u_bar CC!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C IMPORTANT NOTE: TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR HEAT AND C MOISTURE, SUBROUTINE SFCDIF BELOW CAN C C A) BE CALLED HERE FROM THE DRIVER TO SFLX (THUS CH IS INPUT TO SFLX) C (AS IN THE COUPLED ATMOSPHERIC MODEL) OR C C B) BE CALLED LATER IN ROUTINE SFLX (THUS CH IS OUTPUT FROM SFLX), C RIGHT AFTER THE SFLX CALL TO REDPRM AND BEFORE THE SFLX CALL PENMAN. C (AS IN AN UNCOUPLED, OFF-LINE LSM MODE) C C TO ENABLE THIS FLEXIBILITY, THE ADDITIONAL ARGUMENTS CZIL, Z0, SFCSPD C ARE PASSED TO SFLX, IN CASE SFCDIF IS CALLED IN SFLX C C IN THE UNCOUPLED, OFF-LINE LSM REPRESENTED HEREIN BY THIS DRIVER, C WE CALL SFCDIF LATER IN ROUTINE SFLX. THIS ENABLES EXTERNAL USERS C OF THE UNCOUPLED LSM TO MORE EASILY CONSTRUCT THEIR OWN MAIN DRIVER C C THE ROUTINE SFCDIF REPRESENTS THE SO-CALLED "SURFACE LAYER" OR THE C "CONSTANT FLUX LAYER" (THE LOWEST 20-100 M OF THE ATMOSPHERE). C HENCE ROUTINE SFCDIF EMBODIES THE "ATMOSPHERIC RESISTANCE" TO SFC FLUXES C CC!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CC CCC...CALL SFCDIF ( Z, Z0, T1V, TH2V, SFCSPD,CZIL, CM, CH ) C C C DRIVER STEP 9 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 9. CALL LAND-SURFACE PHYSICS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C CALL SFLX(ICE,LSATURATED,DT,Z,NSOIL,NROOT,SLDPTH, I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2, I IVEGTYP,ISOILTP,ISLOPTYP, I SHDFAC,XLAI,PTU,REFKDT,TBOT,ALB,SNOALB, I Z0,SFCSPD,CZIL, 2 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV, O CH,ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX,ALBEDO, O SOILW,SOILM) C AET = ETA C CALCULATE UPWARD LONGWAVE RAD USING UPDATED SKIN TEMPERATURE T14 = T1 * T1 * T1 * T1 FUP = 5.67E-8 * T14 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE RESIDUAL OF ALL SURFACE ENERGY BALANCE EQN TERMS. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC S = -S F = SOLDN*(1.0-ALBEDO) + LWDN C RES = F - H - S - AET - FUP - FLX1 - FLX2 - FLX3 RESTOT = RESTOT + RES CC.....S = S + RES IF ((INDI .LT. 200) .OR. (MOD(INDI,50) .EQ. 0)) THEN C Debug PRINT*,' CALCULATE RESIDUAL OF ALL SURFACE ENERGY' PRINT*,' --------------------------------------' C...PRINT*,'S=',S,' RES=',RES,' RESTOT=',RESTOT,' TIMESTEP=',INDI PRINT*,' RES=',RES,' TIMESTEP=',INDI PRINT*,' RES/(S+ETA) (in %)=',100*RES/(S+ETA),'%' PRINT*,' --------------------------------------' ENDIF C C WRITE HYDROLOGICAL VARIABLES IN GRD-FILES C C Debug PRINT*,' ---------------------------------------' C Debug PRINT*,' IIDAY=',IIDAY,' TIMESTEP=',INDI C Debug PRINT*,' ---------------------------------------' C ACCUMULATE DAILY VALUES FOR WATER BALANCE PCPSUM = PCPSUM + PRCP*DT ETSUM = ETSUM + ETA*DT/2.5E6 ETSUM_O = ETSUM_O + ETA_OBS*DT/2.5E6 RUNOFFSUM = RUNOFFSUM + RUNOFF1*DT*1000. + RUNOFF2*DT*1000. IF (ITIME .EQ. 0) THEN C Remember to initialize SMCMM, SMSCDIF, ... C after reading controlfile. SMCNW = 0.0 SMCNW_O =0.0 do ij = 1,3 SMCNW = SLDPTH(ij)*1000.*SMC(ij) + SMCNW SMCNW_O = SLDPTH(ij)*1000.*SMC_OBS(ij) + SMCNW_O end do SMCDIF = SMCNW - SMCMM SMCMM = SMCNW SMCDIF_O = SMCNW_O - SMCMM_O SMCMM_O = SMCNW_O RUNOFFDAY = RUNOFFSUM PCPDAY = PCPSUM ETADAY = -ETSUM ETADAY_O = -ETSUM_O PCPSUM = 0.0 ETSUM = 0.0 ETSUM_O = 0.0 RUNOFFSUM = 0.0 C C DRIVER STEP 10 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C 10. Write to model output files C EXTENSION .GRS FILES - GrADS-readable format output (unformatted) C EXTENSION .TXT FILES - ASCII formatted. C C WRITE DAILY ACUMMULATION VARIABLES, DAILY.GRS, GrADS FORMAT C or DAILY.TXT, ASCII formatted C CALL PRTDAILY(INDI,NDAILY,NSOIL,jday, & IIDAY,ITIME,ETADAY,ETADAY_O, & PCPDAY,SMCDIF,SMCDIF_O,RUNOFFDAY, & DT,IRECD) ENDIF C C WRITE VARIABLES IN WATER AMMOUNT UNITS (HYDRO.GRS or .TXT) C CALL PRTHYDF(INDI,NOUT1,NSOIL,jday, & IIDAY,ITIME,ETA, & ETP,PRCP,SH2O,SMC,ALBEDO,ALB,SNOALB, & SNEQV, SNMAX, SNOWH, SOILM, SOILW, & RUNOFF1, RUNOFF2, RUNOFF3, & DT,EDIR,EC,ETT,CMC,IREC1) C C WRITE VARIABLES IN THERMODYNAMIC ENERGY UNITS (THERMO.GRS or .TXT) C CALL PRTHMF(INDI,NOUT3,NSOIL,jday, & IIDAY,ITIME,CH,CM,Z,F,S,CMC, & SMCMAX,SFCTMP,T1,Q1,SFCPRS,SFCSPD,ETA,ETP,PRCP, & STC,DT,RNET,H,AET,RES,IREC3) C C WRITE ATNOSPHERIC INPUT FORCING DATA AND INDEPENDENT VALIDATION DATA C IN FILE (BND_DATA.GRS or .TXT) C CALL PRTBND(INDI,NOUT5,jday, & IIDAY, ITIME, SFCTMP, Q2, SFCPRS, Rg, * LW_in, Par_in, Par_out, RNET, S_OBS, PRCP, wet, T1_OBS, * T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm, sm_5cm, * sm_20cm, sm_60cm, w_dir, u_bar, eddyuw, * uprim2, vprim2, wprim2, H_OBS, ETA_OBS, * IMONTH, IDAY, IREC5) C C ============================================================= END DO C | C END OF TIME LOOP | C | C ---------------------------------------------------------------------| C ---------------------------------------------------------------------| C Debug PRINT*,' CALCULATE RESIDUAL OF ALL SURFACE ENERGY' PRINT*,' --------------------------------------' C...PRINT*,'S=',S,' RES=',RES,' RESTOT=',RESTOT,' TIMESTEP=',INDI PRINT*,'SUM OF ALL TIMESTEPS ENERGY RESIDUAL=',RESTOT, &' TIMESTEPS=',INDI PRINT*,' AVG RES PER TIMESTEP =',RESTOT/INDI PRINT*,' --------------------------------------' STOP 0 END C SUBROUTINE CANRES(SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, & SMCWLT,SMCREF,RCMIN,RC,PC,NROOT,Q2SAT,DQSDT2, & TOPT,RSMAX,RGL,HS,XLAI) IMPLICIT NONE 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 ********************************************************************** INTEGER NSOLD, NROOT, NSOIL, K PARAMETER (NSOLD = 4) REAL SIGMA, RD, CP, SLV REAL SOLAR, CH, SFCTMP, Q2, SFCPRS REAL SMC(NSOIL), ZSOIL(NSOIL), PART(NSOLD) REAL SMCWLT, SMCREF, RCMIN, RC, PC, Q2SAT, DQSDT2 REAL TOPT, RSMAX, RGL, HS, XLAI, RCS, RCT, RCQ, RCSOIL, FF REAL P, QS, GX, TAIR4, ST1, SLVCP, RR, DELTA PARAMETER (SIGMA=5.67E-8, RD=287.0, CP=1004.5, SLV=2.5E6) RCS = 0.0 RCT = 0.0 RCQ = 0.0 RCSOIL = 0.0 RC = 0.0 C ---------------------------------------------------------------------- C CONTRIBUTION DUE TO INCOMING SOLAR RADIATION. C ---------------------------------------------------------------------- CC/98/01/05/ FF = 0.55*2.0*SOLAR/RGL FF = 0.55*2.0*SOLAR/(RGL*XLAI) RCS = (FF + RCMIN/RSMAX) / (1.0 + FF) RCS = MAX(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 = MAX(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 = MAX(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. C**** USING SOIL DEPTH AS WEIGHTING FACTOR PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX C**** USING ROOT DISTRIBUTION AS WEIGHTING FACTOR CC PART(1) = RTDIS(1) * GX DO 10 K = 2, NROOT GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT) IF (GX .GT. 1.) GX = 1. IF (GX .LT. 0.) GX = 0. C**** USING SOIL DEPTH AS WEIGHTING FACTOR PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX C**** USING ROOT DISTRIBUTION AS WEIGHTING FACTOR CC PART(K) = RTDIS(K) * GX 10 CONTINUE DO 20 K = 1, NROOT RCSOIL = RCSOIL+PART(K) 20 CONTINUE RCSOIL = MAX(RCSOIL,0.0001) C ---------------------------------------------------------------------- C DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. C CONVERT CANOPY RESISTANCE (RC) TO PLANT COEFFICIENT (PC). C ---------------------------------------------------------------------- CC/98/01/05/ RC = RCMIN/(RCS*RCT*RCQ*RCSOIL) RC = RCMIN/(XLAI*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 FUNCTION CSNOW ( DSNOW ) IMPLICIT NONE REAL C REAL DSNOW REAL CSNOW REAL UNIT PARAMETER ( UNIT=0.11631 ) C **** SIMULATION OF TERMAL SNOW CONDUCTIVITY C **** SIMULATION UNITS OF CSNOW IS CAL/(CM*HR* C) C **** AND IT WILL BE RETURND IN W/(M* C) C **** BASIC VERSION IS DYACHKOVA EQUATION C ***** DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4 C=0.328*10**(2.25*DSNOW) CSNOW=UNIT*C C ***** DE VAUX EQUATION (1933), IN RANGE 0.1-0.6 C CSNOW=0.0293*(1.+100.*DSNOW**2) C ***** E. ANDERSEN FROM FLERCHINGER C CSNOW=0.021+2.51*DSNOW**2 RETURN END FUNCTION DEVAP ( ETP1, SMC, ZSOIL, SHDFAC, SMCMAX, B, & DKSAT, DWSAT, SMCDRY, SMCREF, SMCWLT) IMPLICIT NONE 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 C REAL WCND C REAL WDF REAL ZSOIL REAL SMCREF REAL SMCWLT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GET SOIL H2O DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY GIVEN SMC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FX CALL WDFCND ( WDF, WCND, SMC, SMCMAX, B, DKSAT, DWSAT ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FX REPRESENTS SOIL MOISTURE FLUX / ATMOSPHERIC DEMAND CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FX FX = ( -2. * WDF * ( SMC - SMCDRY ) / ZSOIL - WCND ) / ETP1 C*** TEST OF BETA METHODE, F. CHEN 7/11/96 *** C Change 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! formulation to reduce C excessive evaporation under little green vegetation: C FX = (SMC - SMCWLT) / (SMCREF - SMCWLT) C FX = (SMC - SMCWLT) / (SMCMAX - SMCWLT) FX = (SMC - SMCDRY) / (SMCMAX - SMCDRY) 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 FUNCTION DQS (T) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE VALUES OF VAPOR PRESSURE (E) CC AND P * DQS/DT (P TIMES CHG IN SAT MXG RATIO WITH RESPECT CC TO THE CHG IN TEMP) IN SUBSTITUTION TO THE LOOK-UP TABLES. CC CC SUBSTITUTES LOOK-UP TABLES ASSOCIATED WITH THE DATA CC BLOCK /CHMXR/ . CC CC FORMULAS AND CONSTANTS FROM ROGERS AND YAU, 1989. CC CC ADDED BY PABLO J. GRUNMANN, 6/30/97. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REAL DESDT REAL DQS CK REAL ESD REAL LW REAL T REAL ES C CK REAL CP CK REAL CV CK REAL CVV REAL CPV REAL RV REAL CW REAL EPS REAL ESO REAL TO C CK PARAMETER (CP = 1005.) CK PARAMETER (CV = 718.) CK PARAMETER (CVV = 1410.) PARAMETER (CPV = 1870.) PARAMETER (RV = 461.5) PARAMETER (CW = 4187.) PARAMETER (EPS = 0.622) PARAMETER (ESO = 611.2) PARAMETER (TO = 273.15) C C ABOUT THE PARAMETERS: C C EPS ---------- WATER - DRY AIR MOLECULAR MASS RATIO, EPSILON C C VALUES FOR SPECIFIC HEAT CAPACITY AND INDIVIDUAL GAS CONSTANTS C IN [JOULES/(KG*KELVIN)] UNITS. C C DRY AIR: C CP, CV C WATER VAPOR: C CVV = 1410. C CPV = 1870. C RV = 461.5 C LIQUID WATER: C CW = 4187. C C ESO = ES(T=273.15 K) = SAT. VAPOR PRESSURE (IN PASCAL) AT T=TO C TO = 273.15 C C SAT. MIXING RATIO: QS ~= EPS*ES/P C CLAUSIUS-CLAPEYRON: DES/DT = L*ES/(RV*T**2) C @QS/@T = (EPS/P)*DES/DT C LW = 2.501000E6 - ( CW - CPV ) * ( T - TO ) ES = ESO*EXP (LW*(1/TO - 1/T)/RV) DESDT = LW*ES/(RV*T*T) C C FOR INSERTION IN DQSDT FUNCTION: C DQSDT = DQS/P , WHERE DQS = EPS*DESDT C DQS = EPS*DESDT C RETURN END FUNCTION DQSDT ( SFCTMP, SFCPRS ) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO RETRIEVE THE APPROPRIATE VALUE OF DQSDT (THE CHANGE CC ======= OF THE SATURATION MIXING RATIO WITH RESPECT TO THE CC CHANGE IN TEMPERATURE) FROM: CC CC FORMULAS INTRODUCED IN NEW FUNCTION DQS CC (MODIFIED BY PABLO GRUNMANN, 7/9/97). CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL SFCTMP REAL SFCPRS REAL DQS REAL DQSDT IF ((SFCTMP .GE. 173.0) .AND. (SFCTMP .LE. 373.0)) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF THE INPUT SFC AIR TEMP IS BTWN 173 K AND 373 K, USE C FUNCTION DQS TO DETERMINE THE SLOPE OF SAT.MIX RATIO FUNCTION C -ADAPTED TO USE NEW DQS, 7/9/97. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DQSDT = DQS (SFCTMP) / SFCPRS ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C OTHERWISE, SET DQSDT EQUAL TO ZERO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DQSDT = 0.0 END IF RETURN END FUNCTION E (T) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE VALUES OF SAT. VAPOR PRESSURE (E) CC SUBSTITUTES LOOK-UP TABLES ASSOCIATED WITH THE DATA CC BLOCK /VAPPRS/ . CC FORMULAS AND CONSTANTS FROM ROGERS AND YAU, 1989. CC CC ADDED BY PABLO J. GRUNMANN, 7/9/97. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REAL LW REAL T REAL E C CK REAL EPS CK REAL CP CK REAL CV CK REAL CVV REAL CPV REAL RV REAL CW REAL ESO REAL TO C CK PARAMETER (EPS = 0.622) CK PARAMETER (CP = 1005.) CK PARAMETER (CV = 718.) CK PARAMETER (CVV = 1410.) PARAMETER (CPV = 1870.) PARAMETER (RV = 461.5) PARAMETER (CW = 4187.) PARAMETER (ESO = 611.2) PARAMETER (TO = 273.15) C C ABOUT THE PARAMETERS: C C EPS ---------- WATER - DRY AIR MOLECULAR MASS RATIO, EPSILON C C VALUES FOR SPECIFIC HEAT CAPACITY AND INDIVIDUAL GAS CONSTANTS C IN [JOULES/(KG*KELVIN)] UNITS. C C DRY AIR: C CP, CV C WATER VAPOR: C CVV = 1410. C CPV = 1870. C RV = 461.5 C LIQUID WATER: C CW = 4187. C C ESO = ES(TO) = SAT. VAPOR PRESSURE (IN PASCAL) AT T=TO C TO = 273.15 C_______________________________________________________________________ C C CLAUSIUS-CLAPEYRON: DES/DT = L*ES/(RV*T**2) C C LW = 2.501000D6 - ( CW - CPV ) * ( T - TO ) E = ESO*EXP (LW*(1/TO - 1/T)/RV) C RETURN END FUNCTION FRH2O(T,SMC,SH2O,SMCMAX,B,PSIS) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: CALCULATE FREE WATER CONTENT IF TEMPERATURE CC ======= IS BELOW 273.16K (T0). CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FRH2O, HLICE, GS, CK, DICE, DH2O, T0, ERROR, T REAL SMC, SH2O, SMCMAX, B, PSIS, BX, TR, SWL REAL SMCM, SMCC, DF, DENOM, SWLK, DSWL INTEGER N PARAMETER (HLICE=3.335 E5) PARAMETER (GS = 9.81) PARAMETER (CK=8.) PARAMETER (DICE=920.0) PARAMETER (DH2O=1000.0) PARAMETER (T0=273.16) PARAMETER (ERROR=0.005) C *** LIMITS ON PARAMETER B: B < 5.5 **** C *** SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT **** C *** IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES **** C****************************************************************** C BX = B IF ( B .GT. 5.5 ) BX = 5.5 C------------------------------------------------------------------ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC ASSUMED THAT FREE WATER CONTENT EQUALS ITS VALUE AT TEMPERATURE CC 273.11K IN RANGE OF TEMPERATURES 273.11K - 273.16K TO REDUCE CC ITERRATINS PROBLEM IN THIS RANGE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC TR=T IF(T .GT. 273.11 .AND. T .LE. T0 ) TR=273.11 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF TEMPERATURE ABOVE 273.16K SH2O = SMC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(TR .GT. T0) THEN FRH2O=SMC GOTO 77 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT SOIL MOISTURE VARIABLES INTO % CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SWL=(SMC-SH2O)*100. SMCM=SMCMAX*100. SMCC=SMC*100. N=0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C START OF ITERRATIONS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 10 DF = ( -PSIS*GS/HLICE ) * ( ( 1.+CK*SWL*0.01 )**2. ) * + ( SMCM/(SMCC-SWL ) )**BX DENOM = 2. * CK / ( 1.+CK*SWL*0.01 ) + BX / ( SMCC - SWL ) SWLK = SWL - (1.-(TR-T0)/(TR*DF))/DENOM DSWL=ABS(SWLK-SWL) SWL=SWLK IF(SWL.GE.SMCC) THEN C WRITE(*,*) 'WARN2 ',N,SWL,SMCC SWL=SMCC*0.8 SMCC=SMCC*0.97 ENDIF N=N+1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC EXIT IF MORE THAN 900 ITERRATIONS. FREE WATER CONTENT CC WILL NOT BE CHANGED CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(N .GT. 900 ) THEN IF(TR . LT. 271.16) THEN SWL = ( SMC-0.02 ) * 100. ELSE SWL=0. ENDIF GOTO 20 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC END OF ITERRATIONS IF LAST ESTIMATE DIFFERENCE IS LESS THAN CC AN ERROR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( DSWL .GT. ERROR ) GOTO 10 20 IF(SWL .LT. 0.) THEN C WRITE(*,*) ' SWL < 0. ',SWL SWL=0. ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC RECALCULATE ICE CONTENT IN % INTO VOLUMETRIC FREE WATER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c this eq was using MAX1(), which yields an integer result c changed to MAX() - Pablo Grunmann 12-07-98. - !!!!!!!!!!!! FRH2O = MAX ( (SMC - SWL*0.01 ) , 0.02 ) C *** FLERCHINGER EQUATION ******* CC FK=(((-HLICE/(GS*PSIS))*((TR-T0)/TR))**(-1/B))*SMCMAX CC FRH2O = MIN ( FK, SMC ) CC IF(FRH2O .LT. 0.) FRH2O = 0.02 C----------------------------------------------------- 77 CONTINUE RETURN END SUBROUTINE HRT ( RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, + SMCWLT, PSISAT, SH2O, DT, B, F1, DF1, QUARTZ ) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOLD PARAMETER ( NSOLD = 4 ) INTEGER NSOIL INTEGER K 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 c REAL MXSMC c REAL MXICE REAL QUARTZ REAL RHSTS ( NSOIL ) REAL S REAL SMC ( NSOIL ) REAL SH2O ( NSOIL ) REAL SMCMAX REAL SMCWLT REAL STC ( NSOIL ) REAL TBOT REAL YY REAL ZBOT REAL ZSOIL ( NSOIL ) REAL ZZ1 REAL CICE, T0, TSURF, PSISAT, DT, B, SICE, TBK, TSNSR, TBK1 REAL TBND, SNKSRC INTEGER I CMIC$ TASKCOMMON ABCI COMMON /ABCI/ AI, BI, CI PARAMETER(CAIR=1004.0,CH2O=4.2 E6,CSOIL=1.26E6,ZBOT=-3.0) C THIS COMMON STATEMENT IS NEED ONLY TO PRINT ESTIMATED SOIL C SURFACE TEMPERATURE, AND CAN BE OMITTED COMMON /PRN_TSOIL/ TSURF C ---------- FROZEN GROUND VERSION ---------------------------- C CICE IS ICE HEAT CAPACITY C PARAMETER ( CICE = 2.106 E6, T0 = 273.16 ) C ------------------------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE THERMAL DIFFUSIVITY FOR THE WETTER OF THE TOP SOIL C LAYER OR THE NEXT LOWER SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ---------- FROZEN GROUND VERSION ---------------------------- C SOIL THERMAL CONDUCTIVITY ADJUSTMENT IF IT IS FROZEN C SICE = 0. DO I = 1, NSOIL SICE = SICE + SMC(I) - SH2O(I) END DO C C C C C IF( SMC(1) .GE. SMC(2)) THEN C MXSMC = SMC(1) C MXICE = SMC(1) - SH2O(1) C ELSE C MXSMC = SMC(2) C MXICE = SMC(2) - SH2O(2) c ENDIF CALL TDFCND ( DFMAX1, SMC(1), B, F1,QUARTZ,SMCMAX,SH2O(1)) C__________________________________________________________________ C This IF-statement is no longer used if Peters-Lidard thermal C conductivity is used for all conditions (frozen,unfrozen): C IF ( MXICE .GT. 0.0 ) CALL FRZCND ( DFMAX1, MXICE ) C______________________________________________P.Grunmann_01/99____ C ------------------------------------------------------------ C MXSMC = MAX ( SMC(1), SMC(2) ) C CALL TDFCND ( DFMAX1, MXSMC, B, F1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC HCPCT = SH2O(1)*CH2O + (1.0-SMCMAX)*CSOIL + (SMCMAX-SMC(1))*CAIR C VK correction: SMC to SH2O, first term on R.H. side.-P.Grunmann 1/99 C ---------- FROZEN GROUND VERSION ---------------------------- C ICE HEAT CAPACUTY IS INCLUDED C + + ( SMC(1) - SH2O(1) )*CICE C ------------------------------------------------------------------- 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 C ------------ FROZEN GROUND VERSION -------------------------- C SOIL SURFACE TEMPERATURE ADJUSTMENT BY ZZ1 ( THE SAME AS T1 ) C TSURF = ( YY + ( ZZ1 - 1 ) * STC(1) ) / ZZ1 C ------------------------------------------------------------------- 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 ) C ----------- FROZEN GROUND VERSION -------------------------- C CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR. C AND SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR. C IF ( SICE .NE. 0. .OR. TSURF .LT. T0 ) THEN TBK = TBND ( STC(1), STC(2), ZSOIL, ZBOT, 1, NSOIL ) TSNSR = SNKSRC ( TSURF, STC(1),TBK, SMC(1), SH2O(1), HCPCT, + ZSOIL,SMCMAX, PSISAT, B, F1, DT, 1, QUARTZ) RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT ) ENDIF C ------------------------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE DDZ2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DDZ2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO K = 2, NSOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC THIS SOIL LAYER'S HEAT CAPACITY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC HCPCT = SH2O(K)*CH2O +(1.0-SMCMAX)*CSOIL +(SMCMAX-SMC(K))*CAIR C VK correction: SMC to SH2O, first term on R.H. side.-P.Grunmann 1/99 C ---------- FROZEN GROUND VERSION ---------------------------- C ICE HEAT CAPACUTY IS INCLUDED C + + ( SMC(K) - SH2O(K) )*CICE C ------------------------------------------------------------------- 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 C MXSMC = MAX ( SMC(K), SMC(K+1) ) C ------------ FROZEN GROUND VERSION -------------------------- C SOIL THERMAL DIFFUSIVITY ADJUSTMENT IF SOIL IS FROZEN C C IF ( SMC(K) .GE. SMC(K+1) ) THEN C MXSMC = SMC(K) C MXICE = SMC(K) - SH2O(K) C ELSE C MXSMC = SMC(K+1) C MXICE = SMC(K+1) - SH2O(K+1) C ENDIF CC CALL TDFCND ( DFMAXN, MXSMC, B, F1,QUARTZ,SMCMAX,SH2O(K)) CALL TDFCND ( DFMAXN, SMC(K), B, F1,QUARTZ,SMCMAX,SH2O(K)) CC__________________________________________________________________ C This IF-statement is no longer used if Peters-Lidard thermal C conductivity is used for all conditions (frozen,unfrozen): C IF ( MXICE .GT. 0.0 ) CALL FRZCND ( DFMAXN, MXICE ) C ------------------------------------------------------------------- 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) C ----------- FROZEN GROUND VERSION -------------------------- C CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR. C TBK1 = TBND ( STC(K), STC(K+1), ZSOIL, ZBOT, K, NSOIL ) C ------------------------------------------------------------------- ELSE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE THERMAL DIFFUSIVITY FOR THE LOWEST SOIL LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ----------- FROZEN GROUND VERSION -------------------------- C SOIL THERMAL DIFFUSIVITY ADJUSTMENT C AND CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR. C C MXICE = SMC(K) - SH2O(K) CALL TDFCND ( DFMAXN, SMC(K), B, F1,QUARTZ,SMCMAX,SH2O(K)) C__________________________________________________________________ C This IF-statement is no longer used if Peters-Lidard thermal C conductivity is used for all conditions (frozen,unfrozen): C IF ( MXICE .GT. 0.0 ) CALL FRZCND ( DFMAXN, MXICE ) TBK1 = TBND ( STC(K), TBOT, ZSOIL, ZBOT, K, NSOIL ) C ------------------------------------------------------------------- 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 C ----------- FROZEN GROUND VERSION -------------------------- C SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR. C IF ( SICE .NE. 0. .OR. TBK .LT. 0. ) THEN TSNSR = SNKSRC ( TBK, STC(K),TBK1, SMC(K), SH2O(K), HCPCT, + ZSOIL,SMCMAX, PSISAT, B, F1, DT, K, QUARTZ) RHSTS(K) = RHSTS(K) - TSNSR / DENOM ENDIF C ------------------------------------------------------------------- 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 C ----------- FROZEN GROUND VERSION -------------------------- C REMEMBER UPPER LAYER TEMPERATURE C TBK = TBK1 C ------------------------------------------------------------------- DFMAX1 = DFMAXN DTSDZ = DTSDZ2 DDZ = DDZ2 C 40 CONTINUE END DO RETURN END SUBROUTINE HRTICE (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, & TBOT, B, F1, DF1 ) IMPLICIT NONE 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 INTEGER NSOLD 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 B REAL F1 REAL HCPCT INTEGER K CC 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 CMIC$ TASKCOMMON ABCI COMMON /ABCI/ AI, BI, CI PARAMETER(CAIR=1004.0,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 ( STCOUT, STCIN, RHSTS, DT, NSOIL ) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOLD PARAMETER ( NSOLD=4 ) INTEGER K INTEGER NSOIL REAL AI ( NSOLD ) REAL BI ( NSOLD ) REAL CI ( NSOLD ) REAL DT REAL RHSTS ( NSOIL ) C ----------- FROZEN GROUND VERSION ------------------------- C REAL STCOUT ( NSOIL ) REAL STCIN ( NSOIL ) C ------------------------------------------------------------------ CMIC$ TASKCOMMON 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 C ----------- FROZEN GROUND MODIFICATION ---------------------- C NEW TEMPERATURES IN STCOUT FROM OLD STCIN C STCOUT(K) = STCIN(K) + CI(K) C ------------------------------------------------------------------- C STC(K) = STC(K) + CI(K) 20 CONTINUE RETURN END SUBROUTINE JULDATE(JULD,IMONTH,IDAY,JULM) IMPLICIT NONE CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NAME: JULIAN DAYS TO MM/DD DATE CC CC CC PURPOSE: TO CONVERT FROM JULIAN DAY NUMBER OF THE YEAR TO CC CALENDAR DATE (MONTH, DAY). CC PABLO J. GRUNMANN, 05/98 CC CC VARIABLES: CC ========= CC CC LABEL ............DESCRIPTION............... CC CC IDMONTH(IM) NUMBER OF DAYS IN MONTH NUMBER "IM", FOR IM=1,12 CC JULD JULIAN DAY CC ILDM,ILDMP JULIAN DAY OF THE LAST DAY OF A PARTICULAR CC MONTH AND ITS CONSECUTIVE (ILDMP). CC IMONTH,IDAY MONTH, DAY - CALENDAR DATE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INTEGER IM INTEGER IDMONTH(12) INTEGER JULM(13) INTEGER JULD INTEGER IDAY INTEGER ILDM INTEGER ILDMP INTEGER IMONTH C................................................................. C DEFINE LAST DAY OF FEBRUARY AND MODIFY DATA IDMONTH C IF NECESSARY: C FEB 28 C FEB 29 (LEAP YEAR) C..................JAN,FEB,MAR,APR,MAY,JUN,JUL,AUG.SEP,OCT,NOV,DEC DATA IDMONTH/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ C DATA JULM/0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, C & 334, 365/ C................................................................. ILDMP = 0 IMONTH = 0 C Print*,JULM JULM(1)= 0 DO IM = 1,12 JULM(IM+1) = JULM(IM) + IDMONTH(IM) ILDM = JULM(IM) ILDMP = JULM(IM+1) IF ((JULD .LE. ILDMP) .AND. (JULD .GT. ILDM)) THEN IMONTH = IM IDAY = JULD - JULM(IM) ENDIF END DO RETURN END SUBROUTINE MONTH_D(JDAY,XMON,XDAY) IMPLICIT NONE C INTERPOLATE TO DAY OF YEAR. C THE MONTHLY FIELDS ASSUMED VALID AT 15TH OF MONTH. C C READ FROM A FILE THAT HAS 13 RECORDS, ONE PER MONTH WITH C JANUARY OCCURRING BOTH AS RECORD 1 AND AGAIN AS RECORD 13, C THE LATTER TO SIMPLIFY TIME INTERPOLATION FOR DAYS C BETWEEN DEC 16 AND JAN 15. WE TREAT JAN 1 TO JAN 15 C AS JULIAN DAYS 366 TO 380 BELOW, I.E WRAP AROUND YEAR. C C BASED ON SUBROUTINE CNSTS.f WHICH C INCLUDED REVISION BY F. CHEN 7/96 TO REFLECT A NEW NESDIS C VEGETATION FRACTION PRODUCT (FIVE-YEAR CLIMATOLOGY WITH C 0.144 DEGREE RESOLUTION FROM 89.928S, 180W TO 89.928N, 180E) C C This is PABLO J. GRUNMANN's version for use with the uncoupled, C column mode Land-Surface Model C C DATA JULM/0,31,59,90,120,151,181,212,243,273,304,334,365/ C.......getting JULM from subrooutine JULDATE to avoid C.......problems with leap years. C C **** DO TIME INTERPOLATION **** C INTEGER IDAY INTEGER IMON1 INTEGER IMON2 INTEGER IMONTH INTEGER JDAY INTEGER JULD INTEGER JULM(13) REAL DAY1 REAL DAY2 REAL RDAY REAL WGHT1 REAL WGHT2 REAL XDAY REAL XMON(13) JULD = JDAY CALL JULDATE(JULD,IMONTH,IDAY,JULM) IF(JULD.LE.15) JULD=JULD+365 IMON2=IMONTH IF (IDAY .GT. 15) IMON2 = IMON2 + 1 IF (IMON2 .EQ. 1) IMON2 = 13 IMON1 = IMON2 - 1 C **** ASSUME DATA VALID AT 15TH OF MONTH DAY2 = REAL(JULM(IMON2)+15) DAY1 = REAL(JULM(IMON1)+15) RDAY = REAL(JULD) WGHT1 = (DAY2-RDAY)/(DAY2-DAY1) WGHT2 = (RDAY-DAY1)/(DAY2-DAY1) C XDAY=WGHT1*XMON(IMON1)+WGHT2*XMON(IMON2) 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, & SH2O, SLOPE, KDT, FRZFACT, PSISAT, & ZSOIL, DKSAT, DWSAT, TBOT, RUNOFF1, RUNOFF2, & RUNOFF3, EDIR1, EC1, ETT1, NROOT, ICE,RTDIS, & QUARTZ) IMPLICIT NONE C ----------- FROZEN GROUND VERSION --------------------------- C NEW STATES ADDED: SH2O, AND FROZEN GROUND CORRECTION FACTOR 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 KDT REAL PC REAL PRCP REAL PRCP1 REAL Q2 REAL RCH REAL RIB REAL RR REAL RTDIS (NSOIL) REAL RUNOFF,RUNOXX3 REAL S REAL SFCTMP REAL SHDFAC REAL SIGMA REAL SMC ( NSOIL ) REAL SH2O ( 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 REAL Q1, SLOPE, FRZFACT, PSISAT, RUNOFF1, RUNOFF2, RUNOFF3 REAL EDIR1, EC1, ETT1, QUARTZ, SICE1 INTEGER NROOT, ICE CMIC$ TASKCOMMON RITE COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOFF, & DEW,RIB,RUNOXX3 PARAMETER(CP=1004.5, 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 C print*,'NOPAC calling smflx' C print*,'smc=' , smc C print*,'sh2o=', sh2o IF ( ETP .GT. 0.0 ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT PRCP FROM KG M-2 S-1 TO M S-1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ------------ FROZEN GROUND VERSION -------------------------- C (+) NEW STATES ADDED: SH2O, AND FROZEN GROUND CORRECTION FACTOR C CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, + SH2O, SLOPE, KDT, FRZFACT, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC, & CMCMAX,SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, & EDIR1, EC1, ETT1, SFCTMP,Q2,NROOT,RTDIS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ETA = ETA1 * 1000.0 C print*,'NOPAC, etp>0, AFTER smflx' C print*,'smc=' , smc C print*,'sh2o=', sh2o 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 C ------------ FROZEN GROUND VERSION -------------------------- C (+) NEW STATES ADDED: SH2O, AND FROZEN GROUND CORRECTION FACTOR C CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, + SH2O, SLOPE, KDT, FRZFACT, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC, & CMCMAX,SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, & EDIR1, EC1, ETT1, SFCTMP, Q2, NROOT,RTDIS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ETA = ETA1 * 1000.0 C print*,'NOPAC, etp<0, AFTER smflx' CC print*,'smc=' , smc C print*,'sh2o=', sh2o 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 CCC....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,QUARTZ,SMCMAX,SH2O ) C -------------- FROZEN GROUND VERSION ---------------------- C THERMAL CONDUCTIVITY ADJUSTMENT C SICE1 = SMC(1) - SH2O(1) C__________________________________________________________________ C This IF-statement is no longer used if Peters-Lidard thermal C conductivity is used for all conditions (frozen,unfrozen): C IF ( SICE1 .GT. 0.0 ) CALL FRZCND ( DF1, SICE1 ) YYNUM = F - SIGMA * T24 YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0 C ------------ FROZEN GROUND VERSION -------------------------- C (+) NEW STATES ADDED: SH2O & PARAMETER SMCWLT & PSISAT C CALL SHFLX ( S,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,TBOT, + SMCWLT, PSISAT, SH2O, & B,F1,DF1, ICE, & QUARTZ ) 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) IMPLICIT NONE 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,RUNOXX3 REAL S REAL SFCPRS REAL SFCTMP REAL SIGMA REAL T24 REAL T2V REAL TH2 REAL DQSDT2 CMIC$ TASKCOMMON RITE COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOFF, & DEW,RIB,RUNOXX3 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... CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FLX2 = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DELTA = ELCP * DQSDT2 T24 = SFCTMP * SFCTMP * SFCTMP * 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C !!!!!!!!!!!! TEST DRAINAGE -> ETP = 0.0 !!!!!!!!!!!!!!!! C ETP = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC RETURN END SUBROUTINE PRMSOI(SOLTYP,SBB,SDRYSMC,SF11,SMAXSMC,SREFSMC, + FRZFACT,SSATPSI,SSATDK,SSATDW,SWLTSMC,QTZ) IMPLICIT NONE C C SET-UP SOIL PARAMETERS FOR GIVEN SOIL TYPE C C INPUT: C SOLTYP: SOIL TYPE (INTEGER INDEX) C OUPUT: C C SOIL PARAMETERS: C SMAXSMC: MAX SOIL MOISTURE CONTENT (POROSITY) C SREFSMC: REFERENCE SOIL MOISTURE (ONSET OF SOIL MOISTURE STRESS IN TRANSPIRATION) C SWLTSMC: WILTING PT SOIL MOISTURE CONTENTS C SDRYSMC: AIR DRY SOIL MOIST CONTENT LIMITS C SSATPSI: SATURATED SOIL POTENTIAL C SSATDK: SATURATED SOIL HYDRAULIC CONDUCTIVITY C SBB: THE 'B' PARAMETER C SSATDW: SATURATED SOIL DIFFUSIVITY C SF11: USED TO COMPUTE SOIL DIFFUSIVITY/CONDUCTIVITY C QUARTZ: SOIL QUARTZ CONTENT C ************************************************************************* C C SOIL TYPES ZOBLER (1986) COSBY ET AL (1984) (quartz cont.(1)) C 1 COARSE LOAMY SAND (0.82) C 2 MEDIUM SILTY CLAY LOAM (0.10) C 3 FINE LIGHT CLAY (0.25) C 4 COARSE-MEDIUM SANDY LOAM (0.60) C 5 COARSE-FINE SANDY CLAY (0.52) C 6 MEDIUM-FINE CLAY LOAM (0.35) C 7 COARSE-MED-FINE SANDY CLAY LOAM (0.60) C 8 ORGANIC LOAM (0.40) C 9 LAND ICE LOAMY SAND (NA using 0.82) C C*************************************************************************** C INTEGER SOLTYP REAL SBB, SDRYSMC, SF11, SMAXSMC, SREFSMC, FRZFACT, SSATPSI REAL SSATDK, SSATDW, SWLTSMC, QTZ REAL BB(9),DRYSMC(9),F11(9),MAXSMC(9),REFSMC(9),SATPSI(9), & SATDK(9),SATDW(9),WLTSMC(9), QUARTZ(9) DATA MAXSMC/0.421, 0.464, 0.468, 0.434, 0.406, 0.465, 0.404, & 0.439, 0.421/ 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/ c DATA DRYSMC/0.029, 0.119, 0.139, 0.047, 0.020, 0.103, 0.069, c & 0.066, 0.029/ C DATA DRYSMC/0.07, 0.14, 0.22, 0.08, 0.18, 0.16, 0.12, C & 0.10, 0.07/ DATA DRYSMC/0.07, 0.24, 0.22, 0.08, 0.18, 0.16, 0.12, & 0.10, 0.07/ 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/ DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52, & 0.35, 0.60, 0.40, 0.82/ 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) QTZ = QUARTZ(SOLTYP) C ------------- FROZEN GROUND VERSION --------------------- C ADJUSTMENT FACTOR OF FROZEN GROUND PARAMETER USING SOIL TYPE 3 C AS REFERENCE TYPE C FRZFACT = ( SMAXSMC / SREFSMC ) * ( REFSMC(3) / MAXSMC(3) ) C ------------------------------------------------------------------ C RETURN END SUBROUTINE PRMVEG(VEGTYP,SHDFAC,RCMIN,RGL,HS,SNUP, O RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL) IMPLICIT NONE C NEW PARAMETER: SNOW COVER THRESHOLD, SNUP C C SET-UP VEGETATION PARAMETERS FOR A GIVEN VEGETAION TYPE C C INPUT: C VEGTYP: VEGETATION TYPE (INTEGER INDEX) C OUPUT: C VEGETATION PARAMETERS: C SHDFAC: VEGETATION GREENNESS FRACTION C RCMIN: MIMIMUM STOMATAL RESISTANCE C RGL: PARAMETER USED IN SOLAR RAD TERM OF CANOPY RESISTANCE FUNCTION C HS: PARAMETER USED IN VAPOR PRESSURE DEFICIT TERM OF CANOPY RESISTANCE FUNCTION C SNUP: THRESHOLD SNOW DEPTH (IN WATER EQUIVALENT M) THAT IMPLIES 100% SNOW COVER C C ********************************************************************** C C SSIB VEGETATION TYPES (DORMAN AND SELLERS, 1989; JAM) C C 1: BROADLEAF-EVERGREEN TREES (TROPICAL FOREST) C 2: BROADLEAF-DECIDUOUS TREES C 3: BROADLEAF AND NEEDLELEAF TREES (MIXED FOREST) C 4: NEEDLELEAF-EVERGREEN TREES C 5: NEEDLELEAF-DECIDUOUS TREES (LARCH) C 6: BROADLEAF TREES 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 AS FOR TYPE 7) C 13: GLACIAL C C*********************************************************************** INTEGER VEGTYP INTEGER NSOIL INTEGER NROOT INTEGER I REAL RTDIS(NSOIL) REAL SLDPTH(NSOIL) REAL ZSOIL(NSOIL) REAL SHDFAC, RCMIN, RGL, HS, SNUP REAL RSMTBL(13),RGLTBL(13),HSTBL(13) REAL SNUPX (13) 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/ C THE FOLLOWING IS THRESHOLD WATER-EQUIV SNOWDEPTH (IN M) C IMPLYING 100% SNOW COVERAGE DATA SNUPX/0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.04, * 0.04, 0.04, 0.04, 0.025, 0.04, 0.025/ SNUP = SNUPX ( VEGTYP ) RCMIN = RSMTBL(VEGTYP) RGL = RGLTBL(VEGTYP) HS = HSTBL(VEGTYP) IF(VEGTYP .EQ. 11) SHDFAC=0.0 C C CALCULATE ROOT DISTRIBUTION. C PRESENT VERSION ASSUMES UNIFORM DISTRIBUTION BASED ON SOIL LAYERS. C DO I=1,NROOT RTDIS(I) = -SLDPTH(I)/ZSOIL(NROOT) END DO RETURN END SUBROUTINE PRTBND(INDI,NOUT,jday, & IIDAY, IITIME, SFCTMP, Q2, & SFCPRS, Rg, * LW_in, Par_in, Par_out, rnet, & S_OBS, PRCP, wet, T1_OBS, * T_2cm, T_4cm, T_8cm, T_16cm, & T_32cm, T_64cm, sm_5cm, * sm_20cm, sm_60cm, w_dir, u_bar, eddyuw, * uprim2, vprim2, wprim2, H_OBS, ETA_OBS, * IMONTH, IDAY, IREC) IMPLICIT NONE CCC CCC CC CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NAME: PRINT BONDVILLE, IL (BND_) DATA: CC CC PURPOSE: TO WRITE OUT OBSERVATIONAL DATA THAT CAN BE CC ======= USED FOR VALIDATION. CC CC VARIABLES: CC ========= C OBSERVED SFLX STATE VARIABLES C C (T1_OBS) T1: SKIN (GRND SFC) TEMPERATURE (K) C (sm_XXcm) SMC(3 layers): SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION) C FOR LAYER CENTERED AT XX cm DEPTH. C (STC_OBS) STC(3 layers): SOIL TEMP (K) C C OBSERVED SFLX OUTPUT C C (ETA_OBS) ETA: ACTUAL LATENT HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM SURFACE) C (H_OBS) H: SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM SURFACE) C (S_OBS) S: SOIL HEAT FLUX (W M-2: NEGATIVE, IF DOWNWARD FROM SURFACE) C Q1: EFFECTIVE MIXING RATIO AT GRND SFC ( KG KG-1) C CC ORIGINAL BND_ DATA VARIABLE NAMES C C jday Julian Day C time LST, half hour ending (using IITIME) C Ta air temperature (C), at 3 m C RH relative humidity (%) at 3 m C (converted to specific humidity Q2) C Pres surface pressure in mb C Rg incoming short wave radiation (W/m2) C Par_in incoming visible radiation (0.4-0.7 um) in uE/m2/s C Par_out outgoing or reflected visible light C Rnet net radiation (W/m2) C GHF (as S_OBS) soil or ground heat flux (W/m2) C C RAIN total rain for half hour (inches) C wet wetness sensor (volts), high values indicate wetness. C IRT (as T1_OBS) surface or skin temp (CONVERTED TO K) C T_2cm soil temp at 2 cm (C) C T_4cm soil tmep at 4 cm C T_8cm soil temp at 8 cm C T_16cm soil temp at 16 cm C T_32cm soil temp at 32 cm C T_64cm soil temp at 64 cm c C STC5 linear interpolations, using the soil C STC20 temperatures above, to match the model C STC60 depths (5,20 and 60 cm). C C sm_5cm soil volumetric water content at 5 cm zone C sm_20cm soil volumetric water content at 20 cm zone C sm_60cm soil volumetric water content at 60 cm zone C (as SMC_OBS(3 layers)) C C w_dir wind direction C u_bar average wind vector speed (m/s), at 6m C eddyuw (u'w') kinematic shear stress (m2/s2) C uprim2 (u'2) streamwise velocity variance (m2/s2) C vprim2 (v'2) crosswind velocity variance (m2/s2) C wprim2 (w'2) vertical velocity variance (m2/s2) C H (as H_OBS) sensible heat flux (W/m2) C C LE (as ETA_OBS) latent heat flux (W/m2) CC The eddy covariance sensors are located at 6 m AGL CC The bulk density of the soil is 1.4 gm/cm3 CC The site is currently in corn stubble (like it would look CC after combining) CC CC The units uE/m2/s refer to micro Einsteins per square meter per CC second. A uE is 6.02 x 10 (17) photons. CC______________________________________________________________________ REAL T0 PARAMETER (T0=273.16) INTEGER INDI, NOUT, NASCII, IREC INTEGER IIDAY, IITIME INTEGER IMONTH, IDAY, jday REAL Q2 REAL SFCPRS REAL SDATE REAL DDTIME REAL Rg REAL LW_in REAL Par_in REAL Par_out REAL Rnet REAL S_OBS REAL PRCP REAL wet REAL T1_OBS REAL T_2cm REAL T_4cm REAL T_8cm REAL T_16cm REAL T_32cm REAL T_64cm REAL STC5 REAL STC20 REAL STC60 REAL sm_5cm REAL sm_20cm REAL sm_60cm REAL w_dir REAL u_bar REAL eddyuw REAL uprim2 REAL vprim2 REAL wprim2 REAL H_OBS REAL ETA_OBS REAL RMONTH REAL RIDAY REAL SFCTMP C ____________________________________________________________ C C For now, print soil temperature observations as they are, C T_2cm, T_4cm, T_8cm, T_16cm, T_32cm and T_64cm. C C In the future, may use a cubic spline or C a simple linear scheme to interpolate C to STC_OBS( 5cm ), STC_OBS( 20cm ) and STC_OBS( 60cm ). C ____________________________________________________________ C C DATA CONVERSIONS TO MATCH MODEL'S VARIABLES C Soil Temperature in K T_2cm = T_2cm + T0 T_4cm = T_4cm + T0 T_8cm = T_8cm + T0 T_16cm = T_16cm + T0 T_32cm = T_32cm + T0 T_64cm = T_64cm + T0 C SOIL TEMPERATURE AT SOIL MOISTURE MEASUREMENT C LEVELS (LINEAR INTERP) STC5 = T_4cm + (T_8cm - T_4cm)/4. STC20 = T_16cm + 4*(T_32cm - T_16cm)/16. STC60 = T_64cm + 4*(T_32cm - T_64cm)/32. C ____________________________________________________________ C C ***** WRITE UNFORMATTED FOR GRADS OUTPUT ******** C IF (NOUT .GT. 0) THEN c convert to real for binary output SDATE = FLOAT(IIDAY) DDTIME = FLOAT(IITIME) RMONTH = FLOAT(IMONTH) RIDAY = FLOAT(IDAY ) WRITE(NOUT,REC=IREC) SDATE IREC = IREC + 1 WRITE(NOUT,REC=IREC) DDTIME IREC = IREC + 1 WRITE(NOUT,REC=IREC) SFCTMP IREC = IREC + 1 WRITE(NOUT,REC=IREC) Q2 IREC = IREC + 1 WRITE(NOUT,REC=IREC) SFCPRS IREC = IREC + 1 WRITE(NOUT,REC=IREC) Rg IREC = IREC + 1 WRITE(NOUT,REC=IREC) LW_in IREC = IREC + 1 WRITE(NOUT,REC=IREC) Par_in IREC = IREC + 1 WRITE(NOUT,REC=IREC) Par_out IREC = IREC + 1 WRITE(NOUT,REC=IREC) rnet IREC = IREC + 1 WRITE(NOUT,REC=IREC) S_OBS IREC = IREC + 1 WRITE(NOUT,REC=IREC) PRCP IREC = IREC + 1 WRITE(NOUT,REC=IREC) wet IREC = IREC + 1 WRITE(NOUT,REC=IREC) T1_OBS IREC = IREC + 1 WRITE(NOUT,REC=IREC) T_2cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) T_4cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) T_8cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) T_16cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) T_32cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) T_64cm IREC = IREC + 1 C WRITE(NOUT,REC=IREC) STC5 IREC = IREC + 1 WRITE(NOUT,REC=IREC) STC20 IREC = IREC + 1 WRITE(NOUT,REC=IREC) STC60 IREC = IREC + 1 C WRITE(NOUT,REC=IREC) sm_5cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) sm_20cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) sm_60cm IREC = IREC + 1 WRITE(NOUT,REC=IREC) w_dir IREC = IREC + 1 WRITE(NOUT,REC=IREC) u_bar IREC = IREC + 1 WRITE(NOUT,REC=IREC) eddyuw IREC = IREC + 1 WRITE(NOUT,REC=IREC) uprim2 IREC = IREC + 1 WRITE(NOUT,REC=IREC) vprim2 IREC = IREC + 1 WRITE(NOUT,REC=IREC) wprim2 IREC = IREC + 1 WRITE(NOUT,REC=IREC) H_OBS IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETA_OBS IREC = IREC + 1 WRITE(NOUT,REC=IREC) RMONTH IREC = IREC + 1 WRITE(NOUT,REC=IREC) RIDAY IREC = IREC + 1 ELSE NASCII = -NOUT WRITE(NASCII,200) jday, IITIME, SFCTMP, Q2, & SFCPRS, Rg, LW_in, Par_in, Par_out, & rnet, S_OBS, PRCP, wet, T1_OBS, & T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm, & STC5, STC20, STC60, sm_5cm, sm_20cm, sm_60cm, & w_dir, u_bar, eddyuw, uprim2, vprim2, wprim2, & H_OBS, ETA_OBS, IMONTH, IDAY END IF 200 FORMAT(2(I6,1X),32(F15.4,1X),2(I6,1X)) RETURN END SUBROUTINE PRTDAILY(INDI,NOUT,NSOIL,jday, & IIDAY,ITIME,ETADAY,ETADAY_O, & PCPDAY,SMCDIF,SMCDIF_O,RUNOFFDAY, & DT,IREC) IMPLICIT NONE CCCC CCCC C From MAIN: C C C C CALL PRTDAILY(INDI,NDAILY,NSOIL,IIDAY,ITIME,ETADAY,ETADAY_O, C C & PCPDAY,SMCDIF,SMCDIF_O,RUNOFFDAY, C C & DT,IRECD) C C C CCCC CCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NAME: PRINT DAILY ACCUMULATED HYDROLOGICAL VARIABLES CC CC VARIABLES: CC ========= CC CC LABEL I/O TYPE ............DESCRIPTION............... CC CC AET LOC RS ACTUAL EVAPOTRANSPIRATIVE ENERGY CC (J M-2 S-1) CC AET2 LOC RS ACTUAL EVAPOTRANSPIRATION (MM) CC CH IOC RS DRAG COEF FOR HEAT/MOISTURE CC CM IOC RS DRAG COEFFICIENT FOR MOMENTUM CC CMC IOC RS CANOPY MOISTURE CONTENT (M) CC DEW IB RS DEWFALL AMOUNT (M S-1) CC DRIP IB RS EXCESS CANOPY MOISTURE (M) CC EC(EC1) IB RS CANOPY EVAPORATION (M S-1) CC EDIR(EDIR1)IB RS DIRECT SOIL EVAPORATION (M S-1) CC ETA IOC RA FINAL ACTUAL EVAPOTRANSP (KG M-2 S-1) CC ETAS IOC RA FINAL ACTUAL EVAPOTRANSP (MM) CC ETP IOC RA FINAL POTNTL EVAPOTRANSP (KG M-2 S-1) CC ETPS IOC RA FINAL POTNTL EVAPOTRANSP (MM) CC ETT(ETT1) IB RS ACCUM PLANT TRANSPIRATION (M S-1) CC F IOC RS NET FLUX (TOT DOWNWARD RADIATION) CC FLX1 IB RS 1ST FLUX VALUE (W M-2) CC FLX2 IB RS 2ND FLUX VALUE (W M-2) CC FLX3 IB RS 3RD FLUX VALUE (W M-2) CC FOG IC LS INTERMEDIATE HRLY FOG FLAG CC FUP IB RS UPWARD GRND LW RADIATION (W M-2) CC H LOC RS SENSIBLE HEAT FLUX (W M-2) CC HEAT LOC RS SENSIBLE HEAT FLUX SUB-PRODUCT CC HEMI IC IS CURRENT HEMISPHERE (1=N, 2=S) CC I IC IS 1/8 MESH I COORDINATE CC ICLAMT IC IA INTERMEDIATE HRLY CLD AMOUNT CC ICLTYP IC IA INTERMEDIATE HRLY CLD TYPE CC J IC IS 1/8 MESH J COORDINATE CC K LOC IS LOOP INDEX CC NSOIL IC IS SOIL LAYER NUMBER CC PET LOC RS POTENTIAL EVAPOTRANSPIRATION (MM) CC PRCP IOC RA HALF HOURLY PRECIP AMT (KG M-2 S-1) CC Q2 IOC RS MIXING RATIO AT 1ST MDL LVL ABV SKIN CC Q2SAT IOC RS SAT MXNG RATIO AT 1ST MDL LVL ABV SKIN CC RES LOC RS ENERGY BALANCE EQN RESIDUAL (W M-2) CC RIB IB RS BULK RICHARDSON NUMBER CC RLDOWN IC RS DOWNWARD LONGWAVE RADIATION (W M-2) CC RR LOC RS SENSIBLE HEAT SUB-PRODUCT CC RSOLIN IC RS SOLAR RADIATION (W M-2) CC RUNOFF1 IB RS GRND SFC RUNOFF (M ) CC RUNOFF2 IB RS UNDERGROUND RUNOFF (M ) CC RUNOFF3 IB RS RUNOFF WITHIN SOIL LAYERS (M ) CC SMC IOC RA SOIL MOISTURE CONTENT (VOLUMETRIC) CC ZSOIL IOC RA SOIL LAYER DEPTH ( M ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER JDAY INTEGER IIDAY INTEGER ITIME INTEGER INDI INTEGER IREC INTEGER NOUT INTEGER NASCII INTEGER NSOIL REAL CMCS REAL DT REAL EC1S REAL EDIR1S REAL ETT1S REAL RTIME REAL SDATE REAL ETADAY, ETADAY_O, PCPDAY REAL SMCDIF, SMCDIF_O, RUNOFFDAY IF (NOUT .GT. 0) THEN SDATE=FLOAT(IIDAY) RTIME=FLOAT(ITIME) WRITE(NOUT,REC=IREC) SDATE IREC = IREC + 1 WRITE(NOUT,REC=IREC) RTIME IREC = IREC + 1 WRITE(NOUT,REC=IREC) PCPDAY IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETADAY IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETADAY_O IREC = IREC + 1 WRITE(NOUT,REC=IREC) RUNOFFDAY IREC = IREC + 1 WRITE(NOUT,REC=IREC) EDIR1S IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETT1S IREC = IREC + 1 WRITE(NOUT,REC=IREC) EC1S IREC = IREC + 1 WRITE(NOUT,REC=IREC) CMCS IREC = IREC + 1 WRITE(NOUT,REC=IREC) SMCDIF IREC = IREC + 1 WRITE(NOUT,REC=IREC) SMCDIF_O IREC = IREC + 1 ELSE NASCII = -NOUT WRITE(NASCII,200) & jday, & ITIME, & PCPDAY, & ETADAY, & ETADAY_O, & RUNOFFDAY, & EDIR1S, & ETT1S, & EC1S, & CMCS, & SMCDIF, & SMCDIF_O END IF 200 FORMAT(I6,1X,I6,10(1x,F15.4)) RETURN END SUBROUTINE PRTHMF(INDI,NOUT,NSOIL,jday, & IIDAY,IITIME,CH,CM, & Z,F,S,CMC, & SMCMAX,SFCTMP,T1,Q1,SFCPRS, & SFCSPD,ETA,ETP,PRCP, & STC,DT,RNET,H,AET,RES,IREC) IMPLICIT NONE CCC CCC C From MAIN: C C CALL PRTHMF(INDI,NOUT3,NSOIL,IIDAY,itime,CH,CM,Z,F,S,CMC, C C & SMCMAX,SFCTMP,T1,Q1,SFCPRS,SFCSPD,ETA,ETP,PRCP, C C & STC,DT,RNET,H,AET,IREC3) C CCC CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NAME: PRINT THERMODYNAMIC VARIABLES CC CC CPC: N/A CC === CC CC PURPOSE: TO CALC AND WRITE OUT PARAMETERS TO AID IN SCIENTIFIC CC ======= TUNING AND DEBUGGING OF THE OSU SOIL HYDROLOGY MODEL. CC THIS ROUTINE IS ONLY INVOKED FOR POINTS SELECTED BY CC SPECIFICATION IN THE CONTROL FILE. CC CC METHOD: 1. CALC NEEDED PARTIAL PRODUCTS AND SUMS. CC ====== 2. WRITE OUT RESULTS TO PRINT$ FILE. CC CC REFERENCES: FUNCTIONAL DESCRIPTION, SUBSYSTEM SPECIFICATION CC CC VARIABLES: CC ========= CC CC LABEL I/O TYPE ............DESCRIPTION............... CC CC AET LOC RS ACTUAL EVAPOTRANSPIRATIVE ENERGY CC (J M-2 S-1) CC AET2 LOC RS ACTUAL EVAPOTRANSPIRATION (MM) CC CH IOC RS DRAG COEF FOR HEAT/MOISTURE CC CM IOC RS DRAG COEFFICIENT FOR MOMENTUM CC CMC IOC RS CANOPY MOISTURE CONTENT (M) CC DEW IB RS DEWFALL AMOUNT (M S-1) CC DRIP IB RS EXCESS CANOPY MOISTURE (M) CC EC IB RS CANOPY EVAPORATION (M S-1) CC EDIR IB RS DIRECT SOIL EVAPORATION (M S-1) CC ETA IOC RA FINAL ACTUAL EVAPOTRANSP (KG M-2 S-1) CC ETP IOC RA FINAL POTNTL EVAPOTRANSP (KG M-2 S-1) CC ETPS IOC RA FINAL POTNTL EVAPOTRANSP (J M-2 S-1) CC ETT IB RS ACCUM PLANT TRANSPIRATION (M S-1) CC F IOC RS NET FLUX (TOT DOWNWARD RADIATION) CC FLX1 IB RS 1ST FLUX VALUE (W M-2) CC FLX2 IB RS 2ND FLUX VALUE (W M-2) CC FLX3 IB RS 3RD FLUX VALUE (W M-2) CC FUP IB RS UPWARD GRND LW RADIATION (W M-2) CC H LOC RS SENSIBLE HEAT FLUX (W M-2) CC HEAT LOC RS SENSIBLE HEAT FLUX SUB-PRODUCT CC I IC IS 1/8 MESH I COORDINATE CC ICLAMT IC IA INTERMEDIATE HRLY CLD AMOUNT CC ICLTYP IC IA INTERMEDIATE HRLY CLD TYPE CC J IC IS 1/8 MESH J COORDINATE CC K LOC IS LOOP INDEX CC NSOIL IC IS SOIL LAYER NUMBER CC PET LOC RS POTENTIAL EVAPOTRANSPIRATION (MM) CC PRCP IOC RA HALF HOURLY PRECIP AMT (KG M-2 S-1) CC Q2 IOC RS MIXING RATIO AT 1ST MDL LVL ABV SKIN CC Q2SAT IOC RS SAT MXNG RATIO AT 1ST MDL LVL ABV SKIN CC RES LOC RS ENERGY BALANCE EQN RESIDUAL (W M-2) CC RIB IB RS BULK RICHARDSON NUMBER CC RLDOWN IC RS DOWNWARD LONGWAVE RADIATION (W M-2) CC RR LOC RS SENSIBLE HEAT SUB-PRODUCT CC RSOLIN IC RS SOLAR RADIATION (W M-2) CC RUNOFF IB RS GRND SFC RUNOFF (M S-1) CC S IC RS GRND SFC FLUX (W M-2) CC SFCPRS IOC RS SFC PRESSURE (PASCALS) CC SFCSPD IC RS SFC WIND SPEED (M S-1) CC SFCTMP IOC RS SFC TEMPERATURE (K) CC SMC IOC RA SOIL MOISTURE CONTENT (VOLUMETRIC) CC SMCMAX IOC RS MAXIMUM SOIL MOISTURE CONTENT LIMIT CC SNODEP IOC RA SNOW DEPTH ( M ) CC STC IOC RA SOIL TEMPERATURE ( 5 CM AND 95 CM ) CC T1 IOC RS SKIN (GRND SFC) TEMPERATURE (K) CC T14 LOC RS GRND SFC TEMP TO THE 4TH POWER (K+4) CC TIME IOC IS 1 HRLY TIME (LOOP INDEX) CC Z IOC RS HT ABOVE GRND LVL (M) CC ZSOIL IOC RA SOIL LAYER DEPTH ( M ) INTEGER IIDAY INTEGER jday INTEGER IITIME INTEGER INDI INTEGER IREC INTEGER NOUT INTEGER NASCII INTEGER NSOIL REAL RNET REAL AET REAL AET2 REAL BETA REAL CH REAL CM REAL CMC REAL DEW REAL DRIP REAL DT REAL EC REAL EDIR REAL ETA REAL ETP REAL ETPS REAL ETT REAL F REAL FLX1 REAL FLX2 REAL FLX3 REAL FUP REAL H REAL PRCP REAL Q1 REAL RES REAL RIB REAL RUNOFF,RUNOXX3 REAL S REAL SDATE REAL RTIME REAL SFCPRS REAL SFCSPD REAL SFCTMP REAL SMCMAX REAL STC ( 10 ) REAL T1 REAL T14 REAL Z C COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOFF, & DEW,RIB,RUNOXX3 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT ACTUAL EVAPOTRANSPIRATION VALUE FROM KG M-2 S-1 TO C W M-2 (BY MULTIPLYING BY 2.5E+6) FOR USE IN ENERGY BALANCE. C ALSO CHANGE ACTUAL AND POTENTIAL EVAPORTRANSPIRATION VALUES C FROM KG M-2 S-1 TO MM (BY MULTIPLYING BY DT). CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC T14 = T1 * T1 * T1 * T1 AET = ETA ETPS = 2.56E+6 * ETP AET2 = DT * ETA FUP = 5.67E-8 * T14 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RESIDUAL OF ALL SURFACE ENERGY BALANCE EQN TERMS: C RES = F - H - S - AET - FUP - FLX1 - FLX2 - FLX3 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C NET RADIATION C RNET = F - FUP C IF (NOUT .GT. 0) THEN SDATE=FLOAT(IIDAY) RTIME=FLOAT(IITIME) WRITE(NOUT,REC=IREC) SDATE IREC = IREC + 1 WRITE(NOUT,REC=IREC) RTIME IREC = IREC + 1 WRITE(NOUT,REC=IREC) F IREC = IREC + 1 WRITE(NOUT,REC=IREC) RNET IREC = IREC + 1 WRITE(NOUT,REC=IREC) CH IREC = IREC + 1 WRITE(NOUT,REC=IREC) H IREC = IREC + 1 WRITE(NOUT,REC=IREC) S IREC = IREC + 1 WRITE(NOUT,REC=IREC) AET IREC = IREC + 1 WRITE(NOUT,REC=IREC) RES IREC = IREC + 1 WRITE(NOUT,REC=IREC) T1 IREC = IREC + 1 WRITE(NOUT,REC=IREC) Q1 IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETPS IREC = IREC + 1 WRITE(NOUT,REC=IREC) PRCP IREC = IREC + 1 WRITE(NOUT,REC=IREC) STC(1) IREC = IREC + 1 WRITE(NOUT,REC=IREC) STC(2) IREC = IREC + 1 WRITE(NOUT,REC=IREC) STC(3) IREC = IREC + 1 WRITE(NOUT,REC=IREC) STC(4) IREC = IREC + 1 ELSE NASCII = -NOUT WRITE(NASCII,200) & jday, & IITIME, & F, & RNET, & CH, & H, & S, & AET, & RES, & T1, & Q1, & ETPS, & PRCP, & STC(1), & STC(2), & STC(3), & STC(4) END IF 200 FORMAT(I6,1X,I6,15(1x,F15.4)) RETURN END C SUBROUTINE PRTHYDF(INDI,NOUT,NSOIL,jday, & IIDAY,DDTIME,ETA, & ETP,PRCP,SH2O,SMC,ALBEDO,ALB,SNOALB, & SNEQV, SNMAX, SNOWH, SOILM, SOILW, & RUNOFF1, & RUNOFF2, & RUNOFF3, & DT,EDIR1,EC1,ETT1,CMC,IREC) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCC CCC From MAIN: CCCCC C C C CALL PRTHYDF(INDI,NOUT1,NSOIL,jday, C C & IIDAY,ITIME,ETA, C C & ETP,PRCP,SH2O,SMC,ALBEDO,ALB,SNOALB, C C & RUNOFF1,RUNOFF2,RUNOFF3, C C & SNEQV,SNMAX,SNOWH,SOILM,SOILW, C C & DT,EDIR ,EC ,ETT ,CMC, IREC1) C C C CCCCC CCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C NAME: PRINT HYDROLOGICAL VARIABLES C C VARIABLES: C ========= C C LABEL I/O TYPE ............DESCRIPTION.............. C C ALB SNOW-FREE ALBEDO (FRACTION) C AET LOC RS ACTUAL EVAPOTRANSPIRATIVE ENERGY C (J M-2 S-1) C AET2 LOC RS ACTUAL EVAPOTRANSPIRATION (MM) C CH IOC RS DRAG COEF FOR HEAT/MOISTURE C CM IOC RS DRAG COEFFICIENT FOR MOMENTUM C CMC IOC RS CANOPY MOISTURE CONTENT (M) C DEW IB RS DEWFALL AMOUNT (M S-1) C DRIP IB RS EXCESS CANOPY MOISTURE (M) C EC(EC1) IB RS CANOPY EVAPORATION (M S-1) C EDIR(EDIR1)IB RS DIRECT SOIL EVAPORATION (M S-1) C ETA IOC RA FINAL ACTUAL EVAPOTRANSP (KG M-2 S-1) C ETAS IOC RA FINAL ACTUAL EVAPOTRANSP (MM) C ETP IOC RA FINAL POTNTL EVAPOTRANSP (KG M-2 S-1) C ETPS IOC RA FINAL POTNTL EVAPOTRANSP (MM) C ETT(ETT1) IB RS ACCUM PLANT TRANSPIRATION (M S-1) C F IOC RS NET FLUX (TOT DOWNWARD RADIATION) C FLX1 IB RS 1ST FLUX VALUE (W M-2) C FLX2 IB RS 2ND FLUX VALUE (W M-2) C FLX3 IB RS 3RD FLUX VALUE (W M-2) C FOG IC LS INTERMEDIATE HRLY FOG FLAG C FUP IB RS UPWARD GRND LW RADIATION (W M-2) C H LOC RS SENSIBLE HEAT FLUX (W M-2) C HEAT LOC RS SENSIBLE HEAT FLUX SUB-PRODUCT C HEMI IC IS CURRENT HEMISPHERE (1=N, 2=S) C I IC IS 1/8 MESH I COORDINATE C ICLAMT IC IA INTERMEDIATE HRLY CLD AMOUNT C ICLTYP IC IA INTERMEDIATE HRLY CLD TYPE C J IC IS 1/8 MESH J COORDINATE C K LOC IS LOOP INDEX C NSOIL IC IS SOIL LAYER NUMBER C PET LOC RS POTENTIAL EVAPOTRANSPIRATION (MM) C PRCP IOC RA HALF HOURLY PRECIP AMT (KG M-2 S-1) C Q2 IOC RS MIXING RATIO AT 1ST MDL LVL ABV SKIN C Q2SAT IOC RS SAT MXNG RATIO AT 1ST MDL LVL ABV SKIN C RES LOC RS ENERGY BALANCE EQN RESIDUAL (W M-2) C RIB IB RS BULK RICHARDSON NUMBER C RLDOWN IC RS DOWNWARD LONGWAVE RADIATION (W M-2) C RR LOC RS SENSIBLE HEAT SUB-PRODUCT C RSOLIN IC RS SOLAR RADIATION (W M-2) C RUNOFF1 IB RS GRND SFC RUNOFF (M ) C RUNOFF2 IB RS UNDERGROUND RUNOFF (M ) C RUNOFF3 IB RS RUNOFF WITHIN SOIL LAYERS (M ) C SMC IOC RA SOIL MOISTURE CONTENT (VOLUMETRIC) C SH2O(NSOIL): UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION) C SNEQV: WATER EQUIVALENT SNOW DEPTH (M) (formerly called SNODEP) C SNMAX: SNOW MELT (M) (WATER EQUIVALENT) C SNOWH: SNOW DEPTH (M) C SNOALB: MAX ALBEDO OVER DEEP SNOW (FRACTION) C SOILM: TOTAL SOIL COLUMN WATER CONTENT (M) C SOILW: AVAILABLE SOIL MOISTURE (UNITLESS FRACTION) C ZSOIL IOC RA SOIL LAYER DEPTH ( M ) INTEGER IIDAY INTEGER jday INTEGER DDTIME INTEGER INDI INTEGER IREC INTEGER NASCII INTEGER NOUT INTEGER NSOIL C REAL ALB REAL ALBEDO REAL CMC REAL CMCS REAL DT REAL EC1 REAL EC1S REAL EDIR1 REAL EDIR1S REAL ETA REAL ETP REAL ETT1 REAL ETT1S REAL PRCP REAL RTIME REAL RUNOFF1 REAL RUNOFF11 REAL RUNOFF2 REAL RUNOFF22 REAL RUNOFF3 REAL RUNOFF33 REAL SDATE REAL SMC(10) REAL SH2O(10) REAL SNEQV REAL SNMAX REAL SNOALB REAL SNOWH REAL SOILM REAL SOILW REAL SNEQVMM REAL SNMELT REAL SOILMM REAL TRUNOFF CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC EC1S=EC1*1000.*DT EDIR1S=EDIR1*1000.*DT ETT1S = ETT1*1000.*DT CMCS=CMC*1000. RUNOFF11 = RUNOFF1*DT*1000. RUNOFF22 = RUNOFF2*DT*1000. RUNOFF33 = RUNOFF3*1000. TRUNOFF = RUNOFF11 + RUNOFF22 + RUNOFF33 SNEQVMM = SNEQV*1000. SNMELT = SNMAX*1000. SOILMM = SOILM*1000. IF (NOUT .GT. 0) THEN SDATE=FLOAT(IIDAY) RTIME=FLOAT(DDTIME) WRITE(NOUT,REC=IREC) SDATE IREC = IREC + 1 WRITE(NOUT,REC=IREC) RTIME IREC = IREC + 1 WRITE(NOUT,REC=IREC) PRCP IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETP IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETA IREC = IREC + 1 WRITE(NOUT,REC=IREC) RUNOFF11 IREC = IREC + 1 WRITE(NOUT,REC=IREC) RUNOFF22 IREC = IREC + 1 WRITE(NOUT,REC=IREC) RUNOFF33 IREC = IREC + 1 WRITE(NOUT,REC=IREC) TRUNOFF IREC = IREC + 1 WRITE(NOUT,REC=IREC) EDIR1S IREC = IREC + 1 WRITE(NOUT,REC=IREC) ETT1S IREC = IREC + 1 WRITE(NOUT,REC=IREC) EC1S IREC = IREC + 1 WRITE(NOUT,REC=IREC) CMCS IREC = IREC + 1 WRITE(NOUT,REC=IREC) SH2O(1) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SH2O(2) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SH2O(3) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SH2O(4) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SMC(1) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SMC(2) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SMC(3) IREC = IREC + 1 WRITE(NOUT,REC=IREC) SMC(4) IREC = IREC + 1 WRITE(NOUT,REC=IREC) ALBEDO IREC = IREC + 1 WRITE(NOUT,REC=IREC) ALB IREC = IREC + 1 WRITE(NOUT,REC=IREC) SNOALB IREC = IREC + 1 WRITE(NOUT,REC=IREC) SNOWH IREC = IREC + 1 WRITE(NOUT,REC=IREC) SOILW IREC = IREC + 1 WRITE(NOUT,REC=IREC) SNEQVMM IREC = IREC + 1 WRITE(NOUT,REC=IREC) SNMELT IREC = IREC + 1 WRITE(NOUT,REC=IREC) SOILMM IREC = IREC + 1 ELSE NASCII = -NOUT WRITE(NASCII,200) jday, DDTIME, PRCP, ETP, & ETA, & RUNOFF11, RUNOFF22, RUNOFF33, TRUNOFF, & EDIR1S, ETT1S, EC1S, CMCS, & SH2O(1), SH2O(2), SH2O(3), SH2O(4), & SMC(1), SMC(2), SMC(3), SMC(4), & ALBEDO,ALB,SNOALB, & SNOWH, & SOILW, & SNEQVMM, & SNMELT, & SOILMM END IF 200 FORMAT(I6,1X,I6,15(1x,F15.4)) RETURN END C SUBROUTINE QDATAP (T,P,RH,QD,QS,ES) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: OBTAIN SPECIFIC HUMIDITY (q) FROM RELATIVE HUMIDITY CC AND GIVEN PRESSURE AND TEMPERATURE. CC CC CC FORMULAS AND CONSTANTS FROM ROGERS AND YAU, 1989: 'A CC SHORT COURSE IN CLOUD PHYSICS', PERGAMON PRESS, 3rd ED. CC CC Pablo J. Grunmann, 3/6/98. CC Updated to eliminate subroutine SVP, 6/24/98. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C---------------------------------------- C In: C T Temperature (K) C P Pressure (Pa) C RH Relative humidity (%) C----------------------------------------- C Out: C QD Specific humidity (Kg/Kg) C QS Saturation Specific humidity (Kg/Kg) C ES Saturation vapor pressure for water (Pa) C---------------------------------------- REAL T REAL P REAL RH REAL RHF REAL QD REAL QS REAL ES REAL EP REAL EPS REAL E PARAMETER (eps=0.622 ) C C ABOUT THE PARAMETER: C C eps ---------- (Water)/(dry air) molecular mass ratio, epsilon C _____________________________________________________________________ C C function E(T) = Sat. vapor pressure (in Pascal) at C temperature T (uses Clausius-Clapeyron). Es = E(T) C CONVERT REL. HUMIDITY (%) TO THE FRACTIONAL VALUE RHF = RH/100. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CALCULATE SATURATION MIXING RATIO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C QS = 0.622 * ES /P was substituted by a more precise C formula: -PABLO J. GRUNMANN, 05/28/98. QS = 0.622 * ES /(P - (1.-0.622)*ES) C C CONVERSION FROM REL. HUMIDITY: C (Rogers, pg. 17) C EP = (P*Es*RHF)/(P - Es*(1. - RHF)) QD = eps*EP/(P - (1. - eps)*EP) C RETURN END SUBROUTINE READBND(jday, itime, SFCTMP, RH, SFCPRS, Rg, * Par_in, Par_out, rnet, LW_in, GHF, PRCP, wet, SKN_IRT, * T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm, sm_5cm, * sm_20cm, sm_60cm, w_dir, u_bar, eddyuw, * uprim2, vprim2, wprim2, H, LE, * IMONTH,IDAY,NREAD) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NAME: READBND CC CC CC PURPOSE: CC THIS IS A MODIFICATION OF V.KOREN'S CC SUBROUTINE READPILPS(ICALB,NREAD,NREAD2,IYEAR,IMONTH,IDAY, CC IHOUR,P,T,Q,U,V,LWDN,RAIN,SOLDN,SOUT, CC TSKNEST) CC TO READ FORCING DATA FROM BONDVILLE, IL ATDD SITE. CC CC PABLO J. GRUNMANN, 05/98 CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ARGUMENT LIST IN THE CALL READBND C PREPARED BY PABLO GRUNMANN IN 05/98 C 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C INTEGER jday,itime,IMONTH,IDAY,NREAD,I INTEGER JULM(13) REAL SFCTMP REAL T REAL RH REAL SFCPRS REAL P REAL Rg REAL Par_in REAL Par_out REAL rnet REAL LW_in REAL GHF REAL PRCP REAL RAIN REAL wet REAL SKN_IRT REAL T_2cm REAL T_4cm REAL T_8cm REAL T_16cm REAL T_32cm REAL T_64cm REAL sm_5cm REAL sm_20cm REAL sm_60cm REAL w_dir REAL u_bar REAL eddyuw REAL uprim2 REAL vprim2 REAL wprim2 REAL w_speed REAL CO2 REAL H REAL LE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c Save previous step required forcing data c to make up in case of missing. REAL XOLD(32) do i=1,32 XOLD(i)=0.0 end do XOLD(2)= SFCTMP - 273.16 XOLD(3)= RH XOLD(5)= SFCPRS/1.E2 XOLD(6)= Rg XOLD(7)= LW_in XOLD(12)= 1800.*PRCP/25.4 XOLD(14)= SKN_IRT XOLD(25)= u_bar CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C READ BND_ DATA C ATTENTION! C REMEMBER TO REMOVE TITLE FROM BND_ FILE (1ST LINE), THEN READ DATA READ(NREAD,*) jday, itime, w_speed, w_dir, * T, RH, P, Rg, Par_in, Par_out, * rnet, GHF, RAIN, wet, SKN_IRT, * T_2cm, T_4cm, T_8cm, T_16cm, T_32cm, T_64cm, * u_bar, eddyuw, uprim2, vprim2, wprim2, * H, LE, CO2, LW_in, * sm_5cm, sm_20cm, sm_60cm C----------------------------------------------------------------------- C TILDEN MEYERS' BND_ DATA: C C jday Julian Day C itime LST, half hour ending C Ta air temperature (C), at 3 m C RH relative humidity (%) at 3 m C Pres surface pressure in mb C Rg incoming short wave radiation (W/m2) C Par_in incoming visible radiation (0.4-0.7 um) in uE/m2/s C Par_out outgoing or reflected visible light C Rnet net radiation (W/m2) C GHF soil or ground heat flux (W/m2) C RAIN total rain for half hour (inches) C wet wetness sensor (voltage - higher values indicatE wetness) C SKN_IRT surface or skin temp (C) C T_2_cm soil temp at 2 cm (C) C T_4_cm soil tmep at 4 cm C T_8_cm soil temp at 8 cm C T_16_cm soil temp at 16 cm C T_32_cm soil temp at 32 cm C T_64_cm soil temp at 64 cm C sm_5_cm soil volumetric water content at 5 cm zone C sm_20_cm soil volumetric water content at 20 cm zone C sm_60_cm soil volumetric water content at 60 cm zone C w_dir wind direction C u_bar average wind vector speed (m/s), at 6m C u'w' kinematic shear stress (m2/s2) C u'2 streamwise velocity variance (m2/s2) C v'2 crosswind velocity variance (m2/s2) C w'2 vertical velocity variance (m2/s2) C H sensible heat flux (W/m2) C LE latent energy flux (W/m2) C C The eddy covariance sensors are located at 6 m AGL C The bulk density of the soil is 1.4 gm/cm3 C The site is currently in corn stubble (like it would look C after combining) C C The units uE/m2/s refer to micro Einsteins per square meter per C second. A uE is 6.02 x 10 (17) photons. C C C----------------------------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C THIS PART STILL REQUIRES ATTENTION: C C C C IN LACK OF FORCING DATA, CURRENT PROCEDURE HAS BEEN C C (AND IT IS) TO USE PREVIOUS TIME-STEP RECORDED VALUES. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF (T .LT. -73.) THEN T = XOLD(2) ENDIF IF (RH .LT. 0. ) THEN RH = XOLD(3) ENDIF IF (P .LT. 850. ) THEN P = XOLD(5) ENDIF IF (Rg .LT. -100. ) THEN Rg = XOLD(6) ENDIF IF (LW_in .LT. 0.) THEN LW_in = XOLD(7) ENDIF IF (RAIN .LT. 0. ) THEN RAIN = XOLD(12) ENDIF IF (u_bar .LT. 0. ) THEN C MISSING u_bar, USE W_SPEED WITH LINEAR REGRESSION CONVERSION u_bar = 0.85*w_speed C BUT, IF EVEN W_SPEED IS NOT AVAILABLE, USE PREVIOUS STEP VALUE IF (u_bar .LT. 0. ) THEN u_bar = XOLD(25) ENDIF ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT VARIABLES C C OBTAIN DAY AND MONTH: JULIAN DATE SUBROUTINE CALL JULDATE(JDAY,IMONTH,IDAY,JULM) C C CONVERT RAIN FROM INCHES IN 30MIN TO KG M^-2 S^-1 C (SAME AS PRECIP.RATE IN MM/SEC) C PRCP = RAIN*25.4/1800. C C AIR TEMPERATURE IN KELVIN C SFCTMP = T + 273.16 C C SFC PRESSURE IN PASCAL C SFCPRS = P*1.E2 C C ---- DONE CONVERTING VARIABLES ---------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC RETURN END SUBROUTINE READCNTL(CNTRFL,NRUN,DT,NSOIL,NROOT,Z,SLDPTH, & SOILTP,VEGTYP,SLOPTYP,ALBEDOM,SHDFACM,MLAI,SNOALB,ICE, & TBOT,Z0,CZIL,REFKDT,T1, & STC,SMC,SH2O,CMC,SNODEP,SNEQV,FILENAME, & LATITUDE, & LONGITUDE, & JDAY, & TIME) IMPLICIT NONE C READS A CONTROL FILE WHICH INCLUDES INPUT PARAMETERS FOR OSU MODEL - INTEGER NRUN INTEGER NSOIL INTEGER NROOT INTEGER SOILTP INTEGER VEGTYP INTEGER SLOPTYP INTEGER ICE REAL SLDPTH(4) REAL STC(4) REAL SMC(4) REAL SH2O(4) REAL DT REAL Z REAL ALBEDOM(13) REAL SHDFACM(13) REAL MLAI(13) REAL SNOALB REAL TBOT REAL Z0 REAL T1 REAL CMC REAL SNODEP REAL SNEQV REAL CZIL REAL REFKDT REAL LATITUDE REAL LONGITUDE INTEGER JDAY INTEGER TIME INTEGER NREAD, I, K CHARACTER*72 CNTRFL, LINEFEED CHARACTER*72 FILENAME NREAD = 21 OPEN (UNIT=NREAD, FILE=CNTRFL, STATUS='OLD', & FORM='FORMATTED') C JUMP-READ 2 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED) DO I=1,2 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C Model Configuration: READ (NREAD, * ) LATITUDE READ (NREAD, * ) LONGITUDE READ (NREAD, 100) JDAY READ (NREAD, 100) TIME C READ (NREAD, 150) NRUN READ (NREAD, * ) DT READ (NREAD, 100) NSOIL READ (NREAD, 100) NROOT READ (NREAD, *) Z READ (NREAD, *) (SLDPTH(K), K=1,NSOIL) C NEW: NROOT, SLDPTH(K), K=1,NSOIL C REMOVED: ZSOIL(K), K=1,NSOIL C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED WRITE(*,'(A)') LINEFEED END DO C ATMOSPHERIC DATA FILE TO BE USED FOR FORCING: READ(NREAD,'(A)') FILENAME PRINT*,' FILENAME = ', FILENAME C C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C ------------------------ C Land surface characteristics: C ------------------------ READ (NREAD, 100) SOILTP READ (NREAD, 100) VEGTYP READ (NREAD, 100) SLOPTYP C ...........SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C ........... READ (NREAD, *) (ALBEDOM(K), K=1,12) ALBEDOM(13) = ALBEDOM(1) C ...........SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C ........... READ (NREAD, *) (SHDFACM(K), K=1,12) SHDFACM(13) = SHDFACM(1) C ...........SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C ........... READ (NREAD, *) (MLAI(K), K=1,12) MLAI(13) = MLAI(1) C ...........SKIP COMMENTS (LINE FEED) READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED C ........... READ (NREAD, *) SNOALB READ (NREAD, 100) ICE C NEW: ALBEDO, SHDFAC, ICE, C REMOVED: LAND, CFACTR C C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C ------------------ C PHYSICAL PARAMETERS: C ------------------ READ (NREAD, *) TBOT READ (NREAD, *) Z0 READ (NREAD, *) CZIL READ (NREAD, *) REFKDT C C READ 3 LINES (CHARACTER) TO SKIP COMMENTS (LINE FEED) DO I=1,3 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C ------------------------ C INITIAL STATE VARIABLES: C ------------------------ READ (NREAD, *) T1 READ (NREAD, *) (STC(K), K=1,NSOIL) READ (NREAD, *) (SMC(K), K=1,NSOIL) READ (NREAD, *) (SH2O(K), K=1,NSOIL) READ (NREAD, *) CMC READ (NREAD, *) SNODEP READ (NREAD, *) SNEQV C READ LAST 2 LINES (CHARACTER) TO SHOW THAT CONTROLFILE WAS C COMPLETELY READ DO I=1,2 READ(NREAD,'(A)') LINEFEED C writ WRITE(*,'(A)') LINEFEED END DO C NEW: SNEQV C REMOVED: C ESD C KDT C DKSAT C B C DWSAT C SMCMAX C SMCREF C PSISAT C SMCWLT C DMAX CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C V. KOREN 5/21/97 C TWO PARAMETERS TO RUN FROZEN GROUND CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C READ (NREAD, *) CVFRZ C C READ (NREAD, *) FRZK C C C C READ (NREAD, *) SLOPE C C C C CVFRZ = 3.0 C C FRZCL = 0.50 C C SLOPE = 0.0005 C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 100 FORMAT(I6) 150 FORMAT(I7) C 200 FORMAT(F14.8) C 400 FORMAT(10(2X, F14.8)) CLOSE(NREAD) C RETURN END C SUBROUTINE REDPRM(VEGTYP,SOITYP, SLOPTYP, PTU, + CFACTR, CMCMAX, RSMAX, TOPT,REFKDT,KDT, O Z0,CZIL,SHDFAC,RCMIN,RGL,HS, + ZBOT, FRZX, PSISAT, SLOPE, SNUP,SALP, O B,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF, + SMCDRY,F1,QUARTZ, O RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL) IMPLICIT NONE C --------------- FROZEN GROUND VERSION ------------------ C FROZEN GROUND CORRECTION FACTOR IS ADDED C 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, SLOPTYP INTEGER NSOIL REAL Z0, CZIL, SHDFAC, PTU, RCMIN, RGL, HS, FRZFACT, PSISAT REAL SLOPE, SNUP, SALP REAL B, DKSAT, DWSAT, SMCMAX, SMCWLT, SMCREF REAL SMCDRY, F1, QUARTZ REAL REFDK REAL REFKDT REAL KDT REAL FRZX REAL FRZK REAL RTDIS(NSOIL) REAL SLDPTH(NSOIL) REAL ZSOIL(NSOIL) C *** NEW CLASS PARAMETER 'SLOPTYP' WAS INCLUDED TO ESTIMATE C *** LINEAR RESERVOIR COEFFICIENT 'SLOPE' TO THE BASEFLOW RUNOFF C *** OUT OF THE BOTTOM LAYER. LOWEST CLASS (SLOPTYP=0)MEANS C *** HIGHEST SLOPE PARAMETER= 1 C *** DEFINITION OF SLOPTYP FROM 'ZOBLER' SLOPE TYPE C *** SLOPE CLASS PERCENT SLOPE C 1 0-8 C 2 8-30 C 3 > 30 C 4 0-30 C 5 0-8 & > 30 C 6 8-30 & > 30 C 7 0-8, 8-30, > 30 C 9 GLACIAL ICE C BLANK OCEAN/SEE C NOTE: CLASS 9 FROM 'ZOBLER' FILE SHOULD BE REPLACED BY 8 C AND 'BLANK' 9 C REAL CFACTR REAL CMCMAX REAL RSMAX REAL TOPT REAL ZBOT REAL SLOPES (9) DATA SLOPES/0.1, 0.6, 1.0, 0.35, 0.55, 0.8, 0.63, 0., 0./ C Specify depth[m] of lower boundary soil temperature ZBOT = 3.0 C C SPECIFY SNOW DISTRIBUTION SHAPE PARAMETER C SALP - SHAPE PARAMETER OF DISTRIBUTION FUNCTION C OF SNOW COVER. FROM ANDERSON'S DATA (HYDRO-17) C BEST FIT IS WHEN SALP = 2.6 SALP = 2.6 C SPECIFY THE VEG/SOIL PARAMETERS C Set two canopy water parameters CFACTR = 0.5 CMCMAX = 0.5E-3 C Set max. stomatal resistance RSMAX = 5000.0 C Set optimum transpiration air temperature TOPT = 298.0 C Set drainage parameter based on slope index SLOPE = SLOPES ( SLOPTYP) C SET-UP VEGETATION PARAMETERS CALL PRMVEG(VEGTYP,SHDFAC,RCMIN,RGL,HS,SNUP, O RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL) C SET-UP SOIL PARAMETERS CALL PRMSOI(SOITYP,B,SMCDRY,F1,SMCMAX,SMCREF, FRZFACT, & PSISAT,DKSAT,DWSAT,SMCWLT,QUARTZ) C ***** KDT IS DEFINED BY REFERENCE REFKDT AND DKSAT ***** C REFDK=2.E-6 IS THE SAT. DK. VALUE FOR THE SOIL TYPE 2 C REFDK = 2.0E-6 KDT = REFKDT * DKSAT/REFDK C ----------------------------------------------------------------------- C FROZEN GROUND PARAMETER, FRZK, DEFINITION C FRZK IS ICE CONTENT THRESHOLD ABOVE WHICH FROZEN SOIL IS IMPERMEABLE C REFERENCE VALUE OF THIS PARAMETER FOR THE LIGHT CLAY SOIL (TYPE=3) C FRZK = 0.15 M C FRZK = 0.15 C TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT C FRZFACT IS CALCULATED IN REDPAR SUBROUTINE C FRZX = FRZK * FRZFACT C ----------------------------------------------------------------------- RETURN END C SUBROUTINE ROSR12 ( P, A, B, C, D, DELTA, NSOIL ) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN CC ======= BELOW: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER K INTEGER KK INTEGER NSOIL INTEGER NSOLD PARAMETER ( NSOLD = 4 ) REAL P(NSOLD) REAL A(NSOLD) REAL B(NSOLD) REAL C(NSOLD) REAL D(NSOLD) REAL DELTA(NSOLD) 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 SFCDIF(ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS) IMPLICIT NONE C ****************************************************************** C * * C * SURFACE LAYER * C * * C ****************************************************************** C----------------------------------------------------------------------- REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, SQVISC REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, PSLHS REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, RLMA C REAL ZTFC INTEGER ITRMX, ILECH, ITR PARAMETER &(WWST=1.2,WWST2=WWST*WWST,G=9.8,VKRM=0.40,EXCM=0.001 &,BETA=1./270.,BTG=BETA*G,ELFC=VKRM*BTG &,WOLD=.15,WNEW=1.-WOLD,ITRMX=05,PIHF=3.14159265/2. C----------------------------------------------------------------------- &,EPSU2=1.E-4,EPSUST=0.07,EPSIT=1.E-4,EPSA=1.E-8 &,ZTMIN=-5.,ZTMAX=1.,HPBL=1000.0 &,SQVISC=258.2) PARAMETER &(RIC=0.183,RRIC=1.0/RIC,FHNEU=0.8,RFC=0.191 &,RFAC=RIC/(FHNEU*RFC*RFC)) C----------------------------------------------------------------------- C **** LECH'S SURFACE FUNCTIONS **** PSLMU(ZZ)=-0.96*log(1.0-4.5*ZZ) PSLMS(ZZ)=ZZ*RRIC-2.076*(1.-1./(ZZ+1.)) PSLHU(ZZ)=-0.96*log(1.0-4.5*ZZ) PSLHS(ZZ)=ZZ*RFAC-2.076*(1.-1./(ZZ+1.)) C----------------------------------------------------------------------- C **** PAULSON'S SURFACE FUNCTIONS ***** PSPMU(XX)=-2.*log((XX+1.)*0.5)-log((XX*XX+1.)*0.5)+2.*ATAN(XX) & -PIHF PSPMS(YY)=5.*YY PSPHU(XX)=-2.*log((XX*XX+1.)*0.5) PSPHS(YY)=5.*YY C ILECH=0 C*********************************************************************** C----------------------------------------------------------------------- C ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1 CC......ZTFC=0.1 C CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT C Commented out to allow argument-passed value C CZIL=0.2 ZILFC=-CZIL*VKRM*SQVISC C----------------------------------------------------------------------- ZU=Z0 CC.......ZT=Z0*ZTFC C----------------------------------------------------------------------- RDZ=1./ZLM CXCH=EXCM*RDZ C----------------------------------------------------------------------- DTHV=THLM-THZ0 C DU2=MAX(SFCSPD*SFCSPD,EPSU2) C--------------BELJARS CORRECTION OF USTAR------------------------------ BTGH=BTG*HPBL WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.) USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST) C--------------ZILITINKEVITCH APPROACH FOR ZT-------------------------------- ZT=EXP(ZILFC*SQRT(USTAR*Z0))*Z0 C----------------------------------------------------------------------- ZSLU=ZLM+ZU ZSLT=ZLM+ZT C RLOGU=log(ZSLU/ZU) RLOGT=log(ZSLT/ZT) C RLMO=ELFC*AKHS*DTHV/USTAR**3 C----------------------------------------------------------------------- DO 100 ITR=1,ITRMX C--------------1./MONIN-OBUKKHOV LENGTH-SCALE--------------------------- ZETALT=MAX(ZSLT*RLMO,ZTMIN) RLMO=ZETALT/ZSLT ZETALU=ZSLU*RLMO C ZETAU=ZU*RLMO ZETAT=ZT*RLMO C----------------------------------------------------------------------- IF(ILECH.EQ.0) THEN IF(RLMO.LT.0.)THEN XLU4=1.-16.*ZETALU XLT4=1.-16.*ZETALT XU4 =1.-16.*ZETAU XT4 =1.-16.*ZETAT C XLU=SQRT(SQRT(XLU4)) XLT=SQRT(SQRT(XLT4)) XU =SQRT(SQRT(XU4)) XT =SQRT(SQRT(XT4)) C PSMZ=PSPMU(XU) SIMM=PSPMU(XLU)-PSMZ+RLOGU PSHZ=PSPHU(XT) SIMH=PSPHU(XLT)-PSHZ+RLOGT ELSE ZETALU=MIN(ZETALU,ZTMAX) ZETALT=MIN(ZETALT,ZTMAX) C PSMZ=PSPMS(ZETAU) SIMM=PSPMS(ZETALU)-PSMZ+RLOGU PSHZ=PSPHS(ZETAT) SIMH=PSPHS(ZETALT)-PSHZ+RLOGT ENDIF ELSE C----------------------------------------------------------------------- C ***** LECH'S FUNCTIONS **** IF(RLMO.LT.0.)THEN PSMZ=PSLMU(ZETAU) SIMM=PSLMU(ZETALU)-PSMZ+RLOGU PSHZ=PSLHU(ZETAT) SIMH=PSLHU(ZETALT)-PSHZ+RLOGT ELSE ZETALU=MIN(ZETALU,ZTMAX) ZETALT=MIN(ZETALT,ZTMAX) C PSMZ=PSLMS(ZETAU) SIMM=PSLMS(ZETALU)-PSMZ+RLOGU PSHZ=PSLHS(ZETAT) SIMH=PSLHS(ZETALT)-PSHZ+RLOGT ENDIF ENDIF C--------------BELJAARS CORRECTION FOR USTAR---------------------------- USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST) C--------------ZILITINKEVITCH FIX FOR ZT-------------------------------- ZT=EXP(ZILFC*SQRT(USTAR*Z0))*Z0 ZSLT=ZLM+ZT RLOGT=log(ZSLT/ZT) C----------------------------------------------------------------------- USTARK=USTAR*VKRM AKMS=MAX(USTARK/SIMM,CXCH) AKHS=MAX(USTARK/SIMH,CXCH) C----------------------------------------------------------------------- WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.) RLMN=ELFC*AKHS*DTHV/USTAR**3 C----------------------------------------------------------------------- RLMA=RLMO*WOLD+RLMN*WNEW C----------------------------------------------------------------------- C IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110 C----------------------------------------------------------------------- RLMO=RLMA C----------------------------------------------------------------------- 100 CONTINUE C----------------------------------------------------------------------- 110 CONTINUE C----------------------------------------------------------------------- RETURN END C SUBROUTINE SFLX(ICE,SATURATED,DT,Z,NSOIL,NROOT,SLDPTH, I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2, I VEGTYP,SOITYP,SLOPTYP, I SHDFAC,XLAI,PTU,REFKDT,TBOT,ALB,SNOALB, I Z0,SFCSPD,CZIL, 2 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV, O CH,ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX,ALBEDO, O SOILW,SOILM) C IMPLICIT NONE CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: SUB-DRIVER FOR FAMILY OF PHYSICS SUBROUTINES OF A CC SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL CC MOISTURE, SOIL ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, CC SNOWPACK WATER CONTENT,SNOWDEPTH, AND ALL TERMS CC OF THE SURFACE ENERGY BALANCE AND SURFACE WATER CC BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF CC DOWNWARD RADIATION AND PRECIP) CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC C ------------ FROZEN GROUND VERSION ------------------------------- C NEW STATES: SH2O(NSOIL) - UNFROZEN SOIL MOISTURE C SNOWH - SNOW DEPTH C C ------------------------------------------------------------------------- C C NOTE ON SNOW STATE VARIABLES: C SNOWH = actual physical snow depth in m C SNEQV = liquid water-equivalent snow depth in m C (time-dependent snow density is obtained from SNEQV/SNOWH) C C NOTE ON ALBEDO FRACTIONS: C Input: C ALB = BASELINE SNOW-FREE ALBEDO, FOR JULIAN DAY OF YEAR C (USUALLY FROM TEMPORAL INTERPOLATION OF MONTHLY MEAN VALUES) C (CALLING PROG MAY OR MAY NOT INCLUDE DIURNAL SUN ANGLE EFFECT) C SNOALB = UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW C (E.G. FROM ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.) C Output: C ALBEDO = COMPUTED ALBEDO WITH SNOWCOVER EFFECTS C (COMPUTED USING ALB, SNOALB, SNEQV, AND SHDFAC->>veg greenness) C C ARGUMENT LIST IN THE CALL SFLX C C 1. CALLING STATEMENT C C CALL SFLX(ICE,LSATURATED,DT,Z,NSOIL,NROOT,SLDPTH, C I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,TH2,Q2,Q2SAT,DQSDT2, C I IVEGTYP,ISOILTP,ISLOPTYP, C I SHDFAC,XLAI,PTU,REFKDT,TBOT,ALB,SNOALB, C I Z0,SFCSPD,CZIL, C 2 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV, C O CH,ETP,ETA,H,S,RUNOFF1,RUNOFF2,Q1,SNMAX,ALBEDO, C O SOILW,SOILM) C C 2. INPUT C *** GENERAL PARAMETERS *** C C ICE: SEA-ICE FLAG (=1: SEA-ICE, =0: LAND) C SATURATED: ATMOSPHERIC SATURATION FLAG C (true: SATURATION -> ETP=0.0 IN SFLX C false: OTHERWISE) C DT: TIMESTEP (SEC) C Z: HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES C NSOIL: NUMBER OF SOIL LAYERS C NROOT: NUMBER OF ROOT-ZONE LAYERS C SLDPTH: THE THICKNESS OF EACH SOIL LAYER (M) C C *** ATMOSPHERIC VARIABLES *** C C LWDN: LW DOWNWARD RADIATION (W M-2; POSITIVE, IF DOWNWARD) C SOLDN: SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, IF DOWNWARD) C SFCPRS: PRESSURE AT HEIGHT Z ABOVE GROUND(PASCALS) C PRCP: PRECIP RATE (KG M-2 S-1) C SFCTMP: AIR TEMPERATURE (K) AT HEIGHT Z ABOVE GROUND C TH2: AIR POTENTIAL TEMPERATURE (K) AT GROUND SURFACE C Q2: MIXING RATIO AT HEIGHT Z ABOVE GROUND (KG KG-1) C Q2SAT: SAT MIXING RATIO AT HEIGHT Z ABOVE GROUND (KG KG-1) C DQSDT2: SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP (KG KG-1 K-1) C C *** CANOPY/SOIL CHARACTERISTICS *** C C VEGTYP: VEGETATION TYPE (INTEGER INDEX) C SOITYP: SOIL TYPE (INTEGER INDEX) C SLOPTYP: CLASS OF BOTTOM LAYER SLOPE (INTEGER INDEX) C SHDFAC: AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION C XLAI: LEAF-AREA INDEX (UNITLESS, IN THE RANGE 0.-5.) C PTU: PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS) C REFKDT: BASELINE VALUE FOR RUNOFF PARAMETER KDT C TBOT: BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN AIR TEMPERATURE) C ALB: BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION) C SNOALB: ALBEDO UPPER BOUND OVER DEEP SNOW (FRACTION) C C 3. INPUT AND OUTPUT STATE VARIABLES C C !!! *********** STATE VARIABLES ************** !!! C C CMC: CANOPY MOISTURE CONTENT (M) C T1: GROUND (AND/OR SNOW) SKIN TEMPERATURE (K) C STC(NSOIL): SOIL TEMP (K) C SMC(NSOIL): TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION) C SH2O(NSOIL): UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION) C SNOWH: SNOW DEPTH (M) C SNEQV: WATER-EQUIVALENT SNOW DEPTH (M) C ALBEDO: SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION) C *** NOTE: TSTATE VARIABLES ARE INPUT AND OUTPUT*** C C 4. OUTPUT C CH: SFC EXCH COEF FOR HEAT AND MOISTURE (M S-1) C ETP: POTENTIAL EVAPORATION (W M-2) C ETA: ACTUAL LATENT 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: SURFACE RUNOFF (M S-1) C RUNOFF2: SUBSURFACE RUNOFF (M S-1) C Q1: EFFECTIVE MIXING RATIO AT GRND SFC ( KG KG-1) C SNMAX: SNOW MELT (M) (WATER EQUIVALENT) C ALBEDO: SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION) C SOILM: TOTAL SOIL COLUMN WATER CONTENT (M) C SOILW: AVAILABLE SOIL MOISTURE (UNITLESS FRACTION BETWEEN C SOIL SATURATION AND WILTING POINT) INTEGER NSOLD PARAMETER (NSOLD=4) C LOGICAL SNOWNG LOGICAL FRZGRA LOGICAL SATURATED C INTEGER K INTEGER KZ INTEGER ICE INTEGER NSOIL,VEGTYP,SOITYP,NROOT INTEGER SLOPTYP C REAL R REAL CP C PARAMETER ( R = 287.04, CP = 1004.5 ) REAL B REAL CFACTR C CH IS SFC EXCHANGE COEF FOR HEAT/MOIST C CM IS SFC MOMENTUM DRAG (NOT NEEDED IN SFLX) REAL CH REAL CM C REAL CMC REAL CMCMAX REAL CZIL REAL DF1 REAL DKSAT REAL DT REAL DWSAT REAL EPSCA REAL SNEQV REAL ETA REAL ETP REAL F REAL F1 REAL KDT REAL LWDN REAL PC REAL PRCP REAL PTU REAL Q1 REAL Q2 REAL Q2SAT REAL RCH REAL REFKDT REAL RR REAL RTDIS (NSOLD) REAL S REAL SFCPRS REAL SFCSPD REAL SFCTMP REAL SHDFAC REAL SMC(NSOIL),STC(NSOIL),SLDPTH(NSOIL),SH2O(NSOIL) REAL SMCDRY REAL SMCMAX REAL SMCREF REAL SMCWLT CC.. REAL SNEQV REAL SNOWH REAL T1 REAL T1V REAL T1MIN REAL T24 REAL T2V REAL TBOT REAL TH2 REAL TH2V REAL Z REAL Z0 REAL ZSOIL ( NSOLD ) REAL ZBOT REAL TOPT, TFREEZ, SOLDN, DQSDT2, XLAI, RUNOFF1, RUNOFF2 REAL SNMAX, SOILW, RGL, SOILM REAL BETA, DRIP, EC, EDIR, ETT, FLX1, FLX2, FLX3, RUNOF, DEW REAL RIB, RUNOFF3, RSMAX, ALBEDO, RCMIN, HS, PSISAT REAL SLOPE, SNUP,SALP,SNOALB, ALB REAL QUARTZ, FRZX, SNDENS, SNCOND REAL CSNOW, RSNOW, SNOFAC REAL SN_NEW, PRCP1, SICE, RC, EDIR1, EC1, ETT1, SOILWM REAL SOILWW REAL H PARAMETER (TFREEZ = 273.16) CMIC$ TASKCOMMON RITE COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOF, & DEW,RIB,RUNOFF3 C C INITIALIZATION RUNOFF1 = 0.0 RUNOFF2 = 0.0 RUNOFF3 = 0.0 SNMAX = 0.0 IF(ICE .EQ. 1) THEN C SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS DO KZ = 1, NSOIL ZSOIL(KZ)=-3.*FLOAT(KZ)/FLOAT(NSOIL) END DO ELSE C CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO C BOTTOM OF EACH SOIL LAYER ZSOIL(1)=-SLDPTH(1) DO KZ = 2, NSOIL ZSOIL(KZ)=-SLDPTH(KZ)+ZSOIL(KZ-1) END DO ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, CC INCLUDING SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS. CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CALL REDPRM(VEGTYP,SOITYP, SLOPTYP,PTU, + CFACTR, CMCMAX, RSMAX, TOPT,REFKDT,KDT, O Z0,CZIL,SHDFAC,RCMIN,RGL,HS, + ZBOT, FRZX, PSISAT, SLOPE, SNUP, SALP, O B,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF, + SMCDRY,F1,QUARTZ, O RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC NEXT CALL ROUTINE SFCDIF TO CALCULATE CC THE SFC EXCHANGE COEF (CH) FOR HEAT AND MOISTURE CC CC NOTE NOTE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CC CC COMMENT OUT CALL SFCDIF, IF SFCDIF ALREADY CALLED CC IN CALLING PROGRAM (SUCH AS IN COUPLED ATMOSPHERIC MODEL) CC CC NOTE !! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, CC IN CASE ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND CC ZILINTINKEVICH COEF (CZIL) ARE SET THERE VIA NAMELIST I/O CC CC NOTE !! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD CC TIMES THE "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. CC HENCE THE CH RETURNED FROM SFCDIF HAS UNITS OF M/S. CC THE IMPORTANT COMPANION COEFFICIENT OF CH, CARRIED HERE AS "RCH", CC IS THE CH FROM SFCDIF TIMES AIR DENSITY AND PARAMETER "CP". CC "RCH" IS COMPUTED IN "CALL PENMAN" AND CARRIED IN COMMON/RITE/. CC RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATE IN EQNS. CC CC NOTE !! THE COMPANION COEFFICIENT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC VIRTUAL TEMPS AND POTENTIAL TEMPS FROM AIR TEMP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C TH2 is passed in C.....TH2 = SFCTMP + ( 0.0098 * Z ) T2V = SFCTMP * (1.0 + 0.61 * Q2 ) T1V = T1 * (1.0 + 0.61 * Q2 ) TH2V = TH2 * (1.0 + 0.61 * Q2 ) C CALL SFCDIF ( Z, Z0, T1V, TH2V, SFCSPD, CZIL, CM, CH ) C C C prt PRINT*,' 1ST DT=', DT C prt PRINT*,' *********** INPUT PARAMETERS *************' C prt PRINT*,' ' C prt PRINT*,'VEGTYP= ',VEGTYP C prt PRINT*,'SOITYP= ',SOITYP C prt PRINT*,'SLOPTYP= ',SLOPTYP C prt PRINT*,'CFACTR= ',CFACTR C prt PRINT*,'CMCMAX= ',CMCMAX C prt PRINT*,'RSMAX= ',RSMAX C prt PRINT*,'TOPT= ',TOPT C prt PRINT*,'REFKDT= ',REFKDT C prt PRINT*,'KDT= ',KDT C prt PRINT*,'Z0= ',Z0 C prt PRINT*,'CZIL= ',CZIL C prt PRINT*,'CH= ',CH C prt PRINT*,'SHDFAC= ',SHDFAC C prt PRINT*,'RCMIN= ',RCMIN C prt PRINT*,'RGL= ',RGL C prt PRINT*,'HS= ',HS C prt PRINT*,'ZBOT= ',ZBOT C prt PRINT*,'FRZX= ',FRZX C prt PRINT*,'PSISAT= ',PSISAT C prt PRINT*,'SLOPE= ',SLOPE C prt PRINT*,'SNUP= ',SNUP C prt PRINT*,'SALP= ',SALP C prt PRINT*,'B= ',B C prt PRINT*,'DKSAT= ',DKSAT C prt PRINT*,'DWSAT= ',DWSAT C prt PRINT*,'SMCMAX= ',SMCMAX C prt PRINT*,'SMCWLT= ',SMCWLT C prt PRINT*,'SMCREF= ',SMCREF C prt PRINT*,'SMCDRY= ',SMCDRY C prt PRINT*,'F1= ',F1 C prt PRINT*,'QUARTZ= ',QUARTZ C prt PRINT*,'RTDIS= ',RTDIS C prt PRINT*,'SLDPTH= ',SLDPTH C prt PRINT*,'ZSOIL= ',ZSOIL C prt PRINT*,'NROOT= ',NROOT C prt PRINT*,'NSOIL= ',NSOIL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INITIALIZE MISC VARIABLES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC S = 0.0 SNOWNG = .FALSE. FRZGRA = .FALSE. C IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP IF(ICE .EQ. 1) THEN SNEQV = 0.01 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DETERMINE IF IT'S PRECIPITATING 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. TFREEZ ) THEN SNOWNG = .TRUE. ELSE IF ( T1 .LE. TFREEZ ) 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(SNEQV .EQ. 0.0) THEN SNDENS = 0.0 SNOWH = 0.0 SNCOND = 1.0 ELSE SNDENS=SNEQV/SNOWH SNCOND = CSNOW (SNDENS) ENDIF IF ( ( SNOWNG ) .OR. ( FRZGRA ) ) THEN C ---------- FROZEN GROUND MODIFICATION ---------------------- C SNOW DENSITY & CONDUCTIVITY ESTIMATION C C ...1ST: NEW SNOWFALL AMOUNT SN_NEW = PRCP * DT * 0.001 C ...2ND: UPDATE SNOW DENSITY USING C OLD AND NEW SNOW CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS) C C ** ....CSNOW IS FUNCTION TO CALCULATE SNOW CONDUCTIVITY SNCOND = CSNOW (SNDENS) C C ---------------------------------------------------------------- SNEQV = SNEQV + PRCP * DT * 0.001 PRCP1 = 0.0 ELSE PRCP1 = PRCP ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CHECK: IF SNEQV < 1.E-6 M THEN SNOW ZEROED AND RESIDUAL IS SNOWMELT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(SNEQV.LT.1.E-6) THEN SNMAX=SNEQV SNEQV = 0.0 SNOWH = 0.0 SNOWNG = .FALSE. ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALC GRND HEAT FLUX FOR SNOW PACK CASE (WILL EFFECTIVELY C CALC THE HEAT FLUX BELOW THE SNOW PACK). REMEMBER, IF C T1 > 273.16 K, THERE WILL BE SNOW MELT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( SNEQV .GT. 0.0 ) THEN T1MIN = MIN ( T1, TFREEZ ) C ---- FROZEN GROUND MODIFICATION USING COMPUTED SNOW DENSITY --- S = SNCOND * (T1MIN - STC(1) ) /SNOWH C ------------------------------------------------------------------ CC........S = 0.35 * ( T1MIN - STC(1) ) / ( 10.0 * SNEQV ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C BOUND THE HEAT FLUX MAGNITUDE UNDER SNOW PACK TO 100 W/M-2 C SET THE MINIMUM FLUX WITHIN SNOW TO -200 W M-2 C FOR NUMERICAL REASONS RELATED TO VERY SHALLOW SNOW DEPTH 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,QUARTZ,SMCMAX,SH2O ) C --------- FROZEN GROUND VERSION --------------------------- C SOIL THERMAL DIFFUSIVITY CORRECTION C SICE = SMC(1) - SH2O(1) C__________________________________________________________________ C This IF-statement is no longer used if Peters-Lidard thermal C conductivity is used for all conditions (frozen,unfrozen): C IF (SICE .NE. 0.0) CALL FRZCND ( DF1, SICE ) C ---------------------------------------------------------------- S = DF1 * ( STC(1) - T1 ) / ( 0.5 * ZSOIL(1) ) ENDIF C C MODIFY ALBEDO IF SNOWCOVER C IF ( (SNEQV .EQ. 0.0) .OR. (ALB .GE. SNOALB) ) THEN ALBEDO = ALB ELSE C SNOW EFFECT IS SNOWDEPTH DEPENDENT IF (SNEQV .LT. SNUP) THEN RSNOW = SNEQV/SNUP SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP)) ELSE SNOFAC = 1.0 ENDIF ALBEDO = ALB + (1.0-SHDFAC)*SNOFAC*(SNOALB-ALB) ENDIF C CALCULATE DOWNWARD RADIATION F = SOLDN*(1.0-ALBEDO) + LWDN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP) C (AND OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE/ FOR C LATER CALCULATIONS) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL PENMAN ( SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,F,T24,S,Q2, & Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,DQSDT2) IF(SATURATED) 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 C ---- FROZEN GROUND EXTENSION: SMC WAS REPLACED BY SH2O ----- C CALL CANRES(SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, & SMCWLT,SMCREF,RCMIN,RC,PC,NROOT,Q2SAT,DQSDT2, & TOPT,RSMAX,RGL,HS,XLAI) 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 ( SNEQV .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, + SH2O, SLOPE, KDT, FRZX, PSISAT, & ZSOIL, DKSAT, DWSAT, TBOT,RUNOFF1,RUNOFF2, & RUNOFF3, EDIR1, EC1, ETT1,NROOT,ICE,RTDIS,QUARTZ) 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, SALP, SNEQV, + SNOWH, SH2O, SLOPE, KDT, FRZX, PSISAT, SNUP, & ZSOIL, DWSAT, DKSAT, TBOT, SHDFAC,RUNOFF1, & RUNOFF2,RUNOFF3,EDIR1,EC1,ETT1,NROOT,SNMAX,ICE, & RTDIS,QUARTZ) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C MODIFY ALBEDO WITH UPDATED SNOW (FOR RETURN TO PARENT MODEL RADIATION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( (SNEQV .EQ. 0.0) .OR. (ALB .GE. SNOALB) ) THEN ALBEDO = ALB ELSE C SNOW EFFECT IS SNOWDEPTH DEPENDENT IF (SNEQV .LT. SNUP) THEN RSNOW = SNEQV/SNUP SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP)) ELSE SNOFAC = 1.0 ENDIF ALBEDO = ALB + (1.0-SHDFAC)*SNOFAC*(SNOALB-ALB) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PREPARE SENSIBLE HEAT (H) FOR RETURN TO PARENT MODEL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C T2V = sfc air virtual temperature C RHO = sfc air density (kg m-3) C..... RHO = SFCPRS/(R * T2V) C.....CHKFF = (CH * CP * RHO ) C..... H = -(CHKFF ) * ( TH2 - T1 ) C H = -(CH * CP * SFCPRS)/(R * T2V) * ( TH2 - T1 ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT UNITS AND/OR SIGN OF TOTAL EVAP (ETA), POTENTIAL EVAP (ETP), C SUBSURFACE HEAT FLUX (S), AND RUNOFFS FOR C WHAT PARENT MODEL EXPECTS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1 C AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW C RUNOFF3 = RUNOFF3/DT RUNOFF2 = RUNOFF2+RUNOFF3 C C TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE C SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION SOILM=-1.0*SMC(1)*ZSOIL(1) DO K = 2, NSOIL SOILM=SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K)) END DO 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)) END DO SOILW=SOILWW/SOILWM C RETURN END C SUBROUTINE SHFLX(S,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,TBOT, + SMCWLT, PSISAT, SH2O, B,F1,DF1,ICE,QUARTZ ) IMPLICIT NONE C ------------ FROZEN GROUND VERSION -------------------------- C NEW STATES ADDED: SH2O & PARAMETER SMCWLT & PSISAT C C ------------------------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SOIL HEAT FLUX. SOIL TEMPERATURE CONTENT CC ======= (AN ARRAY CONTAINING SOIL TEMPERATURES) IS ALSO UPDATED. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOLD PARAMETER ( NSOLD = 4 ) INTEGER NSOIL REAL B REAL DF1 REAL DT REAL F1 REAL RHSTS ( NSOLD ) REAL S REAL SMC ( NSOIL ) REAL SH2O ( NSOIL ) REAL SMCMAX REAL SMCWLT REAL STC ( NSOIL ) REAL STCOUT ( NSOLD ) REAL T1 REAL TBOT REAL YY REAL ZSOIL ( NSOIL ) REAL ZZ1 REAL PSISAT, QUARTZ INTEGER ICE, IFRZ, I 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) CALL HSTEP ( STC, STC, RHSTS,DT,NSOIL ) ELSE CALL HRT(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, C ------------- FROZEN GROUND VERSION ----------------------- C + SMCWLT, PSISAT, SH2O, DT, C ----------------------------------------------------------------- + B,F1,DF1,QUARTZ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C HSTEP ROUTINE CALCS/UPDATES SOIL TEMPS BASED ON RHSTS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL HSTEP ( STCOUT,STC,RHSTS,DT,NSOIL ) C ---------------- FROZEN GROUND VERSION -------------------- C DOUBLE CALLING HRT & HSTEP TO REDUCE NOISE FLUCTUATION, THE SAME C AS SOIL MOISTURE CALCULATION. DOUBLE CALL WILL BE ONLY IF THERE C IS A ICE ON AT LEAST ONE SOIL LAYER C IFRZ = 0 DO I = 1,NSOIL IF( (SMC(I) - SH2O(I) ) .GT. 0.0 ) IFRZ = 1 END DO IF ( IFRZ .EQ. 1 ) THEN DO I = 1,NSOIL STC(I) = 0.5 * ( STCOUT(I) + STC(I) ) END DO CALL HRT(RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, C ------------- FROZEN GROUND VERSION ----------------------- C + SMCWLT, PSISAT, SH2O, DT, C ----------------------------------------------------------------- + B,F1,DF1,QUARTZ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C HSTEP ROUTINE CALCS/UPDATES SOIL TEMPS BASED ON RHSTS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL HSTEP ( STC,STC,RHSTS,DT,NSOIL ) ELSE DO I = 1,NSOIL STC(I) = STCOUT(I) END DO ENDIF C ------------------------------------------------------------------ ENDIF 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, & SH2O, SLOPE, KDT, FRZFACT, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT,SMCREF,SHDFAC,CMCMAX, & SMCDRY,CFACTR, RUNOFF1,RUNOFF2, RUNOFF3, EDIR1, EC1, & ETT1, SFCTMP,Q2,NROOT,RTDIS) IMPLICIT NONE C ------------ FROZEN GROUND VERSION -------------------------- C NEW STATES ADDED: SH2O, AND FROZEN GROUD CORRECTION FACTOR, FRZFACT C AND PARAMETER SLOPE C 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 INTEGER NSOLD 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 KDT REAL PC REAL PCPDRP REAL PRCP1 REAL RHSCT REAL RHSTT ( NSOLD ) REAL RIB REAL RTDIS (NSOIL) REAL RUNOF REAL RUNOFF,RUNOXX3 REAL SHDFAC REAL SMC ( NSOIL ) C --------------- FROZEN GROUND VERSION --------------------- REAL SH2O ( NSOIL ) REAL SICE ( NSOLD ) REAL SH2OA ( NSOLD) REAL SH2OFG ( NSOLD ) C ------------------------------------------------------------------- REAL SMCDRY REAL SMCMAX REAL SMCREF REAL SMCWLT REAL TRHSCT REAL ZSOIL ( NSOIL ) REAL SLOPE, FRZFACT, RUNOFF1, RUNOFF2, RUNOFF3, EDIR1, EC1 REAL ETT1, SFCTMP, Q2, DUMMY, CMC2MS, DEVAP INTEGER NROOT, I CMIC$ TASKCOMMON RITE COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOF, & DEW,RIB,RUNOXX3 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 IF ( ETP1 .GT. 0.0 ) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C -------------- FROZEN GROUND VERSION --------------------- C SMC STATES WERE REPLACED BY SH2O STATES C EDIR = DEVAP ( ETP1, SH2O(1), ZSOIL(1), SHDFAC, SMCMAX, C ------------------------------------------------------------------ C 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 C -------------- FROZEN GROUND VERSION --------------------- C SMC STATES WERE REPLACED BY SH2O STATES C CALL TRANSP ( ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, C ------------------------------------------------------------------ C CALL TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS) 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 C --------------- FROZEN GROUND VERSION -------------------- C STORE ISE CONTENT AT EACH LAYER BEFORE CALLING SRT & SSTEP C DO I = 1,NSOIL SICE(I) = SMC(I) - SH2O(I) END DO C ------------------------------------------------------------------ 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 C print*,'SMFLX, BEFORE srt, sstep' C print*,'smc=' , smc C print*,'sh2o=', sh2o IF ( PCPDRP .GT. 0.0 ) THEN C --------------- FROZEN GROUND VERSION --------------------- C SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. C SH2O & SICE STATES INCLUDED IN SSTEP SUBR. C FROZEN GROUND CORRECTION FACTOR, FRZFACT, ADDED C ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER C CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, & DWSAT,DKSAT,SMCMAX, B, RUNOFF1, + RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT, SICE) C print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER srt - 1' C print*,'smc=' , smc C print*,'sh2o=', sh2o CALL SSTEP ( SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, & CMCMAX, RUNOFF3, ZSOIL, SMC, SICE ) C print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER sstep - 1' C print*,'smc=' , smc C print*,'sh2o=', sh2o DO 30 K = 1, NSOIL SH2OA(K) = ( SH2O(K) + SH2OFG(K) ) * 0.5 30 CONTINUE CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, & DWSAT,DKSAT,SMCMAX, B, RUNOFF1, + RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT, SICE) C print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER srt - 2' C print*,'smc=' , smc C print*,'sh2o=', sh2o CALL SSTEP ( SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, & CMCMAX, RUNOFF3, ZSOIL,SMC,SICE) C print*,'SMFLX, (PCPDRP .GT. 0.0) AFTER sstep - 2' C print*,'smc=' , smc C print*,'sh2o=', sh2o ELSE CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, & DWSAT,DKSAT,SMCMAX, B, RUNOFF1, + RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT, SICE) C print*,'SMFLX, AFTER srt ' C print*,'smc=' , smc C print*,'sh2o=', sh2o CALL SSTEP ( SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, & CMCMAX, RUNOFF3, ZSOIL,SMC,SICE) C print*,'SMFLX, AFTER sstep ' C print*,'smc=' , smc C print*,'sh2o=', sh2o C -------------------------------------------------------------------- C ***** BELOW IS OLD CODE ************************************* C CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMC,NSOIL,PCPDRP,ZSOIL, C & DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT) C CALL SSTEP ( SMCFG,SMC,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, C & CMCMAX, RUNOFF3, ZSOIL) C C DO 30 K = 1, NSOIL C SMCA(K) = ( SMC(K) + SMCFG(K) ) * 0.5 C 30 CONTINUE C C CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMCA,NSOIL,PCPDRP,ZSOIL, C & DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT) C CALL SSTEP ( SMC,SMC,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, C & CMCMAX, RUNOFF3, ZSOIL) C ELSE C C CALL SRT ( RHSTT,RUNOFF,EDIR,ET,SMC,SMC,NSOIL,PCPDRP,ZSOIL, C & DWSAT, DKSAT, SMCMAX, B, RUNOFF1, RUNOFF2, DT, SMCWLT) C CALL SSTEP ( SMC,SMC,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, C & CMCMAX, RUNOFF3, ZSOIL) C ******************************************************************** ENDIF RUNOF = RUNOFF RETURN END FUNCTION SNKSRC ( TUP,TM,TDN, SMC, SH2O, HCPCT, ZSOIL, + SMCMAX, PSISAT, B, F1, DT, K, QUARTZ) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION CC ======= EQUATION. (SH2O) IS AVAILABLE LIQUED WATER. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER K REAL B REAL DF REAL DFH2O REAL DFICE REAL DH2O REAL DT REAL DZ REAL DZH REAL F1 REAL FREE REAL FRH2O REAL HLICE REAL PSISAT REAL QDN REAL QTOT REAL QUARTZ REAL QUP REAL SH2O REAL SICE REAL SMC REAL SMCMAX REAL SNKSRC REAL T0 REAL TAVG REAL TDN REAL TM REAL TUP REAL TZ REAL X0 REAL XDN REAL XH2O REAL XUP REAL HCPCT REAL ZSOIL (4) PARAMETER (HLICE=3.3350 E5) PARAMETER (DH2O =1.0000 E3) PARAMETER ( T0 =2.7316 E2) C ************************************************************************** C *** OPTION TO RUN NON-FROZEN VERSION: UNCOMMENT NEXT TWO STATEMENTS **** C ******************* WARNING !!!!! ************************** C *** TO RUN FROZEN GROUND VERSION, NEXT TWO STATEMENTS SHOULD BE **** C *** COMMENTED **** C CC SNKSRC = 0.0 CC GOTO 77 C ************************************************************************** SICE=SMC-SH2O IF(K.EQ.1) THEN DZ=-ZSOIL(1) ELSE DZ=ZSOIL(K-1)-ZSOIL(K) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE DIFFUSIVITIES OF THAWED AND FROZEN GROUND CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C WARNING !!!!!, by Pablo Grunmann, 10/98 C ATTENTION TO THE CALLING ARGUMENT: C C 1. USE SMC instead of SH2O to use VKOREN's Frozen DF (call C and IF statement commented out with C VK) C 2. USE SH2O in calling argument list for Petrs-Lidard (call C commented out with C P-L) C !!! !!!! C P-L CALL TDFCND ( DFH2O, SMC, B, F1,QUARTZ,SMCMAX,SH2O) C VK CALL TDFCND ( DFH2O, SMC, B, F1,QUARTZ,SMCMAX, SMC) CALL TDFCND ( DFH2O, SMC, B, F1,QUARTZ,SMCMAX,SH2O) DFICE = DFH2O C Frozen conductivity for original VK approach C VK IF ( SICE .GT. 0.0 ) CALL FRZCND ( DFICE, SICE ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C HEAT FLUX FROM THE TOP BOUNDARY OF LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(TUP .LE. T0) THEN DF=DFICE ELSE DF=DFH2O ENDIF QUP=-DF*(T0-TUP)/(0.5*DZ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C HEAT FLUX FROM THE BOTTOM BONDARY OF LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(TDN .LE. T0) THEN DF=DFICE ELSE DF=DFH2O ENDIF QDN=-DF*(T0-TDN)/(0.5*DZ) QTOT=QUP+QDN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE POTENTIAL REDUCTION OF LIQUED WATER CONTENT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC XH2O=QTOT*DT/(DH2O*HLICE*DZ) + SH2O CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ESTIMATE UNFROZEN WATER AT TEMPERATURE TAVG, C AND CHECK IF CALCULATED WATER CONTENT IS REASONABLE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C **** NEW CALCULATION OF AVERAGE TEMPERATURE (TAVG) ********** C **** IN FREEZING/THAWING LAYER USING UP, DOWN, AND MIDDLE *** C **** LAYER TEMPERATURES (TUP, TDN, TM) ********** DZH=DZ*0.5 IF(TUP .LT. 0.) THEN IF(TM .LT. 0.) THEN IF(TDN .LT. 0.) THEN C1 ****** TUP, TM, TDN < 0. CASE ************ TAVG=(TUP+2.0*TM+TDN)/4. C C111 **** END TDN_111 THEN ELSE C2 ****** TUP & TM < 0., TDN > 0. CASE ***** X0=(T0-TM)*DZH/(TDN-TM) TAVG=0.5*(TUP*DZH+TM*(DZH+X0)+T0*(2.*DZH-X0))/DZ C C112 **** END TDN_111 ELSE **** ENDIF C C11 *********** END TM_1 THEN *********** ELSE IF(TDN .LT. 0.) THEN C3 ******* TUP < 0. TM > 0. TDN < 0. CASE **** XUP=(T0-TUP)*DZH/(TM-TUP) XDN=DZH-(T0-TZ)*DZH/(TDN-TM) TAVG=0.5*(TUP*XUP+T0*(2.*DZ-XUP-XDN)+TDN*XDN)/DZ C C121 **** END TDN_121 THEN ***** ELSE C4 ****** TUP < 0 TM > 0 TDN > 0 CASE ****** XUP=(T0-TUP)*DZH/(TM-TUP) TAVG=0.5*(TUP*XUP+T0*(2.*DZ-XUP))/DZ C C122 *** END TDN_121 ELSE *** ENDIF C C12 *********** END TM_1 ELSE ********** ENDIF C1 ******************** END TUP THEN *********** ELSE IF(TM .LT. 0.) THEN IF(TDN .LT. 0.) THEN C5 ***** TUP > 0 TM < 0 TDN < 0 CASE ********* XUP=DZH-(T0-TUP)*DZH/(TM-TUP) TAVG=0.5*(T0*(DZ-XUP)+TM*(DZH+XUP)+TDN*DZH)/DZ C C211 **** END TDN_211 THEN ***** ELSE C6 ***** TUP > 0 TM < 0 TDN > 0 CASE ********* XUP=DZH-(T0-TUP)*DZH/(TM-TUP) XDN=(T0-TM)*DZH/(TDN-TM) TAVG=0.5*(T0*(2.*DZ-XUP+XDN)+TM*(XUP+XDN))/DZ C C212 *** END TDN_211 ELSE *** ENDIF C C21 ************ END TM_2 THEN ************ ELSE IF(TDN .LT. 0.) THEN C7 ***** TUP > 0 TM > 0 TDN < 0 CASE ******** XDN=DZH-(T0-TM)*DZH/(TDN-TM) TAVG=(T0*(DZ-XDN)+0.5*(T0+TDN)*XDN)/DZ C C221 **** END TDN_221 THEN ***** ELSE C8 ***** TUP > 0 TM > 0 TDN > 0 CASE ******** TAVG=(TUP+2.*TM+TDN)/4. C C222 *** END TDN_221 ELSE *** ENDIF C C22 ************ END TM_2 ELSE ************ ENDIF C ********************* END TUP ELSE ************** ENDIF C 777 100 FREE=FRH2O(TAVG, SMC, SH2O, SMCMAX, B, PSISAT ) IF ( XH2O .LT. SH2O .AND. XH2O .LT. FREE) THEN IF ( FREE .GT. SH2O ) THEN XH2O = SH2O ELSE XH2O = FREE ENDIF ENDIF IF ( XH2O .GT. SH2O .AND. XH2O .GT. FREE ) THEN IF ( FREE .LT. SH2O ) THEN XH2O = SH2O ELSE XH2O = FREE ENDIF ENDIF IF(XH2O .LT. 0. ) XH2O=0. 10 IF(XH2O .GT. SMC) XH2O=SMC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE SINK/SOURCE TERM AND REPLACE PREVIOUS WATER CONTENT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SNKSRC=-DH2O*HLICE*DZ*(XH2O-SH2O)/DT IF(XH2O .LT. 0.) THEN XH2O=SH2O ENDIF SH2O=XH2O 77 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, SALP, ESD, + SNOWH, SH2O, SLOPE, KDT, FRZFACT, PSISAT,SNUP, & ZSOIL, DWSAT, DKSAT, TBOT, SHDFAC, RUNOFF1, & RUNOFF2,RUNOFF3,EDIR1,EC1,ETT1,NROOT,SNMAX,ICE, & RTDIS,QUARTZ) IMPLICIT NONE 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 SNOWH 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 KDT REAL LSUBF REAL LSUBC REAL LSUBS REAL PC REAL PRCP REAL PRCP1 REAL Q1 REAL Q2 REAL RCH REAL RIB REAL RR REAL RTDIS (NSOIL) REAL RUNOFF REAL S REAL S1 REAL SFCTMP REAL SHDFAC REAL SIGMA REAL SMC ( NSOIL ) REAL SH2O ( 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 REAL TFREEZ, SALP, SFCPRS, SLOPE, FRZFACT, PSISAT, SNUP REAL RUNOFF1, RUNOFF2, RUNOFF3,RUNOXX3 REAL EDIR1, EC1, ETT1, QUARTZ REAL SNDENS, SNCOND, RSNOW, SNCOVER, QSAT, ETP3, SEH, T14 REAL SICE1, CSNOW INTEGER NROOT, ICE CMIC$ TASKCOMMON RITE COMMON/RITE/ BETA,DRIP,EC,EDIR,ETT,FLX1,FLX2,FLX3,RUNOFF, & DEW,RIB,RUNOXX3 PARAMETER(CP=1004.5,CPH2O=4.218E+3,CPICE=2.106E+3, & LSUBF=3.335E+5,LSUBC=2.5E+6,LSUBS=2.83E+6,SIGMA=5.67E-8) C ------------ FROZEN GROUND VERSION ------------------------ C C PARAMETER ( TFREEZ = 273.16) C------------------------------------------------------------------- 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. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC T12A = ((F - SIGMA * T24) / RCH+TH2-SFCTMP-BETA*EPSCA) / RR C --------------- FROZEN GROUND VERSION ----------------------- C VARIABLE SNOW DENSITY & CONDUCTIVITY SNDENS = ESD / SNOWH SNCOND = CSNOW ( SNDENS ) T12B = SNCOND * STC(1) / ( SNOWH * RR * RCH ) DENOM = 1.0 + SNCOND / ( SNOWH * RR * RCH ) C ------------------------------------------------------------------- C T12B = 0.35 * STC(1) / ( 10.0 * ESD * RR * RCH ) C 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. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC.....S = 0.35 * ( T12 - STC(1) ) / ( 10.0 * ESD ) C --------------- FROZEN GROUND VERSION ----------------------- C VARIABLE SNOW DENSITY & CONDUCTIVITY C S = SNCOND * ( T12 - STC(1) ) / ( SNOWH ) C ------------------------------------------------------------------- S = MIN ( MAX ( S,-100.0 ), 100.0 ) IF ( T12 .LE. TFREEZ ) 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 = MAX ( 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 = MIN ( ESD / DT, EX ) C ESD = MAX ( ESD - SNMAX, 0.0 ) C T1 = MAX ( 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) C T1 = 273.16 C ------------ FROZEN GROUND VERSION -------------------------- C ADJUSTMENT T1 TO ACCOUNT FOR SNOW PATCHES C IF ( SNUP .GT. 0.0 .AND. ESD .LT. SNUP ) THEN RSNOW = ESD / SNUP SNCOVER = 1.-( EXP(-SALP*RSNOW)-RSNOW*EXP(-SALP)) C SNCOVER = MIN ( ESD / SNUP, 1.) ELSE SNCOVER = 1. ENDIF T1 = TFREEZ * SNCOVER + T12 * ( 1.0 - SNCOVER ) C ------------------------------------------------------------------- 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 C ------------ FROZEN GROUND VERSION -------------------------- C S = SNCOND * ( T1 - STC(1) ) / ( SNOWH ) C ------------------------------------------------------------------- C 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 C ------------ FROZEN GROUND VERSION -------------------------- C SNOWMELT REDUCTION DEPENDING ON SNOW COVER C IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION C IF ( SNCOVER .GT. 0.05) EX = EX * SNCOVER C ------------------------------------------------------------------- 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 C ------------ FROZEN GROUND VERSION -------------------------- C (+) NEW STATES ADDED: SH2O, AND FROZEN GROUND CORRECTION FACTOR C CALL SMFLX ( ETA1,SMC,NSOIL,CMC,ETP1,DT,PRCP1,ZSOIL, + SH2O, SLOPE, KDT, FRZFACT, & SMCMAX,B,PC,SMCWLT,DKSAT,DWSAT, & SMCREF,SHDFAC,CMCMAX,SMCDRY,CFACTR,RUNOFF1,RUNOFF2, & RUNOFF3, EDIR1, EC1, ETT1,SFCTMP,Q2,NROOT,RTDIS) 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,QUARTZ,SMCMAX,SH2O ) C ------------ FROZEN GROUND VERSION -------------------------- C SOIL CONDUCTIVITY ADJUSTMENT C SICE1 = SMC(1) - SH2O(1) C__________________________________________________________________ C This IF-statement is no longer used if Peters-Lidard thermal C conductivity is used for all conditions (frozen,unfrozen): C IF ( SICE1 .GT. 0.0 ) CALL FRZCND ( DF1, SICE1 ) C ------------------------------------------------------------------- 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 C ------------ FROZEN GROUND VERSION -------------------------- C SNOW DEPTH & DENSITY ADJUSTMENT BASED ON SNOW COMPACTION C CALL SNOWPACK ( ESD, DT, SNOWH, SNDENS, T1, YY ) C ------------------------------------------------------------------- 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 C ------------ FROZEN GROUND VERSION -------------------------- C (+) NEW STATES ADDED: SH2O & PARAMETER SMCWLT & PSISAT C CALL SHFLX ( S1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL,TBOT, + SMCWLT, PSISAT, SH2O, & B,F1,DF1,ICE, & QUARTZ ) 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 SNOWPACK ( W,DTS,HC,DS,TSNOW,TSOIL ) IMPLICIT NONE C ************************************************************** C ** SUBROUTINE TO CALCULATE SNOW COMPACTION *** C ** EQUATION OF INCREASING OF SNOW DENSITY WAS OBTAINED *** C AS AN APPROXIMATE SOLUTION OF E. ANDERSON DIFFERENTIAL *** C EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, *** C BY VICTOR KOREN 03/25/95 *** C ************************************************************** C ************************************************************** C W IS A WATER EQUIVALENT OF SNOW, IN M *** C DTS IS A TIME STEP, IN SEC *** C HC IS A SNOW DEPTH, IN M *** C DS IS A SNOW DENSITY, IN G/CM3 *** C TSNOW IS A SNOW SURFACE TEMPERATURE, K *** C TSOIL IS A SOIL SURFACE TEMPERATURE, K *** C SUBROUTINE WILL RETURN NEW VALUES OF H AND DS *** C ************************************************************** REAL C1, C2, HC, W, DTS, DS, TSNOW, TSOIL, H, WX REAL DT, TSNOWX, TSOILX, TAVG, B, DSX, DW INTEGER IPOL, J PARAMETER (C1=0.01, C2=21.0) REAL PEXP C ** CONVERSION INTO SIMULATION UNITS ************************* H=HC*100. WX=W*100. DT=DTS/3600. TSNOWX=TSNOW-273.16 TSOILX=TSOIL-273.16 C ** CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK *** TAVG=0.5*(TSNOWX+TSOILX) C ** CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION C DS=DS0*(EXP(B*W)-1.)/(B*W) C B=DT*C1*EXP(0.08*TAVG-C2*DS0) C ** C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR)) C ** C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G IF(WX .GT. 1.E-2) THEN B=DT*C1*EXP(0.08*TAVG-C2*DS) C C 9/03/98 This line was changed to increase an accuracy of exponent C DSX=DS*((DEXP(B*WX)-1.)/(B*WX)) C-------------------------------------------------------------------- C The above line was causing numerical difficulties because of the C denominator, that can go to zero (the function C has a well defined limit, equal to 1. when B*WX C goes to zero, though) C C Suggestion to substitute the line above by C Pablo Grunmann 03-09-98 C-------------------------------------------------------------------- C Remember to declare PEXP C Precision C desired: ipol greater than 9 only makes a difference on double C precision. C ipol=9, for rel.error <~ 1.6 e-6 % (8 significant digits) C ipol=8, for rel.error <~ 1.8 e-5 % (7 significant digits) C ipol=7, for rel.error <~ 1.8 e-4 % ... ipol = 9 PEXP = 0. do j = ipol,1,-1 PEXP = (1. + PEXP)*B*WX/real(j+1) end do PEXP = PEXP + 1. C C This result imitates the original formula C DSX=DS*(PEXP) C-------------------------------------------------------------------- C End of substitution C-------------------------------------------------------------------- IF(DSX .GT. 0.40) DSX=0.40 IF(DSX .LT. 0.05) DSX=DS DS=DSX H=WX/DS C ** CORRECTION OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER C ** DURING SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED C ** IN SNOW PER DAY DURING SNOWMELT TILL SNOW DENSITY 0.40 IF((TSNOWX .GE. 0.) .AND. (H .NE. 0.)) THEN DW=0.13*DT/24. DS=DS*(1.-DW)+DW IF(DS .GT. 0.40) DS=0.40 H=WX/DS ENDIF ELSE IF(WX .LE. 0.) THEN WX=0. H=0. DS=0. ELSE H=WX/DS ENDIF ENDIF C ** REPLACE PREVIOUS LIQUID WATER AND CHANGE DEPTH TO M C 7 HC=H*0.01 RETURN END SUBROUTINE SNOW_NEW ( T,P,HC,DS ) IMPLICIT NONE C CALCULATING SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL C T - AIR TEMPERATURE, K C P - NEW SNOWFALL, M C HC - SNOW DEPTH, M C DS - SNOW DENSITY C NEW VALUES OF SNOW DEPTH & DENSITY WILL BE RETURNED REAL HC REAL T REAL P REAL DS REAL H REAL PX REAL TX REAL DS0 REAL HNEW C ** CONVERSION INTO SIMULATION UNITS ************************* H=HC*100. PX=P*100. TX=T-273.16 C ---------------------------------------------------------------------- C *** CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE ** C *** EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED C *** AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE, C *** VEMADOLEN, SWEDEN, 1980, 172-177PP. C----------------------------------------------------------------------- IF(TX .LE. -15.) THEN DS0=0.05 ELSE DS0=0.05+0.0017*(TX+15.)**1.5 ENDIF C ** ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL HNEW=PX/DS0 IF ((H+HNEW) .GT. 0.0) THEN DS=(H*DS+HNEW*DS0)/(H+HNEW) H=H+HNEW ELSE DS=0.0 H=0.0 END IF HC=H*0.01 RETURN END SUBROUTINE SRT (RHSTT,RUNOFF,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, & ZSOIL,DWSAT,DKSAT,SMCMAX,B, RUNOFF1, + RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOLD PARAMETER ( NSOLD = 4 ) INTEGER K INTEGER NSOIL INTEGER CVFRZ 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 KDT REAL MXSMC REAL MXSMC2 REAL NUMER REAL PCPDRP REAL PDDUM REAL RHSTT ( NSOIL ) REAL RUNOFF REAL SH2O ( NSOIL ) REAL SH2OA ( NSOIL ) REAL SICE ( NSOIL ) REAL SMCMAX REAL WCND REAL WCND2 REAL WDF REAL WDF2 REAL ZSOIL ( NSOIL ) REAL RUNOFF1, RUNOFF2, DT, SMCWLT, SLOPE, FRZX, DT1 REAL SMCAV, DICE, DD, VAL, DDT, PX, FCR, ACRT, SUM REAL SSTT, SLOPX INTEGER IOHINF, KS, IALP1, J, JJ CMIC$ TASKCOMMON ABCI COMMON /ABCI/ AI, BI, CI C ----------- FROZEN GROUND VERSION ------------------------- C REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF C AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV. C CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. C BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT C CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. C THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}) C PARAMETER ( CVFRZ = 3 ) C ------------------------------------------------------------------ C print*,'in SRT, Declaration -----------------------' C print*,'NSOIL=' , NSOIL C print*,'NSOLD=' , NSOLD 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 DT1 = DT/86400. SMCAV = SMCMAX - SMCWLT DMAX(1)=-ZSOIL(1)*SMCAV C ----------- FROZEN GROUND VERSION ------------------------ C DICE = -ZSOIL(1) * SICE(1) C------------------------------------------------------------------- DMAX(1)=DMAX(1)*(1.0 - (SH2OA(1)+SICE(1)-SMCWLT)/SMCAV) DD=DMAX(1) DO KS=2,NSOIL C ----------- FROZEN GROUND VERSION ------------------------ C DICE = DICE + ( ZSOIL(KS-1) - ZSOIL(KS) ) * SICE(KS) C------------------------------------------------------------------- DMAX(KS)=(ZSOIL(KS-1)-ZSOIL(KS))*SMCAV DMAX(KS)=DMAX(KS)*(1.0 - (SH2OA(KS)+SICE(KS)-SMCWLT)/SMCAV) DD=DD+DMAX(KS) END DO 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 C ----------- FROZEN GROUND VERSION -------------------------- C REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS C FCR = 1. IF ( DICE .GT. 1.E-2) THEN ACRT = CVFRZ * FRZX / DICE SUM = 1. IALP1 = CVFRZ - 1 DO J = 1,IALP1 K = 1 DO JJ = J+1, IALP1 K = K * JJ END DO SUM = SUM + (ACRT ** ( CVFRZ-J)) / FLOAT (K) END DO FCR = 1. - EXP(-ACRT) * SUM END IF INFMAX = INFMAX * FCR C ------------------------------------------------------------------- C ************ CORRECTION OF INFILTRATION LIMITATION ********** C IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE C VALUE OF HYDROLIC CONDUCTIVITY C C MXSMC = MAX ( SH2OA(1), SH2OA(2) ) MXSMC = SH2OA(1) C print*,'SRT, BEFORE WDFCND - 1 ------------------------------' C print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX C print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT CALL WDFCND ( WDF,WCND,MXSMC,SMCMAX,B,DKSAT,DWSAT ) INFMAX = MAX(INFMAX, WCND) INFMAX= MIN(INFMAX,PX) C print*,'SRT, AFTER WDFCND - 1 ------------------------------' C print*,'WDF,WCND=' , WDF,WCND C print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX C print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT C ****************************************************************** C 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 C_________________________________________________________________ C Comments on August 98 changes: C C This approach disregarded the actual soil moisture in C the upper layer if it was less that that of the next, for C hydraulic conductivity and diffussivity purposes (supposedly C to solve problems with the numerical discretization on a very C coarse vertical resolution setting - two layer model). C In our experience, This would actually allow the upper C layer to go very dry at a fast rate while the next layer C keeps unchanged (except for diffusivity effect) C C Pablo J. Grunmann C C MXSMC = MAX( SH2OA(1), SH2OA(2) ) C___________________________________________________________ C Test a less drastic version of MXSMC: !!!!!!!!!!!!!!!!!!!! C MXSMC = MAX( SH2OA(1), (SH2OA(2)+SH2OA(1))/2.0 ) C___________________________________________________________ C_________________________________________________________________ C New approach: dif/conductivity parameters on upper layer C depend on upper layer only MXSMC = SH2OA(1) C print*,'SRT, BEFORE WDFCND - 2' C print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX C print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT CALL WDFCND ( WDF,WCND,MXSMC,SMCMAX,B,DKSAT,DWSAT ) C print*,'SRT, AFTER WDFCND - 2' C print*,'WDF,WCND=' , WDF,WCND C print*,'MXSMC,SMCMAX=' , MXSMC,SMCMAX C print*,'B,DKSAT,DWSAT=' , 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 = ( SH2O(1) - SH2O(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 SLOPX = 1. 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) ) CC MXSMC2 = MAX ( SH2OA(K), SH2OA(K+1) ) C ************ PILPS TEST ************* MXSMC2 = SH2OA(K) C print*,'SRT, BEFORE WDFCND - 3' C print*,'MXSMC2,SMCMAX=' , MXSMC2,SMCMAX C print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT C print*,'K=' , K CALL WDFCND ( WDF2,WCND2,MXSMC2,SMCMAX,B,DKSAT,DWSAT ) C print*,'SRT, AFTER WDFCND - 3' C print*,'WDF2,WCND2=' , WDF2,WCND2 C print*,'MXSMC2,SMCMAX=' , MXSMC2,SMCMAX C print*,'B,DKSAT,DWSAT=' , 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 = ( SH2O(K) - SH2O(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 C SLOPE OF BOTTOM LAYER IS INTRODUCED ************ C SLOPX = SLOPE C-------------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC C CONDUCTIVITY FOR THIS LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C print*,'SRT, BEFORE WDFCND - 4' C print*,'SH2OA(NSOIL),SMCMAX=' , SH2OA(NSOIL),SMCMAX C print*,'B,DKSAT,DWSAT=' , B,DKSAT,DWSAT C print*,'K=' , K CALL WDFCND( WDF2,WCND2,SH2OA(NSOIL),SMCMAX,B,DKSAT,DWSAT ) C print*,'SRT, AFTER WDFCND - 4' C print*,'WDF2,WCND2=' , WDF2,WCND2 C print*,'SH2OA(NSOIL),SMCMAX=' , SH2OA(NSOIL),SMCMAX C print*,'B,DKSAT,DWSAT=' , 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) + SLOPX * 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 = SLOPX * WCND2 ENDIF IF ( K .NE. NSOIL ) THEN WDF = WDF2 WCND = WCND2 DSMDZ = DSMDZ2 DDZ = DDZ2 END IF 10 CONTINUE C print*,'SRT, final Runoff' C print*,'RUNOFF1=' , RUNOFF1 C print*,'RUNOFF2=' , RUNOFF2 RETURN END SUBROUTINE SSTEP ( SH2OOUT, SH2OIN, CMC, RHSTT, RHSCT, DT, & NSOIL, SMCMAX, CMCMAX, RUNOFF3, ZSOIL,SMC,SICE ) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE/UPDATE THE SOIL MOISTURE CONTENT VALUES CC ======= AND THE CANOPY MOISTURE CONTENT VALUES. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER NSOLD, I PARAMETER ( NSOLD = 4 ) INTEGER K INTEGER KK11 INTEGER NSOIL REAL AI ( NSOLD ) REAL BI ( NSOLD ) REAL CI ( NSOLD ) REAL CMC REAL CMCMAX REAL DT REAL RHSCT REAL RHSTT ( NSOIL ) C -------- FROZEN GROUND CORRECTIONS ------------------ REAL SH2OIN ( NSOIL ) REAL SH2OOUT ( NSOIL ) REAL SICE ( NSOIL ) REAL SMC ( NSOIL ) C -------------------------------------------------------------- REAL SMCMAX REAL ZSOIL(NSOIL) REAL RUNOFF3, RUNOFS, WPLUS, DDZ, STOT, WFREE, DPLUS cmic$ TASKCOMMON 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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 ******* C *********** WATER BALANCE CALCULATIONS WERE WRONG ******** C IT WAS CORRECTED DURING DURING INCLUDING FROZEN GROUND **** C RUNOFS = 0.0 WPLUS = 0.0 RUNOFF3 = 0. DDZ = - ZSOIL(1) DO 20 K = 1 , NSOIL IF ( K .NE. 1 ) DDZ = ZSOIL(K - 1) - ZSOIL(K) SH2OOUT(K) = SH2OIN(K) + CI(K) + WPLUS / DDZ C print*,'IN sstep' C print*,'SH2OOUT=', SH2OOUT STOT = SH2OOUT(K) + SICE(K) IF ( STOT .GT. SMCMAX ) THEN IF ( K .EQ. 1 ) THEN DDZ = -ZSOIL(1) ELSE KK11 = K - 1 DDZ = -ZSOIL(K) + ZSOIL(KK11) END IF WPLUS = ( STOT - SMCMAX ) * DDZ ELSE WPLUS = 0. END IF SMC(K) = MAX ( MIN( STOT, SMCMAX ), 0.02 ) SH2OOUT(K) = MAX ( (SMC(K) - SICE(K)), 0.0 ) 20 CONTINUE C *** V. KOREN 9/01/98 ****** C WATER BALANCE CHECKING UPWARD C IF(WPLUS .GT. 0.) THEN DO I=NSOIL-1,1,-1 IF(I .EQ. 1) THEN DDZ=-ZSOIL(1) ELSE DDZ=-ZSOIL(I)+ZSOIL(I-1) ENDIF WFREE=(SMCMAX-SH2OOUT(I)-SICE(I))*DDZ DPLUS=WFREE-WPLUS IF(DPLUS .GE. 0.) THEN SH2OOUT(I)=SH2OOUT(I)+WPLUS/DDZ SMC(I)=SH2OOUT(I)+SICE(I) WPLUS=0. ELSE SH2OOUT(I)=SH2OOUT(I)+WFREE/DDZ SMC(I)=SH2OOUT(I)+SICE(I) WPLUS=-DPLUS ENDIF END DO 30 RUNOFF3=WPLUS ENDIF 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 (on original source) C IF (CMC .LT. 1.E-20) CMC=0.0 CMC = MIN(CMC,CMCMAX) RETURN END FUNCTION TBND (TU, TB, ZSOIL, ZBOT, K, NSOIL) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER CC ======= BY INTERPOLATION OF THE MIDDLE LAYER TEMPERATURES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL TBND, T0, TU, TB, ZBOT, ZUP, ZB INTEGER NSOIL, K PARAMETER (T0=273.16) REAL ZSOIL (4) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(K .EQ. 1) THEN ZUP=0. ELSE ZUP=ZSOIL(K-1) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE CC TEMPERATURE INTO THE LAST LAYER BOUNDARY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(K .EQ. NSOIL) THEN ZB=2.*ZBOT-ZSOIL(K) ELSE ZB=ZSOIL(K+1) ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC TBND=TU+(TB-TU)*(ZUP-ZSOIL(K))/(ZUP-ZB) RETURN END SUBROUTINE TDFCND ( DF, SMC, B, F1, Q, SMCMAX, SH2O) IMPLICIT NONE C INCLUDED NEW VARS HERE: ,QUARTZ,SMCMAX,SH2O C______________________________________________________________ C NEW CALLING ARGUMENTS in SFLX, HRT 10/98 C (remove CALL TDFCND (DF, SMC, BB(ISOIL), F11(ISOIL)) C C CALL TDFCND (DF, SMC, BB(ISOIL), F11(ISOIL), C & QUARTZ(ISOIL), SMCMAX(ISOIL), SH2O) C-------------------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF CC ======= THE SOIL FOR A GIVEN POINT AND TIME. CC CC VERSION: PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998) CC ======= CC Added by PABLO J. GRUNMANN, 08/1998 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL B REAL DF REAL F1 REAL GAMMD REAL THKDRY REAL AKE REAL THKICE REAL THKO REAL THKQTZ REAL THKSAT REAL THKS REAL THKW REAL Q REAL SATRATIO REAL SH2O REAL SMC REAL SMCMAX REAL XU REAL XUNFROZ C IN PRMSOI: C DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52, C & 0.35, 0.60, 0.40, 0.82/ 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975). C C PABLO GRUNMANN, 08/17/98 C REFS.: C FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK C AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP. C JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS, C UNIVERSITY OF TRONDHEIM, C PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL C CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES C AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES, C VOL. 55, PP. 1209-1224. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THKW WATER THERMAL CONDUCTIVITY C THKQTZ THERMAL CONDUCTIVITY FOR QUARTZ C THKO THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS C THKS THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ + OTHER) C THKICE ICE THERMAL CONDUCTIVITY C SMCMAX POROSITY (= SMCMAX) C Q QUARTZ CONTENT(SOIL TYPE) C C NEEDS PARAMETERS C POROSITY(SOIL TYPE): C POROS = SMCMAX C SATURATION RATIO: SATRATIO = SMC/SMCMAX C print *, 'SATRATIO=',SATRATIO C PARAMETERS W/(M.K) THKICE = 2.2 THKW = 0.57 THKO = 2.0 C IF (Q .LE. 0.2) THKO = 3.0 THKQTZ = 7.7 C SOLIDS' CONDUCTIVITY THKS = (THKQTZ**Q)*(THKO**(1.- Q)) C print *, 'THKS = ',THKS C UNFROZEN FRACTION (FROM 1.0, I.E., 100%LIQUID, TO 0.0 (100% FROZEN)) XUNFROZ=(SH2O + 1.E-9)/(SMC + 1.E-9) C print *, ' ' C print *, 'XUNFROZ = ',XUNFROZ C print *, ' ' C UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ) XU=XUNFROZ*SMCMAX C SATURATED THERMAL CONDUCTIVITY THKSAT = THKS**(1.-SMCMAX)*THKICE**(SMCMAX-XU)*THKW**(XU) C print *, 'THKSAT = ',THKSAT C DRY DENSITY IN KG/M3 GAMMD = (1. - SMCMAX)*2700. C print *, 'GAMMD = ',GAMMD C DRY THERMAL CONDUCTIVITY IN W.M-1.K-1 THKDRY = (0.135*GAMMD + 64.7)/(2700. - 0.947*GAMMD) C print *, 'THKDRY = ',THKDRY C RANGE OF VALIDITY FOR THE KERSTEN NUMBER IF ( SATRATIO .GT. 0.1 ) THEN C KERSTEN NUMBER (FINE FORMULA, AT LEAST 5% OF PARTICLES<(2.E-6)M) IF ( (XUNFROZ + 0.0005) .LT. SMC ) THEN C FROZEN AKE = SATRATIO ELSE C UNFROZEN AKE = LOG10(SATRATIO) + 1.0 ENDIF ELSE C USE K = KDRY AKE = 0.0 C print *, 'AKE (ELSE) = ',AKE ENDIF C THERMAL CONDUCTIVITY DF = AKE*(THKSAT - THKDRY) + THKDRY RETURN END SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC PURPOSE: TO CALCULATE TRANSPIRATION FROM THE VEGTYP FOR THIS PT. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER I, K INTEGER NSOIL INTEGER NROOT REAL CFACTR REAL CMC REAL CMCMAX REAL ET ( NSOIL ) REAL ETP1 REAL ETP1A REAL GX (7) C REAL PART (4) REAL PC REAL RTDIS (NSOIL) REAL SHDFAC REAL SMC ( NSOIL ) REAL SMCREF REAL SMCWLT REAL ZSOIL ( NSOIL ) REAL SFCTMP, Q2, SGX, DENOM, RTX 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) SGX = 0.0 DO I = 1, NROOT GX(I) = ( SMC(I) - SMCWLT ) / ( SMCREF - SMCWLT ) GX(I) = MAX ( MIN ( GX(I), 1. ), 0. ) SGX = SGX + GX (I) END DO SGX = SGX / NROOT DENOM = 0. DO I = 1,NROOT RTX = RTDIS(I) + GX(I) - SGX GX(I) = GX(I) * MAX ( RTX, 0. ) DENOM = DENOM + GX(I) END DO IF ( DENOM .LE. 0.0) DENOM = 1. DO I = 1, NROOT ET(I) = ETP1A * GX(I) / DENOM END DO C ************ BELOW IS OLD CODE ************ C ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A C C ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A C*** USING ROOT DISTRIBUTION AS WEIGHTING FACTOR C C ET(1) = RTDIS(1) * ETP1A C C ET(1) = ETP1A*PART(1) 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 C C 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 C C ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A C*** USING ROOT DISTRIBUTION AS WEIGHTING FACTOR C C ET(K) = RTDIS(K) * ETP1A C C ET(K) = ETP1A*PART(K) C C 10 CONTINUE RETURN END SUBROUTINE WDFCND ( WDF,WCND,SMC,SMCMAX,B,DKSAT,DWSAT ) IMPLICIT NONE 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 C print*,'------------ in WDFCND -------------------------------' C print*,'BEFORE WDFCND' C print*,'B=',B C print*,'DKSAT=',DKSAT C print*,'DWSAT=',DWSAT C print*,'EXPON=',EXPON C print*,'FACTR2=',FACTR2 C print*,'SMC=',SMC C print*,'SMCMAX=',SMCMAX C print*,'WCND=',WCND C print*,'WDF=',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 C Drainage Tests !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C WDF = 0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EXPON = ( 2.0 * B ) + 3.0 WCND = DKSAT * FACTR2 ** EXPON C print*,' WDFCND Results --------------------------------' C print*,'B=',B C print*,'DKSAT=',DKSAT C print*,'DWSAT=',DWSAT C print*,'EXPON=',EXPON C print*,'FACTR2=',FACTR2 C print*,'SMC=',SMC C print*,'SMCMAX=',SMCMAX C print*,'WCND=',WCND C print*,'WDF=',WDF C PRINT*,' SMC WDF WCND B' C PRINT*,SMC,WDF,WCND,B RETURN END