Context Navigation

turbdiff_mod.F90 @ 3436

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

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

Download in other formats:

Original Format