Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_gfdl.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:

Property svn:executable set to *

File size: 66.5 KB

Line
1	!-- XLV latent heat of vaporization for water (J/kg)
2	!
3	MODULE module_sf_gfdl
4
5	!real, dimension(-100:2000,-100:2000), save :: z00
6
7	CONTAINS
8
9	!-------------------------------------------------------------------
10	SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, &
11	CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, &
12	DT, SMOIS,num_soil_layers,ISLTYP,ZNT,UST,PSIM,PSIH, &
13	XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, & ! gopal's doing for Ocean coupling
14	QGH,QSFC,U10,V10, &
15	GZ1OZ0,WSPD,BR,ISFFLX, &
16	EP1,EP2,KARMAN,NTSFLG,SFENTH, &
17	ids,ide, jds,jde, kds,kde, &
18	ims,ime, jms,jme, kms,kme, &
19	its,ite, jts,jte, kts,kte )
20	!-------------------------------------------------------------------
21	USE MODULE_GFS_MACHINE, ONLY : kind_phys
22	USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
23	USE MODULE_GFS_PHYSCONS, grav => con_g
24	!-------------------------------------------------------------------
25	IMPLICIT NONE
26	!-------------------------------------------------------------------
27	!-- U3D 3D u-velocity interpolated to theta points (m/s)
28	!-- V3D 3D v-velocity interpolated to theta points (m/s)
29	!-- T3D temperature (K)
30	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
31	!-- P3D 3D pressure (Pa)
32	!-- DT time step (second)
33	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
34	!-- ROVCP R/CP
35	!-- R gas constant for dry air (J/kg/K)
36	!-- XLV latent heat of vaporization for water (J/kg)
37	!-- PSFC surface pressure (Pa)
38	!-- ZNT roughness length (m)
39	!-- MAVAIL surface moisture availability (between 0 and 1)
40	!-- UST u* in similarity theory (m/s)
41	!-- PSIM similarity stability function for momentum
42	!-- PSIH similarity stability function for heat
43	!-- XLAND land mask (1 for land, 2 for water)
44	!-- HFX upward heat flux at the surface (W/m^2)
45	!-- QFX upward moisture flux at the surface (kg/m^2/s)
46	!-- TAUX RHOU*2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
47	!-- TAUY RHOU*2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
48	!-- LH net upward latent heat flux at surface (W/m^2)
49	!-- GSW downward short wave flux at ground surface (W/m^2)
50	!-- GLW downward long wave flux at ground surface (W/m^2)
51	!-- TSK surface temperature (K)
52	!-- FLHC exchange coefficient for heat (m/s)
53	!-- FLQC exchange coefficient for moisture (m/s)
54	!-- QGH lowest-level saturated mixing ratio
55	!-- U10 diagnostic 10m u wind
56	!-- V10 diagnostic 10m v wind
57	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
58	!-- WSPD wind speed at lowest model level (m/s)
59	!-- BR bulk Richardson number in surface layer
60	!-- ISFFLX isfflx=1 for surface heat and moisture fluxes
61	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
62	!-- KARMAN Von Karman constant
63	!-- SFENTH enthalpy flux factor 0 zot via charnock ..>0 zot enhanced>15m/s
64	!-- ids start index for i in domain
65	!-- ide end index for i in domain
66	!-- jds start index for j in domain
67	!-- jde end index for j in domain
68	!-- kds start index for k in domain
69	!-- kde end index for k in domain
70	!-- ims start index for i in memory
71	!-- ime end index for i in memory
72	!-- jms start index for j in memory
73	!-- jme end index for j in memory
74	!-- kms start index for k in memory
75	!-- kme end index for k in memory
76	!-- its start index for i in tile
77	!-- ite end index for i in tile
78	!-- jts start index for j in tile
79	!-- jte end index for j in tile
80	!-- kts start index for k in tile
81	!-- kte end index for k in tile
82	!-------------------------------------------------------------------
83
84	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
85	ims,ime, jms,jme, kms,kme, &
86	its,ite, jts,jte, kts,kte, &
87	ISFFLX,NUM_SOIL_LAYERS,NTSFLG
88
89	REAL, INTENT(IN) :: &
90	CP, &
91	EP1, &
92	EP2, &
93	KARMAN, &
94	R, &
95	ROVCP, &
96	DT, &
97	SFENTH, &
98	XLV
99
100	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
101	P3D, &
102	QV3D, &
103	T3D, &
104	U3D, &
105	V3D
106	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP
107	REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS
108
109	REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
110	PSFC, &
111	GLW, &
112	GSW, &
113	XLAND
114	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
115	TSK, &
116	BR, &
117	CHS, &
118	CHS2, &
119	CPM, &
120	CQS2, &
121	FLHC, &
122	FLQC, &
123	GZ1OZ0, &
124	HFX, &
125	LH, &
126	PSIM, &
127	PSIH, &
128	QFX, &
129	QGH, &
130	QSFC, &
131	UST, &
132	ZNT, &
133	WSPD, &
134	TAUX, & ! gopal's doing for Ocean coupling
135	TAUY
136
137	REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
138	U10, &
139	V10
140
141
142	!--------------------------- LOCAL VARS ------------------------------
143
144	REAL :: ESAT, &
145	cpcgs, &
146	smc, &
147	smcdry, &
148	smcmax
149
150	REAL (kind=kind_phys) :: &
151	RHOX
152	REAL, DIMENSION(1:30) :: MAXSMC, &
153	DRYSMC
154	REAL (kind=kind_phys), DIMENSION(its:ite) :: &
155	CH, &
156	CM, &
157	DDVEL, &
158	DRAIN, &
159	EP, &
160	EVAP, &
161	FH, &
162	FH2, &
163	FM, &
164	HFLX, &
165	PH, &
166	PM, &
167	PRSL1, &
168	PRSLKI, &
169	PS, &
170	Q1, &
171	Q2M, &
172	QSS, &
173	QSURF, &
174	RB, &
175	RCL, &
176	RHO1, &
177	SLIMSK, &
178	STRESS, &
179	T1, &
180	T2M, &
181	THGB, &
182	THX, &
183	TSKIN, &
184	SHELEG, &
185	U1, &
186	U10M, &
187	USTAR, &
188	V1, &
189	V10M, &
190	WIND, &
191	Z0RL, &
192	Z1
193
194	REAL, DIMENSION(kms:kme, ims:ime) :: &
195	rpc, &
196	tpc, &
197	upc, &
198	vpc
199
200	REAL, DIMENSION(ims:ime) :: &
201	pspc, &
202	pkmax, &
203	tstrc, &
204	zoc, &
205	wetc, &
206	slwdc, &
207	rib, &
208	zkmax, &
209	tkmax, &
210	fxmx, &
211	fxmy, &
212	cdm, &
213	fxh, &
214	fxe, &
215	xxfh, &
216	xxfh2, &
217	wind10, &
218	tjloc
219
220	INTEGER :: &
221	I, &
222	II, &
223	IGPVS, &
224	IM, &
225	J, &
226	K, &
227	KM
228
229
230	DATA MAXSMC/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
231	0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
232	0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
233	0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
234	0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
235
236	DATA DRYSMC/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
237	0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
238	0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
239	0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
240	0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
241	DATA IGPVS/0/
242	save igpvs
243
244
245	if(igpvs.eq.0) then
246	! call readzo(glat,glon,6,ims,ime,jms,jme,its,ite,jts,jte,z00)
247	endif
248	igpvs=1
249
250	IM=ITE-ITS+1
251	KM=KTE-KTS+1
252
253	WRITE(0,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB 2010 UPGRADS',NTSFLG
254
255	DO J=jts,jte
256
257	DO i=its,ite
258	DDVEL(I)=0.
259	RCL(i)=1.
260	PRSL1(i)=P3D(i,kts,j)*.001
261	wetc(i)=1.0
262	if(xland(i,j).lt.1.99) then
263	smc=smois(i,1,j)
264	smcdry=drysmc(isltyp(i,j))
265	smcmax=maxsmc(isltyp(i,j))
266	wetc(i)=(smc-smcdry)/(smcmax-smcdry)
267	wetc(i)=amin1(1.,amax1(wetc(i),0.))
268	endif
269	! convert from Pa to cgs...
270	pspc(i)=PSFC(i,j)*10.
271	pkmax(i)=P3D(i,kts,j)*10.
272	PS(i)=PSFC(i,j)*.001
273	Q1(I) = QV3D(i,kts,j)
274	rpc(kts,i)=QV3D(i,kts,j)
275	! QSURF(I)=QSFC(I,J)
276	QSURF(I)=0.
277	SHELEG(I)=0.
278	SLIMSK(i)=ABS(XLAND(i,j)-2.)
279	TSKIN(i)=TSK(i,j)
280	tstrc(i)=TSK(i,j)
281	T1(I) = T3D(i,kts,j)
282	tpc(kts,i)=T3D(i,kts,j)
283	U1(I) = U3D(i,kts,j)
284	upc(kts,i)=U3D(i,kts,j) * 100.
285	USTAR(I) = UST(i,j)
286	V1(I) = V3D(i,kts,j)
287	vpc(kts,i)=v3D(i,kts,j) * 100.
288	Z0RL(I) = ZNT(i,j)*100.
289	zoc(i)=ZNT(i,j)*100.
290	if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i)
291	! Z0RL(I) = z00(i,j)*100.
292	! slwdc... GFDL downward net flux in units of cal/(cm**2/min)
293	! also divide by 104 to convert from /m2 to /cm**2
294	slwdc(i)=gsw(i,j)+glw(i,j)
295	slwdc(i)=0.23960.slwdc(i)*1.e-4
296	tjloc(i)=float(j)
297
298	ENDDO
299
300	DO i=its,ite
301	PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
302	THGB(I)=TSKIN(i)(100./PS(I))*ROVCP
303	THX(I)=T1(i)(100./PRSL1(I))*ROVCP
304	RHO1(I)=PRSL1(I)1000./(RT1(I)(1.+EP1Q1(I)))
305	Q1(I)=Q1(I)/(1.+Q1(I))
306	ENDDO
307
308	! if(j==2)then
309	! write(0,*)'--------------------------------------------'
310	! write(0,*) 'u, v, t, r, pkmax, pspc,wetc, tjloc,zoc,tstr'
311	! write(0,*)'--------------------------------------------'
312	! endif
313
314	! do i = its,ite
315	! WRITE(0,1010)i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), &
316	! pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i)
317	! enddo
318
319	CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,tstrc, &
320	pspc,pkmax,wetc,slwdc,tjloc, &
321	upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, &
322	ids,ide, jds,jde, kds,kde, &
323	ims,ime, jms,jme, kms,kme, &
324	its,ite, jts,jte, kts,kte )
325
326	! if(j==2)then
327	! write(0,*)'--------------------------------------------'
328	! write(0,*) 'fxh, fxe, fxmx, fxmy, cdm, xxfh zoc,tstrc'
329	! write(0,*)'--------------------------------------------'
330	! endif
331
332	! do i = its,ite
333	! WRITE(0,1010)i,j,fxh(i),fxe(i),fxmx(i),fxmy(i),cdm(i),rib(i),xxfh(i),zoc(i),tstrc(i)
334	! enddo
335
336	1010 format(2I4,9F11.6)
337
338
339
340	!GFDL CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
341	!GFDL SHELEG,TSKIN,QSURF, &
342	!WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, &
343	!WRF SLRAD,SNOWMT,DELT, &
344	!GFDL Z0RL, &
345	!WRF TG3,GFLUX,F10M, &
346	!GFDL U10M,V10M,T2M,Q2M, &
347	!WRF ZSOIL, &
348	!GFDL CM,CH,RB, &
349	!WRF RHSCNPY,RHSMC,AIM,BIM,CIM, &
350	!GFDL RCL,PRSL1,PRSLKI,SLIMSK, &
351	!GFDL DRAIN,EVAP,HFLX,STRESS,EP, &
352	!GFDL FM,FH,USTAR,WIND,DDVEL, &
353	!GFDL PM,PH,FH2,QSS,Z1 )
354
355
356	DO i=its,ite
357	! update skin temp only when using GFDL slab...
358
359	IF(NTSFLG==1) then
360	tsk(i,j) = tstrc(i) ! gopal's doing
361	!bugfix 4
362	! bob's doing patch tsk with neigboring values... are grid boundaries
363	if(j.eq.jde) then
364	tsk(i,j) = tsk(i,j-1)
365	endif
366
367	if(j.eq.jds) then
368	tsk(i,j) = tsk(i,j+1)
369	endif
370
371	if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
372	if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
373	endif
374
375
376	znt(i,j)= 0.01*abs(zoc(i))
377	wspd(i,j) = SQRT(upc(kts,i)upc(kts,i) + vpc(kts,i)vpc(kts,i))
378	wspd(i,j) = amax1(wspd(i,j) ,100.)/100.
379	u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100.
380	v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100.
381	! br =0.0001zfull(i,kmax)dthv/
382	! & (gmkstheta(i,kmax)wspd *wspd )
383	! zkmax = rgastpc(kmax,i)qqlog(kmax)*og
384	zkmax(i) = -Rtpc(kts,i)alog(pkmax(i)/pspc(i))/grav
385	!------------------------------------------------------------------------
386
387	gz1oz0(i,j)=alog(zkmax(i)/znt(i,j))
388	ustar (i)= 0.01sqrt(cdm(i) &
389	(upc(kts,i)upc(kts,i) + vpc(kts,i)vpc(kts,i)))
390	! convert from g/(cmcmsec) to kg/(mmsec)
391	qfx (i,j)=-10.fxe(i) ! BOB: qfx (i,1)=-10.fxe(i)
392	! cpcgs = 1.00464e7
393	! convert from ergs/gram/K to J/kg/K cpmks=1004
394	! hfx (i,1)=-0.001cpcgsfxh(i)
395	hfx (i,j)= -10.CPfxh(i) ! Bob: hfx (i,1)= -10.CPfxh(i)
396	taux (i,j)= fxmx(i)/10. ! gopal's doing for Ocean coupling
397	tauy (i,j)= fxmy(i)/10. ! gopal's doing for Ocean coupling
398	fm(i) = karman/sqrt(cdm(i))
399	fh(i) = karman*xxfh(i)
400	PSIM(i,j)=GZ1OZ0(i,j)-FM(i)
401	PSIH(i,j)=GZ1OZ0(i,j)-FH(i)
402	fh2(i) = karman*xxfh2(i)
403	ch(i) = karmankarman/(fm(i) fh(i))
404	cm(i) = cdm(i)
405
406	U10(i,j)=U10M(i)
407	V10(i,j)=V10M(i)
408	BR(i,j)=rib(i)
409	CHS(I,J)=CH(I)*wspd (i,j)
410	CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
411	CPM(I,J)=CP(1.+0.8QV3D(i,kts,j))
412	esat = fpvs(t1(i))
413	QGH(I,J)=ep2esat/(1000.ps(i)-esat)
414	esat = fpvs(tskin(i))
415	qss(i) = ep2esat/(1000.ps(i)-esat)
416	QSFC(I,J)=qss(i)
417	! PSIH(i,j)=PH(i)
418	! PSIM(i,j)=PM(i)
419	UST(i,j)=ustar(i)
420	! wspd (i,j) = SQRT(upc(kts,i)upc(kts,i) + vpc(kts,i)vpc(kts,i))
421	! wspd (i,j) = amax1(wspd (i,j) ,100.)/100.
422	! WSPD(i,j)=WIND(i)
423	! ZNT(i,j)=Z0RL(i)*.01
424	ENDDO
425
426	! write(0,*)'fm,fh,cm,ch(125)', fm(125),fh(125),cm(125),ch(125)
427
428	DO i=its,ite
429	FLHC(i,j)=CPM(I,J)RHO1(I)CHS(I,J)
430	FLQC(i,j)=RHO1(I)*CHS(I,J)
431	! GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
432	CQS2(i,j)=CHS2(I,J)
433	ENDDO
434
435	IF (ISFFLX.EQ.0) THEN
436	DO i=its,ite
437	HFX(i,j)=0.
438	LH(i,j)=0.
439	QFX(i,j)=0.
440	ENDDO
441	ELSE
442	DO i=its,ite
443	IF(XLAND(I,J)-1.5.GT.0.)THEN
444	! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
445	! cpcgs = 1.00464e7
446	! convert from ergs/gram/K to J/kg/K cpmks=1004
447	! hfx (i,j)=-0.001cpcgsfxh(i)
448	hfx (i,j)= -10.CPfxh(i) ! Bob: hfx (i,1)= -10.CPfxh(i)
449	ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
450	! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
451	! cpcgs = 1.00464e7
452	! convert from ergs/gram/K to J/kg/K cpmks=1004
453	! hfx (i,j)=-0.001cpcgsfxh(i)
454	hfx (i,j)= -10.CPfxh(i) ! Bob: hfx (i,j)= -10.CPfxh(i)
455	HFX(I,J)=AMAX1(HFX(I,J),-250.)
456	ENDIF
457	! QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
458	! convert from g/(cmcmsec) to kg/(mmsec)
459	qfx(i,j)=-10.*fxe(i)
460	QFX(I,J)=AMAX1(QFX(I,J),0.)
461	LH(I,J)=XLV*QFX(I,J)
462	ENDDO
463	ENDIF
464	! if(j.eq.2) write(0,*) 'u3d ,ustar,cdm at end of gfdlsfcmod'
465	! write(0,*) j,(u3d(ii,1,j),ii=70,90)
466	! write(0,*) j,(ustar(ii),ii=70,90)
467	! write(0,*) j,(cdm(ii),ii=70,90)
468	if(j.eq.jds.or.j.eq.jde) then
469
470	write(0,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde
471	write(0,*) "TSFC", (TSK(i,j),i=its,ite)
472	endif
473
474	ENDDO ! FOR THE J LOOP I PRESUME
475	! if(100.ge.its.and.100.le.ite.and.100.ge.jts.and.100.le.jte) then
476	! write(0,*) 'output vars of sf_gfdl at i,j=100'
477	! write(0,*) 'TSK',TSK(100,100)
478	! write(0,*) 'PSFC',PSFC(100,100)
479	! write(0,*) 'GLW',GLW(100,100)
480	! write(0,*) 'GSW',GSW(100,100)
481	! write(0,*) 'XLAND',XLAND(100,100)
482	! write(0,*) 'BR',BR(100,100)
483	! write(0,*) 'CHS',CHS(100,100)
484	! write(0,*) 'CHS2',CHS2(100,100)
485	! write(0,*) 'CPM',CPM(100,100)
486	! write(0,*) 'FLHC',FLHC(100,100)
487	! write(0,*) 'FLQC',FLQC(100,100)
488	! write(0,*) 'GZ1OZ0',GZ1OZ0(100,100)
489	! write(0,*) 'HFX',HFX(100,100)
490	! write(0,*) 'LH',LH(100,100)
491	! write(0,*) 'PSIM',PSIM(100,100)
492	! write(0,*) 'PSIH',PSIH(100,100)
493	! write(0,*) 'QFX',QFX(100,100)
494	! write(0,*) 'QGH',QGH(100,100)
495	! write(0,*) 'QSFC',QSFC(100,100)
496	! write(0,*) 'UST',UST(100,100)
497	! write(0,*) 'ZNT',ZNT(100,100)
498	! write(0,*) 'wet',wet(100)
499	! write(0,*) 'smois',smois(100,1,100)
500	! write(0,*) 'WSPD',WSPD(100,100)
501	! write(0,*) 'U10',U10(100,100)
502	! write(0,*) 'V10',V10(100,100)
503	! endif
504
505
506	END SUBROUTINE SF_GFDL
507
508	!-------------------------------------------------------------------
509	SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,tstrc, &
510	pspc,pkmax,wetc,slwdc,tjloc, &
511	upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, &
512	ids,ide, jds,jde, kds,kde, &
513	ims,ime, jms,jme, kms,kme, &
514	its,ite, jts,jte, kts,kte )
515
516	!------------------------------------------------------------------------
517	!
518	! MFLUX2 computes surface fluxes of momentum, heat,and moisture
519	! using monin-obukhov. the roughness length "z0" is prescribed
520	! over land and over ocean "z0" is computed using charnocks formula.
521	! the universal functions (from similarity theory approach) are
522	! those of hicks. This is Bob's doing.
523	!
524	!------------------------------------------------------------------------
525
526	IMPLICIT NONE
527	! use allocate_mod
528	! use module_TLDATA , ONLY : tab,table,cp,g,rgas,og
529
530	! include 'RESOLUTION.h'
531	! include 'PARAMETERS.h'
532	! include 'STDUNITS.h' stdout
533
534	! include 'FLAGS.h'
535	! include 'BKINFO.h' nstep
536	! include 'CONSLEV.h'
537	! include 'CONMLEV.h'
538	! include 'ESTAB.h'
539	! include 'FILEC.h'
540	! include 'FILEPC.h'
541	! include 'GDINFO.h' ngd
542	! include 'LIMIT.h'
543	! include 'QLOGS.h'
544	! include 'TIME.h' dt(nnst)
545	! include 'WINDD.h'
546	! include 'ZLDATA.h' old MOBFLX?
547
548	!-----------------------------------------------------------------------
549	! user interface variables
550	!-----------------------------------------------------------------------
551	integer,intent(in) :: ids,ide, jds,jde, kds,kde
552	integer,intent(in) :: ims,ime, jms,jme, kms,kme
553	integer,intent(in) :: its,ite, jts,jte, kts,kte
554	integer,intent(in) :: jfix,ntsflg
555
556	real, intent (out), dimension (ims :ime ) :: fxh
557	real, intent (out), dimension (ims :ime ) :: fxe
558	real, intent (out), dimension (ims :ime ) :: fxmx
559	real, intent (out), dimension (ims :ime ) :: fxmy
560	real, intent (out), dimension (ims :ime ) :: cdm
561	! real, intent (out), dimension (ims :ime ) :: cdm2
562	real, intent (out), dimension (ims :ime ) :: rib
563	real, intent (out), dimension (ims :ime ) :: xxfh
564	real, intent (out), dimension (ims :ime ) :: xxfh2
565	real, intent (out), dimension (ims :ime ) :: wind10
566
567	real, intent ( inout), dimension (ims :ime ) :: zoc
568	real, intent ( inout), dimension (ims :ime ) :: tstrc
569
570	real, intent ( in) :: dt
571	real, intent ( in) :: sfenth
572	real, intent ( in), dimension (ims :ime ) :: pspc
573	real, intent ( in), dimension (ims :ime ) :: pkmax
574	real, intent ( in), dimension (ims :ime ) :: wetc
575	real, intent ( in), dimension (ims :ime ) :: slwdc
576	real, intent ( in), dimension (ims :ime ) :: tjloc
577
578	real, intent ( in), dimension (kms:kme, ims :ime ) :: upc
579	real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc
580	real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc
581	real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc
582
583	!-----------------------------------------------------------------------
584	! internal variables
585	!-----------------------------------------------------------------------
586
587	integer, parameter :: icntx = 30
588
589	integer, dimension(1 :ime) :: ifz
590	integer, dimension(1 :ime) :: indx
591	integer, dimension(1 :ime) :: istb
592	integer, dimension(1 :ime) :: it
593	integer, dimension(1 :ime) :: iutb
594
595	real, dimension(1 :ime) :: aap
596	real, dimension(1 :ime) :: bq1
597	real, dimension(1 :ime) :: bq1p
598	real, dimension(1 :ime) :: delsrad
599	real, dimension(1 :ime) :: ecof
600	real, dimension(1 :ime) :: ecofp
601	real, dimension(1 :ime) :: estso
602	real, dimension(1 :ime) :: estsop
603	real, dimension(1 :ime) :: fmz1
604	real, dimension(1 :ime) :: fmz10
605	real, dimension(1 :ime) :: fmz2
606	real, dimension(1 :ime) :: fmzo1
607	real, dimension(1 :ime) :: foft
608	real, dimension(1 :ime) :: foftm
609	real, dimension(1 :ime) :: frac
610	real, dimension(1 :ime) :: land
611	real, dimension(1 :ime) :: pssp
612	real, dimension(1 :ime) :: qf
613	real, dimension(1 :ime) :: rdiff
614	real, dimension(1 :ime) :: rho
615	real, dimension(1 :ime) :: rkmaxp
616	real, dimension(1 :ime) :: rstso
617	real, dimension(1 :ime) :: rstsop
618	real, dimension(1 :ime) :: sf10
619	real, dimension(1 :ime) :: sf2
620	real, dimension(1 :ime) :: sfm
621	real, dimension(1 :ime) :: sfzo
622	real, dimension(1 :ime) :: sgzm
623	real, dimension(1 :ime) :: slwa
624	real, dimension(1 :ime) :: szeta
625	real, dimension(1 :ime) :: szetam
626	real, dimension(1 :ime) :: t1
627	real, dimension(1 :ime) :: t2
628	real, dimension(1 :ime) :: tab1
629	real, dimension(1 :ime) :: tab2
630	real, dimension(1 :ime) :: tempa1
631	real, dimension(1 :ime) :: tempa2
632	real, dimension(1 :ime) :: theta
633	real, dimension(1 :ime) :: thetap
634	real, dimension(1 :ime) :: tsg
635	real, dimension(1 :ime) :: tsm
636	real, dimension(1 :ime) :: tsp
637	real, dimension(1 :ime) :: tss
638	real, dimension(1 :ime) :: ucom
639	real, dimension(1 :ime) :: uf10
640	real, dimension(1 :ime) :: uf2
641	real, dimension(1 :ime) :: ufh
642	real, dimension(1 :ime) :: ufm
643	real, dimension(1 :ime) :: ufzo
644	real, dimension(1 :ime) :: ugzm
645	real, dimension(1 :ime) :: uzeta
646	real, dimension(1 :ime) :: uzetam
647	real, dimension(1 :ime) :: vcom
648	real, dimension(1 :ime) :: vrtkx
649	real, dimension(1 :ime) :: vrts
650	real, dimension(1 :ime) :: wind
651	real, dimension(1 :ime) :: windp
652	! real, dimension(1 :ime) :: xxfh
653	real, dimension(1 :ime) :: xxfm
654	real, dimension(1 :ime) :: xxsh
655	real, dimension(1 :ime) :: z10
656	real, dimension(1 :ime) :: z2
657	real, dimension(1 :ime) :: zeta
658	real, dimension(1 :ime) :: zkmax
659
660	real, dimension(1 :ime) :: pss
661	real, dimension(1 :ime) :: tstar
662	real, dimension(1 :ime) :: ukmax
663	real, dimension(1 :ime) :: vkmax
664	real, dimension(1 :ime) :: tkmax
665	real, dimension(1 :ime) :: rkmax
666	real, dimension(1 :ime) :: zot
667	real, dimension(1 :ime) :: fhzo1
668	real, dimension(1 :ime) :: sfh
669
670	real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
671	real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
672	real :: szet2, zal2,ugz2
673	real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
674	real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
675	real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
676	real :: shfx,sigt4,reflect
677	real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
678	real :: wndm,ckg
679	real :: yz,y1,y2,y3,y4,windmks,znott,znotm
680	integer:: i,j,ii,iq,nnest,icnt,ngd,ip
681	data amask/ -98.0/
682
683	!-----------------------------------------------------------------------
684	! internal variables
685	!-----------------------------------------------------------------------
686
687	real, dimension (223) :: tab
688	real, dimension (223) :: table
689	real, dimension (101) :: tab11
690	real, dimension (41) :: table4
691	real, dimension (42) :: tab3
692	real, dimension (54) :: table2
693	real, dimension (54) :: table3
694	real, dimension (74) :: table1
695	real, dimension (80) :: tab22
696
697	equivalence (tab(1),tab11(1))
698	equivalence (tab(102),tab22(1))
699	equivalence (tab(182),tab3(1))
700	equivalence (table(1),table1(1))
701	equivalence (table(75),table2(1))
702	equivalence (table(129),table3(1))
703	equivalence (table(183),table4(1))
704
705	!-----------------------------------------------------------------------
706	! tables used to obtain the vapor pressures or saturated vapor
707	! pressure
708	!-----------------------------------------------------------------------
709
710	data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
711	&.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
712	&.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
713	&1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
714	&8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
715	&34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
716	&128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
717	&421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
718
719	data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
720	&,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
721	&5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
722	&10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
723	&18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
724	&31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
725	&53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
726	&86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
727	&129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
728	&190220.,199260./
729
730	data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
731	&285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
732	&401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
733	&555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
734	&756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
735	&1013250.,1049940.,1087740./
736
737	data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
738	& .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
739	&.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
740	&.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
741	&.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
742	&.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
743	&.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
744	&.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
745	&.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
746	&.7220e+01/
747
748	data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
749	&.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
750	&.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
751	&.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
752	&.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
753	&.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
754	&.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
755	&.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
756	&.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
757	&.7090e+03/
758
759	data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
760	&.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
761	&.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
762	&.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
763	&.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
764	&.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
765	&.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
766	&.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
767	&.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
768	&.9780e+04/
769
770	data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
771	&.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
772	&.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
773	&.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
774	&.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
775	&.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
776	&.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
777	!
778	! spcify constants needed by MFLUX2
779	!
780	real,parameter :: cp = 1.00464e7
781	real,parameter :: g = 980.6
782	real,parameter :: rgas = 2.87e6
783	real,parameter :: og = 1./g
784
785	!
786	! character*10 routine
787	! routine = 'mflux2'
788	!
789	!------------------------------------------------------------------------
790	! set water availability constant "ecof" and land mask "land".
791	! limit minimum wind speed to 100 cm/s
792	!------------------------------------------------------------------------
793	j = IFIX(tjloc(its))
794	! constants for 10 m winds (correction for knots
795	!
796	cor1 = .120
797	cor2 = 720.
798	yz= -0.0001344
799	y1= 3.015e-05
800	y2= 1.517e-06
801	y3= -3.567e-08
802	y4= 2.046e-10
803
804	do i = its,ite
805	z10(i) = 1000.
806	z2 (i) = 200.
807	pss(i) = pspc(i)
808	tstar(i) = tstrc(i)
809	ukmax(i) = upc(1,i)
810	vkmax(i) = vpc(1,i)
811	tkmax(i) = tpc(1,i)
812	rkmax(i) = rpc(1,i)
813	enddo
814
815	! write(0,*)'z10,pss,tstar,u...rkmax(ite)', &
816	! z10(ite), pss(ite),tstar(ite),ukmax(ite), &
817	! vkmax(ite),tkmax(ite),rkmax(ite)
818
819	do i = its,ite
820	windp(i) = SQRT(ukmax(i)ukmax(i) + vkmax(i)vkmax(i))
821	wind (i) = amax1(windp(i),100.)
822	if (zoc(i) .LT. amask) zoc(i) = -0.01850.001wind(i)wind(i)og
823	if (zoc(i) .GT. 0.0) then
824	ecof(i) = wetc(i)
825	land(i) = 1.0
826	zot (i) = zoc(i)
827	else
828	ecof(i) = wetc(i)
829	land(i) = 0.0
830	#ifdef HWRF
831	zot (i) = zoc(i)
832	! now use 2 regime fit for znot thermal
833	windmks=wind(i)*.01
834	znott=0.2375exp(-0.5250windmks) + 0.0025exp(-0.0211windmks)
835	znott=0.01*znott
836	! go back to moon et al for below 7m/s
837	if(windmks.le. 7.) &
838	znott = (0.0185/9.8(7.59e-8wind(i)**2+ &
839	2.46e-4wind(i))*2)
840	! back to cgs
841	zot (i) = 100.*znott
842	#else
843	! zot (i) = zoc(i)
844	! now use 2 regime fit for znot thermal
845	windmks=wind(i)*.01
846	znott=1.9551e-5 - 2.6338e-7 * windmks
847	if(windmks.le.10.) znott=0.0025542 * windmks **(-1.8023)
848	znott=amax1(1.e-6,znott)
849	! go back to moon et al for below 7m/s
850	if(windmks.le. 7.) &
851	znott = (0.0185/9.8(7.59e-8wind(i)**2+ &
852	2.46e-4wind(i))*2)
853	! back to cgs
854	zot (i) = 100.*znott
855	#endif
856	! end of kwon correction....
857	! in hwrf, thermal znot(zot) is passed as argument zoc
858	! in hwrf, momentum znot is recalculated internally
859	zoc(i)=-(0.0185/9.8(7.59e-8wind(i)**2+ &
860	2.46e-4wind(i))2)100.
861	if(wind(i).ge.1250.0) &
862	zoc(i)=-(.000739793 * wind(i) -0.58)/10
863	if(wind(i).ge.3000.) then
864	windmks=wind(i)*.01
865	! kwon znotm
866	znotm = yz +windmksy1 +windmks2y2 +windmks*3y3 +windmks*4y4 !powell 2003
867	! back to cgs
868	zoc(i) = 100.*znotm
869	endif
870	endif
871
872	!------------------------------------------------------------------------
873	! where necessary modify zo values over ocean.
874	!------------------------------------------------------------------------
875
876	enddo
877
878	!------------------------------------------------------------------------
879	! define constants:
880	! a and c = constants used in evaluating universal function for
881	! stable case
882	! ca = karmen constant
883	! cm01 = constant part of vertical integral of universal
884	! function; stable case ( 0.5 < zeta < or = 10.0)
885	! cm02 = constant part of vertical integral of universal
886	! function; stable case ( zeta > 10.0)
887	!------------------------------------------------------------------------
888
889	en = 2.
890	c = .76
891	a = 5.
892	ca = .4
893	cmo1 = .5*a - 1.648
894	cmo2 = 17.193 + .5a - 10.c
895	boycon = .61
896	rovcp=rgas/cp
897	! write(0,*)'rgas,cp,rovcp ', rgas,cp,rovcp
898
899	! write(0,*)'--------------------------------------------------'
900	! write(0,*)'pkmax, pspc, theta, zkmax, zoc'
901	! write(0,*)'--------------------------------------------------'
902
903	do i = its,ite
904	! theta(i) = tkmax(i)*rqc9
905	theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
906	vrtkx(i) = 1.0 + boycon*rkmax(i)
907	! zkmax(i) = rgastkmax(i)qqlog(kmax)*og
908	zkmax(i) = -rgastkmax(i)alog(pkmax(i)/pspc(i))*og
909	! IF(I==78)write(0,*)I,JFIX,pkmax(i),pspc(i),theta(i),zkmax(i),zoc(i)
910	enddo
911
912	! write(0,*)'pkmax,pspc ', pkmax,pspc
913	! write(0,*)'theta, zkmax, zoc ', theta, zkmax, zoc
914
915	!------------------------------------------------------------------------
916	! get saturation mixing ratios at surface
917	!------------------------------------------------------------------------
918
919	do i = its,ite
920	tsg (i) = tstar(i)
921	tab1 (i) = tstar(i) - 153.16
922	it (i) = IFIX(tab1(i))
923	tab2 (i) = tab1(i) - FLOAT(it(i))
924	t1 (i) = tab(it(i) + 1)
925	t2 (i) = table(it(i) + 1)
926	estso(i) = t1(i) + tab2(i)*t2(i)
927	psps1 = (pss(i) - estso(i))
928	if(psps1 .EQ. 0.0)then
929	psps1 = .1
930	endif
931	rstso(i) = 0.622*estso(i)/psps1
932	vrts (i) = 1. + boyconecof(i)rstso(i)
933	enddo
934
935	!------------------------------------------------------------------------
936	! check if consideration of virtual temperature changes stability.
937	! if so, set "dthetav" to near neutral value (1.0e-4). also check
938	! for very small lapse rates; if ABS(tempa1) <1.0e-4 then
939	! tempa1=1.0e-4
940	!------------------------------------------------------------------------
941
942	do i = its,ite
943	tempa1(i) = theta(i)vrtkx(i) - tstar(i)vrts(i)
944	tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
945	if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
946	tab1(i) = ABS(tempa1(i))
947	if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
948	!------------------------------------------------------------------------
949	! compute bulk richardson number "rib" at each point. if "rib"
950	! exceeds 95% of critical richardson number "tab1" then "rib = tab1"
951	!------------------------------------------------------------------------
952
953	rib (i) = gzkmax(i)tempa1(i)/ &
954	(tkmax(i)vrtkx(i)wind(i)*wind(i))
955	tab2(i) = ABS(zoc(i))
956	tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
957	if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
958	enddo
959
960	do i = its,ite
961	zeta(i) = ca*rib(i)/0.03
962	enddo
963
964	! write(0,*)'rib,zeta,vrtkx,vrts(ite) ', rib(ite),zeta(ite), &
965	! vrtkx(ite),vrts(ite)
966	!------------------------------------------------------------------------
967	! begin looping through points on line, solving wegsteins iteration
968	! for zeta at each point, and using hicks functions
969	!------------------------------------------------------------------------
970
971	!------------------------------------------------------------------------
972	! set initial guess of zeta=non - dimensional height "szeta" for
973	! stable points
974	!------------------------------------------------------------------------
975
976	rca = 1./ca
977	enrca = en*rca
978	! turn off interfacial layer by zeroing out enrca
979	enrca = 0.0
980	zog = .0185*og
981
982	!------------------------------------------------------------------------
983	! stable points
984	!------------------------------------------------------------------------
985
986	ip = 0
987	do i = its,ite
988	if (zeta(i) .GE. 0.0) then
989	ip = ip + 1
990	istb(ip) = i
991	endif
992	enddo
993
994	if (ip .EQ. 0) go to 170
995	do i = 1,ip
996	szetam(i) = 1.0e+30
997	sgzm(i) = 0.0e+00
998	szeta(i) = zeta(istb(i))
999	ifz(i) = 1
1000	enddo
1001
1002	!------------------------------------------------------------------------
1003	! begin wegstein iteration for "zeta" at stable points using
1004	! hicks(1976)
1005	!------------------------------------------------------------------------
1006
1007	do icnt = 1,icntx
1008	do i = 1,ip
1009	if (ifz(i) .EQ. 0) go to 80
1010	zal1g = ALOG(szeta(i))
1011	if (szeta(i) .LE. 0.5) then
1012	fmz1(i) = (zal1g + aszeta(i))rca
1013	else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
1014	rzeta1 = 1./szeta(i)
1015	fmz1(i) = (8.zal1g + 4.25rzeta1 - &
1016	0.5rzeta1rzeta1 + cmo1)*rca
1017	else if (szeta(i) .GT. 10.) then
1018	fmz1(i) = (cszeta(i) + cmo2)rca
1019	endif
1020	szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
1021	zal2g = ALOG(szetao)
1022	if (szetao .LE. 0.5) then
1023	fmzo1(i) = (zal2g + aszetao)rca
1024	sfzo (i) = 1. + a*szetao
1025	else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
1026	rzeta2 = 1./szetao
1027	fmzo1(i) = (8.zal2g + 4.25rzeta2 - &
1028	0.5rzeta2rzeta2 + cmo1)*rca
1029	sfzo (i) = 8.0 - 4.25rzeta2 + rzeta2rzeta2
1030	else if (szetao .GT. 10.) then
1031	fmzo1(i) = (cszetao + cmo2)rca
1032	sfzo (i) = c*szetao
1033	endif
1034
1035
1036	! compute heat & moisture parts of zot.. for calculation of sfh
1037
1038	szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i)
1039	zal2gh = ALOG(szetho)
1040	if (szetho .LE. 0.5) then
1041	fhzo1(i) = (zal2gh + aszetho)rca
1042	sfzo (i) = 1. + a*szetho
1043	else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
1044	rzeta2 = 1./szetho
1045	fhzo1(i) = (8.zal2gh + 4.25rzeta2 - &
1046	0.5rzeta2rzeta2 + cmo1)*rca
1047	sfzo (i) = 8.0 - 4.25rzeta2 + rzeta2rzeta2
1048	else if (szetho .GT. 10.) then
1049	fhzo1(i) = (cszetho + cmo2)rca
1050	sfzo (i) = c*szetho
1051	endif
1052
1053	!------------------------------------------------------------------------
1054	! compute universal function at 10 meters for diagnostic purposes
1055	!------------------------------------------------------------------------
1056
1057	!!!! if (ngd .EQ. nNEST) then
1058	szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i)
1059	zal10 = ALOG(szet10)
1060	if (szet10 .LE. 0.5) then
1061	fmz10(i) = (zal10 + aszet10)rca
1062	else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
1063	rzeta2 = 1./szet10
1064	fmz10(i) = (8.zal10 + 4.25rzeta2 - &
1065	0.5rzeta2rzeta2 + cmo1)*rca
1066	else if (szet10 .GT. 10.) then
1067	fmz10(i) = (cszet10 + cmo2)rca
1068	endif
1069	sf10(i) = fmz10(i) - fmzo1(i)
1070	! compute 2m values for diagnostics in HWRF
1071	szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i)
1072	zal2 = ALOG(szet2 )
1073	if (szet2 .LE. 0.5) then
1074	fmz2 (i) = (zal2 + aszet2 )rca
1075	else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
1076	rzeta2 = 1./szet2
1077	fmz2 (i) = (8.zal2 + 4.25rzeta2 - &
1078	0.5rzeta2rzeta2 + cmo1)*rca
1079	else if (szet2 .GT. 2.) then
1080	fmz2 (i) = (cszet2 + cmo2)rca
1081	endif
1082	sf2 (i) = fmz2 (i) - fmzo1(i)
1083
1084	!!!! endif
1085	sfm(i) = fmz1(i) - fmzo1(i)
1086	sfh(i) = fmz1(i) - fhzo1(i)
1087	sgz = carib(istb(i))sfm(i)*sfm(i)/ &
1088	(sfh(i) + enrca*sfzo(i))
1089	fmz = (sgz - szeta(i))/szeta(i)
1090	fmzo = ABS(fmz)
1091	if (fmzo .GE. 5.0e-5) then
1092	sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
1093	if(sq .EQ. 1) then
1094	write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (STABLE CASE)'
1095	write(0,*)'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
1096	endif
1097	szetam(i) = szeta(i)
1098	szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
1099	sgzm (i) = sgz
1100	else
1101	ifz(i) = 0
1102	endif
1103	80 continue
1104	enddo
1105	enddo
1106
1107	do i = 1,ip
1108	if (ifz(i) .GE. 1) go to 110
1109	enddo
1110
1111	go to 130
1112
1113	110 continue
1114	write(6,120)
1115	120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
1116	! call MPI_CLOSE(1,routine)
1117
1118	!------------------------------------------------------------------------
1119	! update "zo" for ocean points. "zo"cannot be updated within the
1120	! wegsteins iteration as the scheme (for the near neutral case)
1121	! can become unstable
1122	!------------------------------------------------------------------------
1123
1124	130 continue
1125	do i = 1,ip
1126	szo = zoc(istb(i))
1127	if (szo .LT. 0.0) then
1128	wndm=wind(istb(i))*0.01
1129	if(wndm.lt.15.0) then
1130	ckg=0.0185*og
1131	else
1132	!! ckg=(0.000308wndm+0.00925)og
1133	!! ckg=(0.000616wndm)og
1134	ckg=(sfenth(40.000308wndm) + (1.-sfenth)0.0185 )*og
1135	endif
1136
1137	szo = - ckgwind(istb(i))wind(istb(i))/ &
1138	(sfm(i)*sfm(i))
1139	cons_p000001 = .000001
1140	cons_7 = 7.
1141	vis = 1.4E-1
1142
1143	ustar = sqrt( -szo / zog)
1144	restar = -ustar * szo / vis
1145	restar = max(restar,cons_p000001)
1146	! Rat taken from Zeng, Zhao and Dickinson 1997
1147	rat = 2.67 * restar ** .25 - 2.57
1148	rat = min(rat ,cons_7) !constant
1149	rat=0.
1150	zot(istb(i)) = szo * exp(-rat)
1151	else
1152	zot(istb(i)) = zoc(istb(i))
1153	endif
1154
1155	! in hwrf thermal znot is loaded back into the zoc array for next step
1156	zoc(istb(i)) = szo
1157	enddo
1158
1159	do i = 1,ip
1160	xxfm(istb(i)) = sfm(i)
1161	xxfh(istb(i)) = sfh(i)
1162	xxfh2(istb(i)) = sf2 (i)
1163	xxsh(istb(i)) = sfzo(i)
1164	enddo
1165
1166	!------------------------------------------------------------------------
1167	! obtain wind at 10 meters for diagnostic purposes
1168	!------------------------------------------------------------------------
1169
1170	!!! if (ngd .EQ. nNEST) then
1171	do i = 1,ip
1172	wind10(istb(i)) = sf10(i)*wind(istb(i))/sfm(i)
1173	wind10(istb(i)) = wind10(istb(i)) * 1.944
1174	if(wind10(istb(i)) .GT. 6000.0) then
1175	wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
1176	- cor2
1177	endif
1178	! the above correction done by GFDL in centi-kts!!!-change back
1179	wind10(istb(i)) = wind10(istb(i)) / 1.944
1180	enddo
1181	!!! endif
1182	!!! if (ngd .EQ. nNEST-1 .AND. llwe .EQ. 1 ) then
1183	!!! do i = 1,ip
1184	!!! wind10c(istb(i),j) = sf10(i)*wind(istb(i))/sfm(i)
1185	!!! wind10c(istb(i),j) = wind10c(istb(i),j) * 1.944
1186	!!! if(wind10c(istb(i),j) .GT. 6000.0) then
1187	!!! wind10c(istb(i),j)=wind10c(istb(i),j)+wind10c(istb(i),j)*cor1
1188	!!! * - cor2
1189	!!! endif
1190	!!! enddo
1191	!!! endif
1192
1193	!------------------------------------------------------------------------
1194	! unstable points
1195	!------------------------------------------------------------------------
1196
1197	170 continue
1198
1199	iq = 0
1200	do i = its,ite
1201	if (zeta(i) .LT. 0.0) then
1202	iq = iq + 1
1203	iutb(iq) = i
1204	endif
1205	enddo
1206
1207	if (iq .EQ. 0) go to 290
1208	do i = 1,iq
1209	uzeta (i) = zeta(iutb(i))
1210	ifz (i) = 1
1211	uzetam(i) = 1.0e+30
1212	ugzm (i) = 0.0e+00
1213	enddo
1214
1215	!------------------------------------------------------------------------
1216	! begin wegstein iteration for "zeta" at unstable points using
1217	! hicks functions
1218	!------------------------------------------------------------------------
1219
1220	do icnt = 1,icntx
1221	do i = 1,iq
1222	if (ifz(i) .EQ. 0) go to 200
1223	ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i))))
1224	uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1225	ux11 = 1. - 16.*uzeta(i)
1226	ux12 = 1. - 16.*uzetao
1227	y = SQRT(ux11)
1228	yo = SQRT(ux12)
1229	ufzo(i) = 1./yo
1230	ux13 = (1. + y)/(1. + yo)
1231	ux21 = ALOG(ux13)
1232	ufh(i) = (ugzzo - 2.ux21)rca
1233	! recompute scalers for ufm in terms of mom znot... zoc
1234	ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i))))
1235	uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1236	ux11 = 1. - 16.*uzeta(i)
1237	ux12 = 1. - 16.*uzetao
1238	y = SQRT(ux11)
1239	yo = SQRT(ux12)
1240	ux13 = (1. + y)/(1. + yo)
1241	ux21 = ALOG(ux13)
1242	! ufzo(i) = 1./yo
1243	x = SQRT(y)
1244	xo = SQRT(yo)
1245	xnum = (x*2 + 1.)((x + 1.)**2)
1246	xden = (xo*2 + 1.)((xo + 1.)**2)
1247	xtan = ATAN(x) - ATAN(xo)
1248	ux3 = ALOG(xnum/xden)
1249	ufm(i) = (ugzzo - ux3 + 2.xtan)rca
1250	!!!! if (ngd .EQ. nNEST) then
1251
1252	!------------------------------------------------------------------------
1253	! obtain ten meter winds for diagnostic purposes
1254	!------------------------------------------------------------------------
1255
1256	ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i))))
1257	uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1258	uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1259	ux11 = 1. - 16.*uzet1o
1260	ux12 = 1. - 16.*uzetao
1261	y = SQRT(ux11)
1262	y10 = SQRT(ux12)
1263	ux13 = (1. + y)/(1. + y10)
1264	ux21 = ALOG(ux13)
1265	x = SQRT(y)
1266	x10 = SQRT(y10)
1267	xnum = (x*2 + 1.)((x + 1.)**2)
1268	xden = (x10*2 + 1.)((x10 + 1.)**2)
1269	xtan = ATAN(x) - ATAN(x10)
1270	ux3 = ALOG(xnum/xden)
1271	uf10(i) = (ugz10 - ux3 + 2.xtan)rca
1272
1273	! obtain 2m values for diagnostics...
1274
1275
1276	ugz2 = ALOG(z2 (iutb(i))/ABS(zoc(iutb(i))))
1277	uzet1o = ABS(z2 (iutb(i)))/zkmax(iutb(i))*uzeta(i)
1278	uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1279	ux11 = 1. - 16.*uzet1o
1280	ux12 = 1. - 16.*uzetao
1281	y = SQRT(ux11)
1282	yo = SQRT(ux12)
1283	ux13 = (1. + y)/(1. + yo)
1284	ux21 = ALOG(ux13)
1285	uf2 (i) = (ugzzo - 2.ux21)rca
1286
1287
1288	!!! endif
1289	ugz = carib(iutb(i))ufm(i)ufm(i)/(ufh(i) + enrcaufzo(i))
1290	ux1 = (ugz - uzeta(i))/uzeta(i)
1291	ux2 = ABS(ux1)
1292	if (ux2 .GE. 5.0e-5) then
1293	uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
1294	uzetam(i) = uzeta(i)
1295	if(uq .EQ. 1) then
1296	write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (UNSTABLE CASE)'
1297	write(0,*)'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1298	endif
1299	uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
1300	ugzm (i) = ugz
1301	else
1302	ifz(i) = 0
1303	endif
1304	200 continue
1305	enddo
1306	enddo
1307
1308
1309	do i = 1,iq
1310	if (ifz(i) .GE. 1) go to 230
1311	enddo
1312
1313	go to 250
1314
1315	230 continue
1316	write(6,240)
1317	240 format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ')
1318	! call MPI_CLOSE(1,routine)
1319
1320	!------------------------------------------------------------------------
1321	! gather unstable values
1322	!------------------------------------------------------------------------
1323
1324	250 continue
1325
1326	!------------------------------------------------------------------------
1327	! update "zo" for ocean points. zo cannot be updated within the
1328	! wegsteins iteration as the scheme (for the near neutral case)
1329	! can become unstable.
1330	!------------------------------------------------------------------------
1331
1332	do i = 1,iq
1333	uzo = zoc(iutb(i))
1334	if (zoc(iutb(i)) .LT. 0.0) then
1335	wndm=wind(iutb(i))*0.01
1336	if(wndm.lt.15.0) then
1337	ckg=0.0185*og
1338	else
1339	!! ckg=(0.000308wndm+0.00925)og <
1340	!! ckg=(0.000616wndm)og <
1341	ckg=(40.000308wndm)*og
1342	ckg=(sfenth(40.000308wndm) + (1.-sfenth)0.0185 )*og
1343	endif
1344	uzo =-ckgwind(iutb(i))wind(iutb(i))/(ufm(i)*ufm(i))
1345	cons_p000001 = .000001
1346	cons_7 = 7.
1347	vis = 1.4E-1
1348
1349	ustar = sqrt( -uzo / zog)
1350	restar = -ustar * uzo / vis
1351	restar = max(restar,cons_p000001)
1352	! Rat taken from Zeng, Zhao and Dickinson 1997
1353	rat = 2.67 * restar ** .25 - 2.57
1354	rat = min(rat ,cons_7) !constant
1355	rat=0.0
1356	zot(iutb(i)) = uzo * exp(-rat)
1357	else
1358	zot(iutb(i)) = zoc(iutb(i))
1359	endif
1360	! in hwrf thermal znot is loaded back into the zoc array for next step
1361	zoc(iutb(i)) = uzo
1362	enddo
1363
1364	!------------------------------------------------------------------------
1365	! obtain wind at ten meters for diagnostic purposes
1366	!------------------------------------------------------------------------
1367
1368	!!! if (ngd .EQ. nNEST) then
1369	do i = 1,iq
1370	wind10(iutb(i)) = uf10(i)*wind(iutb(i))/ufm(i)
1371	wind10(iutb(i)) = wind10(iutb(i)) * 1.944
1372	if(wind10(iutb(i)) .GT. 6000.0) then
1373	wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
1374	- cor2
1375	endif
1376	! the above correction done by GFDL in centi-kts!!!-change back
1377	wind10(iutb(i)) = wind10(iutb(i)) / 1.944
1378	enddo
1379	!!! endif
1380	!!! if (ngd .EQ. nNEST-1) then
1381	!!! do i = 1,iq
1382	!!! wind10c(iutb(i),j) = uf10(i)*wind(iutb(i))/ufm(i)
1383	!!! wind10c(iutb(i),j) = wind10c(iutb(i),j) * 1.944
1384	!!! if(wind10c(iutb(i),j) .GT. 6000.0) then
1385	!!! wind10c(iutb(i),j)=wind10c(iutb(i),j)+wind10c(iutb(i),j)*cor1
1386	!!! * - cor2
1387	!!! endif
1388	!!! enddo
1389	!!! endif
1390
1391	do i = 1,iq
1392	xxfm(iutb(i)) = ufm(i)
1393	xxfh(iutb(i)) = ufh(i)
1394	xxfh2(iutb(i)) = uf2 (i)
1395	xxsh(iutb(i)) = ufzo(i)
1396	enddo
1397
1398	290 continue
1399
1400	do i = its,ite
1401	ucom(i) = ukmax(i)
1402	vcom(i) = vkmax(i)
1403	if (windp(i) .EQ. 0.0) then
1404	windp(i) = 100.0
1405	ucom (i) = 100.0/SQRT(2.0)
1406	vcom (i) = 100.0/SQRT(2.0)
1407	endif
1408	rho(i) = pss(i)/(rgas(tsg(i) + enrca(theta(i) - &
1409	tsg(i))xxsh(i)/(xxfh(i) + enrcaxxsh(i))))
1410	bq1(i) = wind(i)rho(i)/(xxfm(i)(xxfh(i) + enrca*xxsh(i)))
1411	enddo
1412
1413	! do land sfc temperature prediction if ntsflg=1
1414	! ntsflg = 1 ! gopal's doing
1415
1416	if (ntsflg .EQ. 0) go to 370
1417	alll = 600.
1418	xks = 0.01
1419	hcap = .5/2.39e-8
1420	pith = SQRT(4.*ATAN(1.0))
1421	alfus = alll/2.39e-8
1422	teps = 0.1
1423	! slwdc... in units of cal/min ????
1424	! slwa... in units of ergs/sec/cm*2
1425	! 1 erg=2.39e-8 cal
1426	!------------------------------------------------------------------------
1427	! pack land and sea ice points
1428	!------------------------------------------------------------------------
1429
1430	ip = 0
1431	do i = its,ite
1432	if (land(i) .EQ. 1) then
1433	ip = ip + 1
1434	indx (ip) = i
1435	! slwa is defined as positive down....
1436	slwa (ip) = slwdc(i)/(2.39e-8*60.)
1437	tss (ip) = tstar(i)
1438	thetap (ip) = theta(i)
1439	rkmaxp (ip) = rkmax(i)
1440	aap (ip) = 5.673e-5
1441	pssp (ip) = pss(i)
1442	ecofp (ip) = ecof(i)
1443	estsop (ip) = estso(i)
1444	rstsop (ip) = rstso(i)
1445	bq1p (ip) = bq1(i)
1446	bq1p (ip) = amax1(bq1p(ip),0.1e-3)
1447	delsrad(ip) = dt pith/(hcapSQRT(3600.24.xks))
1448	endif
1449	enddo
1450
1451	!------------------------------------------------------------------------
1452	! initialize variables for first pass of iteration
1453	!------------------------------------------------------------------------
1454
1455	do i = 1,ip
1456	ifz (i) = 1
1457	tsm (i) = tss(i)
1458	rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1459	!!! if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. &
1460	!!! nstep .EQ. -99 .AND. ngd .EQ. 1) then
1461
1462	!!! if (j .EQ. 1 .AND. i .EQ. 1) write(6,300)
1463	300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ')
1464
1465	!! foftm(i) = thetap(i) + 1./(cpbq1p(i))(slwa(i) - aap(i)* &
1466	!! tsm(i)*4 + ecofp(i)alfusbq1p(i)rdiff(i))
1467	!! else
1468
1469	foftm(i) = tss(i) + delsrad(i)(slwa(i) - aap(i)tsm(i)**4 - &
1470	cpbq1p(i)(tsm(i) - thetap(i)) + ecofp(i)alfusbq1p(i)* &
1471	rdiff(i))
1472	!! endif
1473	tsp(i) = foftm(i)
1474	enddo
1475
1476	!------------------------------------------------------------------------
1477	! do iteration to determine "tstar" at new time level
1478	!------------------------------------------------------------------------
1479
1480	do icnt = 1,icntx
1481	do i = 1,ip
1482	if (ifz(i) .EQ. 0) go to 330
1483	tab1 (i) = tsp(i) - 153.16
1484	it (i) = IFIX(tab1(i))
1485	tab2 (i) = tab1(i) - FLOAT(it(i))
1486	t1 (i) = tab(it(i) + 1)
1487	t2 (i) = table(it(i) + 1)
1488	estsop(i) = t1(i) + tab2(i)*t2(i)
1489	psps2 = (pssp(i) - estsop(i))
1490	if(psps2 .EQ. 0.0)then
1491	psps2 = .1
1492	endif
1493	rstsop(i) = 0.622*estsop(i)/psps2
1494	rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1495	!!! if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. &
1496	!!! nstep .EQ. -99 .AND. ngd .EQ. 1) then
1497	!!! foft(i) = thetap(i) + (1./(cpbq1p(i)))(slwa(i) - aap(i)* &
1498	!!! tsp(i)*4 + ecofp(i)alfusbq1p(i)rdiff(i))
1499	!!! else
1500	foft(i) = tss(i) + delsrad(i)(slwa(i) - aap(i)tsp(i)**4 - &
1501	cpbq1p(i)(tsp(i) - thetap(i)) + ecofp(i)alfusbq1p(i)* &
1502	rdiff(i))
1503	!!! endif
1504	!!! if (ngd .EQ. 1 .AND. j .EQ. 48 .AND. i .EQ. 19) then
1505	!!! reflect = slwa(i)
1506	!!! sigt4 = -aap(i)tsp(i)*4
1507	!!! shfx = -cpbq1p(i)(tsp(i) - thetap(i))
1508	!!! alevp = ecofp(i)alfusbq1p(i)*rdiff(i)
1509	!!! delten = delsrad(i)
1510	!!! diffot = foft(i) - tss(i)
1511	!!! endif
1512	frac(i) = ABS((foft(i) - tsp(i))/tsp(i))
1513
1514	!------------------------------------------------------------------------
1515	! check for convergence of all points use wegstein iteration
1516	!------------------------------------------------------------------------
1517
1518	if (frac(i) .GE. teps) then
1519	qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
1520	tsm (i) = tsp(i)
1521	tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
1522	foftm(i) = foft(i)
1523	else
1524	ifz(i) = 0
1525	endif
1526	330 continue
1527	enddo
1528	enddo
1529
1530	!------------------------------------------------------------------------
1531	! check for convergence of "t star" prediction
1532	!------------------------------------------------------------------------
1533
1534	do i = 1,ip
1535	if (ifz(i) .EQ. 1) then
1536	write(6, 340) tsp(i), i, j
1537	340 format(2X, ' NON-CONVERGENCE OF T* PREDICTED (T*,I,J) = ', E14.8, &
1538	2I5)
1539
1540	write(6,345) indx(i), j, tstar(indx(i)), tsp(i), ip
1541	345 format(2X, ' I, J, OLD T, NEW T, NPTS ', 2I5, 2E14.8, I5)
1542
1543	write(6,350) reflect, sigt4, shfx, alevp, delten, diffot
1544	350 format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', &
1545	6E14.8)
1546
1547	! call MPI_CLOSE(1,routine)
1548	endif
1549	enddo
1550
1551	do i = 1,ip
1552	ii = indx(i)
1553	tstrc(ii) = tsp (i)
1554	enddo
1555
1556	!------------------------------------------------------------------------
1557	! compute fluxes and momentum drag coef
1558	!------------------------------------------------------------------------
1559
1560	370 continue
1561	do i = its,ite
1562	fxh(i) = bq1(i)*(theta(i) - tsg(i))
1563	fxe(i) = ecof(i)bq1(i)(rkmax(i) - rstso(i))
1564	if (fxe(i) .GT. 0.0) fxe(i) = 0.0
1565	fxmx(i) = rho(i)/(xxfm(i)xxfm(i))wind(i)wind(i)ucom(i)/ &
1566	windp(i)
1567	fxmy(i) = rho(i)/(xxfm(i)xxfm(i))wind(i)wind(i)vcom(i)/ &
1568	windp(i)
1569	cdm(i) = 1./(xxfm(i)*xxfm(i))
1570	! print *, 'i,zot,zoc,cdm,cdm2,tsg,wind', &
1571	! i, zot(i),zoc(i), cdm(i),cdm2(i), tsg(i),wind(i)
1572	enddo
1573
1574	return
1575	end subroutine MFLUX2
1576
1577	SUBROUTINE hwrfsfcinit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
1578	SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
1579	ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
1580	TMN, &
1581	num_soil_layers, &
1582	allowed_to_read, &
1583	ids,ide, jds,jde, kds,kde, &
1584	ims,ime, jms,jme, kms,kme, &
1585	its,ite, jts,jte, kts,kte )
1586
1587	IMPLICIT NONE
1588
1589	! Arguments
1590	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1591	ims,ime, jms,jme, kms,kme, &
1592	its,ite, jts,jte, kts,kte
1593
1594	INTEGER, INTENT(IN) :: num_soil_layers
1595
1596	REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
1597
1598	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1599	INTENT(INOUT) :: SMOIS, &
1600	TSLB !STEMP
1601
1602	REAL, DIMENSION( ims:ime, jms:jme ) , &
1603	INTENT(INOUT) :: SNOW, &
1604	SNOWC, &
1605	CANWAT, &
1606	SMSTAV, &
1607	SMSTOT, &
1608	SFCRUNOFF, &
1609	UDRUNOFF, &
1610	SFCEVP, &
1611	GRDFLX, &
1612	ACSNOW, &
1613	XICE, &
1614	VEGFRA, &
1615	TMN, &
1616	ACSNOM
1617
1618	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1619	INTENT(INOUT) :: IVGTYP, &
1620	ISLTYP
1621
1622	!
1623
1624	INTEGER, INTENT(IN) :: isn
1625	LOGICAL, INTENT(IN) :: allowed_to_read
1626	! Local
1627	INTEGER :: iseason
1628	INTEGER :: icm,jcm,itf,jtf
1629	INTEGER :: I,J,L
1630
1631
1632	itf=min0(ite,ide-1)
1633	jtf=min0(jte,jde-1)
1634
1635	icm = ide/2
1636	jcm = jde/2
1637
1638	iseason=isn
1639
1640	DO J=jts,jtf
1641	DO I=its,itf
1642	! SNOW(i,j)=0.
1643	SNOWC(i,j)=0.
1644	! SMSTAV(i,j)=
1645	! SMSTOT(i,j)=
1646	! SFCRUNOFF(i,j)=
1647	! UDRUNOFF(i,j)=
1648	! GRDFLX(i,j)=
1649	! ACSNOW(i,j)=
1650	! ACSNOM(i,j)=
1651	ENDDO
1652	ENDDO
1653
1654	END SUBROUTINE hwrfsfcinit
1655
1656	END MODULE module_sf_gfdl

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

Download in other formats:

Original Format