Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_mynn.F @ 1

Last change on this file since 1 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: 43.2 KB

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

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

Download in other formats:

Original Format