Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_slab.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 33.5 KB

Line
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	tml,t0ml,hml,h0ml,huml,hvml,ust,u_phy,v_phy, &
29	f,g,omlcall,oml_gamma )
30	!----------------------------------------------------------------
31	IMPLICIT NONE
32	!----------------------------------------------------------------
33	!
34	! SUBROUTINE SLAB CALCULATES THE GROUND TEMPERATURE TENDENCY
35	! ACCORDING TO THE RESIDUAL OF THE SURFACE ENERGY BUDGET
36	! (BLACKADAR, 1978B).
37	!
38	! CHANGES:
39	! FOR SOIL SUB-TIMESTEPS UPDATE SURFACE HFX AND QFX AS TG
40	! CHANGES TO PREVENT POSSIBLE INSTABILITY FOR LONG MODEL
41	! STEPS (DT > ~200 SEC).
42	!
43	! PUT SNOW COVER CHECK ON SOIL SUB-TIMESTEPS
44	!
45	! MAKE UPPER LIMIT ON SOIL SUB-STEP LENGTH MORE CONSERVATIVE
46	!
47	!----------------------------------------------------------------
48	!-- T3D temperature (K)
49	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
50	!-- P3D 3D pressure (Pa)
51	!-- FLHC exchange coefficient for heat (m/s)
52	!-- FLQC exchange coefficient for moisture (m/s)
53	!-- PSFC surface pressure (Pa)
54	!-- XLAND land mask (1 for land, 2 for water)
55	!-- TMN soil temperature at lower boundary (K)
56	!-- HFX upward heat flux at the surface (W/m^2)
57	!-- QFX upward moisture flux at the surface (kg/m^2/s)
58	!-- LH latent heat flux at the surface (W/m^2)
59	!-- TSK surface temperature (K)
60	!-- GSW downward short wave flux at ground surface (W/m^2)
61	!-- GLW downward long wave flux at ground surface (W/m^2)
62	!-- CAPG heat capacity for soil (J/K/m^3)
63	!-- THC thermal inertia (Cal/cm/K/s^0.5)
64	!-- SNOWC flag indicating snow coverage (1 for snow cover)
65	!-- EMISS surface emissivity (between 0 and 1)
66	!-- DELTSM time step (second)
67	!-- ROVCP R/CP
68	!-- XLV latent heat of melting (J/kg)
69	!-- DTMIN time step (minute)
70	!-- IFSNOW ifsnow=1 for snow-cover effects
71	!-- SVP1 constant for saturation vapor pressure (kPa)
72	!-- SVP2 constant for saturation vapor pressure (dimensionless)
73	!-- SVP3 constant for saturation vapor pressure (K)
74	!-- SVPT0 constant for saturation vapor pressure (K)
75	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
76	!-- EP2 constant for specific humidity calculation
77	! (R_d/R_v) (dimensionless)
78	!-- KARMAN Von Karman constant
79	!-- EOMEG angular velocity of earth's rotation (rad/s)
80	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
81	!-- TSLB soil temperature in 5-layer model
82	!-- ZS depths of centers of soil layers
83	!-- DZS thicknesses of soil layers
84	!-- TML ocean mixed layer temperature (K)
85	!-- T0ML ocean mixed layer temperature (K) at initial time
86	!-- HML ocean mixed layer depth (m)
87	!-- H0ML ocean mixed layer depth (m) at initial time
88	!-- HUML ocean mixed layer u component of wind
89	!-- HVML ocean mixed layer v component of wind
90	!-- OML_GAMMA deep water lapse rate (K m-1)
91	!-- OMLCALL whether to call oml model
92	!-- num_soil_layers the number of soil layers
93	!-- ids start index for i in domain
94	!-- ide end index for i in domain
95	!-- jds start index for j in domain
96	!-- jde end index for j in domain
97	!-- kds start index for k in domain
98	!-- kde end index for k in domain
99	!-- ims start index for i in memory
100	!-- ime end index for i in memory
101	!-- jms start index for j in memory
102	!-- jme end index for j in memory
103	!-- kms start index for k in memory
104	!-- kme end index for k in memory
105	!-- its start index for i in tile
106	!-- ite end index for i in tile
107	!-- jts start index for j in tile
108	!-- jte end index for j in tile
109	!-- kts start index for k in tile
110	!-- kte end index for k in tile
111	!----------------------------------------------------------------
112	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
113	ims,ime, jms,jme, kms,kme, &
114	its,ite, jts,jte, kts,kte
115
116	INTEGER, INTENT(IN) :: num_soil_layers
117	LOGICAL, INTENT(IN) :: radiation
118
119	INTEGER, INTENT(IN ) :: IFSNOW
120
121	!
122	REAL, INTENT(IN ) :: DTMIN,XLV,ROVCP,DELTSM
123
124	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
125	REAL, INTENT(IN ) :: EP2,KARMAN,EOMEG,STBOLT
126	REAL, INTENT(IN ) :: P1000mb
127
128	REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
129	INTENT(INOUT) :: TSLB
130
131	REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::ZS,DZS
132
133	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
134	INTENT(IN ) :: QV3D, &
135	P3D, &
136	T3D
137	!
138	REAL, DIMENSION( ims:ime, jms:jme ) , &
139	INTENT(IN ) :: SNOWC, &
140	XLAND, &
141	EMISS, &
142	MAVAIL, &
143	TMN, &
144	GSW, &
145	GLW, &
146	THC
147
148	!CHKLOWQ is declared as memory size
149	!
150	REAL, DIMENSION( ims:ime, jms:jme ) , &
151	INTENT(INOUT) :: HFX, &
152	QFX, &
153	LH, &
154	CAPG, &
155	TSK, &
156	QSFC, &
157	CHKLOWQ
158
159	REAL, DIMENSION( ims:ime, jms:jme ) , &
160	INTENT(IN ) :: PSFC
161	!
162	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
163	FLHC, &
164	FLQC
165
166	! Ocean Mixed Layer Vars
167	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) :: &
168	TML,T0ML,HML,H0ML,HUML,HVML
169	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, INTENT(IN ) :: &
170	U_PHY,V_PHY
171	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN ) :: &
172	UST, F
173	REAL, OPTIONAL, INTENT(IN ) :: G
174	REAL, OPTIONAL, INTENT(IN ) :: OML_GAMMA
175	INTEGER, OPTIONAL, INTENT(IN ) :: OMLCALL
176
177	! LOCAL VARS
178
179	REAL, DIMENSION( its:ite ) :: QV1D, &
180	P1D, &
181	T1D
182	INTEGER :: I,J
183
184	DO J=jts,jte
185
186	DO i=its,ite
187	T1D(i) =T3D(i,1,j)
188	QV1D(i)=QV3D(i,1,j)
189	P1D(i) =P3D(i,1,j)
190	ENDDO
191
192	! the indices to the PSFC argument in the following call look
193	! wrong; however, it is correct to call with its (and not ims)
194	! because of the way PSFC is defined in SLAB1D. Whether that
195	! is a good idea or not, this commenter cannot comment. JM
196
197	CALL SLAB1D(J,T1D,QV1D,P1D,FLHC(ims,j),FLQC(ims,j), &
198	PSFC(its,j),XLAND(ims,j),TMN(ims,j),HFX(ims,j), &
199	QFX(ims,j),TSK(ims,j),QSFC(ims,j),CHKLOWQ(ims,j), &
200	LH(ims,j),GSW(ims,j),GLW(ims,j), &
201	CAPG(ims,j),THC(ims,j),SNOWC(ims,j),EMISS(ims,j), &
202	MAVAIL(ims,j),DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
203	SVP1,SVP2,SVP3,SVPT0,EP2,KARMAN,EOMEG,STBOLT, &
204	TSLB(ims,1,j),ZS,DZS,num_soil_layers,radiation, &
205	P1000mb, &
206	ids,ide, jds,jde, kds,kde, &
207	ims,ime, jms,jme, kms,kme, &
208	its,ite, jts,jte, kts,kte )
209
210	IF (OMLCALL .EQ. 1) THEN
211	! CALL wrf_debug( 1, 'Call OCEANML' )
212	IF (PRESENT(tml) .AND. PRESENT(t0ml) &
213	.AND. .TRUE. ) THEN
214	DO i=its,ite
215	IF(XLAND(I,J).GT.1.5)THEN
216	CALL OCEANML(I,J,TML(i,j),T0ML(i,j),HML(i,j),H0ML(i,j),&
217	HUML(i,j),HVML(i,j),TSK(i,j),HFX(i,j), &
218	LH(i,j),GSW(i,j),GLW(i,j), &
219	U_PHY(i,kts,j),V_PHY(i,kts,j),UST(i,j),F(i,j), &
220	EMISS(i,j),STBOLT,G,DELTSM,OML_GAMMA, &
221	ids,ide, jds,jde, kds,kde, &
222	ims,ime, jms,jme, kms,kme, &
223	its,ite, jts,jte, kts,kte )
224	ENDIF
225	ENDDO
226	ELSE
227	CALL wrf_error_fatal('Lacking arguments for OCEANML in slab')
228	ENDIF
229	ENDIF
230
231	ENDDO
232
233	END SUBROUTINE SLAB
234
235	!----------------------------------------------------------------
236	SUBROUTINE SLAB1D(J,T1D,QV1D,P1D,FLHC,FLQC, &
237	PSFCPA,XLAND,TMN,HFX,QFX,TSK,QSFC,CHKLOWQ, &
238	LH,GSW,GLW,CAPG,THC,SNOWC,EMISS,MAVAIL, &
239	DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
240	SVP1,SVP2,SVP3,SVPT0,EP2, &
241	KARMAN,EOMEG,STBOLT, &
242	TSLB2D,ZS,DZS,num_soil_layers,radiation, &
243	P1000mb, &
244	ids,ide, jds,jde, kds,kde, &
245	ims,ime, jms,jme, kms,kme, &
246	its,ite, jts,jte, kts,kte )
247	!----------------------------------------------------------------
248	IMPLICIT NONE
249	!----------------------------------------------------------------
250	!
251	! SUBROUTINE SLAB CALCULATES THE GROUND TEMPERATURE TENDENCY
252	! ACCORDING TO THE RESIDUAL OF THE SURFACE ENERGY BUDGET
253	! (BLACKADAR, 1978B).
254	!
255	! CHANGES:
256	! FOR SOIL SUB-TIMESTEPS UPDATE SURFACE HFX AND QFX AS TG
257	! CHANGES TO PREVENT POSSIBLE INSTABILITY FOR LONG MODEL
258	! STEPS (DT > ~200 SEC).
259	!
260	! PUT SNOW COVER CHECK ON SOIL SUB-TIMESTEPS
261	!
262	! MAKE UPPER LIMIT ON SOIL SUB-STEP LENGTH MORE CONSERVATIVE
263	!
264	!----------------------------------------------------------------
265
266	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
267	ims,ime, jms,jme, kms,kme, &
268	its,ite, jts,jte, kts,kte,J
269
270	INTEGER , INTENT(IN) :: num_soil_layers
271	LOGICAL, INTENT(IN ) :: radiation
272
273	INTEGER, INTENT(IN ) :: IFSNOW
274	!
275	REAL, INTENT(IN ) :: DTMIN,XLV,ROVCP,DELTSM
276
277	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
278	REAL, INTENT(IN ) :: EP2,KARMAN,EOMEG,STBOLT
279	REAL, INTENT(IN ) :: P1000mb
280
281	REAL, DIMENSION( ims:ime , 1:num_soil_layers ), &
282	INTENT(INOUT) :: TSLB2D
283
284	REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::ZS,DZS
285
286	!
287	REAL, DIMENSION( ims:ime ) , &
288	INTENT(INOUT) :: HFX, &
289	QFX, &
290	LH, &
291	CAPG, &
292	TSK, &
293	QSFC, &
294	CHKLOWQ
295	!
296	REAL, DIMENSION( ims:ime ) , &
297	INTENT(IN ) :: SNOWC, &
298	XLAND, &
299	EMISS, &
300	MAVAIL, &
301	TMN, &
302	GSW, &
303	GLW, &
304	THC
305	!
306	REAL, DIMENSION( its:ite ) , &
307	INTENT(IN ) :: QV1D, &
308	P1D, &
309	T1D
310	!
311	REAL, DIMENSION( its:ite ) , &
312	INTENT(IN ) :: PSFCPA
313
314	!
315	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: &
316	FLHC, &
317	FLQC
318	! LOCAL VARS
319
320	REAL, DIMENSION( its:ite ) :: PSFC
321
322	REAL, DIMENSION( its:ite ) :: &
323	THX, &
324	QX, &
325	SCR3
326
327	REAL, DIMENSION( its:ite ) :: DTHGDT, &
328	TG0, &
329	THTMN, &
330	XLD1, &
331	TSCVN, &
332	OLTG, &
333	UPFLUX, &
334	HM, &
335	RNET, &
336	XINET, &
337	QS, &
338	DTSDT
339	!
340	REAL, DIMENSION( its:ite, num_soil_layers ) :: FLUX
341	!
342	INTEGER :: I,K,NSOIL,ITSOIL,L,NK,RADSWTCH
343	REAL :: PS,PS1,XLDCOL,TSKX,RNSOIL,RHOG1,RHOG2,RHOG3,LAMDAG
344	REAL :: THG,ESG,QSG,HFXT,QFXT,CS,CSW,LAMG(4),THCON,PL
345
346	!----------------------------------------------------------------------
347	!-----DETERMINE IF ANY POINTS IN COLUMN ARE LAND (RATHER THAN OCEAN)
348	! POINTS. IF NOT, SKIP DOWN TO THE PRINT STATEMENTS SINCE OCEAN
349	! SURFACE TEMPERATURES ARE NOT ALLOWED TO CHANGE.
350	!
351	! from sfcrad
352	!----------------------------------------------------------------------
353	DATA CSW/4.183E6/
354	DATA LAMG/1.407E-8, -1.455E-5, 6.290E-3, 0.16857/
355
356	DO i=its,ite
357	! in cmb
358	PSFC(I)=PSFCPA(I)/1000.
359	ENDDO
360
361
362	DO I=its,ite
363	! PL cmb
364	PL=P1D(I)/1000.
365	SCR3(I)=T1D(I)
366	! THCON=(100./PL)**ROVCP
367	THCON=(P1000mb0.001/PL)*ROVCP
368	THX(I)=SCR3(I)*THCON
369	QX(I)=0.
370	ENDDO
371
372	! IF(IDRY.EQ.1) GOTO 81
373	DO I=its,ite
374	QX(I)=QV1D(I)
375	ENDDO
376	81 CONTINUE
377
378	!
379	!-----THE SLAB THERMAL CAPACITY CAPG(I) ARE DEPENDENT ON:
380	! THC(I) - SOIL THERMAL INERTIAL, ONLY.
381	!
382	DO I=its,ite
383	CAPG(I)=3.298E6*THC(I)
384	IF(num_soil_layers .gt. 1)THEN
385
386	! CAPG REPRESENTS SOIL HEAT CAPACITY (J/K/M^3) WHEN DIFSL=5.E-7 (M^2/S)
387	! TO GIVE A CORRECT THERMAL INERTIA (=CAPG*DIFSL^0.5)
388
389	CAPG(I)=5.9114E7*THC(I)
390	ENDIF
391	ENDDO
392	!
393	XLDCOL=2.0
394	DO 10 I=its,ite
395	XLDCOL=AMIN1(XLDCOL,XLAND(I))
396	10 CONTINUE
397	!
398	IF(XLDCOL.GT.1.5)GOTO 90
399	!
400	!
401	!-----CONVERT SLAB TEMPERATURE TO POTENTIAL TEMPERATURE AND
402	! SET XLD1(I) = 0. FOR OCEAN POINTS:
403	!
404	!
405	DO 20 I=its,ite
406	IF((XLAND(I)-1.5).GE.0)THEN
407	XLD1(I)=0.
408	ELSE
409	XLD1(I)=1.
410	ENDIF
411	20 CONTINUE
412	!
413	!-----CONVERT 'TSK(THETAG)' TO 'TG' FOR 'IUP' CALCULATION ....
414	! IF WE ARE USING THE BLACKADAR MULTI-LEVEL (HIGH-RESOLUTION)
415	! PBL MODEL
416	!
417	DO 50 I=its,ite
418	IF(XLD1(I).LT.0.5)GOTO 50
419
420	! PS cmb
421	PS=PSFC(I)
422
423	! TSK is Temperature at gound sfc
424	! TG0(I)=TSK(I)(PS0.01)**ROVCP
425	TG0(I)=TSK(I)
426	50 CONTINUE
427	!
428	!-----COMPUTE THE SURFACE ENERGY BUDGET:
429	!
430	! IF(ISOIL.EQ.1)NSOIL=1
431	IF(num_soil_layers .gt. 1)NSOIL=1
432
433
434	IF (radiation) then
435	RADSWTCH=1
436	ELSE
437	RADSWTCH=0
438	ENDIF
439
440	DO 70 I=its,ite
441	IF(XLD1(I).LT.0.5)GOTO 70
442	! OLTG(I)=TSK(I)(100./PSFC(I))*ROVCP
443	OLTG(I)=TSK(I)(P1000mb0.001/PSFC(I))**ROVCP
444	UPFLUX(I)=RADSWTCHSTBOLTTG0(I)**4
445	XINET(I)=EMISS(I)*(GLW(I)-UPFLUX(I))
446	RNET(I)=GSW(I)+XINET(I)
447	HM(I)=1.18EOMEG(TG0(I)-TMN(I))
448	! MOISTURE FLUX CALCULATED HERE (OVERWRITES SFC LAYER VALUE FOR LAND)
449	ESG=SVP1EXP(SVP2(TG0(I)-SVPT0)/(TG0(I)-SVP3))
450	QSG=EP2*ESG/(PSFC(I)-ESG)
451	THG=TSK(I)(100./PSFC(I))*ROVCP
452	HFX(I)=FLHC(I)*(THG-THX(I))
453	QFX(I)=FLQC(I)*(QSG-QX(I))
454	LH(I)=QFX(I)*XLV
455	QS(I)=HFX(I)+QFX(I)*XLV
456	! IF(ISOIL.EQ.0)THEN
457	IF(num_soil_layers .EQ. 1)THEN
458	DTHGDT(I)=(RNET(I)-QS(I))/CAPG(I)-HM(I)
459	ELSE
460	DTHGDT(I)=0.
461	ENDIF
462	70 CONTINUE
463	! IF(ISOIL.EQ.1)THEN
464	IF(num_soil_layers .gt. 1)THEN
465	NSOIL=1+IFIX(SOILFAC4DIFSL/DZS(1)*DELTSM/DZS(1))
466	RNSOIL=1./FLOAT(NSOIL)
467	!
468	! SOIL SUB-TIMESTEP
469	!
470	DO ITSOIL=1,NSOIL
471	DO I=its,ite
472	DO L=1,num_soil_layers-1
473	IF(XLD1(I).LT.0.5)GOTO 75
474	IF(L.EQ.1.AND.ITSOIL.GT.1)THEN
475	! PS1=(PSFC(I)0.01)*ROVCP
476	PS1=(PSFCPA(I)/P1000mb)**ROVCP
477
478	! for rk scheme A and B are the same
479	PS=PSFC(I)
480	THG=TSLB2D(I,1)/PS1
481	ESG=SVP1EXP(SVP2(TSLB2D(I,1)-SVPT0)/(TSLB2D(I,1) &
482	-SVP3))
483	QSG=EP2*ESG/(PS-ESG)
484	! UPDATE FLUXES FOR NEW GROUND TEMPERATURE
485	HFXT=FLHC(I)*(THG-THX(I))
486	QFXT=FLQC(I)*(QSG-QX(I))
487	QS(I)=HFXT+QFXT*XLV
488	! SUM HFX AND QFX OVER SOIL TIMESTEPS
489	HFX(I)=HFX(I)+HFXT
490	QFX(I)=QFX(I)+QFXT
491	ENDIF
492	FLUX(I,1)=RNET(I)-QS(I)
493	FLUX(I,L+1)=-DIFSLCAPG(I)(TSLB2D(I,L+1)-TSLB2D(I,L))/( &
494	ZS(L+1)-ZS(L))
495	DTSDT(I)=-(FLUX(I,L+1)-FLUX(I,L))/(DZS(L)*CAPG(I))
496	TSLB2D(I,L)=TSLB2D(I,L)+DTSDT(I)DELTSMRNSOIL
497	IF(IFSNOW.EQ.1.AND.L.EQ.1)THEN
498	IF((SNOWC(I).GT.0..AND.TSLB2D(I,1).GT.273.16))THEN
499	TSLB2D(I,1)=273.16
500	ENDIF
501	ENDIF
502	IF(L.EQ.1)DTHGDT(I)=DTHGDT(I)+RNSOIL*DTSDT(I)
503	IF(ITSOIL.EQ.NSOIL.AND.L.EQ.1)THEN
504	! AVERAGE HFX AND QFX OVER SOIL TIMESTEPS FOR OUTPUT TO PBL
505	HFX(I)=HFX(I)*RNSOIL
506	QFX(I)=QFX(I)*RNSOIL
507	LH(I)=QFX(I)*XLV
508	ENDIF
509	75 CONTINUE
510	ENDDO
511	ENDDO
512	ENDDO
513	ENDIF
514	!
515	DO 80 I=its,ite
516	IF(XLD1(I).LT.0.5) GOTO 80
517	TSKX=TG0(I)+DELTSM*DTHGDT(I)
518
519	! TSK is temperature
520	! TSK(I)=TSKX(100./PS1)*ROVCP
521	TSK(I)=TSKX
522	80 CONTINUE
523
524	!
525	!-----MODIFY THE THE GROUND TEMPERATURE IF THE SNOW COVER EFFECTS ARE
526	! CONSIDERED: LIMIT THE GROUND TEMPERATURE UNDER 0 C.
527	!
528	IF(IFSNOW.EQ.0)GOTO 90
529	DO 85 I=its,ite
530	IF(XLD1(I).LT.0.5)GOTO 85
531	! PS1=(PSFC(I)0.01)*ROVCP
532	! TSCVN(I)=TSK(I)*PS1
533	TSCVN(I)=TSK(I)
534	IF((SNOWC(I).GT.0..AND.TSCVN(I).GT.273.16))THEN
535	TSCVN(I)=273.16
536	ELSE
537	TSCVN(I)=TSCVN(I)
538	ENDIF
539	! TSK(I)=TSCVN(I)/PS1
540	TSK(I)=TSCVN(I)
541	85 CONTINUE
542	!
543	90 CONTINUE
544	DO I=its,ite
545	! QSFC and CHKLOWQ needed by Eta PBL
546	QSFC(I)=QX(I)+QFX(I)/FLQC(I)
547	CHKLOWQ(I)=MAVAIL(I)
548	ENDDO
549	!
550	140 CONTINUE
551
552	END SUBROUTINE SLAB1D
553
554	!================================================================
555	SUBROUTINE slabinit(TSK,TMN, &
556	TSLB,ZS,DZS,num_soil_layers, &
557	allowed_to_read, start_of_simulation, &
558	ids,ide, jds,jde, kds,kde, &
559	ims,ime, jms,jme, kms,kme, &
560	its,ite, jts,jte, kts,kte, &
561	oml_hml0, omlcall, &
562	tml,t0ml,hml,h0ml,huml,hvml )
563	!----------------------------------------------------------------
564	IMPLICIT NONE
565	!----------------------------------------------------------------
566	LOGICAL , INTENT(IN) :: allowed_to_read
567	LOGICAL , INTENT(IN) :: start_of_simulation
568	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
569	ims,ime, jms,jme, kms,kme, &
570	its,ite, jts,jte, kts,kte
571
572	INTEGER, INTENT(IN ) :: num_soil_layers
573	!
574	REAL, DIMENSION( ims:ime , 1:num_soil_layers , jms:jme ), INTENT(INOUT) :: TSLB
575
576	REAL, DIMENSION(1:num_soil_layers), INTENT(IN) :: ZS,DZS
577
578	REAL, DIMENSION( ims:ime, jms:jme ) , &
579	INTENT(IN) :: TSK, &
580	TMN
581	REAL, DIMENSION( ims:ime, jms:jme ) , OPTIONAL , &
582	INTENT(INOUT) :: TML, T0ML, HML, H0ML, HUML, HVML
583	REAL , OPTIONAL, INTENT(IN ) :: oml_hml0
584	INTEGER, OPTIONAL, INTENT(IN ) :: omlcall
585
586	! LOCAR VAR
587
588	INTEGER :: L,J,I,itf,jtf
589	CHARACTER*1024 message
590
591	!----------------------------------------------------------------
592
593	itf=min0(ite,ide-1)
594	jtf=min0(jte,jde-1)
595
596	IF ( PRESENT(omlcall) .AND. omlcall .EQ. 1) THEN
597	WRITE(message,*)'Initializing OML with HML0 = ', oml_hml0
598	CALL wrf_debug (1, TRIM(message))
599
600	IF(start_of_simulation .AND. &
601	PRESENT(tml) .AND. PRESENT(t0ml) ) THEN
602	DO J=jts,jtf
603	DO I=its,itf
604	TML(I,J)=TSK(I,J)
605	T0ML(I,J)=TSK(I,J)
606	! MAY HAVE INPUT OF HML BUT FOR NOW SET HERE
607	HML(I,J)=oml_hml0
608	H0ML(I,J)=HML(I,J)
609	HUML(I,J)=0.
610	HVML(I,J)=0.
611	ENDDO
612	ENDDO
613	ENDIF
614	ENDIF
615
616	END SUBROUTINE slabinit
617
618	!-------------------------------------------------------------------
619
620	SUBROUTINE OCEANML(I,J,TML,T0ML,H,H0,HUML, &
621	HVML,TSK,HFX, &
622	LH,GSW,GLW, &
623	UAIR,VAIR,UST,F,EMISS,STBOLT,G,DT,OML_GAMMA, &
624	ids,ide, jds,jde, kds,kde, &
625	ims,ime, jms,jme, kms,kme, &
626	its,ite, jts,jte, kts,kte )
627
628	!----------------------------------------------------------------
629	IMPLICIT NONE
630	!----------------------------------------------------------------
631	INTEGER, INTENT(IN ) :: I, J
632	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
633	ims,ime, jms,jme, kms,kme, &
634	its,ite, jts,jte, kts,kte
635
636	REAL, INTENT(INOUT) :: TML, H, H0, HUML, HVML, TSK
637
638	REAL, INTENT(IN ) :: T0ML, HFX, LH, GSW, GLW, &
639	UAIR, VAIR, UST, F, EMISS
640
641	REAL, INTENT(IN) :: STBOLT, G, DT, OML_GAMMA
642
643	! Local
644	REAL :: rhoair, rhowater, Gam, alp, BV2, A1, A2, B2, u, v, wspd, &
645	hu1, hv1, hu2, hv2, taux, tauy, tauxair, tauyair, q, hold, &
646	hsqrd, thp, cwater
647	CHARACTER(LEN=120) :: time_series
648
649	hu1=huml
650	hv1=hvml
651	rhoair=1.
652	rhowater=1000.
653	cwater=4200.
654	! Deep ocean lapse rate (K/m) - from Rich
655	Gam=oml_gamma
656	! if(i.eq.1 .and. j.eq.1 .or. i.eq.105.and.j.eq.105) print *, 'gamma = ', gam
657	! Gam=0.14
658	! Gam=5.6/40.
659	! if(i.eq.1 .and. j.eq.1 ) print *, 'gamma = ', gam
660	! Gam=5./100.
661	! Thermal expansion coeff (/K)
662	! alp=.0002
663	! temp dependence (/K)
664	alp=max((tml-273.15)*1.e-5, 1.e-6)
665	BV2=alpgGam
666	thp=t0ml-Gam*(h-h0)
667	A1=(tml-thp)h - 0.5Gamhh
668	if(h.ne.0.)then
669	u=hu1/h
670	v=hv1/h
671	else
672	u=0.
673	v=0.
674	endif
675
676	! time step
677
678	! q=(-hfx-lh+gsw+glw-stboltemisstmltmltmltml)/(rhowatercwater)
679	! wspd=max(sqrt(uairuair+vairvair),0.1)
680	wspd=sqrt(uairuair+vairvair)
681	if (wspd .lt. 1.e-10 ) then
682	! print *, 'i,j,wspd are ', i,j,wspd
683	wspd = 1.e-10
684	endif
685	tauxair=ustustuair/wspd
686	taux=rhoair/rhowater*tauxair
687	tauyair=ustustvair/wspd
688	tauy=rhoair/rhowater*tauyair
689	! note: forward-backward coriolis force for effective time-centering
690	hu2=hu1+dt( fhv1 + taux)
691	hv2=hv1+dt(-fhu2 + tauy)
692	! A2=A1+q*dt
693	A2=A1
694
695	huml=hu2
696	hvml=hv2
697
698	hold=h
699	B2=hu2hu2+hv2hv2
700	hsqrd=-A2/Gam + sqrt(A2A2/(GamGam) + 2.*B2/BV2)
701	h=sqrt(max(hsqrd,0.0))
702	! limit to positive h change
703	if(h.lt.hold)h=hold
704
705	if(h.ne.0.)then
706	tml=max(t0ml - Gam(h-h0) + 0.5Gam*h + A2/h, 273.15)
707	u=hu2/h
708	v=hv2/h
709	else
710	tml=t0ml
711	u=0.
712	v=0.
713	endif
714	tsk=tml
715	! if(h.gt.100.)print *,i,j,h,tml,' h,tml'
716
717	! ww: output point data
718	! if( (i.eq.190 .and. j.eq.115) .or. (i.eq.170 .and. j.eq.125) ) then
719	! write(jtime,fmt='("TS ",f10.0)') float(itimestep)
720	! CALL wrf_message ( TRIM(jtime) )
721	! write(time_series,fmt='("OML",2I4,2F9.5,2F8.2,2E15.5,F8.3)') &
722	! i,j,u,v,tml,h,taux,tauy,a2
723	! CALL wrf_message ( TRIM(time_series) )
724	! end if
725
726	END SUBROUTINE OCEANML
727	!-------------------------------------------------------------------
728
729	END MODULE module_sf_slab

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format