module_sf_slab.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_sf_slab
4
5 !---SPECIFY CONSTANTS AND LAYERS FOR SOIL MODEL
6 !---SOIL DIFFUSION CONSTANT SET (M^2/S)
7
8 REAL, PARAMETER :: DIFSL=5.e-7
9
10 !---FACTOR TO MAKE SOIL STEP MORE CONSERVATIVE
11
12 REAL , PARAMETER :: SOILFAC=1.25
13
14 CONTAINS
15
16 !----------------------------------------------------------------
17 SUBROUTINE SLAB(T3D,QV3D,P3D,FLHC,FLQC, &
18 PSFC,XLAND,TMN,HFX,QFX,LH,TSK,QSFC,CHKLOWQ, &
19 GSW,GLW,CAPG,THC,SNOWC,EMISS,MAVAIL, &
20 DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
21 SVP1,SVP2,SVP3,SVPT0,EP2, &
22 KARMAN,EOMEG,STBOLT, &
23 TSLB,ZS,DZS,num_soil_layers,radiation, &
24 P1000mb, &
25 ids,ide, jds,jde, kds,kde, &
26 ims,ime, jms,jme, kms,kme, &
27 its,ite, jts,jte, kts,kte )
28 !----------------------------------------------------------------
29 IMPLICIT NONE
30 !----------------------------------------------------------------
31 !
32 ! SUBROUTINE SLAB CALCULATES THE GROUND TEMPERATURE TENDENCY
33 ! ACCORDING TO THE RESIDUAL OF THE SURFACE ENERGY BUDGET
34 ! (BLACKADAR, 1978B).
35 !
36 ! CHANGES:
37 ! FOR SOIL SUB-TIMESTEPS UPDATE SURFACE HFX AND QFX AS TG
38 ! CHANGES TO PREVENT POSSIBLE INSTABILITY FOR LONG MODEL
39 ! STEPS (DT > ~200 SEC).
40 !
41 ! PUT SNOW COVER CHECK ON SOIL SUB-TIMESTEPS
42 !
43 ! MAKE UPPER LIMIT ON SOIL SUB-STEP LENGTH MORE CONSERVATIVE
44 !
45 !----------------------------------------------------------------
46 !-- T3D temperature (K)
47 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
48 !-- P3D 3D pressure (Pa)
49 !-- FLHC exchange coefficient for heat (m/s)
50 !-- FLQC exchange coefficient for moisture (m/s)
51 !-- PSFC surface pressure (Pa)
52 !-- XLAND land mask (1 for land, 2 for water)
53 !-- TMN soil temperature at lower boundary (K)
54 !-- HFX upward heat flux at the surface (W/m^2)
55 !-- QFX upward moisture flux at the surface (kg/m^2/s)
56 !-- LH latent heat flux at the surface (W/m^2)
57 !-- TSK surface temperature (K)
58 !-- GSW downward short wave flux at ground surface (W/m^2)
59 !-- GLW downward long wave flux at ground surface (W/m^2)
60 !-- CAPG heat capacity for soil (J/K/m^3)
61 !-- THC thermal inertia (Cal/cm/K/s^0.5)
62 !-- SNOWC flag indicating snow coverage (1 for snow cover)
63 !-- EMISS surface emissivity (between 0 and 1)
64 !-- DELTSM time step (second)
65 !-- ROVCP R/CP
66 !-- XLV latent heat of melting (J/kg)
67 !-- DTMIN time step (minute)
68 !-- IFSNOW ifsnow=1 for snow-cover effects
69 !-- SVP1 constant for saturation vapor pressure (kPa)
70 !-- SVP2 constant for saturation vapor pressure (dimensionless)
71 !-- SVP3 constant for saturation vapor pressure (K)
72 !-- SVPT0 constant for saturation vapor pressure (K)
73 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
74 !-- EP2 constant for specific humidity calculation
75 ! (R_d/R_v) (dimensionless)
76 !-- KARMAN Von Karman constant
77 !-- EOMEG angular velocity of earth's rotation (rad/s)
78 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
79 !-- TSLB soil temperature in 5-layer model
80 !-- ZS depths of centers of soil layers
81 !-- DZS thicknesses of soil layers
82 !-- num_soil_layers the number of soil layers
83 !-- ids start index for i in domain
84 !-- ide end index for i in domain
85 !-- jds start index for j in domain
86 !-- jde end index for j in domain
87 !-- kds start index for k in domain
88 !-- kde end index for k in domain
89 !-- ims start index for i in memory
90 !-- ime end index for i in memory
91 !-- jms start index for j in memory
92 !-- jme end index for j in memory
93 !-- kms start index for k in memory
94 !-- kme end index for k in memory
95 !-- its start index for i in tile
96 !-- ite end index for i in tile
97 !-- jts start index for j in tile
98 !-- jte end index for j in tile
99 !-- kts start index for k in tile
100 !-- kte end index for k in tile
101 !----------------------------------------------------------------
102 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
103 ims,ime, jms,jme, kms,kme, &
104 its,ite, jts,jte, kts,kte
105
106 INTEGER, INTENT(IN) :: num_soil_layers
107 LOGICAL, INTENT(IN) :: radiation
108
109 INTEGER, INTENT(IN ) :: IFSNOW
110
111 !
112 REAL, INTENT(IN ) :: DTMIN,XLV,ROVCP,DELTSM
113
114 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
115 REAL, INTENT(IN ) :: EP2,KARMAN,EOMEG,STBOLT
116 REAL, INTENT(IN ) :: P1000mb
117
118 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
119 INTENT(INOUT) :: TSLB
120
121 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::ZS,DZS
122
123 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
124 INTENT(IN ) :: QV3D, &
125 P3D, &
126 T3D
127 !
128 REAL, DIMENSION( ims:ime, jms:jme ) , &
129 INTENT(IN ) :: SNOWC, &
130 XLAND, &
131 EMISS, &
132 MAVAIL, &
133 TMN, &
134 GSW, &
135 GLW, &
136 THC
137
138 !CHKLOWQ is declared as memory size
139 !
140 REAL, DIMENSION( ims:ime, jms:jme ) , &
141 INTENT(INOUT) :: HFX, &
142 QFX, &
143 LH, &
144 CAPG, &
145 TSK, &
146 QSFC, &
147 CHKLOWQ
148
149 REAL, DIMENSION( ims:ime, jms:jme ) , &
150 INTENT(IN ) :: PSFC
151 !
152 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
153 FLHC, &
154 FLQC
155
156 ! LOCAL VARS
157
158 REAL, DIMENSION( its:ite ) :: QV1D, &
159 P1D, &
160 T1D
161 INTEGER :: I,J
162
163 DO J=jts,jte
164
165 DO i=its,ite
166 T1D(i) =T3D(i,1,j)
167 QV1D(i)=QV3D(i,1,j)
168 P1D(i) =P3D(i,1,j)
169 ENDDO
170
171 ! the indices to the PSFC argument in the following call look
172 ! wrong; however, it is correct to call with its (and not ims)
173 ! because of the way PSFC is defined in SLAB1D. Whether *that*
174 ! is a good idea or not, this commenter cannot comment. JM
175
176 CALL SLAB1D(J,T1D,QV1D,P1D,FLHC(ims,j),FLQC(ims,j), &
177 PSFC(its,j),XLAND(ims,j),TMN(ims,j),HFX(ims,j), &
178 QFX(ims,j),TSK(ims,j),QSFC(ims,j),CHKLOWQ(ims,j), &
179 LH(ims,j),GSW(ims,j),GLW(ims,j), &
180 CAPG(ims,j),THC(ims,j),SNOWC(ims,j),EMISS(ims,j), &
181 MAVAIL(ims,j),DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
182 SVP1,SVP2,SVP3,SVPT0,EP2,KARMAN,EOMEG,STBOLT, &
183 TSLB(ims,1,j),ZS,DZS,num_soil_layers,radiation, &
184 P1000mb, &
185 ids,ide, jds,jde, kds,kde, &
186 ims,ime, jms,jme, kms,kme, &
187 its,ite, jts,jte, kts,kte )
188
189 ENDDO
190
191 END SUBROUTINE SLAB
192
193 !----------------------------------------------------------------
194 SUBROUTINE SLAB1D(J,T1D,QV1D,P1D,FLHC,FLQC, &
195 PSFCPA,XLAND,TMN,HFX,QFX,TSK,QSFC,CHKLOWQ, &
196 LH,GSW,GLW,CAPG,THC,SNOWC,EMISS,MAVAIL, &
197 DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
198 SVP1,SVP2,SVP3,SVPT0,EP2, &
199 KARMAN,EOMEG,STBOLT, &
200 TSLB2D,ZS,DZS,num_soil_layers,radiation, &
201 P1000mb, &
202 ids,ide, jds,jde, kds,kde, &
203 ims,ime, jms,jme, kms,kme, &
204 its,ite, jts,jte, kts,kte )
205 !----------------------------------------------------------------
206 IMPLICIT NONE
207 !----------------------------------------------------------------
208 !
209 ! SUBROUTINE SLAB CALCULATES THE GROUND TEMPERATURE TENDENCY
210 ! ACCORDING TO THE RESIDUAL OF THE SURFACE ENERGY BUDGET
211 ! (BLACKADAR, 1978B).
212 !
213 ! CHANGES:
214 ! FOR SOIL SUB-TIMESTEPS UPDATE SURFACE HFX AND QFX AS TG
215 ! CHANGES TO PREVENT POSSIBLE INSTABILITY FOR LONG MODEL
216 ! STEPS (DT > ~200 SEC).
217 !
218 ! PUT SNOW COVER CHECK ON SOIL SUB-TIMESTEPS
219 !
220 ! MAKE UPPER LIMIT ON SOIL SUB-STEP LENGTH MORE CONSERVATIVE
221 !
222 !----------------------------------------------------------------
223
224 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
225 ims,ime, jms,jme, kms,kme, &
226 its,ite, jts,jte, kts,kte,J
227
228 INTEGER , INTENT(IN) :: num_soil_layers
229 LOGICAL, INTENT(IN ) :: radiation
230
231 INTEGER, INTENT(IN ) :: IFSNOW
232 !
233 REAL, INTENT(IN ) :: DTMIN,XLV,ROVCP,DELTSM
234
235 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
236 REAL, INTENT(IN ) :: EP2,KARMAN,EOMEG,STBOLT
237 REAL, INTENT(IN ) :: P1000mb
238
239 REAL, DIMENSION( ims:ime , 1:num_soil_layers ), &
240 INTENT(INOUT) :: TSLB2D
241
242 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::ZS,DZS
243
244 !
245 REAL, DIMENSION( ims:ime ) , &
246 INTENT(INOUT) :: HFX, &
247 QFX, &
248 LH, &
249 CAPG, &
250 TSK, &
251 QSFC, &
252 CHKLOWQ
253 !
254 REAL, DIMENSION( ims:ime ) , &
255 INTENT(IN ) :: SNOWC, &
256 XLAND, &
257 EMISS, &
258 MAVAIL, &
259 TMN, &
260 GSW, &
261 GLW, &
262 THC
263 !
264 REAL, DIMENSION( its:ite ) , &
265 INTENT(IN ) :: QV1D, &
266 P1D, &
267 T1D
268 !
269 REAL, DIMENSION( its:ite ) , &
270 INTENT(IN ) :: PSFCPA
271
272 !
273 REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: &
274 FLHC, &
275 FLQC
276 ! LOCAL VARS
277
278 REAL, DIMENSION( its:ite ) :: PSFC
279
280 REAL, DIMENSION( its:ite ) :: &
281 THX, &
282 QX, &
283 SCR3
284
285 REAL, DIMENSION( its:ite ) :: DTHGDT, &
286 TG0, &
287 THTMN, &
288 XLD1, &
289 TSCVN, &
290 OLTG, &
291 UPFLUX, &
292 HM, &
293 RNET, &
294 XINET, &
295 QS, &
296 DTSDT
297 !
298 REAL, DIMENSION( its:ite, num_soil_layers ) :: FLUX
299 !
300 INTEGER :: I,K,NSOIL,ITSOIL,L,NK,RADSWTCH
301 REAL :: PS,PS1,XLDCOL,TSKX,RNSOIL,RHOG1,RHOG2,RHOG3,LAMDAG
302 REAL :: THG,ESG,QSG,HFXT,QFXT,CS,CSW,LAMG(4),THCON,PL
303
304 !----------------------------------------------------------------------
305 !-----DETERMINE IF ANY POINTS IN COLUMN ARE LAND (RATHER THAN OCEAN)
306 ! POINTS. IF NOT, SKIP DOWN TO THE PRINT STATEMENTS SINCE OCEAN
307 ! SURFACE TEMPERATURES ARE NOT ALLOWED TO CHANGE.
308 !
309 ! from sfcrad
310 !----------------------------------------------------------------------
311 DATA CSW/4.183E6/
312 DATA LAMG/1.407E-8, -1.455E-5, 6.290E-3, 0.16857/
313
314 DO i=its,ite
315 ! in cmb
316 PSFC(I)=PSFCPA(I)/1000.
317 ENDDO
318
319
320 DO I=its,ite
321 ! PL cmb
322 PL=P1D(I)/1000.
323 SCR3(I)=T1D(I)
324 ! THCON=(100./PL)**ROVCP
325 THCON=(P1000mb*0.001/PL)**ROVCP
326 THX(I)=SCR3(I)*THCON
327 QX(I)=0.
328 ENDDO
329
330 ! IF(IDRY.EQ.1) GOTO 81
331 DO I=its,ite
332 QX(I)=QV1D(I)
333 ENDDO
334 81 CONTINUE
335
336 !
337 !-----THE SLAB THERMAL CAPACITY CAPG(I) ARE DEPENDENT ON:
338 ! THC(I) - SOIL THERMAL INERTIAL, ONLY.
339 !
340 DO I=its,ite
341 CAPG(I)=3.298E6*THC(I)
342 IF(num_soil_layers .gt. 1)THEN
343
344 ! CAPG REPRESENTS SOIL HEAT CAPACITY (J/K/M^3) WHEN DIFSL=5.E-7 (M^2/S)
345 ! TO GIVE A CORRECT THERMAL INERTIA (=CAPG*DIFSL^0.5)
346
347 CAPG(I)=5.9114E7*THC(I)
348 ENDIF
349 ENDDO
350 !
351 XLDCOL=2.0
352 DO 10 I=its,ite
353 XLDCOL=AMIN1(XLDCOL,XLAND(I))
354 10 CONTINUE
355 !
356 IF(XLDCOL.GT.1.5)GOTO 90
357 !
358 !
359 !-----CONVERT SLAB TEMPERATURE TO POTENTIAL TEMPERATURE AND
360 ! SET XLD1(I) = 0. FOR OCEAN POINTS:
361 !
362 !
363 DO 20 I=its,ite
364 IF((XLAND(I)-1.5).GE.0)THEN
365 XLD1(I)=0.
366 ELSE
367 XLD1(I)=1.
368 ENDIF
369 20 CONTINUE
370 !
371 !-----CONVERT 'TSK(THETAG)' TO 'TG' FOR 'IUP' CALCULATION ....
372 ! IF WE ARE USING THE BLACKADAR MULTI-LEVEL (HIGH-RESOLUTION)
373 ! PBL MODEL
374 !
375 DO 50 I=its,ite
376 IF(XLD1(I).LT.0.5)GOTO 50
377
378 ! PS cmb
379 PS=PSFC(I)
380
381 ! TSK is Temperature at gound sfc
382 ! TG0(I)=TSK(I)*(PS*0.01)**ROVCP
383 TG0(I)=TSK(I)
384 50 CONTINUE
385 !
386 !-----COMPUTE THE SURFACE ENERGY BUDGET:
387 !
388 ! IF(ISOIL.EQ.1)NSOIL=1
389 IF(num_soil_layers .gt. 1)NSOIL=1
390
391
392 IF (radiation) then
393 RADSWTCH=1
394 ELSE
395 RADSWTCH=0
396 ENDIF
397
398 DO 70 I=its,ite
399 IF(XLD1(I).LT.0.5)GOTO 70
400 ! OLTG(I)=TSK(I)*(100./PSFC(I))**ROVCP
401 OLTG(I)=TSK(I)*(P1000mb*0.001/PSFC(I))**ROVCP
402 UPFLUX(I)=RADSWTCH*STBOLT*TG0(I)**4
403 XINET(I)=EMISS(I)*(GLW(I)-UPFLUX(I))
404 RNET(I)=GSW(I)+XINET(I)
405 HM(I)=1.18*EOMEG*(TG0(I)-TMN(I))
406 ! MOISTURE FLUX CALCULATED HERE (OVERWRITES SFC LAYER VALUE FOR LAND)
407 ESG=SVP1*EXP(SVP2*(TG0(I)-SVPT0)/(TG0(I)-SVP3))
408 QSG=EP2*ESG/(PSFC(I)-ESG)
409 THG=TSK(I)*(100./PSFC(I))**ROVCP
410 HFX(I)=FLHC(I)*(THG-THX(I))
411 QFX(I)=FLQC(I)*(QSG-QX(I))
412 LH(I)=QFX(I)*XLV
413 QS(I)=HFX(I)+QFX(I)*XLV
414 ! IF(ISOIL.EQ.0)THEN
415 IF(num_soil_layers .EQ. 1)THEN
416 DTHGDT(I)=(RNET(I)-QS(I))/CAPG(I)-HM(I)
417 ELSE
418 DTHGDT(I)=0.
419 ENDIF
420 70 CONTINUE
421 ! IF(ISOIL.EQ.1)THEN
422 IF(num_soil_layers .gt. 1)THEN
423 NSOIL=1+IFIX(SOILFAC*4*DIFSL/DZS(1)*DELTSM/DZS(1))
424 RNSOIL=1./FLOAT(NSOIL)
425 !
426 ! SOIL SUB-TIMESTEP
427 !
428 DO ITSOIL=1,NSOIL
429 DO I=its,ite
430 DO L=1,num_soil_layers-1
431 IF(XLD1(I).LT.0.5)GOTO 75
432 IF(L.EQ.1.AND.ITSOIL.GT.1)THEN
433 ! PS1=(PSFC(I)*0.01)**ROVCP
434 PS1=(PSFCPA(I)/P1000mb)**ROVCP
435
436 ! for rk scheme A and B are the same
437 PS=PSFC(I)
438 THG=TSLB2D(I,1)/PS1
439 ESG=SVP1*EXP(SVP2*(TSLB2D(I,1)-SVPT0)/(TSLB2D(I,1) &
440 -SVP3))
441 QSG=EP2*ESG/(PS-ESG)
442 ! UPDATE FLUXES FOR NEW GROUND TEMPERATURE
443 HFXT=FLHC(I)*(THG-THX(I))
444 QFXT=FLQC(I)*(QSG-QX(I))
445 QS(I)=HFXT+QFXT*XLV
446 ! SUM HFX AND QFX OVER SOIL TIMESTEPS
447 HFX(I)=HFX(I)+HFXT
448 QFX(I)=QFX(I)+QFXT
449 ENDIF
450 FLUX(I,1)=RNET(I)-QS(I)
451 FLUX(I,L+1)=-DIFSL*CAPG(I)*(TSLB2D(I,L+1)-TSLB2D(I,L))/( &
452 ZS(L+1)-ZS(L))
453 DTSDT(I)=-(FLUX(I,L+1)-FLUX(I,L))/(DZS(L)*CAPG(I))
454 TSLB2D(I,L)=TSLB2D(I,L)+DTSDT(I)*DELTSM*RNSOIL
455 IF(IFSNOW.EQ.1.AND.L.EQ.1)THEN
456 IF((SNOWC(I).GT.0..AND.TSLB2D(I,1).GT.273.16))THEN
457 TSLB2D(I,1)=273.16
458 ENDIF
459 ENDIF
460 IF(L.EQ.1)DTHGDT(I)=DTHGDT(I)+RNSOIL*DTSDT(I)
461 IF(ITSOIL.EQ.NSOIL.AND.L.EQ.1)THEN
462 ! AVERAGE HFX AND QFX OVER SOIL TIMESTEPS FOR OUTPUT TO PBL
463 HFX(I)=HFX(I)*RNSOIL
464 QFX(I)=QFX(I)*RNSOIL
465 LH(I)=QFX(I)*XLV
466 ENDIF
467 75 CONTINUE
468 ENDDO
469 ENDDO
470 ENDDO
471 ENDIF
472 !
473 DO 80 I=its,ite
474 IF(XLD1(I).LT.0.5) GOTO 80
475 TSKX=TG0(I)+DELTSM*DTHGDT(I)
476
477 ! TSK is temperature
478 ! TSK(I)=TSKX*(100./PS1)**ROVCP
479 TSK(I)=TSKX
480 80 CONTINUE
481
482 !
483 !-----MODIFY THE THE GROUND TEMPERATURE IF THE SNOW COVER EFFECTS ARE
484 ! CONSIDERED: LIMIT THE GROUND TEMPERATURE UNDER 0 C.
485 !
486 IF(IFSNOW.EQ.0)GOTO 90
487 DO 85 I=its,ite
488 IF(XLD1(I).LT.0.5)GOTO 85
489 ! PS1=(PSFC(I)*0.01)**ROVCP
490 ! TSCVN(I)=TSK(I)*PS1
491 TSCVN(I)=TSK(I)
492 IF((SNOWC(I).GT.0..AND.TSCVN(I).GT.273.16))THEN
493 TSCVN(I)=273.16
494 ELSE
495 TSCVN(I)=TSCVN(I)
496 ENDIF
497 ! TSK(I)=TSCVN(I)/PS1
498 TSK(I)=TSCVN(I)
499 85 CONTINUE
500 !
501 90 CONTINUE
502 DO I=its,ite
503 ! QSFC and CHKLOWQ needed by Eta PBL
504 QSFC(I)=QX(I)+QFX(I)/FLQC(I)
505 CHKLOWQ(I)=MAVAIL(I)
506 ENDDO
507 !
508 140 CONTINUE
509
510 END SUBROUTINE SLAB1D
511
512 !================================================================
513 SUBROUTINE slabinit(TSK,TMN, &
514 TSLB,ZS,DZS,num_soil_layers, &
515 restart, allowed_to_read, &
516 ids,ide, jds,jde, kds,kde, &
517 ims,ime, jms,jme, kms,kme, &
518 its,ite, jts,jte, kts,kte )
519 !----------------------------------------------------------------
520 IMPLICIT NONE
521 !----------------------------------------------------------------
522 LOGICAL , INTENT(IN) :: restart, allowed_to_read
523 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
524 ims,ime, jms,jme, kms,kme, &
525 its,ite, jts,jte, kts,kte
526
527 INTEGER, INTENT(IN ) :: num_soil_layers
528 !
529 REAL, DIMENSION( ims:ime , 1:num_soil_layers , jms:jme ), INTENT(INOUT) :: TSLB
530
531 REAL, DIMENSION(1:num_soil_layers), INTENT(IN) :: ZS,DZS
532
533 REAL, DIMENSION( ims:ime, jms:jme ) , &
534 INTENT(IN) :: TSK, &
535 TMN
536 ! LOCAR VAR
537
538 INTEGER :: L,J,I,itf,jtf
539 !----------------------------------------------------------------
540
541 itf=min0(ite,ide-1)
542 jtf=min0(jte,jde-1)
543
544 END SUBROUTINE slabinit
545
546 !-------------------------------------------------------------------
547
548 END MODULE module_sf_slab