Context Navigation

turbdiff_mod.F90 @ 4037

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

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

Download in other formats:

Original Format