| blob | 05c7aaf04e9dc9cf0e96ce851a696605057af1d4 |
1 MODULE module_sf_noahlsm
2 USE module_model_constants, only : CP, R_D, XLF, XLV, RHOWATER, STBOLT, KARMAN
4 ! REAL, PARAMETER :: CP = 1004.5
5 REAL, PARAMETER :: RD = 287.04, SIGMA = 5.67E-8, &
6 CPH2O = 4.218E+3,CPICE = 2.106E+3, &
7 LSUBF = 3.335E+5, &
8 EMISSI_S = 0.95
10 ! VEGETATION PARAMETERS
11 INTEGER :: LUCATS , BARE
12 INTEGER :: NATURAL
13 integer, PARAMETER :: NLUS=50
14 CHARACTER(LEN=256) LUTYPE
15 INTEGER, DIMENSION(1:NLUS) :: NROTBL
16 real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, &
17 SHDTBL, MAXALB, &
18 EMISSMINTBL, EMISSMAXTBL, &
19 LAIMINTBL, LAIMAXTBL, &
20 Z0MINTBL, Z0MAXTBL, &
21 ALBEDOMINTBL, ALBEDOMAXTBL
22 REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
24 ! SOIL PARAMETERS
25 INTEGER :: SLCATS
26 INTEGER, PARAMETER :: NSLTYPE=30
27 CHARACTER(LEN=256) SLTYPE
28 REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,F11, &
29 MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
31 ! LSM GENERAL PARAMETERS
32 INTEGER :: SLPCATS
33 INTEGER, PARAMETER :: NSLOPE=30
34 REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
35 REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
36 REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
37 CZIL_DATA
38 REAL :: LVCOEF_DATA
40 CHARACTER*256 :: err_message
42 integer, private :: iloc, jloc
43 !
44 CONTAINS
45 !
47 SUBROUTINE SFLX (IILOC,JJLOC,ICE,FFROZP,ISURBAN,DT,ZLVL,NSOIL,SLDPTH, & !C
48 LOCAL, & !L
49 LLANDUSE, LSOIL, & !CL
50 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD, & !F
51 COSZ,PRCPRAIN, SOLARDIRECT, & !F
52 TH2,Q2SAT,DQSDT2, & !I
53 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX, & !I
54 ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, EMBRD, & !S
55 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, & !H
56 ! ----------------------------------------------------------------------
57 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
58 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
59 ! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
60 ! ----------------------------------------------------------------------
61 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
62 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
63 BETA,ETP,SSOIL, & !O
64 FLX1,FLX2,FLX3, & !O
65 SNOMLT,SNCOVR, & !O
66 RUNOFF1,RUNOFF2,RUNOFF3, & !O
67 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
68 SOILW,SOILM,Q1,SMAV, & !D
69 RDLAI2D,USEMONALB, &
70 SNOTIME1, &
71 RIBB, &
72 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT) !P
73 ! ----------------------------------------------------------------------
74 ! SUBROUTINE SFLX - UNIFIED NOAHLSM VERSION 1.0 JULY 2007
75 ! ----------------------------------------------------------------------
76 ! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
77 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
78 ! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
79 ! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
80 ! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
81 ! RADIATION AND PRECIP)
82 ! ----------------------------------------------------------------------
83 ! SFLX ARGUMENT LIST KEY:
84 ! ----------------------------------------------------------------------
85 ! C CONFIGURATION INFORMATION
86 ! L LOGICAL
87 ! CL 4-string character bearing logical meaning
88 ! F FORCING DATA
89 ! I OTHER (INPUT) FORCING DATA
90 ! S SURFACE CHARACTERISTICS
91 ! H HISTORY (STATE) VARIABLES
92 ! O OUTPUT VARIABLES
93 ! D DIAGNOSTIC OUTPUT
94 ! P Parameters
95 ! Msic Miscellaneous terms passed from gridded driver
96 ! ----------------------------------------------------------------------
97 ! 1. CONFIGURATION INFORMATION (C):
98 ! ----------------------------------------------------------------------
99 ! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND (NO ICE), =-1 LAND-ICE).
100 ! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
101 ! 1800 SECS OR LESS)
102 ! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
103 ! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
104 ! PARAMETER NSOLD SET BELOW)
105 ! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
106 ! ----------------------------------------------------------------------
107 ! 2. LOGICAL:
108 ! ----------------------------------------------------------------------
109 ! LCH Exchange coefficient (Ch) calculation flag (false: using
110 ! ch-routine SFCDIF; true: Ch is brought in)
111 ! LOCAL Flag for local-site simulation (where there is no
112 ! maps for albedo, veg fraction, and roughness
113 ! true: all LSM parameters (inluding albedo, veg fraction and
114 ! roughness length) will be defined by three tables
115 ! LLANDUSE (=USGS, using USGS landuse classification)
116 ! LSOIL (=STAS, using FAO/STATSGO soil texture classification)
117 ! ----------------------------------------------------------------------
118 ! 3. FORCING DATA (F):
119 ! ----------------------------------------------------------------------
120 ! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
121 ! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
122 ! SOLNET NET DOWNWARD SOLAR RADIATION ((W M-2; POSITIVE)
123 ! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
124 ! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
125 ! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
126 ! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
127 ! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
128 ! COSZ Solar zenith angle (not used for now)
129 ! PRCPRAIN Liquid-precipitation rate (KG M-2 S-1) (not used)
130 ! SOLARDIRECT Direct component of downward solar radiation (W M-2) (not used)
131 ! FFROZP FRACTION OF FROZEN PRECIPITATION
132 ! ----------------------------------------------------------------------
133 ! 4. OTHER FORCING (INPUT) DATA (I):
134 ! ----------------------------------------------------------------------
135 ! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
136 ! Q2SAT SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
137 ! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
138 ! (KG KG-1 K-1)
139 ! ----------------------------------------------------------------------
140 ! 5. CANOPY/SOIL CHARACTERISTICS (S):
141 ! ----------------------------------------------------------------------
142 ! VEGTYP VEGETATION TYPE (INTEGER INDEX)
143 ! SOILTYP SOIL TYPE (INTEGER INDEX)
144 ! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
145 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
146 ! (FRACTION= 0.0-1.0)
147 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
148 ! (FRACTION= 0.0-1.0) <= SHDFAC
149 ! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
150 ! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
151 ! VEG PARMS)
152 ! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
153 ! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
154 ! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
155 ! INCLUDE DIURNAL SUN ANGLE EFFECT)
156 ! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
157 ! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
158 ! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
159 ! TEMPERATURE)
160 ! Z0BRD Background fixed roughness length (M)
161 ! Z0 Time varying roughness length (M) as function of snow depth
162 !
163 ! EMBRD Background surface emissivity (between 0 and 1)
164 ! EMISSI Surface emissivity (between 0 and 1)
165 ! ----------------------------------------------------------------------
166 ! 6. HISTORY (STATE) VARIABLES (H):
167 ! ----------------------------------------------------------------------
168 ! CMC CANOPY MOISTURE CONTENT (M)
169 ! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
170 ! STC(NSOIL) SOIL TEMP (K)
171 ! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
172 ! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
173 ! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
174 ! SNOWH ACTUAL SNOW DEPTH (M)
175 ! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
176 ! NOTE: SNOW DENSITY = SNEQV/SNOWH
177 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
178 ! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
179 ! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
180 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
181 ! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
182 ! IT HAS BEEN MULTIPLIED BY WIND SPEED.
183 ! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
184 ! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
185 ! MULTIPLIED BY WIND SPEED.
186 ! ----------------------------------------------------------------------
187 ! 7. OUTPUT (O):
188 ! ----------------------------------------------------------------------
189 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
190 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
191 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
192 ! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
193 ! ETA ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
194 ! SURFACE)
195 ! ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
196 ! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
197 ! SURFACE)
198 ! FDOWN Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
199 ! ----------------------------------------------------------------------
200 ! EC CANOPY WATER EVAPORATION (W m-2)
201 ! EDIR DIRECT SOIL EVAPORATION (W m-2)
202 ! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
203 ! (W m-2)
204 ! ETT TOTAL PLANT TRANSPIRATION (W m-2)
205 ! ESNOW SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
206 ! (W m-2)
207 ! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
208 ! WATER-HOLDING CAPACITY (M)
209 ! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
210 ! ----------------------------------------------------------------------
211 ! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
212 ! ETP POTENTIAL EVAPORATION (W m-2)
213 ! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
214 ! ----------------------------------------------------------------------
215 ! FLX1 PRECIP-SNOW SFC (W M-2)
216 ! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
217 ! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
218 ! ----------------------------------------------------------------------
219 ! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
220 ! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
221 ! ----------------------------------------------------------------------
222 ! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
223 ! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
224 ! SOIL LAYER (BASEFLOW)
225 ! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
226 ! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
227 ! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
228 ! ----------------------------------------------------------------------
229 ! RC CANOPY RESISTANCE (S M-1)
230 ! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
231 ! = ACTUAL TRANSP
232 ! XLAI LEAF AREA INDEX (DIMENSIONLESS)
233 ! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
234 ! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
235 ! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
236 ! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
237 ! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
238 ! ----------------------------------------------------------------------
239 ! 8. DIAGNOSTIC OUTPUT (D):
240 ! ----------------------------------------------------------------------
241 ! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
242 ! BETWEEN SMCWLT AND SMCMAX)
243 ! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
244 ! Q1 Effective mixing ratio at surface (kg kg-1), used for
245 ! diagnosing the mixing ratio at 2 meter for coupled model
246 ! SMAV Soil Moisture Availability for each layer, as a fraction
247 ! between SMCWLT and SMCMAX.
248 ! Documentation for SNOTIME1 and SNOABL2 ?????
249 ! What categories of arguments do these variables fall into ????
250 ! Documentation for RIBB ?????
251 ! What category of argument does RIBB fall into ?????
252 ! ----------------------------------------------------------------------
253 ! 9. PARAMETERS (P):
254 ! ----------------------------------------------------------------------
255 ! SMCWLT WILTING POINT (VOLUMETRIC)
256 ! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
257 ! LAYER ENDS (VOLUMETRIC)
258 ! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
259 ! STRESS (VOLUMETRIC)
260 ! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
261 ! (VOLUMETRIC)
262 ! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
263 ! IN SUBROUTINE REDPRM.
264 ! ----------------------------------------------------------------------
267 IMPLICIT NONE
268 ! ----------------------------------------------------------------------
270 ! DECLARATIONS - LOGICAL AND CHARACTERS
271 ! ----------------------------------------------------------------------
273 INTEGER, INTENT(IN) :: IILOC, JJLOC
274 LOGICAL, INTENT(IN):: LOCAL
275 LOGICAL :: FRZGRA, SNOWNG
276 CHARACTER (LEN=256), INTENT(IN):: LLANDUSE, LSOIL
278 ! ----------------------------------------------------------------------
279 ! 1. CONFIGURATION INFORMATION (C):
280 ! ----------------------------------------------------------------------
281 INTEGER,INTENT(IN) :: ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
282 INTEGER, INTENT(IN) :: ISURBAN
283 INTEGER,INTENT(OUT):: NROOT
284 INTEGER KZ, K, iout
286 ! ----------------------------------------------------------------------
287 ! 2. LOGICAL:
288 ! ----------------------------------------------------------------------
289 LOGICAL, INTENT(IN) :: RDLAI2D
290 LOGICAL, INTENT(IN) :: USEMONALB
292 REAL, INTENT(IN) :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN, &
293 Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB, &
294 SOLDN,SOLNET,TBOT,TH2,ZLVL, &
295 FFROZP
296 REAL, INTENT(OUT) :: EMBRD
297 REAL, INTENT(OUT) :: ALBEDO
298 REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,CH,CM, &
299 CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD, &
300 EMISSI, ALB
301 REAL, INTENT(INOUT):: SNOTIME1
302 REAL, INTENT(INOUT):: RIBB
303 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
304 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
305 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: SMAV
306 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O, SMC, STC
307 REAL,DIMENSION(1:NSOIL):: RTDIS, ZSOIL
309 REAL,INTENT(OUT) :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA, &
310 ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2, &
311 RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL, &
312 SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM, &
313 SOILW,FDOWN,Q1
314 REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT, &
315 DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS, &
316 KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL, &
317 RSMAX, &
318 RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM, &
319 SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF, &
320 ETNS,PTU,LSUBS
321 REAL :: LVCOEF
322 REAL :: INTERP_FRACTION
323 REAL :: LAIMIN, LAIMAX
324 REAL :: ALBEDOMIN, ALBEDOMAX
325 REAL :: EMISSMIN, EMISSMAX
326 REAL :: Z0MIN, Z0MAX
328 ! ----------------------------------------------------------------------
329 ! DECLARATIONS - PARAMETERS
330 ! ----------------------------------------------------------------------
331 PARAMETER (TFREEZ = 273.15)
332 PARAMETER (LVH2O = 2.501E+6)
333 PARAMETER (LSUBS = 2.83E+6)
334 PARAMETER (R = 287.04)
336 ! ----------------------------------------------------------------------
337 ! INITIALIZATION
338 ! ----------------------------------------------------------------------
339 ILOC = IILOC
340 JLOC = JJLOC
342 RUNOFF1 = 0.0
343 RUNOFF2 = 0.0
344 RUNOFF3 = 0.0
345 SNOMLT = 0.0
347 ! ----------------------------------------------------------------------
348 ! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE / LAND-ICE / ICE-FREE LAND
349 ! SEA-ICE CASE, ICE = 1
350 ! NON-GLACIAL LAND, ICE = 0
351 ! GLACIAL ICE, ICE = -1
352 IF (ICE == 1) call wrf_error_fatal("Sea-ice point in Noah-LSM")
353 IF (ICE == -1) SHDFAC = 0.0
355 ! ----------------------------------------------------------------------
356 ! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
357 ! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
358 ! GROUND)
359 ! ----------------------------------------------------------------------
360 ZSOIL (1) = - SLDPTH (1)
361 DO KZ = 2,NSOIL
362 ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
363 END DO
365 ! ----------------------------------------------------------------------
366 ! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
367 ! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
368 ! ----------------------------------------------------------------------
369 CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT, &
370 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
371 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
372 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
373 RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,CZIL, &
374 LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN, &
375 ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE, &
376 LSOIL,LOCAL,LVCOEF)
378 !urban
379 IF(VEGTYP==ISURBAN)THEN
380 SHDFAC=0.05
381 RSMIN=400.0
382 SMCMAX = 0.45
383 SMCREF = 0.42
384 SMCWLT = 0.40
385 SMCDRY = 0.40
386 ENDIF
388 IF ( SHDFAC >= SHDMAX ) THEN
389 EMBRD = EMISSMAX
390 IF (.NOT. RDLAI2D) THEN
391 XLAI = LAIMAX
392 ENDIF
393 IF (.NOT. USEMONALB) THEN
394 ALB = ALBEDOMIN
395 ENDIF
396 Z0BRD = Z0MAX
397 ELSE IF ( SHDFAC <= SHDMIN ) THEN
398 EMBRD = EMISSMIN
399 IF(.NOT. RDLAI2D) THEN
400 XLAI = LAIMIN
401 ENDIF
402 IF(.NOT. USEMONALB) then
403 ALB = ALBEDOMAX
404 ENDIF
405 Z0BRD = Z0MIN
406 ELSE
408 IF ( SHDMAX > SHDMIN ) THEN
410 INTERP_FRACTION = ( SHDFAC - SHDMIN ) / ( SHDMAX - SHDMIN )
411 ! Bound INTERP_FRACTION between 0 and 1
412 INTERP_FRACTION = MIN ( INTERP_FRACTION, 1.0 )
413 INTERP_FRACTION = MAX ( INTERP_FRACTION, 0.0 )
414 ! Scale Emissivity and LAI between EMISSMIN and EMISSMAX by INTERP_FRACTION
415 EMBRD = ( ( 1.0 - INTERP_FRACTION ) * EMISSMIN ) + ( INTERP_FRACTION * EMISSMAX )
416 IF (.NOT. RDLAI2D) THEN
417 XLAI = ( ( 1.0 - INTERP_FRACTION ) * LAIMIN ) + ( INTERP_FRACTION * LAIMAX )
418 ENDIF
419 if (.not. USEMONALB) then
420 ALB = ( ( 1.0 - INTERP_FRACTION ) * ALBEDOMAX ) + ( INTERP_FRACTION * ALBEDOMIN )
421 endif
422 Z0BRD = ( ( 1.0 - INTERP_FRACTION ) * Z0MIN ) + ( INTERP_FRACTION * Z0MAX )
424 ELSE
426 EMBRD = 0.5 * EMISSMIN + 0.5 * EMISSMAX
427 IF (.NOT. RDLAI2D) THEN
428 XLAI = 0.5 * LAIMIN + 0.5 * LAIMAX
429 ENDIF
430 if (.not. USEMONALB) then
431 ALB = 0.5 * ALBEDOMIN + 0.5 * ALBEDOMAX
432 endif
433 Z0BRD = 0.5 * Z0MIN + 0.5 * Z0MAX
435 ENDIF
437 ENDIF
438 ! ----------------------------------------------------------------------
439 ! INITIALIZE PRECIPITATION LOGICALS.
440 ! ----------------------------------------------------------------------
441 SNOWNG = .FALSE.
442 FRZGRA = .FALSE.
444 ! ----------------------------------------------------------------------
445 IF ( ICE == -1 ) THEN
446 !
447 ! FOR GLACIAL ICE, IF S.W.E. (SNEQV) BELOW THRESHOLD LOWER
448 ! BOUND (0.10 M FOR GLACIAL ICE), THEN SET AT LOWER BOUND.
449 !
450 IF ( SNEQV < 0.10 ) THEN
451 SNEQV = 0.10
452 SNOWH = 0.50
453 ENDIF
454 !
455 ! FOR GLACIAL ICE, SET SMC AND SH20 VALUES = 1.0
456 !
457 DO KZ = 1,NSOIL
458 SMC(KZ) = 1.0
459 SH2O(KZ) = 1.0
460 END DO
461 ENDIF
462 ! ----------------------------------------------------------------------
463 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
464 ! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
465 ! SUBROUTINE)
466 ! ----------------------------------------------------------------------
467 IF ( SNEQV <= 1.E-7 ) THEN ! safer IF kmh (2008/03/25)
468 SNEQV = 0.0
469 SNDENS = 0.0
470 SNOWH = 0.0
471 SNCOND = 1.0
472 ELSE
473 SNDENS = SNEQV / SNOWH
474 IF(SNDENS > 1.0) THEN
475 CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
476 ENDIF
477 CALL CSNOW (SNCOND,SNDENS)
478 END IF
479 ! ----------------------------------------------------------------------
480 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
481 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
482 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
483 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
484 ! ----------------------------------------------------------------------
485 IF (PRCP > 0.0) THEN
486 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
487 ! passed in from model microphysics.
488 IF (FFROZP .GT. 0.5) THEN
489 SNOWNG = .TRUE.
490 ELSE
491 IF (T1 <= TFREEZ) FRZGRA = .TRUE.
492 END IF
493 END IF
494 ! ----------------------------------------------------------------------
495 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
496 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
497 ! IT TO THE EXISTING SNOWPACK.
498 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
499 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
500 ! ----------------------------------------------------------------------
501 IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
502 SN_NEW = PRCP * DT * 0.001
503 SNEQV = SNEQV + SN_NEW
504 PRCPF = 0.0
506 ! ----------------------------------------------------------------------
507 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
508 ! UPDATE SNOW THERMAL CONDUCTIVITY
509 ! ----------------------------------------------------------------------
510 CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
511 !
512 ! kmh 09/04/2006 set Snow Density at 0.2 g/cm**3
513 ! for "cold permanent ice" or new "dry" snow
514 !
515 IF ( (ICE == -1) .and. (SNCOVR .GT. 0.99) ) THEN
516 ! if soil temperature less than 268.15 K, treat as typical Antarctic/Greenland snow firn
517 IF ( STC(1) .LT. (TFREEZ - 5.) ) SNDENS = 0.2
518 IF ( SNOWNG .AND. (T1.LT.273.) .AND. (SFCTMP.LT.273.) ) SNDENS=0.2
519 ENDIF
520 !
521 CALL CSNOW (SNCOND,SNDENS)
523 ! ----------------------------------------------------------------------
524 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
525 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
526 ! ANY CANOPY "DRIP" ADDED TO THIS LATER)
527 ! ----------------------------------------------------------------------
528 ELSE
529 PRCPF = PRCP
530 ENDIF
531 ! ----------------------------------------------------------------------
532 ! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
533 ! ----------------------------------------------------------------------
534 ! ----------------------------------------------------------------------
535 ! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
536 ! ----------------------------------------------------------------------
537 IF (SNEQV == 0.0) THEN
538 SNCOVR = 0.0
539 ALBEDO = ALB
540 EMISSI = EMBRD
541 ELSE
542 ! ----------------------------------------------------------------------
543 ! DETERMINE SNOW FRACTIONAL COVERAGE.
544 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
545 ! ----------------------------------------------------------------------
546 CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
548 ! kmh 2008/03/25
549 ! Don't limit snow cover fraction over permanent ice
550 IF ( ICE == 0 ) SNCOVR = MIN(SNCOVR,0.98)
552 CALL ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,T1, &
553 ALBEDO,EMISSI,DT,SNOWNG,SNOTIME1,LVCOEF)
554 ENDIF
555 ! ----------------------------------------------------------------------
556 ! THERMAL CONDUCTIVITY FOR GLACIAL ICE CASE
557 ! ----------------------------------------------------------------------
558 IF ( ICE == -1 ) THEN
559 DF1 = 2.2
560 !
561 ! kmh 09/03/2006
562 ! kmh 03/25/2008 change SNCOVR threshold to 0.97
563 !
564 ! only apply (small) DF1 conductivity for permanent land ice
565 !
566 IF ( SNCOVR .GT. 0.97 ) THEN
567 DF1 = SNCOND
568 ENDIF
569 ELSEIF ( ICE == 0 ) THEN
570 ! ----------------------------------------------------------------------
571 ! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
572 ! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
573 ! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
574 ! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
575 ! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
576 ! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
577 ! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
578 ! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
579 ! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
580 ! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
581 ! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
582 ! ----------------------------------------------------------------------
583 ! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
584 ! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
585 ! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
586 ! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
587 ! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
588 ! ----------------------------------------------------------------------
589 ! ----------------------------------------------------------------------
590 ! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
591 ! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
592 ! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
593 ! ----------------------------------------------------------------------
594 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
596 !urban
597 IF ( VEGTYP == ISURBAN ) DF1=3.24
599 DF1 = DF1 * EXP (SBETA * SHDFAC)
600 !
601 ! kmh 09/03/2006
602 ! kmh 03/25/2008 change SNCOVR threshold to 0.97
603 !
604 IF ( SNCOVR .GT. 0.97 ) THEN
605 DF1 = SNCOND
606 ENDIF
607 !
608 ! ----------------------------------------------------------------------
609 ! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
610 ! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
611 ! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
612 ! ----------------------------------------------------------------------
613 END IF
615 DSOIL = - (0.5 * ZSOIL (1))
616 IF (SNEQV == 0.) THEN
617 SSOIL = DF1 * (T1- STC (1) ) / DSOIL
618 ELSE
619 DTOT = SNOWH + DSOIL
620 FRCSNO = SNOWH / DTOT
622 ! 1. HARMONIC MEAN (SERIES FLOW)
623 ! DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
624 FRCSOI = DSOIL / DTOT
625 ! 2. ARITHMETIC MEAN (PARALLEL FLOW)
626 ! DF1 = FRCSNO*SNCOND + FRCSOI*DF1
627 DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)
629 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
630 ! DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
631 ! weigh DF by snow fraction
632 ! DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
633 ! DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
634 DF1A = FRCSNO * SNCOND+ FRCSOI * DF1
636 ! ----------------------------------------------------------------------
637 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
638 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
639 ! MID-LAYER SOIL TEMPERATURE
640 ! ----------------------------------------------------------------------
641 DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
642 IF ( ICE == -1 ) then
643 ! kmh 12/15/2005 correct for too deep snow layer
644 ! kmh 09/03/2006 adjust DTOT
645 IF ( DTOT .GT. 2.*DSOIL ) then
646 DTOT = 2.*DSOIL
647 ENDIF
648 ENDIF
649 SSOIL = DF1 * (T1- STC (1) ) / DTOT
650 END IF
651 ! ----------------------------------------------------------------------
652 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
653 ! THE PREVIOUS TIMESTEP.
654 ! ----------------------------------------------------------------------
655 IF (SNCOVR > 0. ) THEN
656 CALL SNOWZ0 (SNCOVR,Z0,Z0BRD,SNOWH)
657 ELSE
658 Z0=Z0BRD
659 END IF
660 ! ----------------------------------------------------------------------
661 ! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
662 ! HEAT AND MOISTURE.
664 ! NOTE !!!
665 ! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
666 ! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
667 ! (CZIL) ARE SET THERE VIA NAMELIST I/O.
669 ! NOTE !!!
670 ! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
671 ! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
672 ! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
673 ! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
674 ! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
675 ! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
677 ! NOTE !!!
678 ! ----------------------------------------------------------------------
679 ! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
680 ! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
681 ! for iterative/implicit solution of CH in SFCDIF
682 ! ----------------------------------------------------------------------
683 ! IF(.NOT.LCH) THEN
684 ! T1V = T1 * (1.0+ 0.61 * Q2)
685 ! TH2V = TH2 * (1.0+ 0.61 * Q2)
686 ! CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
687 ! ENDIF
689 ! ----------------------------------------------------------------------
690 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
691 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
692 ! CALCULATIONS.
693 ! ----------------------------------------------------------------------
694 ! ----------------------------------------------------------------------
695 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
696 ! PENMAN EP SUBROUTINE THAT FOLLOWS
697 ! ----------------------------------------------------------------------
698 ! FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
699 FDOWN = SOLNET + LWDN
700 ! ----------------------------------------------------------------------
701 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
702 ! PENMAN.
703 T2V = SFCTMP * (1.0+ 0.61 * Q2 )
705 iout=0
706 if(iout.eq.1) then
707 print*,'before penman'
708 print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V, &
709 'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL, &
710 'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH, &
711 'EPSCA',EPSCA,'RR',RR ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA, &
712 'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV, &
713 ' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT, &
714 ' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1), &
715 'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
716 endif
718 CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
719 Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
720 !
721 ! kmh 01/09/2007 add T1,ICE,SNCOVR to call
722 !
723 DQSDT2,FLX2,EMISSI,SNEQV,T1,ICE,SNCOVR)
724 !
725 ! ----------------------------------------------------------------------
726 ! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
727 ! IF NONZERO GREENNESS FRACTION
728 ! ----------------------------------------------------------------------
730 ! ----------------------------------------------------------------------
731 ! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
732 ! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
733 ! ----------------------------------------------------------------------
734 IF (SHDFAC > 0.) THEN
735 CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
736 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
737 TOPT,RSMAX,RGL,HS,XLAI, &
738 RCS,RCT,RCQ,RCSOIL,EMISSI)
739 ELSE
740 RC = 0.0
741 END IF
742 ! ----------------------------------------------------------------------
743 ! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
744 ! EXISTS OR NOT:
745 ! ----------------------------------------------------------------------
746 ESNOW = 0.0
747 IF (SNEQV == 0.0) THEN
748 CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
749 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
750 SHDFAC, &
751 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
752 SSOIL, &
753 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
754 SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
755 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
756 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
757 QUARTZ,FXEXP,CSOIL, &
758 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
759 ETA_KINEMATIC = ETA
760 ELSE
761 CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
762 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
763 SBETA,DF1, &
764 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
765 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
766 SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT, &
767 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
768 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
769 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
770 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
771 RIBB,SOLDN, &
772 ISURBAN, &
773 VEGTYP)
774 ETA_KINEMATIC = ESNOW + ETNS
775 END IF
777 ! Calculate effective mixing ratio at grnd level (skin)
778 !
779 ! Q1=Q2+ETA*CP/RCH
780 Q1=Q2+ETA_KINEMATIC*CP/RCH
781 !
782 ! ----------------------------------------------------------------------
783 ! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
784 ! ----------------------------------------------------------------------
785 SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
787 ! ----------------------------------------------------------------------
788 ! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
789 ! ----------------------------------------------------------------------
790 EDIR = EDIR * LVH2O
791 EC = EC * LVH2O
792 DO K=1,4
793 ET(K) = ET(K) * LVH2O
794 ENDDO
795 ETT = ETT * LVH2O
796 ESNOW = ESNOW * LSUBS
797 ETP = ETP*((1.-SNCOVR)*LVH2O + SNCOVR*LSUBS)
798 IF (ETP .GT. 0.) THEN
799 ETA = EDIR + EC + ETT + ESNOW
800 ELSE
801 ETA = ETP
802 ENDIF
803 ! ----------------------------------------------------------------------
804 ! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
805 ! ----------------------------------------------------------------------
806 IF (ETP == 0.0) THEN
807 BETA = 0.0
808 ELSE
809 BETA = ETA/ETP
810 ENDIF
812 ! ----------------------------------------------------------------------
813 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
814 ! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
815 ! SSOIL<0: COOL THE SURFACE (DAY TIME)
816 ! ----------------------------------------------------------------------
817 SSOIL = -1.0* SSOIL
819 ! ----------------------------------------------------------------------
820 ! FOR THE CASE OF LAND (BUT NOT GLACIAL ICE):
821 ! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
822 ! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW. RUNOFF2 IS ALREADY
823 ! A RATE AT THIS POINT
824 ! ----------------------------------------------------------------------
825 IF (ICE == 0) THEN
826 RUNOFF3 = RUNOFF3/ DT
827 RUNOFF2 = RUNOFF2+ RUNOFF3
828 SOILM = -1.0* SMC (1)* ZSOIL (1)
829 DO K = 2,NSOIL
830 SOILM = SOILM + SMC (K)* (ZSOIL (K -1) - ZSOIL (K))
831 END DO
832 SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
833 SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
834 !
835 DO K = 1,NSOIL
836 SMAV(K)=(SMC(K) - SMCWLT)/(SMCMAX - SMCWLT)
837 END DO
839 IF (NROOT >= 2) THEN
840 DO K = 2,NROOT
841 SOILWM = SOILWM + (SMCMAX - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
842 SOILWW = SOILWW + (SMC(K) - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
843 END DO
844 END IF
845 IF (SOILWM .LT. 1.E-6) THEN
846 SOILWM = 0.0
847 SOILW = 0.0
848 SOILM = 0.0
849 ELSE
850 SOILW = SOILWW / SOILWM
851 END IF
852 ELSEIF ( ICE == -1 ) THEN
853 ! ----------------------------------------------------------------------
854 ! FOR THE CASE OF GLACIAL ICE (ICE == -1), ADD ANY SNOWMELT DIRECTLY TO
855 ! SURFACE RUNOFF (RUNOFF1) SINCE THERE IS NO SOIL MEDIUM, AND THUS NO
856 ! CALL TO SUBROUTINE SMFLX (FOR SOIL MOISTURE TENDENCY).
857 ! ----------------------------------------------------------------------
858 RUNOFF1 = SNOMLT/DT
859 SOILWM = 0.0
860 SOILW = 0.0
861 SOILM = 0.0
862 DO K = 1,NSOIL
863 SMAV(K)= 1.0
864 END DO
865 END IF
867 ! ----------------------------------------------------------------------
868 END SUBROUTINE SFLX
869 ! ----------------------------------------------------------------------
871 SUBROUTINE ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO,EMISSI, &
872 DT,SNOWNG,SNOTIME1,LVCOEF)
874 ! ----------------------------------------------------------------------
875 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
876 ! ALB SNOWFREE ALBEDO
877 ! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
878 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
879 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
880 ! SNCOVR FRACTIONAL SNOW COVER
881 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
882 ! TSNOW SNOW SURFACE TEMPERATURE (K)
883 ! ----------------------------------------------------------------------
884 IMPLICIT NONE
886 ! ----------------------------------------------------------------------
887 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
888 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
889 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
890 ! (1985, JCAM, VOL 24, 402-411)
891 ! ----------------------------------------------------------------------
892 REAL, INTENT(IN) :: ALB, SNOALB, EMBRD, SHDFAC, SHDMIN, SNCOVR, TSNOW
893 REAL, INTENT(IN) :: DT
894 LOGICAL, INTENT(IN) :: SNOWNG
895 REAL, INTENT(INOUT):: SNOTIME1
896 REAL, INTENT(OUT) :: ALBEDO, EMISSI
897 REAL :: SNOALB2
898 REAL :: TM,SNOALB1
899 REAL, INTENT(IN) :: LVCOEF
900 REAL, PARAMETER :: SNACCA=0.94,SNACCB=0.58,SNTHWA=0.82,SNTHWB=0.46
901 ! turn of vegetation effect
902 ! ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
903 ! ALBEDO = (1.0-SNCOVR)*ALB + SNCOVR*SNOALB !this is equivalent to below
904 ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
905 EMISSI = EMBRD + SNCOVR*(EMISSI_S - EMBRD)
907 ! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
908 ! IF (TSNOW.LE.263.16) THEN
909 ! ALBEDO=SNOALB
910 ! ELSE
911 ! IF (TSNOW.LT.273.16) THEN
912 ! TM=0.1*(TSNOW-263.16)
913 ! SNOALB1=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
914 ! ELSE
915 ! SNOALB1=0.67
916 ! IF(SNCOVR.GT.0.95) SNOALB1= 0.6
917 ! SNOALB1 = ALB + SNCOVR*(SNOALB-ALB)
918 ! ENDIF
919 ! ENDIF
920 ! ALBEDO = ALB + SNCOVR*(SNOALB1-ALB)
922 ! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
923 ! SNOALB1 = SNOALB+COEF*(0.85-SNOALB)
924 ! SNOALB2=SNOALB1
925 !!m LSTSNW=LSTSNW+1
926 ! SNOTIME1 = SNOTIME1 + DT
927 ! IF (SNOWNG) THEN
928 ! SNOALB2=SNOALB
929 !!m LSTSNW=0
930 ! SNOTIME1 = 0.0
931 ! ELSE
932 ! IF (TSNOW.LT.273.16) THEN
933 !! SNOALB2=SNOALB-0.008*LSTSNW*DT/86400
934 !!m SNOALB2=SNOALB-0.008*SNOTIME1/86400
935 ! SNOALB2=(SNOALB2-0.65)*EXP(-0.05*DT/3600)+0.65
936 !! SNOALB2=(ALBEDO-0.65)*EXP(-0.01*DT/3600)+0.65
937 ! ELSE
938 ! SNOALB2=(SNOALB2-0.5)*EXP(-0.0005*DT/3600)+0.5
939 !! SNOALB2=(SNOALB-0.5)*EXP(-0.24*LSTSNW*DT/86400)+0.5
940 !!m SNOALB2=(SNOALB-0.5)*EXP(-0.24*SNOTIME1/86400)+0.5
941 ! ENDIF
942 ! ENDIF
943 !
944 !! print*,'SNOALB2',SNOALB2,'ALBEDO',ALBEDO,'DT',DT
945 ! ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
946 ! IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
947 !!m LSTSNW1=LSTSNW
948 !! SNOTIME = SNOTIME1
950 ! formulation by Livneh
951 ! ----------------------------------------------------------------------
952 ! SNOALB IS CONSIDERED AS THE MAXIMUM SNOW ALBEDO FOR NEW SNOW, AT
953 ! A VALUE OF 85%. SNOW ALBEDO CURVE DEFAULTS ARE FROM BRAS P.263. SHOULD
954 ! NOT BE CHANGED EXCEPT FOR SERIOUS PROBLEMS WITH SNOW MELT.
955 ! TO IMPLEMENT ACCUMULATIN PARAMETERS, SNACCA AND SNACCB, ASSERT THAT IT
956 ! IS INDEED ACCUMULATION SEASON. I.E. THAT SNOW SURFACE TEMP IS BELOW
957 ! ZERO AND THE DATE FALLS BETWEEN OCTOBER AND FEBRUARY
958 ! ----------------------------------------------------------------------
959 SNOALB1 = SNOALB+LVCOEF*(0.85-SNOALB)
960 SNOALB2=SNOALB1
961 ! ---------------- Initial LSTSNW --------------------------------------
962 IF (SNOWNG) THEN
963 SNOTIME1 = 0.
964 ELSE
965 SNOTIME1=SNOTIME1+DT
966 ! IF (TSNOW.LT.273.16) THEN
967 SNOALB2=SNOALB1*(SNACCA**((SNOTIME1/86400.0)**SNACCB))
968 ! ELSE
969 ! SNOALB2 =SNOALB1*(SNTHWA**((SNOTIME1/86400.0)**SNTHWB))
970 ! ENDIF
971 ENDIF
972 !
973 SNOALB2 = MAX ( SNOALB2, ALB )
974 ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
975 IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
977 ! IF (TSNOW.LT.273.16) THEN
978 ! ALBEDO=SNOALB-0.008*DT/86400
979 ! ELSE
980 ! ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
981 ! ENDIF
983 ! IF (ALBEDO > SNOALB) ALBEDO = SNOALB
985 ! ----------------------------------------------------------------------
986 END SUBROUTINE ALCALC
987 ! ----------------------------------------------------------------------
989 SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
990 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
991 TOPT,RSMAX,RGL,HS,XLAI, &
992 RCS,RCT,RCQ,RCSOIL,EMISSI)
994 ! ----------------------------------------------------------------------
995 ! SUBROUTINE CANRES
996 ! ----------------------------------------------------------------------
997 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
998 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
999 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1000 ! MOISTURE RATHER THAN TOTAL)
1001 ! ----------------------------------------------------------------------
1002 ! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1003 ! NOILHAN (1990, BLM)
1004 ! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1005 ! AND TABLE 2 OF SEC. 3.1.2
1006 ! ----------------------------------------------------------------------
1007 ! INPUT:
1008 ! SOLAR INCOMING SOLAR RADIATION
1009 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1010 ! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1011 ! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1012 ! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1013 ! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1014 ! SFCPRS SURFACE PRESSURE
1015 ! SMC VOLUMETRIC SOIL MOISTURE
1016 ! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1017 ! NSOIL NO. OF SOIL LAYERS
1018 ! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1019 ! XLAI LEAF AREA INDEX
1020 ! SMCWLT WILTING POINT
1021 ! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1022 ! SETS IN)
1023 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1024 ! SURBOUTINE REDPRM
1025 ! OUTPUT:
1026 ! PC PLANT COEFFICIENT
1027 ! RC CANOPY RESISTANCE
1028 ! ----------------------------------------------------------------------
1030 IMPLICIT NONE
1031 INTEGER, INTENT(IN) :: NROOT,NSOIL
1032 INTEGER K
1033 REAL, INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX, &
1034 SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
1035 EMISSI
1036 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
1037 REAL, INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
1038 REAL :: DELTA,FF,GX,P,RR
1039 REAL, DIMENSION(1:NSOIL) :: PART
1040 REAL, PARAMETER :: SLV = 2.501000E6
1043 ! ----------------------------------------------------------------------
1044 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1045 ! ----------------------------------------------------------------------
1046 RCS = 0.0
1047 RCT = 0.0
1048 RCQ = 0.0
1049 RCSOIL = 0.0
1051 ! ----------------------------------------------------------------------
1052 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1053 ! ----------------------------------------------------------------------
1054 RC = 0.0
1055 FF = 0.55*2.0* SOLAR / (RGL * XLAI)
1056 RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
1058 ! ----------------------------------------------------------------------
1059 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1060 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1061 ! ----------------------------------------------------------------------
1062 RCS = MAX (RCS,0.0001)
1063 RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)
1065 ! ----------------------------------------------------------------------
1066 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1067 ! RCQ EXPRESSION FROM SSIB
1068 ! ----------------------------------------------------------------------
1069 RCT = MAX (RCT,0.0001)
1070 RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))
1072 ! ----------------------------------------------------------------------
1073 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1074 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1075 ! ----------------------------------------------------------------------
1076 RCQ = MAX (RCQ,0.01)
1077 GX = (SMC (1) - SMCWLT) / (SMCREF - SMCWLT)
1078 IF (GX > 1.) GX = 1.
1079 IF (GX < 0.) GX = 0.
1081 ! ----------------------------------------------------------------------
1082 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1083 ! ----------------------------------------------------------------------
1084 ! ----------------------------------------------------------------------
1085 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1086 ! PART(1) = RTDIS(1) * GX
1087 ! ----------------------------------------------------------------------
1088 PART (1) = (ZSOIL (1)/ ZSOIL (NROOT)) * GX
1089 DO K = 2,NROOT
1090 GX = (SMC (K) - SMCWLT) / (SMCREF - SMCWLT)
1091 IF (GX > 1.) GX = 1.
1092 IF (GX < 0.) GX = 0.
1093 ! ----------------------------------------------------------------------
1094 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1095 ! ----------------------------------------------------------------------
1096 ! ----------------------------------------------------------------------
1097 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1098 ! PART(K) = RTDIS(K) * GX
1099 ! ----------------------------------------------------------------------
1100 PART (K) = ( (ZSOIL (K) - ZSOIL (K -1))/ ZSOIL (NROOT)) * GX
1101 END DO
1102 DO K = 1,NROOT
1103 RCSOIL = RCSOIL + PART (K)
1104 END DO
1106 ! ----------------------------------------------------------------------
1107 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
1108 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1109 ! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
1110 ! PC * LINERIZED PENMAN POTENTIAL EVAP =
1111 ! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1112 ! ----------------------------------------------------------------------
1113 RCSOIL = MAX (RCSOIL,0.0001)
1115 RC = RSMIN / (XLAI * RCS * RCT * RCQ * RCSOIL)
1116 ! RR = (4.* SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) + 1.0
1117 RR = (4.* EMISSI *SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) &
1118 + 1.0
1120 DELTA = (SLV / CP)* DQSDT2
1122 PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
1124 ! ----------------------------------------------------------------------
1125 END SUBROUTINE CANRES
1126 ! ----------------------------------------------------------------------
1128 SUBROUTINE CSNOW (SNCOND,DSNOW)
1130 ! ----------------------------------------------------------------------
1131 ! SUBROUTINE CSNOW
1132 ! FUNCTION CSNOW
1133 ! ----------------------------------------------------------------------
1134 ! CALCULATE SNOW TERMAL CONDUCTIVITY
1135 ! ----------------------------------------------------------------------
1136 IMPLICIT NONE
1137 REAL, INTENT(IN) :: DSNOW
1138 REAL, INTENT(OUT):: SNCOND
1139 REAL :: C
1140 REAL, PARAMETER :: UNIT = 0.11631
1142 ! ----------------------------------------------------------------------
1143 ! SNCOND IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
1144 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
1145 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
1146 ! ----------------------------------------------------------------------
1147 C = 0.328*10** (2.25* DSNOW)
1148 ! CSNOW=UNIT*C
1150 ! ----------------------------------------------------------------------
1151 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
1152 ! ----------------------------------------------------------------------
1153 ! SNCOND=0.0293*(1.+100.*DSNOW**2)
1154 ! CSNOW=0.0293*(1.+100.*DSNOW**2)
1156 ! ----------------------------------------------------------------------
1157 ! E. ANDERSEN FROM FLERCHINGER
1158 ! ----------------------------------------------------------------------
1159 ! SNCOND=0.021+2.51*DSNOW**2
1160 ! CSNOW=0.021+2.51*DSNOW**2
1162 ! SNCOND = UNIT * C
1163 ! double snow thermal conductivity
1164 SNCOND = 2.0 * UNIT * C
1166 ! ----------------------------------------------------------------------
1167 END SUBROUTINE CSNOW
1168 ! ----------------------------------------------------------------------
1170 SUBROUTINE DEVAP (EDIR,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
1171 DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1173 ! ----------------------------------------------------------------------
1174 ! SUBROUTINE DEVAP
1175 ! FUNCTION DEVAP
1176 ! ----------------------------------------------------------------------
1177 ! CALCULATE DIRECT SOIL EVAPORATION
1178 ! ----------------------------------------------------------------------
1179 IMPLICIT NONE
1180 REAL, INTENT(IN) :: ETP1,SMC,BEXP,DKSAT,DWSAT,FXEXP, &
1181 SHDFAC,SMCDRY,SMCMAX,ZSOIL,SMCREF,SMCWLT
1182 REAL, INTENT(OUT):: EDIR
1183 REAL :: FX, SRATIO
1186 ! ----------------------------------------------------------------------
1187 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1188 ! WHEN FXEXP=1.
1189 ! ----------------------------------------------------------------------
1190 ! ----------------------------------------------------------------------
1191 ! FX > 1 REPRESENTS DEMAND CONTROL
1192 ! FX < 1 REPRESENTS FLUX CONTROL
1193 ! ----------------------------------------------------------------------
1195 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1196 IF (SRATIO > 0.) THEN
1197 FX = SRATIO**FXEXP
1198 FX = MAX ( MIN ( FX, 1. ) ,0. )
1199 ELSE
1200 FX = 0.
1201 ENDIF
1203 ! ----------------------------------------------------------------------
1204 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1205 ! ----------------------------------------------------------------------
1206 EDIR = FX * ( 1.0- SHDFAC ) * ETP1
1208 ! ----------------------------------------------------------------------
1209 END SUBROUTINE DEVAP
1210 ! ----------------------------------------------------------------------
1212 SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1213 SH2O, &
1214 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1215 SMCREF,SHDFAC,CMCMAX, &
1216 SMCDRY,CFACTR, &
1217 EDIR,EC,ET,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1219 ! ----------------------------------------------------------------------
1220 ! SUBROUTINE EVAPO
1221 ! ----------------------------------------------------------------------
1222 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
1223 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1224 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1225 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
1226 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1227 ! ----------------------------------------------------------------------
1228 IMPLICIT NONE
1229 INTEGER, INTENT(IN) :: NSOIL, NROOT
1230 INTEGER :: I,K
1231 REAL, INTENT(IN) :: BEXP, CFACTR,CMC,CMCMAX,DKSAT, &
1232 DT,DWSAT,ETP1,FXEXP,PC,Q2,SFCTMP, &
1233 SHDFAC,SMCDRY,SMCMAX,SMCREF,SMCWLT
1234 REAL, INTENT(OUT) :: EC,EDIR,ETA1,ETT
1235 REAL :: CMC2MS
1236 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: RTDIS, SMC, SH2O, ZSOIL
1237 REAL,DIMENSION(1:NSOIL), INTENT(OUT) :: ET
1239 ! ----------------------------------------------------------------------
1240 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1241 ! GREATER THAN ZERO.
1242 ! ----------------------------------------------------------------------
1243 EDIR = 0.
1244 EC = 0.
1245 ETT = 0.
1246 DO K = 1,NSOIL
1247 ET (K) = 0.
1248 END DO
1250 ! ----------------------------------------------------------------------
1251 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
1252 ! ONLY IF VEG COVER NOT COMPLETE.
1253 ! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
1254 ! ----------------------------------------------------------------------
1255 IF (ETP1 > 0.0) THEN
1256 IF (SHDFAC < 1.) THEN
1257 CALL DEVAP (EDIR,ETP1,SMC (1),ZSOIL (1),SHDFAC,SMCMAX, &
1258 BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1259 END IF
1260 ! ----------------------------------------------------------------------
1261 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1262 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1263 ! ----------------------------------------------------------------------
1265 IF (SHDFAC > 0.0) THEN
1266 CALL TRANSP (ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
1267 CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1268 DO K = 1,NSOIL
1269 ETT = ETT + ET ( K )
1270 END DO
1271 ! ----------------------------------------------------------------------
1272 ! CALCULATE CANOPY EVAPORATION.
1273 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1274 ! ----------------------------------------------------------------------
1275 IF (CMC > 0.0) THEN
1276 EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1277 ELSE
1278 EC = 0.0
1279 END IF
1280 ! ----------------------------------------------------------------------
1281 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1282 ! CANOPY. -F.CHEN, 18-OCT-1994
1283 ! ----------------------------------------------------------------------
1284 CMC2MS = CMC / DT
1285 EC = MIN ( CMC2MS, EC )
1286 END IF
1287 END IF
1288 ! ----------------------------------------------------------------------
1289 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1290 ! ----------------------------------------------------------------------
1291 ETA1 = EDIR + ETT + EC
1293 ! ----------------------------------------------------------------------
1294 END SUBROUTINE EVAPO
1295 ! ----------------------------------------------------------------------
1297 SUBROUTINE FAC2MIT(SMCMAX,FLIMIT)
1298 IMPLICIT NONE
1299 REAL, INTENT(IN) :: SMCMAX
1300 REAL, INTENT(OUT) :: FLIMIT
1302 FLIMIT = 0.90
1304 IF ( SMCMAX == 0.395 ) THEN
1305 FLIMIT = 0.59
1306 ELSE IF ( ( SMCMAX == 0.434 ) .OR. ( SMCMAX == 0.404 ) ) THEN
1307 FLIMIT = 0.85
1308 ELSE IF ( ( SMCMAX == 0.465 ) .OR. ( SMCMAX == 0.406 ) ) THEN
1309 FLIMIT = 0.86
1310 ELSE IF ( ( SMCMAX == 0.476 ) .OR. ( SMCMAX == 0.439 ) ) THEN
1311 FLIMIT = 0.74
1312 ELSE IF ( ( SMCMAX == 0.200 ) .OR. ( SMCMAX == 0.464 ) ) THEN
1313 FLIMIT = 0.80
1314 ENDIF
1316 ! ----------------------------------------------------------------------
1317 END SUBROUTINE FAC2MIT
1318 ! ----------------------------------------------------------------------
1320 SUBROUTINE FRH2O (FREE,TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
1322 ! ----------------------------------------------------------------------
1323 ! SUBROUTINE FRH2O
1324 ! ----------------------------------------------------------------------
1325 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
1326 ! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
1327 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
1328 ! (1999, JGR, VOL 104(D16), 19569-19585).
1329 ! ----------------------------------------------------------------------
1330 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
1331 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
1332 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
1333 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
1334 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
1335 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
1336 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
1337 ! ----------------------------------------------------------------------
1338 ! INPUT:
1340 ! TKELV.........TEMPERATURE (Kelvin)
1341 ! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
1342 ! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
1343 ! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
1344 ! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
1345 ! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
1347 ! OUTPUT:
1348 ! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
1349 ! FREE..........SUPERCOOLED LIQUID WATER CONTENT
1350 ! ----------------------------------------------------------------------
1351 IMPLICIT NONE
1352 REAL, INTENT(IN) :: BEXP,PSIS,SH2O,SMC,SMCMAX,TKELV
1353 REAL, INTENT(OUT) :: FREE
1354 REAL :: BX,DENOM,DF,DSWL,FK,SWL,SWLK
1355 INTEGER :: NLOG,KCOUNT
1356 ! PARAMETER(CK = 0.0)
1357 REAL, PARAMETER :: CK = 8.0, BLIM = 5.5, ERROR = 0.005, &
1358 HLICE = 3.335E5, GS = 9.81,DICE = 920.0, &
1359 DH2O = 1000.0, T0 = 273.15
1361 ! ----------------------------------------------------------------------
1362 ! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
1363 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
1364 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
1365 ! ----------------------------------------------------------------------
1366 BX = BEXP
1368 ! ----------------------------------------------------------------------
1369 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
1370 ! ----------------------------------------------------------------------
1371 IF (BEXP > BLIM) BX = BLIM
1372 NLOG = 0
1374 ! ----------------------------------------------------------------------
1375 ! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
1376 ! ----------------------------------------------------------------------
1377 KCOUNT = 0
1378 ! FRH2O = SMC
1379 IF (TKELV > (T0- 1.E-3)) THEN
1380 FREE = SMC
1381 ELSE
1383 ! ----------------------------------------------------------------------
1384 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
1385 ! IN KOREN ET AL, JGR, 1999, EQN 17
1386 ! ----------------------------------------------------------------------
1387 ! INITIAL GUESS FOR SWL (frozen content)
1388 ! ----------------------------------------------------------------------
1389 IF (CK /= 0.0) THEN
1390 SWL = SMC - SH2O
1391 ! ----------------------------------------------------------------------
1392 ! KEEP WITHIN BOUNDS.
1393 ! ----------------------------------------------------------------------
1394 IF (SWL > (SMC -0.02)) SWL = SMC -0.02
1396 ! ----------------------------------------------------------------------
1397 ! START OF ITERATIONS
1398 ! ----------------------------------------------------------------------
1399 IF (SWL < 0.) SWL = 0.
1400 1001 Continue
1401 IF (.NOT.( (NLOG < 10) .AND. (KCOUNT == 0))) goto 1002
1402 NLOG = NLOG +1
1403 DF = ALOG ( ( PSIS * GS / HLICE ) * ( ( 1. + CK * SWL )**2.) * &
1404 ( SMCMAX / (SMC - SWL) )** BX) - ALOG ( - ( &
1405 TKELV - T0)/ TKELV)
1406 DENOM = 2. * CK / ( 1. + CK * SWL ) + BX / ( SMC - SWL )
1407 SWLK = SWL - DF / DENOM
1408 ! ----------------------------------------------------------------------
1409 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
1410 ! ----------------------------------------------------------------------
1411 IF (SWLK > (SMC -0.02)) SWLK = SMC - 0.02
1412 IF (SWLK < 0.) SWLK = 0.
1414 ! ----------------------------------------------------------------------
1415 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
1416 ! ----------------------------------------------------------------------
1417 DSWL = ABS (SWLK - SWL)
1419 ! ----------------------------------------------------------------------
1420 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
1421 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
1422 ! ----------------------------------------------------------------------
1423 SWL = SWLK
1424 IF ( DSWL <= ERROR ) THEN
1425 KCOUNT = KCOUNT +1
1426 END IF
1427 ! ----------------------------------------------------------------------
1428 ! END OF ITERATIONS
1429 ! ----------------------------------------------------------------------
1430 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
1431 ! ----------------------------------------------------------------------
1432 ! FRH2O = SMC - SWL
1433 goto 1001
1434 1002 continue
1435 FREE = SMC - SWL
1436 END IF
1437 ! ----------------------------------------------------------------------
1438 ! END OPTION 1
1439 ! ----------------------------------------------------------------------
1440 ! ----------------------------------------------------------------------
1441 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
1442 ! IN KOREN ET AL., JGR, 1999, EQN 17
1443 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
1444 ! ----------------------------------------------------------------------
1445 IF (KCOUNT == 0) THEN
1446 PRINT *,'Flerchinger USEd in NEW version. Iterations=',NLOG
1447 FK = ( ( (HLICE / (GS * ( - PSIS)))* &
1448 ( (TKELV - T0)/ TKELV))** ( -1/ BX))* SMCMAX
1449 ! FRH2O = MIN (FK, SMC)
1450 IF (FK < 0.02) FK = 0.02
1451 FREE = MIN (FK, SMC)
1452 ! ----------------------------------------------------------------------
1453 ! END OPTION 2
1454 ! ----------------------------------------------------------------------
1455 END IF
1456 END IF
1457 ! ----------------------------------------------------------------------
1458 END SUBROUTINE FRH2O
1459 ! ----------------------------------------------------------------------
1461 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
1462 TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
1463 F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP,ISURBAN)
1465 ! ----------------------------------------------------------------------
1466 ! SUBROUTINE HRT
1467 ! ----------------------------------------------------------------------
1468 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1469 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1470 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1471 ! ----------------------------------------------------------------------
1472 IMPLICIT NONE
1473 LOGICAL :: ITAVG
1474 INTEGER, INTENT(IN) :: NSOIL, VEGTYP
1475 INTEGER, INTENT(IN) :: ISURBAN
1476 INTEGER :: I, K
1478 REAL, INTENT(IN) :: BEXP, CSOIL, DF1, DT,F1,PSISAT,QUARTZ, &
1479 SMCMAX ,TBOT,YY,ZZ1, ZBOT
1480 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,STC,ZSOIL
1481 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SH2O
1482 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
1483 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
1484 REAL :: DDZ, DDZ2, DENOM, DF1N, DF1K, DTSDZ, &
1485 DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK, &
1486 TBK1,TSNSR,TSURF,CSOIL_LOC
1487 REAL, PARAMETER :: T0 = 273.15, CAIR = 1004.0, CICE = 2.106E6,&
1488 CH2O = 4.2E6
1491 !urban
1492 IF( VEGTYP == ISURBAN ) then
1493 CSOIL_LOC=3.0E6
1494 ELSE
1495 CSOIL_LOC=CSOIL
1496 ENDIF
1498 ! ----------------------------------------------------------------------
1499 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1500 ! ----------------------------------------------------------------------
1501 ITAVG = .TRUE.
1502 ! ----------------------------------------------------------------------
1503 ! BEGIN SECTION FOR TOP SOIL LAYER
1504 ! ----------------------------------------------------------------------
1505 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1506 ! ----------------------------------------------------------------------
1507 HCPCT = SH2O (1)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (1))&
1508 * CAIR &
1509 + ( SMC (1) - SH2O (1) )* CICE
1511 ! ----------------------------------------------------------------------
1512 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1513 ! ----------------------------------------------------------------------
1514 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
1515 AI (1) = 0.0
1516 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
1518 ! ----------------------------------------------------------------------
1519 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1520 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1521 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1522 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1523 ! ----------------------------------------------------------------------
1524 BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT * &
1525 ZZ1)
1526 DTSDZ = (STC (1) - STC (2)) / ( -0.5 * ZSOIL (2))
1527 SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
1528 ! RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1529 DENOM = (ZSOIL (1) * HCPCT)
1531 ! ----------------------------------------------------------------------
1532 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1533 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1534 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1535 ! ----------------------------------------------------------------------
1536 ! QTOT = SSOIL - DF1*DTSDZ
1537 RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
1539 ! ----------------------------------------------------------------------
1540 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1541 ! ----------------------------------------------------------------------
1542 QTOT = -1.0* RHSTS (1)* DENOM
1544 ! ----------------------------------------------------------------------
1545 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1546 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1547 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
1548 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
1549 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1550 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
1551 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1552 ! LATER IN FUNCTION SUBROUTINE SNKSRC
1553 ! ----------------------------------------------------------------------
1554 SICE = SMC (1) - SH2O (1)
1555 IF (ITAVG) THEN
1556 TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
1557 ! ----------------------------------------------------------------------
1558 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1559 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1560 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1561 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1562 ! ----------------------------------------------------------------------
1563 CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
1564 IF ( (SICE > 0.) .OR. (STC (1) < T0) .OR. &
1565 (TSURF < T0) .OR. (TBK < T0) ) THEN
1566 ! TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1),
1567 CALL TMPAVG (TAVG,TSURF,STC (1),TBK,ZSOIL,NSOIL,1)
1568 CALL SNKSRC (TSNSR,TAVG,SMC (1),SH2O (1), &
1569 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1570 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1571 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
1572 END IF
1573 ELSE
1574 ! TSNSR = SNKSRC (STC(1),SMC(1),SH2O(1),
1575 IF ( (SICE > 0.) .OR. (STC (1) < T0) ) THEN
1576 CALL SNKSRC (TSNSR,STC (1),SMC (1),SH2O (1), &
1577 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1578 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1579 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
1580 END IF
1581 ! ----------------------------------------------------------------------
1582 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
1583 ! ----------------------------------------------------------------------
1584 END IF
1586 ! INITIALIZE DDZ2
1587 ! ----------------------------------------------------------------------
1589 DDZ2 = 0.0
1590 DF1K = DF1
1592 ! ----------------------------------------------------------------------
1593 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1594 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
1595 ! ----------------------------------------------------------------------
1596 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
1597 ! ----------------------------------------------------------------------
1598 DO K = 2,NSOIL
1599 HCPCT = SH2O (K)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC ( &
1600 K))* CAIR + ( SMC (K) - SH2O (K) )* CICE
1601 ! ----------------------------------------------------------------------
1602 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
1603 ! ----------------------------------------------------------------------
1604 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
1605 ! ----------------------------------------------------------------------
1606 IF (K /= NSOIL) THEN
1608 ! ----------------------------------------------------------------------
1609 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
1610 ! ----------------------------------------------------------------------
1611 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
1613 !urban
1614 IF ( VEGTYP == ISURBAN ) DF1N = 3.24
1616 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
1618 ! ----------------------------------------------------------------------
1619 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
1620 ! ----------------------------------------------------------------------
1621 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
1622 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
1624 ! ----------------------------------------------------------------------
1625 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
1626 ! TEMP AT BOTTOM OF LAYER.
1627 ! ----------------------------------------------------------------------
1628 CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
1629 HCPCT)
1630 IF (ITAVG) THEN
1631 CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
1632 END IF
1634 ELSE
1635 ! ----------------------------------------------------------------------
1636 ! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
1637 ! BOTTOM LAYER.
1638 ! ----------------------------------------------------------------------
1640 ! ----------------------------------------------------------------------
1641 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
1642 ! ----------------------------------------------------------------------
1643 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
1646 !urban
1647 IF ( VEGTYP == ISURBAN ) DF1N = 3.24
1649 DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
1651 ! ----------------------------------------------------------------------
1652 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
1653 ! ----------------------------------------------------------------------
1654 DTSDZ2 = (STC (K) - TBOT) / DENOM
1656 ! ----------------------------------------------------------------------
1657 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
1658 ! TEMP AT BOTTOM OF LAST LAYER.
1659 ! ----------------------------------------------------------------------
1660 CI (K) = 0.
1661 IF (ITAVG) THEN
1662 CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
1663 END IF
1664 ! ----------------------------------------------------------------------
1665 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
1666 END IF
1667 ! ----------------------------------------------------------------------
1668 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
1669 ! ----------------------------------------------------------------------
1670 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
1671 RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
1672 QTOT = -1.0* DENOM * RHSTS (K)
1674 SICE = SMC (K) - SH2O (K)
1675 IF (ITAVG) THEN
1676 CALL TMPAVG (TAVG,TBK,STC (K),TBK1,ZSOIL,NSOIL,K)
1677 ! TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL,
1678 IF ( (SICE > 0.) .OR. (STC (K) < T0) .OR. &
1679 (TBK .lt. T0) .OR. (TBK1 .lt. T0) ) THEN
1680 CALL SNKSRC (TSNSR,TAVG,SMC (K),SH2O (K),ZSOIL,NSOIL, &
1681 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
1682 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
1683 END IF
1684 ELSE
1685 ! TSNSR = SNKSRC(STC(K),SMC(K),SH2O(K),ZSOIL,NSOIL,
1686 IF ( (SICE > 0.) .OR. (STC (K) < T0) ) THEN
1687 CALL SNKSRC (TSNSR,STC (K),SMC (K),SH2O (K),ZSOIL,NSOIL, &
1688 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
1689 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
1690 END IF
1691 END IF
1693 ! ----------------------------------------------------------------------
1694 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
1695 ! ----------------------------------------------------------------------
1696 AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
1698 ! ----------------------------------------------------------------------
1699 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
1700 ! ----------------------------------------------------------------------
1701 BI (K) = - (AI (K) + CI (K))
1702 TBK = TBK1
1703 DF1K = DF1N
1704 DTSDZ = DTSDZ2
1705 DDZ = DDZ2
1706 END DO
1707 ! ----------------------------------------------------------------------
1708 END SUBROUTINE HRT
1709 ! ----------------------------------------------------------------------
1711 SUBROUTINE HRTICE_GLACIAL (RHSTS,STC,TBOT,NSOIL,ZSOIL,YY,ZZ1,DF1, &
1712 AI,BI,CI)
1714 ! ----------------------------------------------------------------------
1715 ! SUBROUTINE HRTICE_GLACIAL
1716 ! ----------------------------------------------------------------------
1717 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1718 ! THERMAL DIFFUSION EQUATION IN THE CASE OF GLACIAL ICE (ICE == -1).
1719 ! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
1720 ! OF THE IMPLICIT TIME SCHEME.
1721 !
1722 ! (NOTE: THIS SUBROUTINE ONLY CALLED FOR GLACIAL ICE (ICE == -1), BUT
1723 ! NOT FOR NON-GLACIAL LAND (ICE == 0).
1724 ! ----------------------------------------------------------------------
1725 IMPLICIT NONE
1728 INTEGER, INTENT(IN) :: NSOIL
1729 INTEGER :: K
1731 REAL, INTENT(IN) :: DF1,YY,ZZ1
1732 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: AI, BI,CI
1733 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: STC, ZSOIL
1734 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: RHSTS
1735 REAL, INTENT(IN) :: TBOT
1736 REAL :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL, &
1737 ZBOT
1738 REAL :: HCPCT
1739 REAL :: DF1K
1740 REAL :: DF1N
1741 REAL :: ZMD
1743 ! SET A NOMINAL UNIVERSAL VALUE OF GLACIAL ICE SPECIFIC HEAT CAPACITY,
1744 ! HCPCT = 2100.0*900.0 = 1.89000E+6 (SOURCE: BOB GRUMBINE, 2005)
1745 ! TBOT PASSED IN AS ARGUMENT, VALUE FROM GLOBAL DATA SET
1747 !
1748 ! A least-squares fit for the four points provided by
1749 ! Keith Hines for the Yen (1981) values for Antarctic
1750 ! snow firn.
1751 !
1752 HCPCT = 1.E6 * (0.8194 - 0.1309*0.5*ZSOIL(1))
1753 DF1K = DF1
1755 ! ----------------------------------------------------------------------
1756 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1757 ! ----------------------------------------------------------------------
1758 ZBOT = -25.0
1759 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
1760 AI (1) = 0.0
1761 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
1763 ! ----------------------------------------------------------------------
1764 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
1765 ! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
1766 ! RHSTS FOR THE TOP SOIL LAYER.
1767 ! ----------------------------------------------------------------------
1768 BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT * &
1769 ZZ1)
1770 DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
1771 SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
1773 ! ----------------------------------------------------------------------
1774 ! INITIALIZE DDZ2
1775 ! ----------------------------------------------------------------------
1776 RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
1778 ! ----------------------------------------------------------------------
1779 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1780 ! ----------------------------------------------------------------------
1781 DDZ2 = 0.0
1782 DF1K = DF1
1783 DF1N = DF1
1784 DO K = 2,NSOIL
1786 ZMD = 0.5 * (ZSOIL(K)+ZSOIL(K-1))
1787 ! For the land-ice case
1788 ! kmh 09/03/2006 use Yen (1981)'s values for Antarctic snow firn
1789 ! IF ( K .eq. 2 ) HCPCT = 0.855108E6
1790 ! IF ( K .eq. 3 ) HCPCT = 0.922906E6
1791 ! IF ( K .eq. 4 ) HCPCT = 1.009986E6
1793 ! Least squares fit to the four points supplied by Keith Hines
1794 ! from Yen (1981) for Antarctic snow firn. Not optimal, but
1795 ! probably better than just a constant.
1796 HCPCT = 1.E6 * ( 0.8194 - 0.1309*ZMD )
1798 ! IF ( K .eq. 2 ) DF1N = 0.345356
1799 ! IF ( K .eq. 3 ) DF1N = 0.398777
1800 ! IF ( K .eq. 4 ) DF1N = 0.472653
1802 ! Least squares fit to the three points supplied by Keith Hines
1803 ! from Yen (1981) for Antarctic snow firn. Not optimal, but
1804 ! probably better than just a constant.
1805 DF1N = 0.32333 - ( 0.10073 * ZMD )
1807 ! ----------------------------------------------------------------------
1808 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
1809 ! ----------------------------------------------------------------------
1810 IF (K /= NSOIL) THEN
1811 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
1813 ! ----------------------------------------------------------------------
1814 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
1815 ! ----------------------------------------------------------------------
1816 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
1817 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
1818 CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
1820 ! ----------------------------------------------------------------------
1821 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
1822 ! ----------------------------------------------------------------------
1823 ELSE
1825 ! ----------------------------------------------------------------------
1826 ! SET MATRIX COEF, CI TO ZERO.
1827 ! ----------------------------------------------------------------------
1828 DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
1829 - ZBOT)
1830 CI (K) = 0.
1831 ! ----------------------------------------------------------------------
1832 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
1833 ! ----------------------------------------------------------------------
1834 END IF
1835 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
1837 ! ----------------------------------------------------------------------
1838 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
1839 ! ----------------------------------------------------------------------
1840 RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
1841 AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
1843 ! ----------------------------------------------------------------------
1844 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
1845 ! ----------------------------------------------------------------------
1846 BI (K) = - (AI (K) + CI (K))
1847 DF1K = DF1N
1848 DTSDZ = DTSDZ2
1849 DDZ = DDZ2
1850 END DO
1851 ! ----------------------------------------------------------------------
1852 END SUBROUTINE HRTICE_GLACIAL
1853 ! ----------------------------------------------------------------------
1855 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
1857 ! ----------------------------------------------------------------------
1858 ! SUBROUTINE HSTEP
1859 ! ----------------------------------------------------------------------
1860 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
1861 ! ----------------------------------------------------------------------
1862 IMPLICIT NONE
1863 INTEGER, INTENT(IN) :: NSOIL
1864 INTEGER :: K
1866 REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
1867 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
1868 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
1869 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
1870 REAL, DIMENSION(1:NSOIL) :: RHSTSin
1871 REAL, DIMENSION(1:NSOIL) :: CIin
1872 REAL :: DT
1874 ! ----------------------------------------------------------------------
1875 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
1876 ! ----------------------------------------------------------------------
1877 DO K = 1,NSOIL
1878 RHSTS (K) = RHSTS (K) * DT
1879 AI (K) = AI (K) * DT
1880 BI (K) = 1. + BI (K) * DT
1881 CI (K) = CI (K) * DT
1882 END DO
1883 ! ----------------------------------------------------------------------
1884 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
1885 ! ----------------------------------------------------------------------
1886 DO K = 1,NSOIL
1887 RHSTSin (K) = RHSTS (K)
1888 END DO
1889 DO K = 1,NSOIL
1890 CIin (K) = CI (K)
1891 END DO
1892 ! ----------------------------------------------------------------------
1893 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
1894 ! ----------------------------------------------------------------------
1895 CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
1896 ! ----------------------------------------------------------------------
1897 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
1898 ! ----------------------------------------------------------------------
1899 DO K = 1,NSOIL
1900 STCOUT (K) = STCIN (K) + CI (K)
1901 END DO
1902 ! ----------------------------------------------------------------------
1903 END SUBROUTINE HSTEP
1904 ! ----------------------------------------------------------------------
1906 SUBROUTINE NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1907 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
1908 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
1909 SSOIL, &
1910 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1911 SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
1912 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1913 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1914 QUARTZ,FXEXP,CSOIL, &
1915 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
1917 ! ----------------------------------------------------------------------
1918 ! SUBROUTINE NOPAC
1919 ! ----------------------------------------------------------------------
1920 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
1921 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
1922 ! PRESENT.
1923 ! ----------------------------------------------------------------------
1924 IMPLICIT NONE
1926 INTEGER, INTENT(IN) :: ICE, NROOT,NSOIL,VEGTYP
1927 INTEGER, INTENT(IN) :: ISURBAN
1928 INTEGER :: K
1930 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DKSAT,DT,DWSAT, &
1931 EPSCA,ETP,FDOWN,F1,FXEXP,FRZFACT,KDT,PC, &
1932 PRCP,PSISAT,Q2,QUARTZ,RCH,RR,SBETA,SFCTMP,&
1933 SHDFAC,SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, &
1934 T24,TBOT,TH2,ZBOT,EMISSI
1935 REAL, INTENT(INOUT) :: CMC,BETA,T1
1936 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR,ETA,ETT,FLX1,FLX3, &
1937 RUNOFF1,RUNOFF2,RUNOFF3,SSOIL
1938 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
1939 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
1940 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
1941 REAL, DIMENSION(1:NSOIL) :: ET1
1942 REAL :: EC1,EDIR1,ETT1,DF1,ETA1,ETP1,PRCP1,YY, &
1943 YYNUM,ZZ1
1945 ! ----------------------------------------------------------------------
1946 ! EXECUTABLE CODE BEGINS HERE:
1947 ! CONVERT ETP Fnd PRCP FROM KG M-2 S-1 TO M S-1 AND INITIALIZE DEW.
1948 ! ----------------------------------------------------------------------
1949 PRCP1 = PRCP * 0.001
1950 ETP1 = ETP * 0.001
1951 DEW = 0.0
1952 ! ----------------------------------------------------------------------
1953 ! INITIALIZE EVAP TERMS.
1954 ! ----------------------------------------------------------------------
1955 EDIR = 0.
1956 EDIR1 = 0.
1957 EC1 = 0.
1958 EC = 0.
1959 DO K = 1,NSOIL
1960 ET(K) = 0.
1961 ET1(K) = 0.
1962 END DO
1963 ETT = 0.
1964 ETT1 = 0.
1966 IF (ETP > 0.0) THEN
1967 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1968 SH2O, &
1969 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1970 SMCREF,SHDFAC,CMCMAX, &
1971 SMCDRY,CFACTR, &
1972 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1973 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
1974 SH2O,SLOPE,KDT,FRZFACT, &
1975 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
1976 SHDFAC,CMCMAX, &
1977 RUNOFF1,RUNOFF2,RUNOFF3, &
1978 EDIR1,EC1,ET1, &
1979 DRIP)
1981 ! ----------------------------------------------------------------------
1982 ! CONVERT MODELED EVAPOTRANSPIRATION FROM M S-1 TO KG M-2 S-1.
1983 ! ----------------------------------------------------------------------
1985 ETA = ETA1 * 1000.0
1987 ! ----------------------------------------------------------------------
1988 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
1989 ! ETP1 TO ZERO).
1990 ! ----------------------------------------------------------------------
1991 ELSE
1992 DEW = - ETP1
1994 ! ----------------------------------------------------------------------
1995 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
1996 ! ----------------------------------------------------------------------
1998 PRCP1 = PRCP1+ DEW
1999 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2000 SH2O,SLOPE,KDT,FRZFACT, &
2001 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2002 SHDFAC,CMCMAX, &
2003 RUNOFF1,RUNOFF2,RUNOFF3, &
2004 EDIR1,EC1,ET1, &
2005 DRIP)
2007 ! ----------------------------------------------------------------------
2008 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2009 ! ----------------------------------------------------------------------
2010 ! ETA = ETA1 * 1000.0
2011 END IF
2013 ! ----------------------------------------------------------------------
2014 ! BASED ON ETP AND E VALUES, DETERMINE BETA
2015 ! ----------------------------------------------------------------------
2017 IF ( ETP <= 0.0 ) THEN
2018 BETA = 0.0
2019 ETA = ETP
2020 IF ( ETP < 0.0 ) THEN
2021 BETA = 1.0
2022 END IF
2023 ELSE
2024 BETA = ETA / ETP
2025 END IF
2027 ! ----------------------------------------------------------------------
2028 ! CONVERT MODELED EVAPOTRANSPIRATION COMPONENTS 'M S-1' TO 'KG M-2 S-1'.
2029 ! ----------------------------------------------------------------------
2030 EDIR = EDIR1*1000.
2031 EC = EC1*1000.
2032 DO K = 1,NSOIL
2033 ET(K) = ET1(K)*1000.
2034 END DO
2035 ETT = ETT1*1000.
2037 ! ----------------------------------------------------------------------
2038 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2039 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2040 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2041 ! ----------------------------------------------------------------------
2043 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
2045 !urban
2046 IF ( VEGTYP == ISURBAN ) DF1=3.24
2047 !
2049 ! ----------------------------------------------------------------------
2050 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2051 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2052 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2053 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2054 ! ROUTINE SFLX)
2055 ! ----------------------------------------------------------------------
2056 DF1 = DF1 * EXP (SBETA * SHDFAC)
2057 ! ----------------------------------------------------------------------
2058 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2059 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2060 ! ----------------------------------------------------------------------
2061 YYNUM = FDOWN - EMISSI*SIGMA * T24
2062 YY = SFCTMP + (YYNUM / RCH + TH2- SFCTMP - BETA * EPSCA) / RR
2064 ZZ1 = DF1 / ( -0.5 * ZSOIL (1) * RCH * RR ) + 1.0
2066 !urban
2067 CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2068 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2069 QUARTZ,CSOIL,VEGTYP,ISURBAN)
2071 ! ----------------------------------------------------------------------
2072 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2073 ! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
2074 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2075 ! ----------------------------------------------------------------------
2076 FLX1 = CPH2O * PRCP * (T1- SFCTMP)
2077 FLX3 = 0.0
2079 ! ----------------------------------------------------------------------
2080 END SUBROUTINE NOPAC
2081 ! ----------------------------------------------------------------------
2083 SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2084 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
2085 & DQSDT2,FLX2,EMISSI_IN,SNEQV,T1,ICE,SNCOVR)
2087 ! ----------------------------------------------------------------------
2088 ! SUBROUTINE PENMAN
2089 ! ----------------------------------------------------------------------
2090 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
2091 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2092 ! CALLING ROUTINE FOR LATER USE.
2093 ! ----------------------------------------------------------------------
2094 IMPLICIT NONE
2095 LOGICAL, INTENT(IN) :: SNOWNG, FRZGRA
2096 REAL, INTENT(IN) :: CH, DQSDT2,FDOWN,PRCP, &
2097 Q2, Q2SAT,SSOIL, SFCPRS, SFCTMP, &
2098 T2V, TH2,EMISSI_IN,SNEQV
2099 REAL, INTENT(IN) :: T1 , SNCOVR
2100 !
2101 ! kmh 09/13/2006
2102 INTEGER, INTENT(IN) :: ICE
2103 ! kmh 09/03/2006
2104 !
2105 REAL, INTENT(OUT) :: EPSCA,ETP,FLX2,RCH,RR,T24
2106 REAL :: A, DELTA, FNET,RAD,RHO,EMISSI,ELCP1,LVS
2108 REAL, PARAMETER :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6,CP = 1004.6
2109 REAL, PARAMETER :: LSUBS = 2.83E+6
2111 ! ----------------------------------------------------------------------
2112 ! EXECUTABLE CODE BEGINS HERE:
2113 ! ----------------------------------------------------------------------
2114 ! ----------------------------------------------------------------------
2115 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2116 ! ----------------------------------------------------------------------
2117 EMISSI=EMISSI_IN
2118 IF ( ICE == 0 ) THEN
2119 ! Non-glacial land
2120 ELCP1 = (1.0-SNCOVR)*ELCP + SNCOVR*ELCP*LSUBS/LSUBC
2121 LVS = (1.0-SNCOVR)*LSUBC + SNCOVR*LSUBS
2122 ELSEIF ( ICE == -1 ) THEN
2123 ! Glacial ice
2124 IF ( T1 > 273.15 ) THEN
2125 ELCP1=ELCP
2126 LVS=LSUBC
2127 ELSE
2128 ELCP1 = ELCP*LSUBS/LSUBC
2129 LVS = LSUBS
2130 ENDIF
2131 ENDIF
2133 FLX2 = 0.0
2134 ! DELTA = ELCP * DQSDT2
2135 DELTA = ELCP1 * DQSDT2
2136 T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2137 ! RR = T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
2138 RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
2139 RHO = SFCPRS / (RD * T2V)
2141 ! ----------------------------------------------------------------------
2142 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2143 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
2144 ! ----------------------------------------------------------------------
2145 RCH = RHO * CP * CH
2146 IF (.NOT. SNOWNG) THEN
2147 IF (PRCP > 0.0) RR = RR + CPH2O * PRCP / RCH
2148 ELSE
2149 RR = RR + CPICE * PRCP / RCH
2150 END IF
2152 ! ----------------------------------------------------------------------
2153 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2154 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2155 ! ----------------------------------------------------------------------
2156 ! FNET = FDOWN - SIGMA * T24- SSOIL
2157 FNET = FDOWN - EMISSI*SIGMA * T24- SSOIL
2158 IF (FRZGRA) THEN
2159 FLX2 = - LSUBF * PRCP
2160 FNET = FNET - FLX2
2161 ! ----------------------------------------------------------------------
2162 ! FINISH PENMAN EQUATION CALCULATIONS.
2163 ! ----------------------------------------------------------------------
2164 END IF
2165 RAD = FNET / RCH + TH2- SFCTMP
2166 ! A = ELCP * (Q2SAT - Q2)
2167 A = ELCP1 * (Q2SAT - Q2)
2168 EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
2169 ! ETP = EPSCA * RCH / LSUBC
2170 ETP = EPSCA * RCH / LVS
2172 ! ----------------------------------------------------------------------
2173 END SUBROUTINE PENMAN
2174 ! ----------------------------------------------------------------------
2176 SUBROUTINE REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX, &
2177 TOPT, &
2178 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
2179 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
2180 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
2181 RTDIS,SLDPTH,ZSOIL, NROOT,NSOIL,CZIL, &
2182 LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN, &
2183 ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE, &
2184 LSOIL, LOCAL,LVCOEF)
2186 IMPLICIT NONE
2187 ! ----------------------------------------------------------------------
2188 ! Internally set (default valuess)
2189 ! all soil and vegetation parameters required for the execusion oF
2190 ! the Noah lsm are defined in VEGPARM.TBL, SOILPARM.TB, and GENPARM.TBL.
2191 ! ----------------------------------------------------------------------
2192 ! Vegetation parameters:
2193 ! ALBBRD: SFC background snow-free albedo
2194 ! CMXTBL: MAX CNPY Capacity
2195 ! Z0BRD: Background roughness length
2196 ! SHDFAC: Green vegetation fraction
2197 ! NROOT: Rooting depth
2198 ! RSMIN: Mimimum stomatal resistance
2199 ! RSMAX: Max. stomatal resistance
2200 ! RGL: Parameters used in radiation stress function
2201 ! HS: Parameter used in vapor pressure deficit functio
2202 ! TOPT: Optimum transpiration air temperature.
2203 ! CMCMAX: Maximum canopy water capacity
2204 ! CFACTR: Parameter used in the canopy inteception calculation
2205 ! SNUP: Threshold snow depth (in water equivalent m) that
2206 ! implies 100 percent snow cover
2207 ! LAI: Leaf area index
2208 !
2209 ! ----------------------------------------------------------------------
2210 ! Soil parameters:
2211 ! SMCMAX: MAX soil moisture content (porosity)
2212 ! SMCREF: Reference soil moisture (field capacity)
2213 ! SMCWLT: Wilting point soil moisture
2214 ! SMCWLT: Air dry soil moist content limits
2215 ! SSATPSI: SAT (saturation) soil potential
2216 ! DKSAT: SAT soil conductivity
2217 ! BEXP: B parameter
2218 ! SSATDW: SAT soil diffusivity
2219 ! F1: Soil thermal diffusivity/conductivity coef.
2220 ! QUARTZ: Soil quartz content
2221 ! Modified by F. Chen (12/22/97) to use the STATSGO soil map
2222 ! Modified By F. Chen (01/22/00) to include PLaya, Lava, and White San
2223 ! Modified By F. Chen (08/05/02) to include additional parameters for the Noah
2224 ! NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC)
2225 ! F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
2226 ! REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 m/s= 0.5 mm
2227 ! REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
2228 ! WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB) (Wetzel and Chang, 198
2229 ! WLTSMC=WLTSMC1-0.5*WLTSMC1
2230 ! Note: the values for playa is set for it to have a thermal conductivit
2231 ! as sand and to have a hydrulic conductivity as clay
2232 !
2233 ! ----------------------------------------------------------------------
2234 ! Class parameter 'SLOPETYP' was included to estimate linear reservoir
2235 ! coefficient 'SLOPE' to the baseflow runoff out of the bottom layer.
2236 ! lowest class (slopetyp=0) means highest slope parameter = 1.
2237 ! definition of slopetyp from 'zobler' slope type:
2238 ! slope class percent slope
2239 ! 1 0-8
2240 ! 2 8-30
2241 ! 3 > 30
2242 ! 4 0-30
2243 ! 5 0-8 & > 30
2244 ! 6 8-30 & > 30
2245 ! 7 0-8, 8-30, > 30
2246 ! 9 GLACIAL ICE
2247 ! BLANK OCEAN/SEA
2248 ! SLOPE_DATA: linear reservoir coefficient
2249 ! SBETA_DATA: parameter used to caluculate vegetation effect on soil heat
2250 ! FXEXP_DAT: soil evaporation exponent used in DEVAP
2251 ! CSOIL_DATA: soil heat capacity [J M-3 K-1]
2252 ! SALP_DATA: shape parameter of distribution function of snow cover
2253 ! REFDK_DATA and REFKDT_DATA: parameters in the surface runoff parameteriz
2254 ! FRZK_DATA: frozen ground parameter
2255 ! ZBOT_DATA: depth[M] of lower boundary soil temperature
2256 ! CZIL_DATA: calculate roughness length of heat
2257 ! SMLOW_DATA and MHIGH_DATA: two soil moisture wilt, soil moisture referen
2258 ! parameters
2259 ! Set maximum number of soil-, veg-, and slopetyp in data statement.
2260 ! ----------------------------------------------------------------------
2261 INTEGER, PARAMETER :: MAX_SLOPETYP=30,MAX_SOILTYP=30,MAX_VEGTYP=30
2262 LOGICAL :: LOCAL
2263 CHARACTER (LEN=256), INTENT(IN):: LLANDUSE, LSOIL
2265 ! Veg parameters
2266 INTEGER, INTENT(IN) :: VEGTYP
2267 INTEGER, INTENT(OUT) :: NROOT
2268 REAL, INTENT(INOUT) :: SHDFAC
2269 REAL, INTENT(OUT) :: HS,RSMIN,RGL,SNUP, &
2270 CMCMAX,RSMAX,TOPT, &
2271 EMISSMIN, EMISSMAX, &
2272 LAIMIN, LAIMAX, &
2273 Z0MIN, Z0MAX, &
2274 ALBEDOMIN, ALBEDOMAX
2275 ! Soil parameters
2276 INTEGER, INTENT(IN) :: SOILTYP
2277 REAL, INTENT(OUT) :: BEXP,DKSAT,DWSAT,F1,QUARTZ,SMCDRY, &
2278 SMCMAX,SMCREF,SMCWLT,PSISAT
2279 ! General parameters
2280 INTEGER, INTENT(IN) :: SLOPETYP,NSOIL
2281 INTEGER :: I
2283 REAL, INTENT(OUT) :: SLOPE,CZIL,SBETA,FXEXP, &
2284 CSOIL,SALP,FRZX,KDT,CFACTR, &
2285 ZBOT,REFKDT,PTU
2286 REAL, INTENT(OUT) :: LVCOEF
2287 REAL,DIMENSION(1:NSOIL),INTENT(IN) :: SLDPTH,ZSOIL
2288 REAL,DIMENSION(1:NSOIL),INTENT(OUT):: RTDIS
2289 REAL :: FRZFACT,FRZK,REFDK
2291 ! SAVE
2292 ! ----------------------------------------------------------------------
2293 !
2294 IF (SOILTYP .gt. SLCATS) THEN
2295 CALL wrf_error_fatal ( 'Warning: too many input soil types' )
2296 END IF
2297 IF (VEGTYP .gt. LUCATS) THEN
2298 CALL wrf_error_fatal ( 'Warning: too many input landuse types' )
2299 END IF
2300 IF (SLOPETYP .gt. SLPCATS) THEN
2301 CALL wrf_error_fatal ( 'Warning: too many input slope types' )
2302 END IF
2304 ! ----------------------------------------------------------------------
2305 ! SET-UP SOIL PARAMETERS
2306 ! ----------------------------------------------------------------------
2307 CSOIL = CSOIL_DATA
2308 BEXP = BB (SOILTYP)
2309 DKSAT = SATDK (SOILTYP)
2310 DWSAT = SATDW (SOILTYP)
2311 F1 = F11 (SOILTYP)
2312 PSISAT = SATPSI (SOILTYP)
2313 QUARTZ = QTZ (SOILTYP)
2314 SMCDRY = DRYSMC (SOILTYP)
2315 SMCMAX = MAXSMC (SOILTYP)
2316 SMCREF = REFSMC (SOILTYP)
2317 SMCWLT = WLTSMC (SOILTYP)
2318 ! ----------------------------------------------------------------------
2319 ! Set-up universal parameters (not dependent on SOILTYP, VEGTYP or
2320 ! SLOPETYP)
2321 ! ----------------------------------------------------------------------
2322 ZBOT = ZBOT_DATA
2323 SALP = SALP_DATA
2324 SBETA = SBETA_DATA
2325 REFDK = REFDK_DATA
2326 FRZK = FRZK_DATA
2327 FXEXP = FXEXP_DATA
2328 REFKDT = REFKDT_DATA
2329 PTU = 0. ! (not used yet) to satisify intent(out)
2330 KDT = REFKDT * DKSAT / REFDK
2331 CZIL = CZIL_DATA
2332 SLOPE = SLOPE_DATA (SLOPETYP)
2333 LVCOEF = LVCOEF_DATA
2335 ! ----------------------------------------------------------------------
2336 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
2337 ! ----------------------------------------------------------------------
2338 FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
2339 FRZX = FRZK * FRZFACT
2341 ! ----------------------------------------------------------------------
2342 ! SET-UP VEGETATION PARAMETERS
2343 ! ----------------------------------------------------------------------
2344 TOPT = TOPT_DATA
2345 CMCMAX = CMCMAX_DATA
2346 CFACTR = CFACTR_DATA
2347 RSMAX = RSMAX_DATA
2348 NROOT = NROTBL (VEGTYP)
2349 SNUP = SNUPTBL (VEGTYP)
2350 RSMIN = RSTBL (VEGTYP)
2351 RGL = RGLTBL (VEGTYP)
2352 HS = HSTBL (VEGTYP)
2353 EMISSMIN = EMISSMINTBL (VEGTYP)
2354 EMISSMAX = EMISSMAXTBL (VEGTYP)
2355 LAIMIN = LAIMINTBL (VEGTYP)
2356 LAIMAX = LAIMAXTBL (VEGTYP)
2357 Z0MIN = Z0MINTBL (VEGTYP)
2358 Z0MAX = Z0MAXTBL (VEGTYP)
2359 ALBEDOMIN = ALBEDOMINTBL (VEGTYP)
2360 ALBEDOMAX = ALBEDOMAXTBL (VEGTYP)
2362 IF (VEGTYP .eq. BARE) SHDFAC = 0.0
2363 IF (NROOT .gt. NSOIL) THEN
2364 WRITE (err_message,*) 'Error: too many root layers ', &
2365 NSOIL,NROOT
2366 CALL wrf_error_fatal ( err_message )
2367 ! ----------------------------------------------------------------------
2368 ! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
2369 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
2370 ! ----------------------------------------------------------------------
2371 END IF
2372 DO I = 1,NROOT
2373 RTDIS (I) = - SLDPTH (I)/ ZSOIL (NROOT)
2374 ! ----------------------------------------------------------------------
2375 ! SET-UP SLOPE PARAMETER
2376 ! ----------------------------------------------------------------------
2377 END DO
2379 ! print*,'end of PRMRED'
2380 ! print*,'VEGTYP',VEGTYP,'SOILTYP',SOILTYP,'SLOPETYP',SLOPETYP, &
2381 ! & 'CFACTR',CFACTR,'CMCMAX',CMCMAX,'RSMAX',RSMAX,'TOPT',TOPT, &
2382 ! & 'REFKDT',REFKDT,'KDT',KDT,'SBETA',SBETA, 'SHDFAC',SHDFAC, &
2383 ! & 'RSMIN',RSMIN,'RGL',RGL,'HS',HS,'ZBOT',ZBOT,'FRZX',FRZX, &
2384 ! & 'PSISAT',PSISAT,'SLOPE',SLOPE,'SNUP',SNUP,'SALP',SALP,'BEXP', &
2385 ! & BEXP, &
2386 ! & 'DKSAT',DKSAT,'DWSAT',DWSAT, &
2387 ! & 'SMCMAX',SMCMAX,'SMCWLT',SMCWLT,'SMCREF',SMCREF,'SMCDRY',SMCDRY, &
2388 ! & 'F1',F1,'QUARTZ',QUARTZ,'FXEXP',FXEXP, &
2389 ! & 'RTDIS',RTDIS,'SLDPTH',SLDPTH,'ZSOIL',ZSOIL, 'NROOT',NROOT, &
2390 ! & 'NSOIL',NSOIL,'Z0',Z0,'CZIL',CZIL,'LAI',LAI, &
2391 ! & 'CSOIL',CSOIL,'PTU',PTU, &
2392 ! & 'LOCAL', LOCAL
2394 END SUBROUTINE REDPRM
2396 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
2398 ! ----------------------------------------------------------------------
2399 ! SUBROUTINE ROSR12
2400 ! ----------------------------------------------------------------------
2401 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
2402 ! ### ### ### ### ### ###
2403 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
2404 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
2405 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
2406 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
2407 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
2408 ! # . . # # . # = # . #
2409 ! # . . # # . # # . #
2410 ! # . . # # . # # . #
2411 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
2412 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
2413 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
2414 ! ### ### ### ### ### ###
2415 ! ----------------------------------------------------------------------
2416 IMPLICIT NONE
2418 INTEGER, INTENT(IN) :: NSOIL
2419 INTEGER :: K, KK
2421 REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
2422 REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
2424 ! ----------------------------------------------------------------------
2425 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
2426 ! ----------------------------------------------------------------------
2427 C (NSOIL) = 0.0
2428 P (1) = - C (1) / B (1)
2429 ! ----------------------------------------------------------------------
2430 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
2431 ! ----------------------------------------------------------------------
2433 ! ----------------------------------------------------------------------
2434 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
2435 ! ----------------------------------------------------------------------
2436 DELTA (1) = D (1) / B (1)
2437 DO K = 2,NSOIL
2438 P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
2439 DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
2440 * P (K -1)))
2441 END DO
2442 ! ----------------------------------------------------------------------
2443 ! SET P TO DELTA FOR LOWEST SOIL LAYER
2444 ! ----------------------------------------------------------------------
2445 P (NSOIL) = DELTA (NSOIL)
2447 ! ----------------------------------------------------------------------
2448 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
2449 ! ----------------------------------------------------------------------
2450 DO K = 2,NSOIL
2451 KK = NSOIL - K + 1
2452 P (KK) = P (KK) * P (KK +1) + DELTA (KK)
2453 END DO
2454 ! ----------------------------------------------------------------------
2455 END SUBROUTINE ROSR12
2456 ! ----------------------------------------------------------------------
2459 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2460 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2461 QUARTZ,CSOIL,VEGTYP,ISURBAN)
2463 ! ----------------------------------------------------------------------
2464 ! SUBROUTINE SHFLX
2465 ! ----------------------------------------------------------------------
2466 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
2467 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
2468 ! ON THE TEMPERATURE.
2469 ! ----------------------------------------------------------------------
2470 IMPLICIT NONE
2472 INTEGER, INTENT(IN) :: ICE, NSOIL, VEGTYP, ISURBAN
2473 INTEGER :: I
2475 REAL, INTENT(IN) :: BEXP,CSOIL,DF1,DT,F1,PSISAT,QUARTZ, &
2476 SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1
2477 REAL, INTENT(INOUT) :: T1
2478 REAL, INTENT(OUT) :: SSOIL
2479 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
2480 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O
2481 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
2482 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
2483 REAL, PARAMETER :: T0 = 273.15
2485 ! ----------------------------------------------------------------------
2486 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
2487 ! ----------------------------------------------------------------------
2489 IF ( ICE == -1 ) THEN
2490 ! Glacial ice case
2492 CALL HRTICE_GLACIAL (RHSTS,STC,TBOT,NSOIL,ZSOIL,YY,ZZ1,DF1,&
2493 AI,BI,CI)
2495 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
2497 ELSEIF ( ICE == 0 ) THEN
2499 ! Non-glacial land case
2501 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
2502 ZBOT,PSISAT,SH2O,DT, &
2503 BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP,ISURBAN)
2505 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
2507 ENDIF
2509 DO I = 1,NSOIL
2510 STC (I) = STCF (I)
2511 ENDDO
2513 ! ----------------------------------------------------------------------
2514 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
2515 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
2516 ! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
2517 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
2518 ! DIFFERENTLY IN ROUTINE SNOPAC)
2519 ! ----------------------------------------------------------------------
2520 ! ----------------------------------------------------------------------
2521 ! CALCULATE SURFACE SOIL HEAT FLUX
2522 ! ----------------------------------------------------------------------
2523 T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
2524 SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
2526 ! ----------------------------------------------------------------------
2527 END SUBROUTINE SHFLX
2528 ! ----------------------------------------------------------------------
2530 SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2531 & SH2O,SLOPE,KDT,FRZFACT, &
2532 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2533 & SHDFAC,CMCMAX, &
2534 & RUNOFF1,RUNOFF2,RUNOFF3, &
2535 & EDIR,EC,ET, &
2536 & DRIP)
2538 ! ----------------------------------------------------------------------
2539 ! SUBROUTINE SMFLX
2540 ! ----------------------------------------------------------------------
2541 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
2542 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
2543 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
2544 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
2545 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
2546 ! ----------------------------------------------------------------------
2547 IMPLICIT NONE
2549 INTEGER, INTENT(IN) :: NSOIL
2550 INTEGER :: I,K
2552 REAL, INTENT(IN) :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR, &
2553 KDT, PRCP1, SHDFAC, SLOPE, SMCMAX, SMCWLT
2554 REAL, INTENT(OUT) :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
2555 REAL, INTENT(INOUT) :: CMC
2556 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET,ZSOIL
2557 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SMC, SH2O
2558 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS, RHSTT, &
2559 SICE, SH2OA, SH2OFG
2560 REAL :: DUMMY, EXCESS,FRZFACT,PCPDRP,RHSCT,TRHSCT
2561 REAL :: FAC2
2562 REAL :: FLIMIT
2564 ! ----------------------------------------------------------------------
2565 ! EXECUTABLE CODE BEGINS HERE.
2566 ! ----------------------------------------------------------------------
2567 ! ----------------------------------------------------------------------
2568 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
2569 ! ----------------------------------------------------------------------
2570 DUMMY = 0.
2572 ! ----------------------------------------------------------------------
2573 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
2574 ! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
2575 ! FALL TO THE GRND.
2576 ! ----------------------------------------------------------------------
2577 RHSCT = SHDFAC * PRCP1- EC
2578 DRIP = 0.
2579 TRHSCT = DT * RHSCT
2580 EXCESS = CMC + TRHSCT
2582 ! ----------------------------------------------------------------------
2583 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
2584 ! SOIL
2585 ! ----------------------------------------------------------------------
2586 IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
2587 PCPDRP = (1. - SHDFAC) * PRCP1+ DRIP / DT
2589 ! ----------------------------------------------------------------------
2590 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT and SSTEP
2591 !
2592 DO I = 1,NSOIL
2593 SICE (I) = SMC (I) - SH2O (I)
2594 END DO
2595 ! ----------------------------------------------------------------------
2596 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
2597 ! TENDENCY EQUATIONS.
2598 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
2599 ! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
2600 ! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
2601 ! THE FIRST SOIL LAYER)
2602 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
2603 ! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
2604 ! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
2605 ! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
2606 ! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
2607 ! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
2608 ! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
2609 ! SOIL MOISTURE STATE
2610 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
2611 ! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
2612 ! OF SECTION 2 OF KALNAY AND KANAMITSU
2613 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
2614 ! ----------------------------------------------------------------------
2615 ! According to Dr. Ken Mitchell's suggestion, add the second contraint
2616 ! to remove numerical instability of runoff and soil moisture
2617 ! FLIMIT is a limit value for FAC2
2618 FAC2=0.0
2619 DO I=1,NSOIL
2620 FAC2=MAX(FAC2,SH2O(I)/SMCMAX)
2621 ENDDO
2622 CALL FAC2MIT(SMCMAX,FLIMIT)
2624 ! ----------------------------------------------------------------------
2625 ! FROZEN GROUND VERSION:
2626 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
2627 ! INC&UDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
2628 ! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
2629 ! ----------------------------------------------------------------------
2630 IF ( ( (PCPDRP * DT) > (0.0001*1000.0* (- ZSOIL (1))* SMCMAX) ) &
2631 .OR. (FAC2 > FLIMIT) ) THEN
2632 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
2633 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
2634 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2635 CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
2636 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2637 DO K = 1,NSOIL
2638 SH2OA (K) = (SH2O (K) + SH2OFG (K)) * 0.5
2639 END DO
2640 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
2641 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
2642 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2643 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
2644 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2646 ELSE
2647 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
2648 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
2649 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2650 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
2651 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2652 ! RUNOF = RUNOFF
2654 END IF
2656 ! ----------------------------------------------------------------------
2657 END SUBROUTINE SMFLX
2658 ! ----------------------------------------------------------------------
2661 SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
2663 ! ----------------------------------------------------------------------
2664 ! SUBROUTINE SNFRAC
2665 ! ----------------------------------------------------------------------
2666 ! CALCULATE SNOW FRACTION (0 -> 1)
2667 ! SNEQV SNOW WATER EQUIVALENT (M)
2668 ! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
2669 ! SALP TUNING PARAMETER
2670 ! SNCOVR FRACTIONAL SNOW COVER
2671 ! ----------------------------------------------------------------------
2672 IMPLICIT NONE
2674 REAL, INTENT(IN) :: SNEQV,SNUP,SALP,SNOWH
2675 REAL, INTENT(OUT) :: SNCOVR
2676 REAL :: RSNOW, Z0N
2678 ! ----------------------------------------------------------------------
2679 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
2680 ! REDPRM) ABOVE WHICH SNOCVR=1.
2681 ! ----------------------------------------------------------------------
2682 IF (SNEQV < SNUP) THEN
2683 RSNOW = SNEQV / SNUP
2684 SNCOVR = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
2685 ELSE
2686 SNCOVR = 1.0
2687 END IF
2689 ! FORMULATION OF DICKINSON ET AL. 1986
2690 ! Z0N = 0.035
2692 ! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
2694 ! FORMULATION OF MARSHALL ET AL. 1994
2695 ! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
2697 ! ----------------------------------------------------------------------
2698 END SUBROUTINE SNFRAC
2699 ! ----------------------------------------------------------------------
2701 SUBROUTINE SNKSRC (TSNSR,TAVG,SMC,SH2O,ZSOIL,NSOIL, &
2702 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2703 ! ----------------------------------------------------------------------
2704 ! SUBROUTINE SNKSRC
2705 ! ----------------------------------------------------------------------
2706 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
2707 ! AVAILABLE LIQUED WATER.
2708 ! ----------------------------------------------------------------------
2709 IMPLICIT NONE
2711 INTEGER, INTENT(IN) :: K,NSOIL
2712 REAL, INTENT(IN) :: BEXP, DT, PSISAT, QTOT, SMC, SMCMAX, &
2713 TAVG
2714 REAL, INTENT(INOUT) :: SH2O
2716 REAL, DIMENSION(1:NSOIL), INTENT(IN):: ZSOIL
2718 REAL :: DF, DZ, DZH, FREE, TSNSR, &
2719 TDN, TM, TUP, TZ, X0, XDN, XH2O, XUP
2721 REAL, PARAMETER :: DH2O = 1.0000E3, HLICE = 3.3350E5, &
2722 T0 = 2.7315E2
2724 IF (K == 1) THEN
2725 DZ = - ZSOIL (1)
2726 ELSE
2727 DZ = ZSOIL (K -1) - ZSOIL (K)
2728 END IF
2729 ! ----------------------------------------------------------------------
2730 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
2731 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
2732 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
2733 ! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
2734 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
2735 ! ----------------------------------------------------------------------
2736 ! FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
2738 ! ----------------------------------------------------------------------
2739 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
2740 ! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
2741 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
2742 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
2743 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
2744 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
2745 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
2746 ! ----------------------------------------------------------------------
2747 CALL FRH2O (FREE,TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
2749 ! ----------------------------------------------------------------------
2750 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
2751 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
2752 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
2753 ! ----------------------------------------------------------------------
2754 XH2O = SH2O + QTOT * DT / (DH2O * HLICE * DZ)
2755 IF ( XH2O < SH2O .AND. XH2O < FREE) THEN
2756 IF ( FREE > SH2O ) THEN
2757 XH2O = SH2O
2758 ELSE
2759 XH2O = FREE
2760 END IF
2761 END IF
2762 ! ----------------------------------------------------------------------
2763 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
2764 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
2765 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
2766 ! ----------------------------------------------------------------------
2767 IF ( XH2O > SH2O .AND. XH2O > FREE ) THEN
2768 IF ( FREE < SH2O ) THEN
2769 XH2O = SH2O
2770 ELSE
2771 XH2O = FREE
2772 END IF
2773 END IF
2775 ! ----------------------------------------------------------------------
2776 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
2777 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
2778 ! ----------------------------------------------------------------------
2779 ! SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
2780 IF (XH2O < 0.) XH2O = 0.
2781 IF (XH2O > SMC) XH2O = SMC
2782 TSNSR = - DH2O * HLICE * DZ * (XH2O - SH2O)/ DT
2783 SH2O = XH2O
2785 ! ----------------------------------------------------------------------
2786 END SUBROUTINE SNKSRC
2787 ! ----------------------------------------------------------------------
2789 SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
2790 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
2791 SBETA,DF1, &
2792 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,&
2793 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
2794 SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT, &
2795 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
2796 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
2797 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
2798 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI,&
2799 RIBB,SOLDN, &
2800 ISURBAN, &
2802 VEGTYP)
2804 ! ----------------------------------------------------------------------
2805 ! SUBROUTINE SNOPAC
2806 ! ----------------------------------------------------------------------
2807 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
2808 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
2809 ! PRESENT.
2810 ! ----------------------------------------------------------------------
2811 IMPLICIT NONE
2813 INTEGER, INTENT(IN) :: ICE, NROOT, NSOIL,VEGTYP
2814 INTEGER, INTENT(IN) :: ISURBAN
2815 INTEGER :: K
2816 !
2817 ! kmh 09/03/2006 add IT16 for surface temperature iteration
2818 !
2819 INTEGER :: IT16
2820 LOGICAL, INTENT(IN) :: SNOWNG
2821 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DF1,DKSAT, &
2822 DT,DWSAT, EPSCA,FDOWN,F1,FXEXP, &
2823 FRZFACT,KDT,PC, PRCP,PSISAT,Q2,QUARTZ, &
2824 RCH,RR,SBETA,SFCPRS, SFCTMP, SHDFAC, &
2825 SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, T24, &
2826 TBOT,TH2,ZBOT,EMISSI,SOLDN
2827 REAL, INTENT(INOUT) :: CMC, BETA, ESD,FLX2,PRCPF,SNOWH,SNCOVR, &
2828 SNDENS, T1, RIBB, ETP
2829 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR, ETNS, ESNOW,ETT, &
2830 FLX1,FLX3, RUNOFF1,RUNOFF2,RUNOFF3, &
2831 SSOIL,SNOMLT
2832 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
2833 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
2834 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
2835 REAL, DIMENSION(1:NSOIL) :: ET1
2836 REAL :: DENOM,DSOIL,DTOT,EC1,EDIR1,ESDFLX,ETA, &
2837 ETT1, ESNOW1, ESNOW2, ETA1,ETP1,ETP2, &
2838 ETP3, ETNS1, ETANRG, ETAX, EX, FLX3X, &
2839 FRCSNO,FRCSOI, PRCP1, QSAT,RSNOW, SEH, &
2840 SNCOND,SSOIL1, T11,T12, T12A, T12AX, &
2841 T12B, T14, YY, ZZ1
2842 ! T12B, T14, YY, ZZ1,EMISSI_S
2843 !
2844 ! kmh 01/11/2007 add T15, T16, and DTOT2 for SFC T iteration and snow heat flux
2845 !
2846 REAL :: T15, T16, DTOT2
2847 REAL, PARAMETER :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6, &
2848 LSUBS = 2.83E+6, TFREEZ = 273.15, &
2849 SNOEXP = 2.0
2851 ! ----------------------------------------------------------------------
2852 ! EXECUTABLE CODE BEGINS HERE:
2853 ! ----------------------------------------------------------------------
2854 ! IF GLACIAL ICE (ICE == -1), SNOWCOVER FRACTION = 1.0,
2855 ! AND SUBLIMATION IS AT THE POTENTIAL RATE.
2856 ! FOR NON-GLACIAL LAND (ICE == 0), IF SNOWCOVER FRACTION < 1.0, TOTAL
2857 ! EVAPORATION < POTENTIAL DUE TO NON-POTENTIAL CONTRIBUTION FROM
2858 ! NON-SNOW COVERED FRACTION.
2859 ! ----------------------------------------------------------------------
2860 ! INITIALIZE EVAP TERMS.
2861 ! ----------------------------------------------------------------------
2862 ! conversions:
2863 ! ESNOW [KG M-2 S-1]
2864 ! ESDFLX [KG M-2 S-1] .le. ESNOW
2865 ! ESNOW1 [M S-1]
2866 ! ESNOW2 [M]
2867 ! ETP [KG M-2 S-1]
2868 ! ETP1 [M S-1]
2869 ! ETP2 [M]
2870 ! ----------------------------------------------------------------------
2871 DEW = 0.
2872 EDIR = 0.
2873 EDIR1 = 0.
2874 EC1 = 0.
2875 EC = 0.
2876 ! EMISSI_S=0.95 ! For snow
2878 DO K = 1,NSOIL
2879 ET (K) = 0.
2880 ET1 (K) = 0.
2881 END DO
2882 ETT = 0.
2883 ETT1 = 0.
2884 ETNS = 0.
2885 ETNS1 = 0.
2886 ESNOW = 0.
2887 ESNOW1 = 0.
2888 ESNOW2 = 0.
2890 ! ----------------------------------------------------------------------
2891 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO ETP1 IN M S-1
2892 ! ----------------------------------------------------------------------
2893 PRCP1 = PRCPF *0.001
2894 ! ----------------------------------------------------------------------
2895 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
2896 ! ----------------------------------------------------------------------
2897 BETA = 1.0
2898 IF (ETP <= 0.0) THEN
2899 IF ( ( RIBB >= 0.1 ) .AND. ( FDOWN > 150.0 ) ) THEN
2900 ETP=(MIN(ETP*(1.0-RIBB),0.)*SNCOVR/0.980 + ETP*(0.980-SNCOVR))/0.980
2901 ENDIF
2902 IF(ETP == 0.) BETA = 0.0
2903 ETP1 = ETP * 0.001
2904 DEW = -ETP1
2905 ESNOW2 = ETP1*DT
2906 ETANRG = ETP*((1.-SNCOVR)*LSUBC + SNCOVR*LSUBS)
2907 ELSE
2908 ETP1 = ETP * 0.001
2909 IF ( ICE == -1 ) THEN
2910 ! GLACIAL ICE CASE
2911 ESNOW = ETP
2912 ESNOW1 = ESNOW*0.001
2913 ESNOW2 = ESNOW1*DT
2914 ETANRG = ESNOW*LSUBS
2915 ELSEIF ( ICE == 0) THEN
2916 ! NON-GLACIAL LAND CASE
2917 IF (SNCOVR < 1.) THEN
2918 CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2919 SH2O, &
2920 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2921 SMCREF,SHDFAC,CMCMAX, &
2922 SMCDRY,CFACTR, &
2923 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS, &
2924 FXEXP)
2925 ! ----------------------------------------------------------------------------
2926 EDIR1 = EDIR1* (1. - SNCOVR)
2927 EC1 = EC1* (1. - SNCOVR)
2928 DO K = 1,NSOIL
2929 ET1 (K) = ET1 (K)* (1. - SNCOVR)
2930 END DO
2931 ETT1 = ETT1*(1.-SNCOVR)
2932 ! ETNS1 = EDIR1+ EC1+ ETT1
2933 ETNS1 = ETNS1*(1.-SNCOVR)
2934 ! ----------------------------------------------------------------------------
2935 EDIR = EDIR1*1000.
2936 EC = EC1*1000.
2937 DO K = 1,NSOIL
2938 ET (K) = ET1 (K)*1000.
2939 END DO
2940 ETT = ETT1*1000.
2941 ETNS = ETNS1*1000.
2942 ! ----------------------------------------------------------------------
2944 ENDIF
2945 ENDIF
2946 ESNOW = ETP*SNCOVR
2947 ESNOW1 = ESNOW*0.001
2948 ESNOW2 = ESNOW1*DT
2949 ETANRG = ESNOW*LSUBS + ETNS*LSUBC
2950 ENDIF
2952 ! ----------------------------------------------------------------------
2953 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
2954 ! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
2955 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
2956 ! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
2957 ! ----------------------------------------------------------------------
2958 FLX1 = 0.0
2959 IF (SNOWNG) THEN
2960 FLX1 = CPICE * PRCP * (T1- SFCTMP)
2961 ELSE
2962 IF (PRCP > 0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
2963 ! ----------------------------------------------------------------------
2964 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
2965 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
2966 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
2967 ! FLUXES. FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
2968 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
2969 ! PENMAN.
2970 ! ----------------------------------------------------------------------
2971 END IF
2972 DSOIL = - (0.5 * ZSOIL (1))
2973 DTOT = SNOWH + DSOIL
2974 DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
2975 ! surface emissivity weighted by snow cover fraction
2976 ! T12A = ( (FDOWN - FLX1 - FLX2 - &
2977 ! & ((SNCOVR*EMISSI_S)+EMISSI*(1.0-SNCOVR))*SIGMA *T24)/RCH &
2978 ! & + TH2 - SFCTMP - ETANRG/RCH ) / RR
2979 T12A = ( (FDOWN - FLX1- FLX2- EMISSI * SIGMA * T24)/ RCH &
2980 + TH2- SFCTMP - ETANRG / RCH ) / RR
2982 T12B = DF1 * STC (1) / (DTOT * RR * RCH)
2984 ! ----------------------------------------------------------------------
2985 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
2986 ! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
2987 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
2988 ! DEPENDING ON SIGN OF ETP.
2989 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
2990 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
2991 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
2992 ! TO ZERO.
2993 ! ----------------------------------------------------------------------
2994 ! SUB-FREEZING BLOCK
2995 ! ----------------------------------------------------------------------
2996 T12 = (SFCTMP + T12A + T12B) / DENOM
2997 IF (T12 <= TFREEZ) THEN
2998 T1 = T12
2999 SSOIL = DF1 * (T1- STC (1)) / DTOT
3000 ! ESD = MAX (0.0, ESD- ETP2)
3001 ESD = MAX(0.0, ESD-ESNOW2)
3002 FLX3 = 0.0
3003 EX = 0.0
3005 SNOMLT = 0.0
3006 ! ----------------------------------------------------------------------
3007 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3008 ! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
3009 ! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3010 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3011 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3012 ! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3013 ! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
3014 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3015 ! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3016 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3017 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3018 ! ----------------------------------------------------------------------
3019 ! ABOVE FREEZING BLOCK
3020 ! ----------------------------------------------------------------------
3021 ELSE
3022 T1 = TFREEZ * SNCOVR ** SNOEXP + T12 * (1.0- SNCOVR ** SNOEXP)
3023 BETA = 1.0
3025 ! ----------------------------------------------------------------------
3026 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3027 ! BETA<1
3028 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3029 ! ----------------------------------------------------------------------
3030 IF ( ICE == -1 ) then
3031 ! kmh 12/15/2005 modify SSOIL
3032 ! kmh 09/03/2006 modify DTOT
3033 IF ( DTOT .GT. 2.0*DSOIL ) THEN
3034 DTOT = 2.0*DSOIL
3035 ENDIF
3036 ENDIF
3037 SSOIL = DF1 * (T1- STC (1)) / DTOT
3038 IF (ESD-ESNOW2 <= ESDMIN) THEN
3039 ESD = 0.0
3040 EX = 0.0
3041 SNOMLT = 0.0
3042 FLX3 = 0.0
3043 ! ----------------------------------------------------------------------
3044 ! SUBLIMATION LESS THAN DEPTH OF SNOWPACK
3045 ! SNOWPACK (ESD) REDUCED BY ESNOW2 (DEPTH OF SUBLIMATED SNOW)
3046 ! ----------------------------------------------------------------------
3047 ELSE
3048 ESD = ESD-ESNOW2
3049 ETP3 = ETP * LSUBC
3050 SEH = RCH * (T1- TH2)
3051 T14 = T1* T1
3052 T14 = T14* T14
3053 ! FLX3 = FDOWN - FLX1 - FLX2 - &
3054 ! ((SNCOVR*EMISSI_S)+EMISSI*(1-SNCOVR))*SIGMA*T14 - &
3055 ! SSOIL - SEH - ETANRG
3056 FLX3 = FDOWN - FLX1- FLX2- EMISSI*SIGMA * T14- SSOIL - SEH - ETANRG
3057 IF (FLX3 <= 0.0) FLX3 = 0.0
3058 ! ----------------------------------------------------------------------
3059 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
3060 ! ----------------------------------------------------------------------
3061 EX = FLX3*0.001/ LSUBF
3063 ! ----------------------------------------------------------------------
3064 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
3065 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
3066 ! ----------------------------------------------------------------------
3067 SNOMLT = EX * DT
3068 IF (ESD- SNOMLT >= ESDMIN) THEN
3069 ESD = ESD- SNOMLT
3070 ! ----------------------------------------------------------------------
3071 ! SNOWMELT EXCEEDS SNOW DEPTH
3072 ! ----------------------------------------------------------------------
3073 ELSE
3074 EX = ESD / DT
3075 FLX3 = EX *1000.0* LSUBF
3076 SNOMLT = ESD
3078 ESD = 0.0
3079 ! ----------------------------------------------------------------------
3080 ! END OF 'ESD .LE. ETP2' IF-BLOCK
3081 ! ----------------------------------------------------------------------
3082 END IF
3083 END IF
3085 ! ----------------------------------------------------------------------
3086 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
3087 ! ----------------------------------------------------------------------
3088 ! ----------------------------------------------------------------------
3089 ! IF NON-GLACIAL LAND, ADD SNOWMELT RATE (EX) TO PRECIP RATE TO BE USED
3090 ! IN SUBROUTINE SMFLX (SOIL MOISTURE EVOLUTION) VIA INFILTRATION.
3091 !
3092 ! FOR GLACIAL ICE, THE SNOWMELT WILL BE ADDED TO SUBSURFACE
3093 ! RUNOFF/BASEFLOW LATER NEAR THE END OF SFLX (AFTER RETURN FROM CALL TO
3094 ! SUBROUTINE SNOPAC)
3095 ! ----------------------------------------------------------------------
3096 IF (ICE == 0) PRCP1 = PRCP1+ EX
3098 ! ----------------------------------------------------------------------
3099 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
3100 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
3101 ! (BELOW).
3102 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES FOR NON-GLACIAL LAND.
3103 ! IF GLACIAL ICE (ICE == -1), SKIP CALL TO SMFLX,
3104 ! SINCE NO SOIL MEDIUM FOR GLACIAL ICE.
3105 ! ----------------------------------------------------------------------
3106 END IF
3107 IF (ICE == 0) THEN
3108 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3109 SH2O,SLOPE,KDT,FRZFACT, &
3110 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3111 SHDFAC,CMCMAX, &
3112 RUNOFF1,RUNOFF2,RUNOFF3, &
3113 EDIR1,EC1,ET1, &
3114 DRIP)
3115 ! ----------------------------------------------------------------------
3116 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
3117 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
3118 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
3119 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
3120 ! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
3121 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
3122 ! ----------------------------------------------------------------------
3123 END IF
3124 ZZ1 = 1.0
3125 YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
3127 ! ----------------------------------------------------------------------
3128 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
3129 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
3130 ! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
3131 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
3132 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
3133 ! ----------------------------------------------------------------------
3134 T11 = T1
3135 CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
3136 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3137 QUARTZ,CSOIL,VEGTYP,ISURBAN)
3139 ! ----------------------------------------------------------------------
3140 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
3141 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
3142 ! ----------------------------------------------------------------------
3143 IF (ICE == 0) THEN
3144 ! NON-GLACIAL LAND
3145 IF (ESD > 0.) THEN
3146 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3147 ELSE
3148 ESD = 0.
3149 SNOWH = 0.
3150 SNDENS = 0.
3151 SNCOND = 1.
3152 SNCOVR = 0.
3153 END IF
3154 ELSEIF (ICE == -1) THEN
3155 ! GLACIAL ICE
3156 IF (ESD .GE. 0.10) THEN
3157 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3158 ELSE
3159 ESD = 0.10
3160 SNOWH = 0.50
3161 !KWM???? SNDENS =
3162 !KWM???? SNCOND =
3163 SNCOVR = 1.0
3164 ENDIF
3165 ENDIF
3166 ! ----------------------------------------------------------------------
3167 END SUBROUTINE SNOPAC
3168 ! ----------------------------------------------------------------------
3171 SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
3173 ! ----------------------------------------------------------------------
3174 ! SUBROUTINE SNOWPACK
3175 ! ----------------------------------------------------------------------
3176 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
3177 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
3178 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
3179 ! KOREN, 03/25/95.
3180 ! ----------------------------------------------------------------------
3181 ! ESD WATER EQUIVALENT OF SNOW (M)
3182 ! DTSEC TIME STEP (SEC)
3183 ! SNOWH SNOW DEPTH (M)
3184 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
3185 ! TSNOW SNOW SURFACE TEMPERATURE (K)
3186 ! TSOIL SOIL SURFACE TEMPERATURE (K)
3188 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
3189 ! ----------------------------------------------------------------------
3190 IMPLICIT NONE
3192 INTEGER :: IPOL, J
3193 REAL, INTENT(IN) :: ESD, DTSEC,TSNOW,TSOIL
3194 REAL, INTENT(INOUT) :: SNOWH, SNDENS
3195 REAL :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP, &
3196 TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
3197 REAL, PARAMETER :: C1 = 0.01, C2 = 21.0, G = 9.81, &
3198 KN = 4000.0
3199 ! ----------------------------------------------------------------------
3200 ! CONVERSION INTO SIMULATION UNITS
3201 ! ----------------------------------------------------------------------
3202 SNOWHC = SNOWH *100.
3203 ESDC = ESD *100.
3204 DTHR = DTSEC /3600.
3205 TSNOWC = TSNOW -273.15
3206 TSOILC = TSOIL -273.15
3208 ! ----------------------------------------------------------------------
3209 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
3210 ! ----------------------------------------------------------------------
3211 ! ----------------------------------------------------------------------
3212 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
3213 ! SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
3214 ! BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
3215 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
3216 ! NUMERICALLY BELOW:
3217 ! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
3218 ! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
3219 ! ----------------------------------------------------------------------
3220 TAVGC = 0.5* (TSNOWC + TSOILC)
3221 IF (ESDC > 1.E-2) THEN
3222 ESDCX = ESDC
3223 ELSE
3224 ESDCX = 1.E-2
3225 END IF
3227 ! DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
3228 ! ----------------------------------------------------------------------
3229 ! THE FUNCTION OF THE FORM (e**x-1)/x EMBEDDED IN ABOVE EXPRESSION
3230 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
3231 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
3232 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
3233 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
3234 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
3235 ! POLYNOMIAL EXPANSION.
3237 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
3238 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
3239 ! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
3240 ! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
3241 ! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
3242 ! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
3243 ! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
3244 ! ----------------------------------------------------------------------
3245 BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
3246 IPOL = 4
3247 PEXP = 0.
3248 ! PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
3249 DO J = IPOL,1, -1
3250 PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
3251 END DO
3253 PEXP = PEXP + 1.
3254 ! ----------------------------------------------------------------------
3255 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
3256 ! ----------------------------------------------------------------------
3257 ! END OF KOREAN FORMULATION
3259 ! BASE FORMULATION (COGLEY ET AL., 1990)
3260 ! CONVERT DENSITY FROM G/CM3 TO KG/M3
3261 ! DSM=SNDENS*1000.0
3263 ! DSX=DSM+DTSEC*0.5*DSM*G*ESD/
3264 ! & (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
3266 ! & CONVERT DENSITY FROM KG/M3 TO G/CM3
3267 ! DSX=DSX/1000.0
3269 ! END OF COGLEY ET AL. FORMULATION
3271 ! ----------------------------------------------------------------------
3272 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
3273 ! ----------------------------------------------------------------------
3274 DSX = SNDENS * (PEXP)
3275 IF (DSX > 0.40) DSX = 0.40
3276 IF (DSX < 0.05) DSX = 0.05
3277 ! ----------------------------------------------------------------------
3278 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
3279 ! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
3280 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
3281 ! ----------------------------------------------------------------------
3282 SNDENS = DSX
3283 IF (TSNOWC >= 0.) THEN
3284 DW = 0.13* DTHR /24.
3285 SNDENS = SNDENS * (1. - DW) + DW
3286 IF (SNDENS >= 0.40) SNDENS = 0.40
3287 ! ----------------------------------------------------------------------
3288 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
3289 ! CHANGE SNOW DEPTH UNITS TO METERS
3290 ! ----------------------------------------------------------------------
3291 END IF
3292 SNOWHC = ESDC / SNDENS
3293 SNOWH = SNOWHC *0.01
3295 ! ----------------------------------------------------------------------
3296 END SUBROUTINE SNOWPACK
3297 ! ----------------------------------------------------------------------
3299 SUBROUTINE SNOWZ0 (SNCOVR,Z0, Z0BRD, SNOWH)
3301 ! ----------------------------------------------------------------------
3302 ! SUBROUTINE SNOWZ0
3303 ! ----------------------------------------------------------------------
3304 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
3305 ! SNCOVR FRACTIONAL SNOW COVER
3306 ! Z0 ROUGHNESS LENGTH (m)
3307 ! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
3308 ! ----------------------------------------------------------------------
3309 IMPLICIT NONE
3310 REAL, INTENT(IN) :: SNCOVR, Z0BRD
3311 REAL, INTENT(OUT) :: Z0
3312 REAL, PARAMETER :: Z0S=0.001
3313 REAL, INTENT(IN) :: SNOWH
3314 REAL :: BURIAL
3315 REAL :: Z0EFF
3317 !m Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0S
3318 BURIAL = 7.0*Z0BRD - SNOWH
3319 IF(BURIAL.LE.0.0007) THEN
3320 Z0EFF = Z0S
3321 ELSE
3322 Z0EFF = BURIAL/7.0
3323 ENDIF
3325 Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0EFF
3327 ! ----------------------------------------------------------------------
3328 END SUBROUTINE SNOWZ0
3329 ! ----------------------------------------------------------------------
3332 SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
3334 ! ----------------------------------------------------------------------
3335 ! SUBROUTINE SNOW_NEW
3336 ! ----------------------------------------------------------------------
3337 ! CALCULATE SNOW DEPTH AND DENSITY TO ACCOUNT FOR THE NEW SNOWFALL.
3338 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
3340 ! TEMP AIR TEMPERATURE (K)
3341 ! NEWSN NEW SNOWFALL (M)
3342 ! SNOWH SNOW DEPTH (M)
3343 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
3344 ! ----------------------------------------------------------------------
3345 IMPLICIT NONE
3346 REAL, INTENT(IN) :: NEWSN, TEMP
3347 REAL, INTENT(INOUT) :: SNDENS, SNOWH
3348 REAL :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
3350 ! ----------------------------------------------------------------------
3351 ! CONVERSION INTO SIMULATION UNITS
3352 ! ----------------------------------------------------------------------
3353 SNOWHC = SNOWH *100.
3354 NEWSNC = NEWSN *100.
3356 ! ----------------------------------------------------------------------
3357 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
3358 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
3359 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
3360 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
3361 !-----------------------------------------------------------------------
3362 TEMPC = TEMP -273.15
3363 IF (TEMPC <= -15.) THEN
3364 DSNEW = 0.05
3365 ELSE
3366 DSNEW = 0.05+0.0017* (TEMPC +15.)**1.5
3367 END IF
3368 ! ----------------------------------------------------------------------
3369 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
3370 ! ----------------------------------------------------------------------
3371 HNEWC = NEWSNC / DSNEW
3372 IF (SNOWHC + HNEWC .LT. 1.0E-3) THEN
3373 SNDENS = MAX(DSNEW,SNDENS)
3374 ELSE
3375 SNDENS = (SNOWHC * SNDENS + HNEWC * DSNEW)/ (SNOWHC + HNEWC)
3376 ENDIF
3377 SNOWHC = SNOWHC + HNEWC
3378 SNOWH = SNOWHC *0.01
3380 ! ----------------------------------------------------------------------
3381 END SUBROUTINE SNOW_NEW
3382 ! ----------------------------------------------------------------------
3384 SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
3385 ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3386 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
3388 ! ----------------------------------------------------------------------
3389 ! SUBROUTINE SRT
3390 ! ----------------------------------------------------------------------
3391 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3392 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3393 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3394 ! ----------------------------------------------------------------------
3395 IMPLICIT NONE
3396 INTEGER, INTENT(IN) :: NSOIL
3397 INTEGER :: IALP1, IOHINF, J, JJ, K, KS
3398 REAL, INTENT(IN) :: BEXP, DKSAT, DT, DWSAT, EDIR, FRZX, &
3399 KDT, PCPDRP, SLOPE, SMCMAX, SMCWLT
3400 REAL, INTENT(OUT) :: RUNOFF1, RUNOFF2
3401 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET, SH2O, SH2OA, SICE, &
3402 ZSOIL
3403 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTT
3404 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI, CI
3405 REAL, DIMENSION(1:NSOIL) :: DMAX
3406 REAL :: ACRT, DD, DDT, DDZ, DDZ2, DENOM, &
3407 DENOM2,DICE, DSMDZ, DSMDZ2, DT1, &
3408 FCR,INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
3409 PX, SICEMAX,SLOPX, SMCAV, SSTT, &
3410 SUM, VAL, WCND, WCND2, WDF, WDF2
3411 INTEGER, PARAMETER :: CVFRZ = 3
3413 ! ----------------------------------------------------------------------
3414 ! FROZEN GROUND VERSION:
3415 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
3416 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
3417 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
3418 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
3419 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
3420 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
3421 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
3422 ! ----------------------------------------------------------------------
3424 ! ----------------------------------------------------------------------
3425 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
3426 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
3427 ! MODIFIED BY Q DUAN
3428 ! ----------------------------------------------------------------------
3429 ! ----------------------------------------------------------------------
3430 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
3431 ! LAYERS.
3432 ! ----------------------------------------------------------------------
3433 IOHINF = 1
3434 SICEMAX = 0.0
3435 DO KS = 1,NSOIL
3436 IF (SICE (KS) > SICEMAX) SICEMAX = SICE (KS)
3437 ! ----------------------------------------------------------------------
3438 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
3439 ! ----------------------------------------------------------------------
3440 END DO
3441 PDDUM = PCPDRP
3442 RUNOFF1 = 0.0
3444 ! ----------------------------------------------------------------------
3445 ! MODIFIED BY Q. DUAN, 5/16/94
3446 ! ----------------------------------------------------------------------
3447 ! IF (IOHINF == 1) THEN
3449 IF (PCPDRP /= 0.0) THEN
3450 DT1 = DT /86400.
3451 SMCAV = SMCMAX - SMCWLT
3453 ! ----------------------------------------------------------------------
3454 ! FROZEN GROUND VERSION:
3455 ! ----------------------------------------------------------------------
3456 DMAX (1)= - ZSOIL (1)* SMCAV
3458 DICE = - ZSOIL (1) * SICE (1)
3459 DMAX (1)= DMAX (1)* (1.0- (SH2OA (1) + SICE (1) - SMCWLT)/ &
3460 SMCAV)
3462 DD = DMAX (1)
3464 ! ----------------------------------------------------------------------
3465 ! FROZEN GROUND VERSION:
3466 ! ----------------------------------------------------------------------
3467 DO KS = 2,NSOIL
3469 DICE = DICE+ ( ZSOIL (KS -1) - ZSOIL (KS) ) * SICE (KS)
3470 DMAX (KS) = (ZSOIL (KS -1) - ZSOIL (KS))* SMCAV
3471 DMAX (KS) = DMAX (KS)* (1.0- (SH2OA (KS) + SICE (KS) &
3472 - SMCWLT)/ SMCAV)
3473 DD = DD+ DMAX (KS)
3474 ! ----------------------------------------------------------------------
3475 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
3476 ! IN BELOW, REMOVE THE SQRT IN ABOVE
3477 ! ----------------------------------------------------------------------
3478 END DO
3479 VAL = (1. - EXP ( - KDT * DT1))
3480 DDT = DD * VAL
3481 PX = PCPDRP * DT
3482 IF (PX < 0.0) PX = 0.0
3484 ! ----------------------------------------------------------------------
3485 ! FROZEN GROUND VERSION:
3486 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
3487 ! ----------------------------------------------------------------------
3488 INFMAX = (PX * (DDT / (PX + DDT)))/ DT
3489 FCR = 1.
3490 IF (DICE > 1.E-2) THEN
3491 ACRT = CVFRZ * FRZX / DICE
3492 SUM = 1.
3493 IALP1 = CVFRZ - 1
3494 DO J = 1,IALP1
3495 K = 1
3496 DO JJ = J +1,IALP1
3497 K = K * JJ
3498 END DO
3499 SUM = SUM + (ACRT ** ( CVFRZ - J)) / FLOAT (K)
3500 END DO
3501 FCR = 1. - EXP ( - ACRT) * SUM
3502 END IF
3504 ! ----------------------------------------------------------------------
3505 ! CORRECTION OF INFILTRATION LIMITATION:
3506 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
3507 ! HYDROLIC CONDUCTIVITY
3508 ! ----------------------------------------------------------------------
3509 ! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
3510 INFMAX = INFMAX * FCR
3512 MXSMC = SH2OA (1)
3513 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
3514 SICEMAX)
3515 INFMAX = MAX (INFMAX,WCND)
3517 INFMAX = MIN (INFMAX,PX/DT)
3518 IF (PCPDRP > INFMAX) THEN
3519 RUNOFF1 = PCPDRP - INFMAX
3520 PDDUM = INFMAX
3521 END IF
3522 ! ----------------------------------------------------------------------
3523 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
3524 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
3525 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
3526 ! ----------------------------------------------------------------------
3527 END IF
3529 MXSMC = SH2OA (1)
3530 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
3531 SICEMAX)
3532 ! ----------------------------------------------------------------------
3533 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
3534 ! ----------------------------------------------------------------------
3535 DDZ = 1. / ( - .5 * ZSOIL (2) )
3536 AI (1) = 0.0
3537 BI (1) = WDF * DDZ / ( - ZSOIL (1) )
3539 ! ----------------------------------------------------------------------
3540 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
3541 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
3542 ! ----------------------------------------------------------------------
3543 CI (1) = - BI (1)
3544 DSMDZ = ( SH2O (1) - SH2O (2) ) / ( - .5 * ZSOIL (2) )
3545 RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET (1))/ ZSOIL (1)
3547 ! ----------------------------------------------------------------------
3548 ! INITIALIZE DDZ2
3549 ! ----------------------------------------------------------------------
3550 SSTT = WDF * DSMDZ + WCND+ EDIR + ET (1)
3552 ! ----------------------------------------------------------------------
3553 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
3554 ! ----------------------------------------------------------------------
3555 DDZ2 = 0.0
3556 DO K = 2,NSOIL
3557 DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
3558 IF (K /= NSOIL) THEN
3560 ! ----------------------------------------------------------------------
3561 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
3562 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
3563 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
3564 ! ----------------------------------------------------------------------
3565 SLOPX = 1.
3567 MXSMC2 = SH2OA (K)
3568 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
3569 SICEMAX)
3570 ! -----------------------------------------------------------------------
3571 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
3572 ! ----------------------------------------------------------------------
3573 DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
3575 ! ----------------------------------------------------------------------
3576 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
3577 ! ----------------------------------------------------------------------
3578 DSMDZ2 = (SH2O (K) - SH2O (K +1)) / (DENOM * 0.5)
3579 DDZ2 = 2.0 / DENOM
3580 CI (K) = - WDF2 * DDZ2 / DENOM2
3582 ELSE
3583 ! ----------------------------------------------------------------------
3584 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
3585 ! ----------------------------------------------------------------------
3587 ! ----------------------------------------------------------------------
3588 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
3589 ! THIS LAYER
3590 ! ----------------------------------------------------------------------
3591 SLOPX = SLOPE
3592 CALL WDFCND (WDF2,WCND2,SH2OA (NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
3593 SICEMAX)
3595 ! ----------------------------------------------------------------------
3596 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
3597 ! ----------------------------------------------------------------------
3599 ! ----------------------------------------------------------------------
3600 ! SET MATRIX COEF CI TO ZERO
3601 ! ----------------------------------------------------------------------
3602 DSMDZ2 = 0.0
3603 CI (K) = 0.0
3604 ! ----------------------------------------------------------------------
3605 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
3606 ! ----------------------------------------------------------------------
3607 END IF
3608 NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2- (WDF * DSMDZ) &
3609 - WCND+ ET (K)
3611 ! ----------------------------------------------------------------------
3612 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
3613 ! ----------------------------------------------------------------------
3614 RHSTT (K) = NUMER / ( - DENOM2)
3615 AI (K) = - WDF * DDZ / DENOM2
3617 ! ----------------------------------------------------------------------
3618 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
3619 ! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
3620 ! ----------------------------------------------------------------------
3621 BI (K) = - ( AI (K) + CI (K) )
3622 IF (K .eq. NSOIL) THEN
3623 RUNOFF2 = SLOPX * WCND2
3624 END IF
3625 IF (K .ne. NSOIL) THEN
3626 WDF = WDF2
3627 WCND = WCND2
3628 DSMDZ = DSMDZ2
3629 DDZ = DDZ2
3630 END IF
3631 END DO
3632 ! ----------------------------------------------------------------------
3633 END SUBROUTINE SRT
3634 ! ----------------------------------------------------------------------
3636 SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
3637 NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
3638 AI,BI,CI)
3640 ! ----------------------------------------------------------------------
3641 ! SUBROUTINE SSTEP
3642 ! ----------------------------------------------------------------------
3643 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
3644 ! CONTENT VALUES.
3645 ! ----------------------------------------------------------------------
3646 IMPLICIT NONE
3647 INTEGER, INTENT(IN) :: NSOIL
3648 INTEGER :: I, K, KK11
3650 REAL, INTENT(IN) :: CMCMAX, DT, SMCMAX
3651 REAL, INTENT(OUT) :: RUNOFF3
3652 REAL, INTENT(INOUT) :: CMC
3653 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2OIN, SICE, ZSOIL
3654 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SH2OOUT
3655 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTT, SMC
3656 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI, BI, CI
3657 REAL, DIMENSION(1:NSOIL) :: RHSTTin
3658 REAL, DIMENSION(1:NSOIL) :: CIin
3659 REAL :: DDZ, RHSCT, STOT, WPLUS
3661 ! ----------------------------------------------------------------------
3662 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
3663 ! TRI-DIAGONAL MATRIX ROUTINE.
3664 ! ----------------------------------------------------------------------
3665 DO K = 1,NSOIL
3666 RHSTT (K) = RHSTT (K) * DT
3667 AI (K) = AI (K) * DT
3668 BI (K) = 1. + BI (K) * DT
3669 CI (K) = CI (K) * DT
3670 END DO
3671 ! ----------------------------------------------------------------------
3672 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
3673 ! ----------------------------------------------------------------------
3674 DO K = 1,NSOIL
3675 RHSTTin (K) = RHSTT (K)
3676 END DO
3677 DO K = 1,NSOIL
3678 CIin (K) = CI (K)
3679 END DO
3680 ! ----------------------------------------------------------------------
3681 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
3682 ! ----------------------------------------------------------------------
3683 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
3684 ! ----------------------------------------------------------------------
3685 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
3686 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
3687 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
3688 ! ----------------------------------------------------------------------
3689 WPLUS = 0.0
3690 RUNOFF3 = 0.
3692 DDZ = - ZSOIL (1)
3693 DO K = 1,NSOIL
3694 IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
3695 SH2OOUT (K) = SH2OIN (K) + CI (K) + WPLUS / DDZ
3696 STOT = SH2OOUT (K) + SICE (K)
3697 IF (STOT > SMCMAX) THEN
3698 IF (K .eq. 1) THEN
3699 DDZ = - ZSOIL (1)
3700 ELSE
3701 KK11 = K - 1
3702 DDZ = - ZSOIL (K) + ZSOIL (KK11)
3703 END IF
3704 WPLUS = (STOT - SMCMAX) * DDZ
3705 ELSE
3706 WPLUS = 0.
3707 END IF
3708 SMC (K) = MAX ( MIN (STOT,SMCMAX),0.02 )
3709 SH2OOUT (K) = MAX ( (SMC (K) - SICE (K)),0.0)
3710 END DO
3712 ! ----------------------------------------------------------------------
3713 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
3714 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
3715 ! ----------------------------------------------------------------------
3716 RUNOFF3 = WPLUS
3717 CMC = CMC + DT * RHSCT
3718 IF (CMC < 1.E-20) CMC = 0.0
3719 CMC = MIN (CMC,CMCMAX)
3721 ! ----------------------------------------------------------------------
3722 END SUBROUTINE SSTEP
3723 ! ----------------------------------------------------------------------
3725 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
3727 ! ----------------------------------------------------------------------
3728 ! SUBROUTINE TBND
3729 ! ----------------------------------------------------------------------
3730 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
3731 ! THE MIDDLE LAYER TEMPERATURES
3732 ! ----------------------------------------------------------------------
3733 IMPLICIT NONE
3734 INTEGER, INTENT(IN) :: NSOIL
3735 INTEGER :: K
3736 REAL, INTENT(IN) :: TB, TU, ZBOT
3737 REAL, INTENT(OUT) :: TBND1
3738 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
3739 REAL :: ZB, ZUP
3740 REAL, PARAMETER :: T0 = 273.15
3742 ! ----------------------------------------------------------------------
3743 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
3744 ! ----------------------------------------------------------------------
3745 IF (K == 1) THEN
3746 ZUP = 0.
3747 ELSE
3748 ZUP = ZSOIL (K -1)
3749 END IF
3750 ! ----------------------------------------------------------------------
3751 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
3752 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
3753 ! ----------------------------------------------------------------------
3754 IF (K == NSOIL) THEN
3755 ZB = 2.* ZBOT - ZSOIL (K)
3756 ELSE
3757 ZB = ZSOIL (K +1)
3758 END IF
3759 ! ----------------------------------------------------------------------
3760 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
3761 ! ----------------------------------------------------------------------
3763 TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
3764 ! ----------------------------------------------------------------------
3765 END SUBROUTINE TBND
3766 ! ----------------------------------------------------------------------
3769 SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
3771 ! ----------------------------------------------------------------------
3772 ! SUBROUTINE TDFCND
3773 ! ----------------------------------------------------------------------
3774 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
3775 ! POINT AND TIME.
3776 ! ----------------------------------------------------------------------
3777 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
3778 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
3779 ! ----------------------------------------------------------------------
3780 IMPLICIT NONE
3781 REAL, INTENT(IN) :: QZ, SMC, SMCMAX, SH2O
3782 REAL, INTENT(OUT) :: DF
3783 REAL :: AKE, GAMMD, THKDRY, THKICE, THKO, &
3784 THKQTZ,THKSAT,THKS,THKW,SATRATIO,XU, &
3785 XUNFROZ
3787 ! ----------------------------------------------------------------------
3788 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
3789 ! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
3790 ! & 0.35, 0.60, 0.40, 0.82/
3791 ! ----------------------------------------------------------------------
3792 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
3793 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
3794 ! ----------------------------------------------------------------------
3795 ! THKW ......WATER THERMAL CONDUCTIVITY
3796 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
3797 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
3798 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
3799 ! THKICE ....ICE THERMAL CONDUCTIVITY
3800 ! SMCMAX ....POROSITY (= SMCMAX)
3801 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
3802 ! ----------------------------------------------------------------------
3803 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
3805 ! PABLO GRUNMANN, 08/17/98
3806 ! REFS.:
3807 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
3808 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
3809 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
3810 ! UNIVERSITY OF TRONDHEIM,
3811 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
3812 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
3813 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
3814 ! VOL. 55, PP. 1209-1224.
3815 ! ----------------------------------------------------------------------
3816 ! NEEDS PARAMETERS
3817 ! POROSITY(SOIL TYPE):
3818 ! POROS = SMCMAX
3819 ! SATURATION RATIO:
3820 ! PARAMETERS W/(M.K)
3821 SATRATIO = SMC / SMCMAX
3822 ! ICE CONDUCTIVITY:
3823 THKICE = 2.2
3824 ! WATER CONDUCTIVITY:
3825 THKW = 0.57
3826 ! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
3827 ! IF (QZ .LE. 0.2) THKO = 3.0
3828 THKO = 2.0
3829 ! QUARTZ' CONDUCTIVITY
3830 THKQTZ = 7.7
3831 ! SOLIDS' CONDUCTIVITY
3832 THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
3834 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
3835 XUNFROZ = SH2O / SMC
3836 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
3837 XU = XUNFROZ * SMCMAX
3839 ! SATURATED THERMAL CONDUCTIVITY
3840 THKSAT = THKS ** (1. - SMCMAX)* THKICE ** (SMCMAX - XU)* THKW ** &
3841 (XU)
3843 ! DRY DENSITY IN KG/M3
3844 GAMMD = (1. - SMCMAX)*2700.
3846 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
3847 THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
3848 ! FROZEN
3849 IF ( (SH2O + 0.0005) < SMC ) THEN
3850 AKE = SATRATIO
3851 ! UNFROZEN
3852 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
3853 ELSE
3855 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
3856 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
3857 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
3859 IF ( SATRATIO > 0.1 ) THEN
3861 AKE = LOG10 (SATRATIO) + 1.0
3863 ! USE K = KDRY
3864 ELSE
3866 AKE = 0.0
3867 END IF
3868 ! THERMAL CONDUCTIVITY
3870 END IF
3872 DF = AKE * (THKSAT - THKDRY) + THKDRY
3873 ! ----------------------------------------------------------------------
3874 END SUBROUTINE TDFCND
3875 ! ----------------------------------------------------------------------
3877 SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
3879 ! ----------------------------------------------------------------------
3880 ! SUBROUTINE TMPAVG
3881 ! ----------------------------------------------------------------------
3882 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
3883 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
3884 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
3885 ! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
3886 ! ----------------------------------------------------------------------
3887 IMPLICIT NONE
3888 INTEGER K
3890 INTEGER NSOIL
3891 REAL DZ
3892 REAL DZH
3893 REAL T0
3894 REAL TAVG
3895 REAL TDN
3896 REAL TM
3897 REAL TUP
3898 REAL X0
3899 REAL XDN
3900 REAL XUP
3902 REAL ZSOIL (NSOIL)
3904 ! ----------------------------------------------------------------------
3905 PARAMETER (T0 = 2.7315E2)
3906 IF (K .eq. 1) THEN
3907 DZ = - ZSOIL (1)
3908 ELSE
3909 DZ = ZSOIL (K -1) - ZSOIL (K)
3910 END IF
3912 DZH = DZ *0.5
3913 IF (TUP .lt. T0) THEN
3914 IF (TM .lt. T0) THEN
3915 ! ----------------------------------------------------------------------
3916 ! TUP, TM, TDN < T0
3917 ! ----------------------------------------------------------------------
3918 IF (TDN .lt. T0) THEN
3919 TAVG = (TUP + 2.0* TM + TDN)/ 4.0
3920 ! ----------------------------------------------------------------------
3921 ! TUP & TM < T0, TDN .ge. T0
3922 ! ----------------------------------------------------------------------
3923 ELSE
3924 X0 = (T0- TM) * DZH / (TDN - TM)
3925 TAVG = 0.5 * (TUP * DZH + TM * (DZH + X0) + T0* ( &
3926 & 2.* DZH - X0)) / DZ
3927 END IF
3928 ELSE
3929 ! ----------------------------------------------------------------------
3930 ! TUP < T0, TM .ge. T0, TDN < T0
3931 ! ----------------------------------------------------------------------
3932 IF (TDN .lt. T0) THEN
3933 XUP = (T0- TUP) * DZH / (TM - TUP)
3934 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
3935 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP - XDN) &
3936 & + TDN * XDN) / DZ
3937 ! ----------------------------------------------------------------------
3938 ! TUP < T0, TM .ge. T0, TDN .ge. T0
3939 ! ----------------------------------------------------------------------
3940 ELSE
3941 XUP = (T0- TUP) * DZH / (TM - TUP)
3942 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP)) / DZ
3943 END IF
3944 END IF
3945 ELSE
3946 IF (TM .lt. T0) THEN
3947 ! ----------------------------------------------------------------------
3948 ! TUP .ge. T0, TM < T0, TDN < T0
3949 ! ----------------------------------------------------------------------
3950 IF (TDN .lt. T0) THEN
3951 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
3952 TAVG = 0.5 * (T0* (DZ - XUP) + TM * (DZH + XUP) &
3953 & + TDN * DZH) / DZ
3954 ! ----------------------------------------------------------------------
3955 ! TUP .ge. T0, TM < T0, TDN .ge. T0
3956 ! ----------------------------------------------------------------------
3957 ELSE
3958 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
3959 XDN = (T0- TM) * DZH / (TDN - TM)
3960 TAVG = 0.5 * (T0* (2.* DZ - XUP - XDN) + TM * &
3961 & (XUP + XDN)) / DZ
3962 END IF
3963 ELSE
3964 ! ----------------------------------------------------------------------
3965 ! TUP .ge. T0, TM .ge. T0, TDN < T0
3966 ! ----------------------------------------------------------------------
3967 IF (TDN .lt. T0) THEN
3968 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
3969 TAVG = (T0* (DZ - XDN) +0.5* (T0+ TDN)* XDN) / DZ
3970 ! ----------------------------------------------------------------------
3971 ! TUP .ge. T0, TM .ge. T0, TDN .ge. T0
3972 ! ----------------------------------------------------------------------
3973 ELSE
3974 TAVG = (TUP + 2.0* TM + TDN) / 4.0
3975 END IF
3976 END IF
3977 END IF
3978 ! ----------------------------------------------------------------------
3979 END SUBROUTINE TMPAVG
3980 ! ----------------------------------------------------------------------
3982 SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
3983 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT, &
3984 & RTDIS)
3986 ! ----------------------------------------------------------------------
3987 ! SUBROUTINE TRANSP
3988 ! ----------------------------------------------------------------------
3989 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
3990 ! ----------------------------------------------------------------------
3991 IMPLICIT NONE
3992 INTEGER I
3993 INTEGER K
3994 INTEGER NSOIL
3996 INTEGER NROOT
3997 REAL CFACTR
3998 REAL CMC
3999 REAL CMCMAX
4000 REAL DENOM
4001 REAL ET (NSOIL)
4002 REAL ETP1
4003 REAL ETP1A
4004 !.....REAL PART(NSOIL)
4005 REAL GX (NROOT)
4006 REAL PC
4007 REAL Q2
4008 REAL RTDIS (NSOIL)
4009 REAL RTX
4010 REAL SFCTMP
4011 REAL SGX
4012 REAL SHDFAC
4013 REAL SMC (NSOIL)
4014 REAL SMCREF
4015 REAL SMCWLT
4017 ! ----------------------------------------------------------------------
4018 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
4019 ! ----------------------------------------------------------------------
4020 REAL ZSOIL (NSOIL)
4021 DO K = 1,NSOIL
4022 ET (K) = 0.
4023 ! ----------------------------------------------------------------------
4024 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
4025 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
4026 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
4027 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
4028 ! TOTAL ETP1A.
4029 ! ----------------------------------------------------------------------
4030 END DO
4031 IF (CMC .ne. 0.0) THEN
4032 ETP1A = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR)
4033 ELSE
4034 ETP1A = SHDFAC * PC * ETP1
4035 END IF
4036 SGX = 0.0
4037 DO I = 1,NROOT
4038 GX (I) = ( SMC (I) - SMCWLT ) / ( SMCREF - SMCWLT )
4039 GX (I) = MAX ( MIN ( GX (I), 1. ), 0. )
4040 SGX = SGX + GX (I)
4041 END DO
4043 SGX = SGX / NROOT
4044 DENOM = 0.
4045 DO I = 1,NROOT
4046 RTX = RTDIS (I) + GX (I) - SGX
4047 GX (I) = GX (I) * MAX ( RTX, 0. )
4048 DENOM = DENOM + GX (I)
4049 END DO
4051 IF (DENOM .le. 0.0) DENOM = 1.
4052 DO I = 1,NROOT
4053 ET (I) = ETP1A * GX (I) / DENOM
4054 ! ----------------------------------------------------------------------
4055 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
4056 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
4057 ! ----------------------------------------------------------------------
4058 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
4059 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
4060 ! ----------------------------------------------------------------------
4061 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
4062 ! ----------------------------------------------------------------------
4063 ! ET(1) = RTDIS(1) * ETP1A
4064 ! ET(1) = ETP1A * PART(1)
4065 ! ----------------------------------------------------------------------
4066 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
4067 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
4068 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
4069 ! ----------------------------------------------------------------------
4070 ! DO K = 2,NROOT
4071 ! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
4072 ! GX = MAX ( MIN ( GX, 1. ), 0. )
4073 ! TEST CANOPY RESISTANCE
4074 ! GX = 1.0
4075 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
4076 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
4077 ! ----------------------------------------------------------------------
4078 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
4079 ! ----------------------------------------------------------------------
4080 ! ET(K) = RTDIS(K) * ETP1A
4081 ! ET(K) = ETP1A*PART(K)
4082 ! END DO
4083 END DO
4084 ! ----------------------------------------------------------------------
4085 END SUBROUTINE TRANSP
4086 ! ----------------------------------------------------------------------
4088 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4089 & SICEMAX)
4091 ! ----------------------------------------------------------------------
4092 ! SUBROUTINE WDFCND
4093 ! ----------------------------------------------------------------------
4094 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
4095 ! ----------------------------------------------------------------------
4096 IMPLICIT NONE
4097 REAL BEXP
4098 REAL DKSAT
4099 REAL DWSAT
4100 REAL EXPON
4101 REAL FACTR1
4102 REAL FACTR2
4103 REAL SICEMAX
4104 REAL SMC
4105 REAL SMCMAX
4106 REAL VKwgt
4107 REAL WCND
4109 ! ----------------------------------------------------------------------
4110 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
4111 ! ----------------------------------------------------------------------
4112 REAL WDF
4113 FACTR1 = 0.05 / SMCMAX
4115 ! ----------------------------------------------------------------------
4116 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
4117 ! ----------------------------------------------------------------------
4118 FACTR2 = SMC / SMCMAX
4119 FACTR1 = MIN(FACTR1,FACTR2)
4120 EXPON = BEXP + 2.0
4122 ! ----------------------------------------------------------------------
4123 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
4124 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
4125 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
4126 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
4127 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
4128 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
4129 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
4130 ! --
4131 ! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
4132 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
4133 ! ----------------------------------------------------------------------
4134 WDF = DWSAT * FACTR2 ** EXPON
4135 IF (SICEMAX .gt. 0.0) THEN
4136 VKWGT = 1./ (1. + (500.* SICEMAX)**3.)
4137 WDF = VKWGT * WDF + (1. - VKWGT)* DWSAT * FACTR1** EXPON
4138 ! ----------------------------------------------------------------------
4139 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
4140 ! ----------------------------------------------------------------------
4141 END IF
4142 EXPON = (2.0 * BEXP) + 3.0
4143 WCND = DKSAT * FACTR2 ** EXPON
4145 ! ----------------------------------------------------------------------
4146 END SUBROUTINE WDFCND
4147 ! ----------------------------------------------------------------------
4149 SUBROUTINE SFCDIF_off (ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)
4151 ! ----------------------------------------------------------------------
4152 ! SUBROUTINE SFCDIF (renamed SFCDIF_off to avoid clash with Eta PBL)
4153 ! ----------------------------------------------------------------------
4154 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
4155 ! SEE CHEN ET AL (1997, BLM)
4156 ! ----------------------------------------------------------------------
4158 IMPLICIT NONE
4159 REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
4160 REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
4161 & SQVISC
4162 REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
4163 & PSLHS
4164 REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
4165 REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
4166 REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
4167 REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
4168 !CC ......REAL ZTFC
4170 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
4171 & RLMA
4173 INTEGER ITRMX, ILECH, ITR
4174 PARAMETER &
4175 & (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
4176 & EXCM = 0.001 &
4177 & ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
4178 & ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
4179 & PIHF = 3.14159265/2.)
4180 PARAMETER &
4181 & (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
4182 & ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
4183 & ,SQVISC = 258.2)
4184 PARAMETER &
4185 & (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
4186 & ,RFAC = RIC / (FHNEU * RFC * RFC))
4188 ! ----------------------------------------------------------------------
4189 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
4190 ! ----------------------------------------------------------------------
4191 ! LECH'S SURFACE FUNCTIONS
4192 ! ----------------------------------------------------------------------
4193 PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
4194 PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
4195 PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
4197 ! ----------------------------------------------------------------------
4198 ! PAULSON'S SURFACE FUNCTIONS
4199 ! ----------------------------------------------------------------------
4200 PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
4201 PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) &
4202 & +2.* ATAN (XX) &
4203 &- PIHF
4204 PSPMS (YY)= 5.* YY
4205 PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
4207 ! ----------------------------------------------------------------------
4208 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
4209 ! OVER SOLID SURFACE (LAND, SEA-ICE).
4210 ! ----------------------------------------------------------------------
4211 PSPHS (YY)= 5.* YY
4213 ! ----------------------------------------------------------------------
4214 ! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
4215 ! C......ZTFC=0.1
4216 ! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
4217 ! ----------------------------------------------------------------------
4218 ILECH = 0
4220 ! ----------------------------------------------------------------------
4221 ZILFC = - CZIL * VKRM * SQVISC
4222 ! C.......ZT=Z0*ZTFC
4223 ZU = Z0
4224 RDZ = 1./ ZLM
4225 CXCH = EXCM * RDZ
4226 DTHV = THLM - THZ0
4228 ! ----------------------------------------------------------------------
4229 ! BELJARS CORRECTION OF USTAR
4230 ! ----------------------------------------------------------------------
4231 DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
4232 !cc If statements to avoid TANGENT LINEAR problems near zero
4233 BTGH = BTG * HPBL
4234 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
4235 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
4236 ELSE
4237 WSTAR2 = 0.0
4238 END IF
4240 ! ----------------------------------------------------------------------
4241 ! ZILITINKEVITCH APPROACH FOR ZT
4242 ! ----------------------------------------------------------------------
4243 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
4245 ! ----------------------------------------------------------------------
4246 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
4247 ZSLU = ZLM + ZU
4248 ! PRINT*,'ZSLT=',ZSLT
4249 ! PRINT*,'ZLM=',ZLM
4250 ! PRINT*,'ZT=',ZT
4252 ZSLT = ZLM + ZT
4253 RLOGU = log (ZSLU / ZU)
4255 RLOGT = log (ZSLT / ZT)
4256 ! PRINT*,'RLMO=',RLMO
4257 ! PRINT*,'ELFC=',ELFC
4258 ! PRINT*,'AKHS=',AKHS
4259 ! PRINT*,'DTHV=',DTHV
4260 ! PRINT*,'USTAR=',USTAR
4262 RLMO = ELFC * AKHS * DTHV / USTAR **3
4263 ! ----------------------------------------------------------------------
4264 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
4265 ! ----------------------------------------------------------------------
4266 DO ITR = 1,ITRMX
4267 ZETALT = MAX (ZSLT * RLMO,ZTMIN)
4268 RLMO = ZETALT / ZSLT
4269 ZETALU = ZSLU * RLMO
4270 ZETAU = ZU * RLMO
4272 ZETAT = ZT * RLMO
4273 IF (ILECH .eq. 0) THEN
4274 IF (RLMO .lt. 0.)THEN
4275 XLU4 = 1. -16.* ZETALU
4276 XLT4 = 1. -16.* ZETALT
4277 XU4 = 1. -16.* ZETAU
4279 XT4 = 1. -16.* ZETAT
4280 XLU = SQRT (SQRT (XLU4))
4281 XLT = SQRT (SQRT (XLT4))
4282 XU = SQRT (SQRT (XU4))
4284 XT = SQRT (SQRT (XT4))
4285 ! PRINT*,'-----------1------------'
4286 ! PRINT*,'PSMZ=',PSMZ
4287 ! PRINT*,'PSPMU(ZETAU)=',PSPMU(ZETAU)
4288 ! PRINT*,'XU=',XU
4289 ! PRINT*,'------------------------'
4290 PSMZ = PSPMU (XU)
4291 SIMM = PSPMU (XLU) - PSMZ + RLOGU
4292 PSHZ = PSPHU (XT)
4293 SIMH = PSPHU (XLT) - PSHZ + RLOGT
4294 ELSE
4295 ZETALU = MIN (ZETALU,ZTMAX)
4296 ZETALT = MIN (ZETALT,ZTMAX)
4297 ! PRINT*,'-----------2------------'
4298 ! PRINT*,'PSMZ=',PSMZ
4299 ! PRINT*,'PSPMS(ZETAU)=',PSPMS(ZETAU)
4300 ! PRINT*,'ZETAU=',ZETAU
4301 ! PRINT*,'------------------------'
4302 PSMZ = PSPMS (ZETAU)
4303 SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
4304 PSHZ = PSPHS (ZETAT)
4305 SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
4306 END IF
4307 ! ----------------------------------------------------------------------
4308 ! LECH'S FUNCTIONS
4309 ! ----------------------------------------------------------------------
4310 ELSE
4311 IF (RLMO .lt. 0.)THEN
4312 ! PRINT*,'-----------3------------'
4313 ! PRINT*,'PSMZ=',PSMZ
4314 ! PRINT*,'PSLMU(ZETAU)=',PSLMU(ZETAU)
4315 ! PRINT*,'ZETAU=',ZETAU
4316 ! PRINT*,'------------------------'
4317 PSMZ = PSLMU (ZETAU)
4318 SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
4319 PSHZ = PSLHU (ZETAT)
4320 SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
4321 ELSE
4322 ZETALU = MIN (ZETALU,ZTMAX)
4324 ZETALT = MIN (ZETALT,ZTMAX)
4325 ! PRINT*,'-----------4------------'
4326 ! PRINT*,'PSMZ=',PSMZ
4327 ! PRINT*,'PSLMS(ZETAU)=',PSLMS(ZETAU)
4328 ! PRINT*,'ZETAU=',ZETAU
4329 ! PRINT*,'------------------------'
4330 PSMZ = PSLMS (ZETAU)
4331 SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
4332 PSHZ = PSLHS (ZETAT)
4333 SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
4334 END IF
4335 ! ----------------------------------------------------------------------
4336 ! BELJAARS CORRECTION FOR USTAR
4337 ! ----------------------------------------------------------------------
4338 END IF
4340 ! ----------------------------------------------------------------------
4341 ! ZILITINKEVITCH FIX FOR ZT
4342 ! ----------------------------------------------------------------------
4343 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
4345 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
4346 ZSLT = ZLM + ZT
4347 !-----------------------------------------------------------------------
4348 RLOGT = log (ZSLT / ZT)
4349 USTARK = USTAR * VKRM
4350 AKMS = MAX (USTARK / SIMM,CXCH)
4351 !-----------------------------------------------------------------------
4352 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
4353 !-----------------------------------------------------------------------
4354 AKHS = MAX (USTARK / SIMH,CXCH)
4355 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
4356 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
4357 ELSE
4358 WSTAR2 = 0.0
4359 END IF
4360 !-----------------------------------------------------------------------
4361 RLMN = ELFC * AKHS * DTHV / USTAR **3
4362 !-----------------------------------------------------------------------
4363 ! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
4364 !-----------------------------------------------------------------------
4365 RLMA = RLMO * WOLD+ RLMN * WNEW
4366 !-----------------------------------------------------------------------
4367 RLMO = RLMA
4368 ! PRINT*,'----------------------------'
4369 ! PRINT*,'SFCDIF OUTPUT ! ! ! ! ! ! ! ! ! ! ! !'
4371 ! PRINT*,'ZLM=',ZLM
4372 ! PRINT*,'Z0=',Z0
4373 ! PRINT*,'THZ0=',THZ0
4374 ! PRINT*,'THLM=',THLM
4375 ! PRINT*,'SFCSPD=',SFCSPD
4376 ! PRINT*,'CZIL=',CZIL
4377 ! PRINT*,'AKMS=',AKMS
4378 ! PRINT*,'AKHS=',AKHS
4379 ! PRINT*,'----------------------------'
4381 END DO
4382 ! ----------------------------------------------------------------------
4383 END SUBROUTINE SFCDIF_off
4384 ! ----------------------------------------------------------------------
4386 END MODULE module_sf_noahlsm