Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_sfclay.F @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 46.0 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3	MODULE module_sf_sfclay
4
5	REAL , PARAMETER :: VCONVC=1.
6	REAL , PARAMETER :: CZO=0.0185
7	REAL , PARAMETER :: OZO=1.59E-5
8
9	REAL, DIMENSION(0:1000 ),SAVE :: PSIMTB,PSIHTB
10
11	CONTAINS
12
13	!-------------------------------------------------------------------
14	SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w, &
15	CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
16	ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
17	XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
18	U10,V10,TH2,T2,Q2, &
19	GZ1OZ0,WSPD,BR,ISFFLX,DX, &
20	SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
21	KARMAN,EOMEG,STBOLT, &
22	P1000mb, &
23	ids,ide, jds,jde, kds,kde, &
24	ims,ime, jms,jme, kms,kme, &
25	its,ite, jts,jte, kts,kte, &
26	ustm,ck,cka,cd,cda,isftcflx,iz0tlnd )
27	!-------------------------------------------------------------------
28	IMPLICIT NONE
29	!-------------------------------------------------------------------
30	!-- U3D 3D u-velocity interpolated to theta points (m/s)
31	!-- V3D 3D v-velocity interpolated to theta points (m/s)
32	!-- T3D temperature (K)
33	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
34	!-- P3D 3D pressure (Pa)
35	!-- dz8w dz between full levels (m)
36	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
37	!-- G acceleration due to gravity (m/s^2)
38	!-- ROVCP R/CP
39	!-- R gas constant for dry air (J/kg/K)
40	!-- XLV latent heat of vaporization for water (J/kg)
41	!-- PSFC surface pressure (Pa)
42	!-- ZNT roughness length (m)
43	!-- UST u* in similarity theory (m/s)
44	!-- USTM u* in similarity theory (m/s) without vconv correction
45	! used to couple with TKE scheme
46	!-- PBLH PBL height from previous time (m)
47	!-- MAVAIL surface moisture availability (between 0 and 1)
48	!-- ZOL z/L height over Monin-Obukhov length
49	!-- MOL T* (similarity theory) (K)
50	!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
51	!-- PSIM similarity stability function for momentum
52	!-- PSIH similarity stability function for heat
53	!-- XLAND land mask (1 for land, 2 for water)
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 net upward latent heat flux at surface (W/m^2)
57	!-- TSK surface temperature (K)
58	!-- FLHC exchange coefficient for heat (W/m^2/K)
59	!-- FLQC exchange coefficient for moisture (kg/m^2/s)
60	!-- CHS heat/moisture exchange coefficient for LSM (m/s)
61	!-- QGH lowest-level saturated mixing ratio
62	!-- QSFC ground saturated mixing ratio
63	!-- U10 diagnostic 10m u wind
64	!-- V10 diagnostic 10m v wind
65	!-- TH2 diagnostic 2m theta (K)
66	!-- T2 diagnostic 2m temperature (K)
67	!-- Q2 diagnostic 2m mixing ratio (kg/kg)
68	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
69	!-- WSPD wind speed at lowest model level (m/s)
70	!-- BR bulk Richardson number in surface layer
71	!-- ISFFLX isfflx=1 for surface heat and moisture fluxes
72	!-- DX horizontal grid size (m)
73	!-- SVP1 constant for saturation vapor pressure (kPa)
74	!-- SVP2 constant for saturation vapor pressure (dimensionless)
75	!-- SVP3 constant for saturation vapor pressure (K)
76	!-- SVPT0 constant for saturation vapor pressure (K)
77	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
78	!-- EP2 constant for specific humidity calculation
79	! (R_d/R_v) (dimensionless)
80	!-- KARMAN Von Karman constant
81	!-- EOMEG angular velocity of earth's rotation (rad/s)
82	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
83	!-- ck enthalpy exchange coeff at 10 meters
84	!-- cd momentum exchange coeff at 10 meters
85	!-- cka enthalpy exchange coeff at the lowest model level
86	!-- cda momentum exchange coeff at the lowest model level
87	!-- isftcflx =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd, =2 Garratt
88	!-- iz0tlnd =0 Carlson-Boland, =1 Czil_new
89	!-- ids start index for i in domain
90	!-- ide end index for i in domain
91	!-- jds start index for j in domain
92	!-- jde end index for j in domain
93	!-- kds start index for k in domain
94	!-- kde end index for k in domain
95	!-- ims start index for i in memory
96	!-- ime end index for i in memory
97	!-- jms start index for j in memory
98	!-- jme end index for j in memory
99	!-- kms start index for k in memory
100	!-- kme end index for k in memory
101	!-- its start index for i in tile
102	!-- ite end index for i in tile
103	!-- jts start index for j in tile
104	!-- jte end index for j in tile
105	!-- kts start index for k in tile
106	!-- kte end index for k in tile
107	!-------------------------------------------------------------------
108	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
109	ims,ime, jms,jme, kms,kme, &
110	its,ite, jts,jte, kts,kte
111	!
112	INTEGER, INTENT(IN ) :: ISFFLX
113	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
114	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
115	REAL, INTENT(IN ) :: P1000mb
116	!
117	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
118	INTENT(IN ) :: dz8w
119
120	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
121	INTENT(IN ) :: QV3D, &
122	P3D, &
123	T3D
124
125	REAL, DIMENSION( ims:ime, jms:jme ) , &
126	INTENT(IN ) :: MAVAIL, &
127	PBLH, &
128	XLAND, &
129	TSK
130	REAL, DIMENSION( ims:ime, jms:jme ) , &
131	INTENT(OUT ) :: U10, &
132	V10, &
133	TH2, &
134	T2, &
135	Q2, &
136	QSFC
137
138	!
139	REAL, DIMENSION( ims:ime, jms:jme ) , &
140	INTENT(INOUT) :: REGIME, &
141	HFX, &
142	QFX, &
143	LH, &
144	MOL,RMOL
145	!m the following 5 are change to memory size
146	!
147	REAL, DIMENSION( ims:ime, jms:jme ) , &
148	INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
149	PSIM,PSIH
150
151	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
152	INTENT(IN ) :: U3D, &
153	V3D
154
155	REAL, DIMENSION( ims:ime, jms:jme ) , &
156	INTENT(IN ) :: PSFC
157
158	REAL, DIMENSION( ims:ime, jms:jme ) , &
159	INTENT(INOUT) :: ZNT, &
160	ZOL, &
161	UST, &
162	CPM, &
163	CHS2, &
164	CQS2, &
165	CHS
166
167	REAL, DIMENSION( ims:ime, jms:jme ) , &
168	INTENT(INOUT) :: FLHC,FLQC
169
170	REAL, DIMENSION( ims:ime, jms:jme ) , &
171	INTENT(INOUT) :: &
172	QGH
173
174
175
176	REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
177
178	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , &
179	INTENT(OUT) :: ck,cka,cd,cda,ustm
180
181	INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
182
183	! LOCAL VARS
184
185	REAL, DIMENSION( its:ite ) :: U1D, &
186	V1D, &
187	QV1D, &
188	P1D, &
189	T1D
190
191	REAL, DIMENSION( its:ite ) :: dz8w1d
192
193	INTEGER :: I,J
194
195	DO J=jts,jte
196	DO i=its,ite
197	dz8w1d(I) = dz8w(i,1,j)
198	ENDDO
199
200	DO i=its,ite
201	U1D(i) =U3D(i,1,j)
202	V1D(i) =V3D(i,1,j)
203	QV1D(i)=QV3D(i,1,j)
204	P1D(i) =P3D(i,1,j)
205	T1D(i) =T3D(i,1,j)
206	ENDDO
207
208	! Sending array starting locations of optional variables may cause
209	! troubles, so we explicitly change the call.
210
211	CALL SFCLAY1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, &
212	CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
213	CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
214	ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
215	MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
216	XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
217	U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), &
218	Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j), &
219	QSFC(ims,j),LH(ims,j), &
220	GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, &
221	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
222	P1000mb, &
223	ids,ide, jds,jde, kds,kde, &
224	ims,ime, jms,jme, kms,kme, &
225	its,ite, jts,jte, kts,kte &
226	#if ( EM_CORE == 1 )
227	,isftcflx,iz0tlnd, &
228	USTM(ims,j),CK(ims,j),CKA(ims,j), &
229	CD(ims,j),CDA(ims,j) &
230	#endif
231	)
232	ENDDO
233
234
235	END SUBROUTINE SFCLAY
236
237
238	!-------------------------------------------------------------------
239	SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d, &
240	CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
241	ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
242	XLAND,HFX,QFX,TSK, &
243	U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH, &
244	QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
245	SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
246	KARMAN,EOMEG,STBOLT, &
247	P1000mb, &
248	ids,ide, jds,jde, kds,kde, &
249	ims,ime, jms,jme, kms,kme, &
250	its,ite, jts,jte, kts,kte, &
251	isftcflx, iz0tlnd, &
252	ustm,ck,cka,cd,cda )
253	!-------------------------------------------------------------------
254	IMPLICIT NONE
255	!-------------------------------------------------------------------
256	REAL, PARAMETER :: XKA=2.4E-5
257	REAL, PARAMETER :: PRT=1.
258
259	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
260	ims,ime, jms,jme, kms,kme, &
261	its,ite, jts,jte, kts,kte, &
262	J
263	!
264	INTEGER, INTENT(IN ) :: ISFFLX
265	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
266	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
267	REAL, INTENT(IN ) :: P1000mb
268
269	!
270	REAL, DIMENSION( ims:ime ) , &
271	INTENT(IN ) :: MAVAIL, &
272	PBLH, &
273	XLAND, &
274	TSK
275	!
276	REAL, DIMENSION( ims:ime ) , &
277	INTENT(IN ) :: PSFCPA
278
279	REAL, DIMENSION( ims:ime ) , &
280	INTENT(INOUT) :: REGIME, &
281	HFX, &
282	QFX, &
283	MOL,RMOL
284	!m the following 5 are changed to memory size---
285	!
286	REAL, DIMENSION( ims:ime ) , &
287	INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
288	PSIM,PSIH
289
290	REAL, DIMENSION( ims:ime ) , &
291	INTENT(INOUT) :: ZNT, &
292	ZOL, &
293	UST, &
294	CPM, &
295	CHS2, &
296	CQS2, &
297	CHS
298
299	REAL, DIMENSION( ims:ime ) , &
300	INTENT(INOUT) :: FLHC,FLQC
301
302	REAL, DIMENSION( ims:ime ) , &
303	INTENT(INOUT) :: &
304	QGH
305
306	REAL, DIMENSION( ims:ime ) , &
307	INTENT(OUT) :: U10,V10, &
308	TH2,T2,Q2,QSFC,LH
309
310
311	REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
312
313	! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
314	REAL, DIMENSION( its:ite ), INTENT(IN ) :: dz8w1d
315
316	REAL, DIMENSION( its:ite ), INTENT(IN ) :: UX, &
317	VX, &
318	QV1D, &
319	P1D, &
320	T1D
321
322	REAL, OPTIONAL, DIMENSION( ims:ime ) , &
323	INTENT(OUT) :: ck,cka,cd,cda,ustm
324
325	INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
326
327	! LOCAL VARS
328
329	REAL, DIMENSION( its:ite ) :: ZA, &
330	THVX,ZQKL, &
331	ZQKLP1, &
332	THX,QX, &
333	PSIH2, &
334	PSIM2, &
335	PSIH10, &
336	PSIM10, &
337	DENOMQ, &
338	DENOMQ2, &
339	DENOMT2, &
340	WSPDI, &
341	GZ2OZ0, &
342	GZ10OZ0
343	!
344	REAL, DIMENSION( its:ite ) :: &
345	RHOX,GOVRTH, &
346	TGDSA
347	!
348	REAL, DIMENSION( its:ite) :: SCR3,SCR4
349	REAL, DIMENSION( its:ite ) :: THGB, PSFC
350	!
351	INTEGER :: KL
352
353	INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
354
355	REAL :: PL,THCON,TVCON,E1
356	REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
357	REAL :: DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
358	REAL :: FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,RESTAR2
359	!-------------------------------------------------------------------
360	KL=kte
361
362	DO i=its,ite
363	! PSFC cb
364	PSFC(I)=PSFCPA(I)/1000.
365	ENDDO
366	!
367	!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
368	!
369	DO 5 I=its,ite
370	TGDSA(I)=TSK(I)
371	! PSFC cb
372	! THGB(I)=TSK(I)(100./PSFC(I))*ROVCP
373	THGB(I)=TSK(I)(P1000mb/PSFCPA(I))*ROVCP
374	5 CONTINUE
375	!
376	!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
377	! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
378	!
379	! * NOTE *
380	! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
381	! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
382	! TENDENCIES.
383	!
384	10 CONTINUE
385
386	! DO 24 I=its,ite
387	! UX(I)=U1D(I)
388	! VX(I)=V1D(I)
389	! 24 CONTINUE
390
391	26 CONTINUE
392
393	!.....SCR3(I,K) STORE TEMPERATURE,
394	! SCR4(I,K) STORE VIRTUAL TEMPERATURE.
395
396	DO 30 I=its,ite
397	! PL cb
398	PL=P1D(I)/1000.
399	SCR3(I)=T1D(I)
400	! THCON=(100./PL)**ROVCP
401	THCON=(P1000mb0.001/PL)*ROVCP
402	THX(I)=SCR3(I)*THCON
403	SCR4(I)=SCR3(I)
404	THVX(I)=THX(I)
405	QX(I)=0.
406	30 CONTINUE
407	!
408	DO I=its,ite
409	QGH(I)=0.
410	FLHC(I)=0.
411	FLQC(I)=0.
412	CPM(I)=CP
413	ENDDO
414	!
415	! IF(IDRY.EQ.1)GOTO 80
416	DO 50 I=its,ite
417	QX(I)=QV1D(I)
418	TVCON=(1.+EP1*QX(I))
419	THVX(I)=THX(I)*TVCON
420	SCR4(I)=SCR3(I)*TVCON
421	50 CONTINUE
422	!
423	DO 60 I=its,ite
424	E1=SVP1EXP(SVP2(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
425	! for land points QSFC can come from previous time step
426	if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1)
427	! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
428	! Q2SAT = QGH IN LSM
429	E1=SVP1EXP(SVP2(T1D(I)-SVPT0)/(T1D(I)-SVP3))
430	PL=P1D(I)/1000.
431	QGH(I)=EP2*E1/(PL-E1)
432	CPM(I)=CP(1.+0.8QX(I))
433	60 CONTINUE
434	80 CONTINUE
435
436	!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
437	! LEVEL, AND THE LAYER THICKNESSES.
438
439	DO 90 I=its,ite
440	ZQKLP1(I)=0.
441	RHOX(I)=PSFC(I)1000./(RSCR4(I))
442	90 CONTINUE
443	!
444	DO 110 I=its,ite
445	ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
446	110 CONTINUE
447	!
448	DO 120 I=its,ite
449	ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
450	120 CONTINUE
451	!
452	DO 160 I=its,ite
453	GOVRTH(I)=G/THX(I)
454	160 CONTINUE
455
456	!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
457	! AKB(1976), EQ(12).
458
459	DO 260 I=its,ite
460	GZ1OZ0(I)=ALOG(ZA(I)/ZNT(I))
461	GZ2OZ0(I)=ALOG(2./ZNT(I))
462	GZ10OZ0(I)=ALOG(10./ZNT(I))
463	IF((XLAND(I)-1.5).GE.0)THEN
464	ZL=ZNT(I)
465	ELSE
466	ZL=0.01
467	ENDIF
468	WSPD(I)=SQRT(UX(I)UX(I)+VX(I)VX(I))
469
470	TSKV=THGB(I)(1.+EP1QSFC(I))
471	DTHVDZ=(THVX(I)-TSKV)
472	! Convective velocity scale Vc and subgrid-scale velocity Vsg
473	! following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
474	! ... HONG Aug. 2001
475	!
476	! VCONV = 0.25sqrt(g/tskvpblh(i)*dthvm)
477	! Use Beljaars over land, old MM5 (Wyngaard) formula over water
478	if (xland(i).lt.1.5) then
479	fluxc = max(hfx(i)/rhox(i)/cp &
480	+ ep1tskvqfx(i)/rhox(i),0.)
481	VCONV = vconvc(g/tgdsa(i)pblh(i)fluxc)*.33
482	else
483	IF(-DTHVDZ.GE.0)THEN
484	DTHVM=-DTHVDZ
485	ELSE
486	DTHVM=0.
487	ENDIF
488	VCONV = 2.*SQRT(DTHVM)
489	endif
490	! Mahrt and Sun low-res correction
491	VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
492	WSPD(I)=SQRT(WSPD(I)WSPD(I)+VCONVVCONV+vsgd*vsgd)
493	WSPD(I)=AMAX1(WSPD(I),0.1)
494	BR(I)=GOVRTH(I)ZA(I)DTHVDZ/(WSPD(I)*WSPD(I))
495	! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
496	IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
497	!jdf
498	RMOL(I)=-GOVRTH(I)DTHVDZZA(I)*KARMAN
499	!jdf
500
501	260 CONTINUE
502
503	!
504	!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
505	!
506	!
507	! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
508	! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
509	!
510	! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
511	!
512	! 1. BR .GE. 0.2;
513	! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
514	!
515	! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
516	! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
517	! (REGIME=2),
518	!
519	! 3. BR .EQ. 0.0
520	! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
521	!
522	! 4. BR .LT. 0.0
523	! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
524	!
525	!CCCCC
526
527	DO 320 I=its,ite
528	!CCCCC
529	!CC REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
530	!CC IF(BR(I).LT.0..AND.HOL(I,J).GT.1.5)GOTO 310
531	IF(BR(I).LT.0.)GOTO 310
532	!
533	!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
534	!
535	IF(BR(I).LT.0.2)GOTO 270
536	REGIME(I)=1.
537	PSIM(I)=-10.*GZ1OZ0(I)
538	! LOWER LIMIT ON PSI IN STABLE CONDITIONS
539	PSIM(I)=AMAX1(PSIM(I),-10.)
540	PSIH(I)=PSIM(I)
541	PSIM10(I)=10./ZA(I)*PSIM(I)
542	PSIM10(I)=AMAX1(PSIM10(I),-10.)
543	PSIH10(I)=PSIM10(I)
544	PSIM2(I)=2./ZA(I)*PSIM(I)
545	PSIM2(I)=AMAX1(PSIM2(I),-10.)
546	PSIH2(I)=PSIM2(I)
547
548	! 1.0 over Monin-Obukhov length
549	IF(UST(I).LT.0.01)THEN
550	RMOL(I)=BR(I)*GZ1OZ0(I) !ZA/L
551	ELSE
552	RMOL(I)=KARMANGOVRTH(I)ZA(I)MOL(I)/(UST(I)UST(I)) !ZA/L
553	ENDIF
554	RMOL(I)=AMIN1(RMOL(I),9.999) ! ZA/L
555	RMOL(I) = RMOL(I)/ZA(I) !1.0/L
556
557	GOTO 320
558	!
559	!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
560	!
561	270 IF(BR(I).EQ.0.0)GOTO 280
562	REGIME(I)=2.
563	PSIM(I)=-5.0BR(I)GZ1OZ0(I)/(1.1-5.0*BR(I))
564	! LOWER LIMIT ON PSI IN STABLE CONDITIONS
565	PSIM(I)=AMAX1(PSIM(I),-10.)
566	!.....AKB(1976), EQ(16).
567	PSIH(I)=PSIM(I)
568	PSIM10(I)=10./ZA(I)*PSIM(I)
569	PSIM10(I)=AMAX1(PSIM10(I),-10.)
570	PSIH10(I)=PSIM10(I)
571	PSIM2(I)=2./ZA(I)*PSIM(I)
572	PSIM2(I)=AMAX1(PSIM2(I),-10.)
573	PSIH2(I)=PSIM2(I)
574
575	! Linear form: PSIM = -0.5*ZA/L; e.g, see eqn 16 of
576	! Blackadar, Modeling the nocturnal boundary layer, Preprints,
577	! Third Symposium on Atmospheric Turbulence Diffusion and Air Quality,
578	! Raleigh, NC, 1976
579	ZOL(I) = BR(I)GZ1OZ0(I)/(1.00001-5.0BR(I))
580
581	if ( ZOL(I) .GT. 0.5 ) then ! linear form ok
582	! Holtslag and de Bruin, J. App. Meteor 27, 689-704, 1988;
583	! see also, Launiainen, Boundary-Layer Meteor 76,165-179, 1995
584	! Eqn (8) of Launiainen, 1995
585	ZOL(I) = ( 1.89GZ1OZ0(I) + 44.2 ) BR(I)*BR(I) &
586	+ ( 1.18GZ1OZ0(I) - 1.37 ) BR(I)
587	ZOL(I)=AMIN1(ZOL(I),9.999)
588	end if
589
590	! 1.0 over Monin-Obukhov length
591	RMOL(I)= ZOL(I)/ZA(I)
592
593	GOTO 320
594	!
595	!-----CLASS 3; FORCED CONVECTION:
596	!
597	280 REGIME(I)=3.
598	PSIM(I)=0.0
599	PSIH(I)=PSIM(I)
600	PSIM10(I)=0.
601	PSIH10(I)=PSIM10(I)
602	PSIM2(I)=0.
603	PSIH2(I)=PSIM2(I)
604
605
606	IF(UST(I).LT.0.01)THEN
607	ZOL(I)=BR(I)*GZ1OZ0(I)
608	ELSE
609	ZOL(I)=KARMANGOVRTH(I)ZA(I)MOL(I)/(UST(I)UST(I))
610	ENDIF
611
612	RMOL(I) = ZOL(I)/ZA(I)
613
614	GOTO 320
615	!
616	!-----CLASS 4; FREE CONVECTION:
617	!
618	310 CONTINUE
619	REGIME(I)=4.
620	IF(UST(I).LT.0.01)THEN
621	ZOL(I)=BR(I)*GZ1OZ0(I)
622	ELSE
623	ZOL(I)=KARMANGOVRTH(I)ZA(I)MOL(I)/(UST(I)UST(I))
624	ENDIF
625	ZOL10=10./ZA(I)*ZOL(I)
626	ZOL2=2./ZA(I)*ZOL(I)
627	ZOL(I)=AMIN1(ZOL(I),0.)
628	ZOL(I)=AMAX1(ZOL(I),-9.9999)
629	ZOL10=AMIN1(ZOL10,0.)
630	ZOL10=AMAX1(ZOL10,-9.9999)
631	ZOL2=AMIN1(ZOL2,0.)
632	ZOL2=AMAX1(ZOL2,-9.9999)
633	NZOL=INT(-ZOL(I)*100.)
634	RZOL=-ZOL(I)*100.-NZOL
635	NZOL10=INT(-ZOL10*100.)
636	RZOL10=-ZOL10*100.-NZOL10
637	NZOL2=INT(-ZOL2*100.)
638	RZOL2=-ZOL2*100.-NZOL2
639	PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
640	PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
641	PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))
642	PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
643	PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))
644	PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))
645
646	!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
647	!--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
648	! PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
649	! PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
650	PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
651	PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
652	PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
653	PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
654	! AHW: mods to compute ck, cd
655	PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
656
657	RMOL(I) = ZOL(I)/ZA(I)
658
659	320 CONTINUE
660	!
661	!-----COMPUTE THE FRICTIONAL VELOCITY:
662	! ZA(1982) EQS(2.60),(2.61).
663	!
664	DO 330 I=its,ite
665	DTG=THX(I)-THGB(I)
666	PSIX=GZ1OZ0(I)-PSIM(I)
667	PSIX10=GZ10OZ0(I)-PSIM10(I)
668	! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
669	! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
670	PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
671
672	IF((XLAND(I)-1.5).GE.0)THEN
673	ZL=ZNT(I)
674	ELSE
675	ZL=0.01
676	ENDIF
677	PSIQ=ALOG(KARMANUST(I)ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)
678	PSIT2=GZ2OZ0(I)-PSIH2(I)
679	PSIQ2=ALOG(KARMANUST(I)2./XKA+2./ZL)-PSIH2(I)
680	! AHW: mods to compute ck, cd
681	PSIQ10=ALOG(KARMANUST(I)10./XKA+10./ZL)-PSIH10(I)
682	IF ( PRESENT(ISFTCFLX) ) THEN
683	IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
684	! v3.1
685	! Z0Q = 1.e-4 + 1.e-3(MAX(0.,UST(I)-1.))*2
686	! hfip1
687	! Z0Q = 0.622.0E-5/UST(I) + 1.E-3(MAX(0.,UST(I)-1.5))**2
688	! v3.2
689	Z0Q = 1.e-4
690	PSIQ=ALOG(ZA(I)/Z0Q)-PSIH(I)
691	PSIT=PSIQ
692	PSIQ2=ALOG(2./Z0Q)-PSIH2(I)
693	PSIQ10=ALOG(10./Z0Q)-PSIH10(I)
694	PSIT2=PSIQ2
695	ENDIF
696	IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN
697	! AHW: Garratt formula: Calculate roughness Reynolds number
698	! Kinematic viscosity of air (linear approc to
699	! temp dependence at sea levle)
700	VISC=(1.32+0.009(SCR3(I)-273.15))1.E-5
701	!! VISC=1.5E-5
702	RESTAR=UST(I)*ZNT(I)/VISC
703	RESTAR2=2.48*SQRT(SQRT(RESTAR))-2.
704	PSIT=GZ1OZ0(I)-PSIH(I)+RESTAR2
705	PSIQ=GZ1OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
706	PSIT2=GZ2OZ0(I)-PSIH2(I)+RESTAR2
707	PSIQ2=GZ2OZ0(I)-PSIH2(I)+2.28*SQRT(SQRT(RESTAR))-2.
708	PSIQ10=GZ10OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
709	ENDIF
710	ENDIF
711	IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
712	Ck(I)=(karman/psix10)*(karman/psiq10)
713	Cd(I)=(karman/psix10)*(karman/psix10)
714	Cka(I)=(karman/psix)*(karman/psiq)
715	Cda(I)=(karman/psix)*(karman/psix)
716	ENDIF
717	IF ( PRESENT(IZ0TLND) ) THEN
718	IF ( IZ0TLND.EQ.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN
719	ZL=ZNT(I)
720	! CZIL RELATED CHANGES FOR LAND
721	VISC=(1.32+0.009(SCR3(I)-273.15))1.E-5
722	RESTAR=UST(I)*ZL/VISC
723	! Modify CZIL according to Chen & Zhang, 2009
724
725	CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) )
726
727	PSIT=GZ1OZ0(I)-PSIH(I)+CZILKARMANSQRT(RESTAR)
728	PSIQ=GZ1OZ0(I)-PSIH(I)+CZILKARMANSQRT(RESTAR)
729	PSIT2=GZ2OZ0(I)-PSIH2(I)+CZILKARMANSQRT(RESTAR)
730	PSIQ2=GZ2OZ0(I)-PSIH2(I)+CZILKARMANSQRT(RESTAR)
731
732	ENDIF
733	ENDIF
734	! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
735	UST(I)=0.5UST(I)+0.5KARMAN*WSPD(I)/PSIX
736	! TKE coupling: compute ust without vconv for use in tke scheme
737	WSPDI(I)=SQRT(UX(I)UX(I)+VX(I)VX(I))
738	IF ( PRESENT(USTM) ) THEN
739	USTM(I)=0.5USTM(I)+0.5KARMAN*WSPDI(I)/PSIX
740	ENDIF
741	U10(I)=UX(I)*PSIX10/PSIX
742	V10(I)=VX(I)*PSIX10/PSIX
743	TH2(I)=THGB(I)+DTG*PSIT2/PSIT
744	Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
745	! T2(I) = TH2(I)(PSFC(I)/100.)*ROVCP
746	T2(I) = TH2(I)(PSFCPA(I)/P1000mb)*ROVCP
747	! LATER Q2 WILL BE OVERWRITTEN FOR LAND POINTS IN SURFCE
748	! QA2(I,J) = Q2(I)
749	! UA10(I,J) = U10(I)
750	! VA10(I,J) = V10(I)
751	! write(,1002)UST(I),KARMANWSPD(I),PSIX,KARMAN*WSPD(I)/PSIX
752	!
753	IF((XLAND(I)-1.5).LT.0.)THEN
754	UST(I)=AMAX1(UST(I),0.1)
755	ENDIF
756	MOL(I)=KARMAN*DTG/PSIT/PRT
757	DENOMQ(I)=PSIQ
758	DENOMQ2(I)=PSIQ2
759	DENOMT2(I)=PSIT2
760	330 CONTINUE
761	!
762	335 CONTINUE
763
764	!-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:
765
766	DO i=its,ite
767	QFX(i)=0.
768	HFX(i)=0.
769	ENDDO
770
771	IF (ISFFLX.EQ.0) GOTO 410
772
773	!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
774
775	DO 360 I=its,ite
776	IF((XLAND(I)-1.5).GE.0)THEN
777	ZNT(I)=CZOUST(I)UST(I)/G+OZO
778	! AHW: change roughness length, and hence the drag coefficients Ck and Cd
779	IF ( PRESENT(ISFTCFLX) ) THEN
780	IF ( ISFTCFLX.NE.0 ) THEN
781	! ZNT(I)=10.exp(-9.UST(I)**(-.3333))
782	ZNT(I)=10.exp(-9.5UST(I)**(-.3333))
783	ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
784	ZNT(I)=MIN(ZNT(I),2.85e-3)
785	ZNT(I)=MAX(ZNT(I),1.27e-7)
786	ENDIF
787	ENDIF
788	ZL = ZNT(I)
789	ELSE
790	ZL = 0.01
791	ENDIF
792	FLQC(I)=RHOX(I)MAVAIL(I)UST(I)*KARMAN/DENOMQ(I)
793	! FLQC(I)=RHOX(I)MAVAIL(I)UST(I)*KARMAN/( &
794	! ALOG(KARMANUST(I)ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
795	DTTHX=ABS(THX(I)-THGB(I))
796	IF(DTTHX.GT.1.E-5)THEN
797	FLHC(I)=CPM(I)RHOX(I)UST(I)*MOL(I)/(THX(I)-THGB(I))
798	! write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
799	1001 format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
800	ELSE
801	FLHC(I)=0.
802	ENDIF
803	360 CONTINUE
804
805	!
806	!-----COMPUTE SURFACE MOIST FLUX:
807	!
808	! IF(IDRY.EQ.1)GOTO 390
809	!
810	DO 370 I=its,ite
811	QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
812	QFX(I)=AMAX1(QFX(I),0.)
813	LH(I)=XLV*QFX(I)
814	370 CONTINUE
815
816	!-----COMPUTE SURFACE HEAT FLUX:
817	!
818	390 CONTINUE
819	DO 400 I=its,ite
820	IF(XLAND(I)-1.5.GT.0.)THEN
821	HFX(I)=FLHC(I)*(THGB(I)-THX(I))
822	IF ( PRESENT(ISFTCFLX) ) THEN
823	IF ( ISFTCFLX.NE.0 ) THEN
824	! AHW: add dissipative heating term
825	HFX(I)=HFX(I)+RHOX(I)USTM(I)USTM(I)*WSPDI(I)
826	ENDIF
827	ENDIF
828	ELSEIF(XLAND(I)-1.5.LT.0.)THEN
829	HFX(I)=FLHC(I)*(THGB(I)-THX(I))
830	HFX(I)=AMAX1(HFX(I),-250.)
831	ENDIF
832	400 CONTINUE
833
834	DO I=its,ite
835	IF((XLAND(I)-1.5).GE.0)THEN
836	ZL=ZNT(I)
837	ELSE
838	ZL=0.01
839	ENDIF
840	CHS(I)=UST(I)*KARMAN/DENOMQ(I)
841	! GZ2OZ0(I)=ALOG(2./ZNT(I))
842	! PSIM2(I)=-10.*GZ2OZ0(I)
843	! PSIM2(I)=AMAX1(PSIM2(I),-10.)
844	! PSIH2(I)=PSIM2(I)
845	CQS2(I)=UST(I)*KARMAN/DENOMQ2(I)
846	CHS2(I)=UST(I)*KARMAN/DENOMT2(I)
847	ENDDO
848
849	410 CONTINUE
850	!jdf
851	! DO I=its,ite
852	! IF(UST(I).GE.0.1) THEN
853	! RMOL(I)=RMOL(I)(-FLHC(I))/(UST(I)UST(I)*UST(I))
854	! ELSE
855	! RMOL(I)=RMOL(I)(-FLHC(I))/(0.10.1*0.1)
856	! ENDIF
857	! ENDDO
858	!jdf
859
860	!
861	END SUBROUTINE SFCLAY1D
862
863	!====================================================================
864	SUBROUTINE sfclayinit( allowed_to_read )
865
866	LOGICAL , INTENT(IN) :: allowed_to_read
867	INTEGER :: N
868	REAL :: ZOLN,X,Y
869
870	DO N=0,1000
871	ZOLN=-FLOAT(N)*0.01
872	X=(1-16.ZOLN)*0.25
873	PSIMTB(N)=2ALOG(0.5(1+X))+ALOG(0.5(1+XX))- &
874	2.ATAN(X)+2.ATAN(1.)
875	Y=(1-16ZOLN)*0.5
876	PSIHTB(N)=2ALOG(0.5(1+Y))
877	ENDDO
878
879	END SUBROUTINE sfclayinit
880
881	!-------------------------------------------------------------------
882
883	END MODULE module_sf_sfclay

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

Download in other formats:

Original Format