Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_sfclay.F @ 3988

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

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

Download in other formats:

Original Format