Context Navigation

turbdiff.F90 @ 1834

Last change on this file since 1834 was 1834, checked in by mlefevre, 7 years ago
MESOSCALE GENERIC. added turb_mod and changed turbdiff to send info on ustar to dynamical core in LES mode. by the way ustar is now calculated.
File size: 27.8 KB

Line
1	subroutine turbdiff(ngrid,nlay,nq,rnat, &
2	ptimestep,pcapcal,lecrit, &
3	pplay,pplev,pzlay,pzlev,pz0, &
4	pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf, &
5	pdtfi,pdqfi,pfluxsrf, &
6	Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2, &
7	pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
8
9	use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
10	use radcommon_h, only : sigma, glat
11	use surfdat_h, only: dryness
12	use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
13	use comcstfi_mod, only: rcp, g, r, cpp
14	use callkeys_mod, only: water,tracer,nosurf
15	#ifdef MESOSCALE
16	use turb_mod, only : ustar
17	use comm_wrf, only : comm_LATENT_HF
18	#endif
19	implicit none
20
21	!==================================================================
22	!
23	! Purpose
24	! -------
25	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
26	!
27	! Implicit scheme
28	! We start by adding to variables x the physical tendencies
29	! already computed. We resolve the equation:
30	!
31	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
32	!
33	! Authors
34	! -------
35	! F. Hourdin, F. Forget, R. Fournier (199X)
36	! R. Wordsworth, B. Charnay (2010)
37	! J. Leconte (2012): To f90
38	! - Rewritten the diffusion scheme to conserve total enthalpy
39	! by accounting for dissipation of turbulent kinetic energy.
40	! - Accounting for lost mean flow kinetic energy should come soon.
41	!
42	!==================================================================
43
44	!-----------------------------------------------------------------------
45	! declarations
46	! ------------
47
48	! arguments
49	! ---------
50	INTEGER,INTENT(IN) :: ngrid
51	INTEGER,INTENT(IN) :: nlay
52	REAL,INTENT(IN) :: ptimestep
53	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
54	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
55	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
56	REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
57	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
58	REAL,INTENT(IN) :: pemis(ngrid)
59	REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
60	REAL,INTENT(IN) :: pfluxsrf(ngrid)
61	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
62	REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
63	REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
64	REAL,INTENT(OUT) :: sensibFlux(ngrid)
65	REAL,INTENT(IN) :: pcapcal(ngrid)
66	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
67	REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
68	REAL,INTENT(IN) :: rnat(ngrid)
69	LOGICAL,INTENT(IN) :: lastcall ! not used
70
71	! Arguments added for condensation
72	logical,intent(in) :: lecrit ! not used.
73	REAL,INTENT(IN) :: pz0
74
75	! Tracers
76	! --------
77	integer,intent(in) :: nq
78	real,intent(in) :: pqsurf(ngrid,nq)
79	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
80	real,intent(out) :: pdqdif(ngrid,nlay,nq)
81	real,intent(out) :: pdqsdif(ngrid,nq)
82
83	! local
84	! -----
85	integer ilev,ig,ilay,nlev
86
87	REAL z4st,zdplanck(ngrid)
88	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
89	REAL zcdv(ngrid),zcdh(ngrid)
90	REAL zcdv_true(ngrid),zcdh_true(ngrid)
91	REAL zu(ngrid,nlay),zv(ngrid,nlay)
92	REAL zh(ngrid,nlay),zt(ngrid,nlay)
93	REAL ztsrf(ngrid)
94	REAL z1(ngrid),z2(ngrid)
95	REAL zmass(ngrid,nlay)
96	REAL zfluxv(ngrid,nlay),zfluxt(ngrid,nlay),zfluxq(ngrid,nlay)
97	REAL zb0(ngrid,nlay)
98	REAL zExner(ngrid,nlay),zovExner(ngrid,nlay)
99	REAL zcv(ngrid,nlay),zdv(ngrid,nlay) !inversion coefficient for winds
100	REAL zct(ngrid,nlay),zdt(ngrid,nlay) !inversion coefficient for temperature
101	REAL zcq(ngrid,nlay),zdq(ngrid,nlay) !inversion coefficient for tracers
102	REAL zcst1
103	REAL zu2!, a
104	REAL zcq0(ngrid),zdq0(ngrid)
105	REAL zx_alf1(ngrid),zx_alf2(ngrid)
106
107	LOGICAL,SAVE :: firstcall=.true.
108	!$OMP THREADPRIVATE(firstcall)
109
110	! Tracers
111	! -------
112	INTEGER iq
113	REAL zq(ngrid,nlay,nq)
114	REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
115	REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
116	REAL zdmassevap(ngrid)
117	REAL rho(ngrid) ! near-surface air density
118	REAL qsat(ngrid),psat(ngrid)
119	REAL kmixmin
120
121	! Variables added for implicit latent heat inclusion
122	! --------------------------------------------------
123	real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
124
125	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
126	!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
127
128	real, parameter :: karman=0.4
129	real cd0, roughratio
130
131	real dqsdif_total(ngrid)
132	real zq0(ngrid)
133
134
135	! Coherence test
136	! --------------
137
138	IF (firstcall) THEN
139
140	if(water)then
141	ivap=igcm_h2o_vap
142	iliq=igcm_h2o_ice
143	iliq_surf=igcm_h2o_vap
144	iice_surf=igcm_h2o_ice ! simply to make the code legible
145	! to be generalised
146	else
147	ivap=0
148	iliq=0
149	iliq_surf=0
150	iice_surf=0 ! simply to make the code legible
151	endif
152	sensibFlux(:)=0.
153
154	firstcall=.false.
155	ENDIF
156
157	!-----------------------------------------------------------------------
158	! 1. Initialisation
159	! -----------------
160
161	nlev=nlay+1
162
163	! Calculate rhodz, (P/Ps)(R/cp) and dtrho/dz=dtrho*2 g/dp
164	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
165	! ---------------------------------------------
166
167	DO ilay=1,nlay
168	DO ig=1,ngrid
169	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
170	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
171	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
172	ENDDO
173	ENDDO
174
175	zcst1=4.gptimestep/(R*R)
176	DO ilev=2,nlev-1
177	DO ig=1,ngrid
178	zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev))
179	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
180	ENDDO
181	ENDDO
182	DO ig=1,ngrid
183	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
184	ENDDO
185	dqsdif_total(:)=0.0
186
187	!-----------------------------------------------------------------------
188	! 2. Add the physical tendencies computed so far
189	! ----------------------------------------------
190
191	DO ilev=1,nlay
192	DO ig=1,ngrid
193	zu(ig,ilev)=pu(ig,ilev)
194	zv(ig,ilev)=pv(ig,ilev)
195	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
196	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
197	ENDDO
198	ENDDO
199	if(tracer) then
200	DO iq =1, nq
201	DO ilev=1,nlay
202	DO ig=1,ngrid
203	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
204	ENDDO
205	ENDDO
206	ENDDO
207	if (water) then
208	DO ilev=1,nlay
209	DO ig=1,ngrid
210	zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
211	ENDDO
212	ENDDO
213	Endif
214	end if
215
216	!-----------------------------------------------------------------------
217	! 3. Turbulence scheme
218	! --------------------
219	!
220	! Source of turbulent kinetic energy at the surface
221	! -------------------------------------------------
222	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
223
224	DO ig=1,ngrid
225	roughratio = 1. + pzlay(ig,1)/pz0
226	cd0 = karman/log(roughratio)
227	cd0 = cd0*cd0
228	zcdv_true(ig) = cd0
229	zcdh_true(ig) = cd0
230	if(nosurf)then
231	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
232	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
233	endif
234	ENDDO
235
236	DO ig=1,ngrid
237	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
238	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
239	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
240	ENDDO
241
242	! Turbulent diffusion coefficients in the boundary layer
243	! ------------------------------------------------------
244
245	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
246
247	! Adding eddy mixing to mimic 3D general circulation in 1D
248	! R. Wordsworth & F. Forget (2010)
249	if ((ngrid.eq.1)) then
250	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
251	do ilev=1,nlay
252	do ig=1,ngrid
253	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
254	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
255	end do
256	end do
257	end if
258
259	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
260	DO ig=1,ngrid
261	zkv(ig,1)=zcdv(ig)
262	ENDDO
263	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
264	!JL12 calculate the flux coefficients (tables multiplied elements by elements)
265	zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
266
267	!-----------------------------------------------------------------------
268	! 4. Implicit inversion of u
269	! --------------------------
270
271	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
272	! avec
273	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
274	! et
275	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
276	! donc les entrees sont /zcdv/ pour la condition a la limite sol
277	! et /zkv/ = Ku
278
279	DO ig=1,ngrid
280	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
281	zcv(ig,nlay)=zmass(ig,nlay)zu(ig,nlay)z1(ig)
282	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
283	ENDDO
284
285	DO ilay=nlay-1,1,-1
286	DO ig=1,ngrid
287	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay) + zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
288	zcv(ig,ilay)=(zmass(ig,ilay)zu(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
289	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
290	ENDDO
291	ENDDO
292
293	DO ig=1,ngrid
294	zu(ig,1)=zcv(ig,1)
295	ENDDO
296	DO ilay=2,nlay
297	DO ig=1,ngrid
298	zu(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zu(ig,ilay-1)
299	ENDDO
300	ENDDO
301
302	!-----------------------------------------------------------------------
303	! 5. Implicit inversion of v
304	! --------------------------
305
306	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
307	! avec
308	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
309	! et
310	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
311	! donc les entrees sont /zcdv/ pour la condition a la limite sol
312	! et /zkv/ = Kv
313
314	DO ig=1,ngrid
315	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
316	zcv(ig,nlay)=zmass(ig,nlay)zv(ig,nlay)z1(ig)
317	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
318	ENDDO
319
320	DO ilay=nlay-1,1,-1
321	DO ig=1,ngrid
322	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay)+zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
323	zcv(ig,ilay)=(zmass(ig,ilay)zv(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
324	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
325	ENDDO
326	ENDDO
327
328	DO ig=1,ngrid
329	zv(ig,1)=zcv(ig,1)
330	ENDDO
331	DO ilay=2,nlay
332	DO ig=1,ngrid
333	zv(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zv(ig,ilay-1)
334	ENDDO
335	ENDDO
336
337	! Calcul of wind stress
338
339	DO ig=1,ngrid
340	flux_u(ig) = zfluxv(ig,1)/ptimestep*zu(ig,1)
341	flux_v(ig) = zfluxv(ig,1)/ptimestep*zv(ig,1)
342	ENDDO
343
344
345	!----------------------------------------------------------------------------
346	! 6. Implicit inversion of h, not forgetting the coupling with the ground
347
348	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
349	! avec
350	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
351	! et
352	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
353	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
354	! et /zkh/ = Kh
355
356	! Using the wind modified by friction for lifting and sublimation
357	! ---------------------------------------------------------------
358	DO ig=1,ngrid
359	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
360	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
361	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
362	zkh(ig,1)= zcdh(ig)
363	ustar(ig)=sqrt(zcdv_true(ig))*sqrt(zu2)
364	ENDDO
365
366
367	! JL12 calculate the flux coefficients (tables multiplied elements by elements)
368	! ---------------------------------------------------------------
369	zfluxq(1:ngrid,1:nlay)=zkh(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay) !JL12 we save zfluxq which doesn't need the Exner factor
370	zfluxt(1:ngrid,1:nlay)=zfluxq(1:ngrid,1:nlay)*zExner(1:ngrid,1:nlay)
371
372	DO ig=1,ngrid
373	z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay))
374	zct(ig,nlay)=zmass(ig,nlay)zt(ig,nlay)z1(ig)
375	zdt(ig,nlay)=zfluxt(ig,nlay)zovExner(ig,nlay-1)z1(ig)
376	ENDDO
377
378	DO ilay=nlay-1,2,-1
379	DO ig=1,ngrid
380	z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ &
381	zfluxt(ig,ilay+1)(zovExner(ig,ilay)-zdt(ig,ilay+1)zovExner(ig,ilay+1)))
382	zct(ig,ilay)=(zmass(ig,ilay)zt(ig,ilay)+zfluxt(ig,ilay+1)zct(ig,ilay+1)zovExner(ig,ilay+1))z1(ig)
383	zdt(ig,ilay)=zfluxt(ig,ilay)z1(ig)zovExner(ig,ilay-1)
384	ENDDO
385	ENDDO
386
387	!JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there
388	DO ig=1,ngrid
389	z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ &
390	zfluxt(ig,2)(zovExner(ig,1)-zdt(ig,2)zovExner(ig,2)))
391	zct(ig,1)=(zmass(ig,1)zt(ig,1)+zfluxt(ig,2)zct(ig,2)zovExner(ig,2))z1(ig)
392	zdt(ig,1)=zfluxt(ig,1)*z1(ig)
393	ENDDO
394
395
396	! Calculate (d Planck / dT) at the interface temperature
397	! ------------------------------------------------------
398
399	z4st=4.0sigmaptimestep
400	DO ig=1,ngrid
401	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
402	ENDDO
403
404	! Calculate temperature tendency at the interface (dry case)
405	! ----------------------------------------------------------
406	! Sum of fluxes at interface at time t + \delta t gives change in T:
407	! radiative fluxes
408	! turbulent convective (sensible) heat flux
409	! flux (if any) from subsurface
410
411	if(.not.water) then
412
413	DO ig=1,ngrid
414	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
415	+ pfluxsrf(ig)ptimestep + zdplanck(ig)ptsrf(ig)
416	z2(ig) = pcapcal(ig)+zdplanck(ig)+cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1))
417	ztsrf(ig) = z1(ig) / z2(ig)
418	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
419	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
420	ENDDO
421	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
422
423
424	! Recalculate temperature to top of atmosphere, starting from ground
425	! ------------------------------------------------------------------
426
427	DO ilay=2,nlay
428	DO ig=1,ngrid
429	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
430	ENDDO
431	ENDDO
432
433	endif ! not water
434
435	!-----------------------------------------------------------------------
436	! TRACERS (no vapour)
437	! -------
438
439	if(tracer) then
440
441	! Calculate vertical flux from the bottom to the first layer (dust)
442	! -----------------------------------------------------------------
443	do ig=1,ngrid
444	rho(ig) = zb0(ig,1) /ptimestep
445	end do
446
447	pdqsdif(:,:)=0.
448
449	! Implicit inversion of q
450	! -----------------------
451	do iq=1,nq
452
453	if (iq.ne.ivap) then
454
455	DO ig=1,ngrid
456	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
457	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
458	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
459	ENDDO
460
461	DO ilay=nlay-1,2,-1
462	DO ig=1,ngrid
463	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
464	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
465	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
466	ENDDO
467	ENDDO
468
469	if ((water).and.(iq.eq.iliq)) then
470	! special case for condensed water tracer: do not include
471	! h2o ice tracer from surface (which is set when handling
472	! h2o vapour case (see further down).
473	! zb(ig,1)=0 if iq ne ivap
474	DO ig=1,ngrid
475	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
476	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
477	ENDDO
478	else ! general case
479	do ig=1,ngrid
480	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
481	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2)+(-pdqsdif(ig,iq))ptimestep)z1(ig)
482	! tracer flux from surface
483	! currently pdqsdif always zero here,
484	! so last line is superfluous
485	enddo
486	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
487
488
489	! Starting upward calculations for simple tracer mixing (e.g., dust)
490	do ig=1,ngrid
491	zq(ig,1,iq)=zcq(ig,1)
492	end do
493
494	do ilay=2,nlay
495	do ig=1,ngrid
496	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
497	end do
498	end do
499
500	endif ! if (iq.ne.ivap)
501
502	! Calculate temperature tendency including latent heat term
503	! and assuming an infinite source of water on the ground
504	! ------------------------------------------------------------------
505
506	if (water.and.(iq.eq.ivap)) then
507
508	! compute evaporation efficiency
509	do ig=1,ngrid
510	if(nint(rnat(ig)).eq.1)then
511	dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf)
512	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
513	dryness(ig)=MAX(0.,dryness(ig))
514	endif
515	enddo
516
517	do ig=1,ngrid
518	! Calculate the value of qsat at the surface (water)
519	call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
520	call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
521	call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
522	dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
523	! calculate dQsat / dT by finite differences
524	! we cannot use the updated temperature value yet...
525	enddo
526
527	! coefficients for q
528
529	do ig=1,ngrid
530	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
531	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
532	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
533	enddo
534
535	do ilay=nlay-1,2,-1
536	do ig=1,ngrid
537	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
538	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
539	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
540	enddo
541	enddo
542
543	do ig=1,ngrid
544	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,1)dryness(ig)+zfluxq(ig,2)(1.-zdq(ig,2)))
545	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
546	zdq(ig,1)=dryness(ig)zfluxq(ig,1)z1(ig)
547	enddo
548
549	do ig=1,ngrid
550	!calculation of surface temperature
551	zdq0(ig) = dqsat(ig)
552	zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
553
554	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
555	+ zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep &
556	+ zfluxq(ig,1)dryness(ig)RLVTT((zdq(ig,1)-1.0)zcq0(ig)+zcq(ig,1))
557
558	z2(ig) = pcapcal(ig) + cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1)) &
559	+ zdplanck(ig)+zfluxq(ig,1)dryness(ig)RLVTTzdq0(ig)(1.0-zdq(ig,1))
560
561	ztsrf(ig) = z1(ig) / z2(ig)
562
563	! calculation of qs and q1
564	zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf(ig)
565	zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
566
567	! calculation of evaporation
568	dqsdif_total(ig)=zfluxq(ig,1)dryness(ig)(zq(ig,1,ivap)-zq0(ig))
569
570	! --------------------------------------------------------
571	! Now check if we've taken too much water from the surface
572	! This can only occur on the continent
573	! If we do, we recompute Tsurf, T1 and q1 accordingly
574	if((-dqsdif_total(ig).gt.(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))).and.rnat(ig).eq.1)then
575	!water flux * ptimestep
576	dqsdif_total(ig)=-(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))
577
578	!recompute surface temperature
579	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzfluxq(ig,1)zct(ig,1)zovExner(ig,1) &
580	+ zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep &
581	+ RLVTT*dqsdif_total(ig)
582	z2(ig) = pcapcal(ig) + cppzfluxq(ig,1)(1.-zovExner(ig,1)*zdt(ig,1)) &
583	+ zdplanck(ig)
584	ztsrf(ig) = z1(ig) / z2(ig)
585
586	!recompute q1 with new water flux from surface
587	zq(ig,1,iq) = (zmass(ig,1)(pq(ig,1,iq)+ptimesteppdqfi(ig,1,iq)) &
588	+zfluxq(ig,2)*zcq(ig,2)-dqsdif_total(ig)) &
589	/ (zmass(ig,1)+(1.-zdq(ig,2))*zfluxq(ig,2))
590	end if
591
592	! calculation surface T tendency and T(1)
593	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
594	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
595	enddo
596
597
598	! recalculate temperature and q(vap) to top of atmosphere, starting from ground
599	do ilay=2,nlay
600	do ig=1,ngrid
601	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
602	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
603	end do
604	end do
605
606
607	do ig=1,ngrid
608	! --------------------------------------------------------------------------
609	! On the ocean, if T > 0 C then the vapour tendency must replace the ice one
610	! The surface vapour tracer is actually liquid. To make things difficult.
611
612	if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
613
614	pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
615	pdqsdif(ig,iice_surf)=0.0
616
617	elseif (nint(rnat(ig)).eq.1) then ! (continent)
618	! If water is evaporating / subliming, we take it from ice before liquid
619	! -- is this valid??
620	if(dqsdif_total(ig).lt.0)then
621	if (-dqsdif_total(ig).gt.pqsurf(ig,iice_surf))then
622	pdqsdif(ig,iice_surf) = -pqsurf(ig,iice_surf)/ptimestep ! removes all the ice!
623	pdqsdif(ig,iliq_surf) = dqsdif_total(ig)/ptimestep- pdqsdif(ig,iice_surf) ! take the remainder from the liquid instead
624	else
625	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
626	pdqsdif(ig,iliq_surf)=0.
627	end if
628	else !dqsdif_total(ig).ge.0
629	!If water vapour is condensing, we must decide whether it forms ice or liquid.
630	if(ztsrf(ig).gt.T_h2O_ice_liq)then
631	pdqsdif(ig,iice_surf)=0.0
632	pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
633	else
634	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
635	pdqsdif(ig,iliq_surf)=0.0
636	endif
637	endif
638
639	elseif (nint(rnat(ig)).eq.2) then ! (continental glaciers)
640	pdqsdif(ig,iliq_surf)=0.0
641	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
642
643	endif !rnat
644	end do ! of DO ig=1,ngrid
645
646	endif ! if (water et iq=ivap)
647	end do ! of do iq=1,nq
648
649	if (water) then ! special case where we recompute water mixing without any evaporation.
650	! The difference with the first calculation then tells us where evaporated water has gone
651
652	DO ig=1,ngrid
653	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
654	zcq(ig,nlay)=zmass(ig,nlay)zqnoevap(ig,nlay)z1(ig)
655	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
656	ENDDO
657
658	DO ilay=nlay-1,2,-1
659	DO ig=1,ngrid
660	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
661	zcq(ig,ilay)=(zmass(ig,ilay)zqnoevap(ig,ilay)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
662	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
663	ENDDO
664	ENDDO
665
666	do ig=1,ngrid
667	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
668	zcq(ig,1)=(zmass(ig,1)zqnoevap(ig,1)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
669	enddo
670
671	! Starting upward calculations for simple tracer mixing (e.g., dust)
672	do ig=1,ngrid
673	zqnoevap(ig,1)=zcq(ig,1)
674	end do
675
676	do ilay=2,nlay
677	do ig=1,ngrid
678	zqnoevap(ig,ilay)=zcq(ig,ilay)+zdq(ig,ilay)*zqnoevap(ig,ilay-1)
679	end do
680	end do
681
682	endif ! if water
683
684
685	endif ! tracer
686
687
688	!-----------------------------------------------------------------------
689	! 8. Final calculation of the vertical diffusion tendencies
690	! -----------------------------------------------------------------
691
692	do ilev = 1, nlay
693	do ig=1,ngrid
694	pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)))/ptimestep
695	pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)))/ptimestep
696	pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
697	enddo
698	enddo
699
700	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
701	sensibFlux(ig)=cppzfluxt(ig,1)/ptimestep(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
702	ENDDO
703
704	if (tracer) then
705	do iq = 1, nq
706	do ilev = 1, nlay
707	do ig=1,ngrid
708	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-(pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
709	enddo
710	enddo
711	enddo
712	if (water) then
713	do ilev = 1, nlay
714	do ig=1,ngrid
715	pdqevap(ig,ilev)=(zq(ig,ilev,ivap)-zqnoevap(ig,ilev))/ptimestep
716	enddo
717	enddo
718	do ig=1,ngrid
719	zdmassevap(ig)=SUM(pdqevap(ig,:)zmass(ig,:))ptimestep
720	end do
721	endif
722	endif
723
724	if(water)then
725	#ifndef MESOSCALE
726	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
727	#endif
728	if (tracer) then
729	#ifndef MESOSCALE
730	call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
731	#else
732	comm_LATENT_HF(:)=0.0
733	comm_LATENT_HF(1:ngrid)=RLVTT*dqsdif_total(1:ngrid)/ptimestep
734	#endif
735	endif
736	endif
737
738	! if(lastcall)then
739	! if(ngrid.eq.1)then
740	! print*,'Saving k.out...'
741	! OPEN(12,file='k.out',form='formatted')
742	! DO ilay=1,nlay
743	! write(12,*) zkh(1,ilay), pplay(1,ilay)
744	! ENDDO
745	! CLOSE(12)
746	! endif
747	! endif
748
749	end

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

Download in other formats:

Original Format