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