Context Navigation

turbdiff_mod.F90 @ 2447

Last change on this file since 2447 was 2427, checked in by emillour, 5 years ago
Generic GCM: Bug fix on call arguments sent from physiq to vdifc (probably not as bad as it sounds, as turbdiff is usually used instead of vdifc). In the process turned vdifc into a module, as well as turbdiff, for better control, and removed unused arguments. EM+YJ
File size: 27.5 KB

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

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

Download in other formats:

Original Format