Context Navigation

turbdiff.F90 @ 263

Last change on this file since 263 was 227, checked in by milmd, 10 years ago
Last LMDZ version (1315) with OpenMP directives and other stuff
File size: 27.7 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	pdufi,pdvfi,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
14	implicit none
15
16	!==================================================================
17	!
18	! Purpose
19	! -------
20	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
21	!
22	! Implicit scheme
23	! We start by adding to variables x the physical tendencies
24	! already computed. We resolve the equation:
25	!
26	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
27	!
28	! Authors
29	! -------
30	! F. Hourdin, F. Forget, R. Fournier (199X)
31	! R. Wordsworth, B. Charnay (2010)
32	! J. Leconte (2012): To f90
33	! - Rewritten the diffusion scheme to conserve total enthalpy
34	! by accounting for dissipation of turbulent kinetic energy.
35	! - Accounting for lost mean flow kinetic energy should come soon.
36	!
37	!==================================================================
38
39	!-----------------------------------------------------------------------
40	! declarations
41	! ------------
42
43	! include "dimensions.h"
44	! include "dimphys.h"
45	include "comcstfi.h"
46	include "callkeys.h"
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) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
60	REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
61	REAL,INTENT(IN) :: pfluxsrf(ngrid)
62	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
63	REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
64	REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
65	REAL,INTENT(OUT) :: sensibFlux(ngrid)
66	REAL,INTENT(IN) :: pcapcal(ngrid)
67	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
68	REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
69	REAL,INTENT(IN) :: rnat(ngrid)
70	LOGICAL,INTENT(IN) :: lastcall ! not used
71
72	! Arguments added for condensation
73	logical,intent(in) :: lecrit ! not used.
74	REAL,INTENT(IN) :: pz0
75
76	! Tracers
77	! --------
78	integer,intent(in) :: nq
79	real,intent(in) :: pqsurf(ngrid,nq)
80	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
81	real,intent(out) :: pdqdif(ngrid,nlay,nq)
82	real,intent(out) :: pdqsdif(ngrid,nq)
83
84	! local
85	! -----
86	integer ilev,ig,ilay,nlev
87
88	REAL z4st,zdplanck(ngrid)
89	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
90	REAL zcdv(ngrid),zcdh(ngrid)
91	REAL zcdv_true(ngrid),zcdh_true(ngrid)
92	REAL zu(ngrid,nlay),zv(ngrid,nlay)
93	REAL zh(ngrid,nlay),zt(ngrid,nlay)
94	REAL ztsrf(ngrid)
95	REAL z1(ngrid),z2(ngrid)
96	REAL zmass(ngrid,nlay)
97	REAL zfluxv(ngrid,nlay),zfluxt(ngrid,nlay),zfluxq(ngrid,nlay)
98	REAL zb0(ngrid,nlay)
99	REAL zExner(ngrid,nlay),zovExner(ngrid,nlay)
100	REAL zcv(ngrid,nlay),zdv(ngrid,nlay) !inversion coefficient for winds
101	REAL zct(ngrid,nlay),zdt(ngrid,nlay) !inversion coefficient for temperature
102	REAL zcq(ngrid,nlay),zdq(ngrid,nlay) !inversion coefficient for tracers
103	REAL zcst1
104	REAL zu2!, a
105	REAL zcq0(ngrid),zdq0(ngrid)
106	REAL zx_alf1(ngrid),zx_alf2(ngrid)
107
108	LOGICAL,SAVE :: firstcall=.true.
109	!$OMP THREADPRIVATE(firstcall)
110
111	! Tracers
112	! -------
113	INTEGER iq
114	REAL zq(ngrid,nlay,nq)
115	REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
116	REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
117	REAL zdmassevap(ngrid)
118	REAL rho(ngrid) ! near-surface air density
119	REAL qsat(ngrid),psat(ngrid)
120	REAL kmixmin
121
122	! Variables added for implicit latent heat inclusion
123	! --------------------------------------------------
124	real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
125
126	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
127	!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
128
129	real, parameter :: karman=0.4
130	real cd0, roughratio
131
132	real dqsdif_total(ngrid)
133	real zq0(ngrid)
134
135
136	! Coherence test
137	! --------------
138
139	IF (firstcall) THEN
140
141	if(water)then
142	ivap=igcm_h2o_vap
143	iliq=igcm_h2o_ice
144	iliq_surf=igcm_h2o_vap
145	iice_surf=igcm_h2o_ice ! simply to make the code legible
146	! to be generalised
147	else
148	ivap=0
149	iliq=0
150	iliq_surf=0
151	iice_surf=0 ! simply to make the code legible
152	endif
153	sensibFlux(:)=0.
154
155	firstcall=.false.
156	ENDIF
157
158	!-----------------------------------------------------------------------
159	! 1. Initialisation
160	! -----------------
161
162	nlev=nlay+1
163
164	! Calculate rhodz, (P/Ps)(R/cp) and dtrho/dz=dtrho*2 g/dp
165	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
166	! ---------------------------------------------
167
168	DO ilay=1,nlay
169	DO ig=1,ngrid
170	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
171	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
172	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
173	ENDDO
174	ENDDO
175
176	zcst1=4.gptimestep/(R*R)
177	DO ilev=2,nlev-1
178	DO ig=1,ngrid
179	zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev))
180	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
181	ENDDO
182	ENDDO
183	DO ig=1,ngrid
184	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
185	ENDDO
186	dqsdif_total(:)=0.0
187
188	!-----------------------------------------------------------------------
189	! 2. Add the physical tendencies computed so far
190	! ----------------------------------------------
191
192	DO ilev=1,nlay
193	DO ig=1,ngrid
194	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
195	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
196	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
197	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
198	ENDDO
199	ENDDO
200	if(tracer) then
201	DO iq =1, nq
202	DO ilev=1,nlay
203	DO ig=1,ngrid
204	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
205	ENDDO
206	ENDDO
207	ENDDO
208	if (water) then
209	DO ilev=1,nlay
210	DO ig=1,ngrid
211	zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
212	ENDDO
213	ENDDO
214	Endif
215	end if
216
217	!-----------------------------------------------------------------------
218	! 3. Turbulence scheme
219	! --------------------
220	!
221	! Source of turbulent kinetic energy at the surface
222	! -------------------------------------------------
223	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
224
225	DO ig=1,ngrid
226	roughratio = 1. + pzlay(ig,1)/pz0
227	cd0 = karman/log(roughratio)
228	cd0 = cd0*cd0
229	zcdv_true(ig) = cd0
230	zcdh_true(ig) = cd0
231	if(nosurf)then
232	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
233	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
234	endif
235	ENDDO
236
237	DO ig=1,ngrid
238	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
239	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
240	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
241	ENDDO
242
243	! Turbulent diffusion coefficients in the boundary layer
244	! ------------------------------------------------------
245
246	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
247
248	! Adding eddy mixing to mimic 3D general circulation in 1D
249	! R. Wordsworth & F. Forget (2010)
250	if ((ngrid.eq.1)) then
251	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
252	do ilev=1,nlay
253	do ig=1,ngrid
254	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
255	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
256	end do
257	end do
258	end if
259
260	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
261	DO ig=1,ngrid
262	zkv(ig,1)=zcdv(ig)
263	ENDDO
264	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
265
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	ENDDO
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)+pdufi(ig,ilev)*ptimestep))/ptimestep
695	pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/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	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
726	if (tracer) then
727	call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
728	endif
729	endif
730
731	! if(lastcall)then
732	! if(ngrid.eq.1)then
733	! print*,'Saving k.out...'
734	! OPEN(12,file='k.out',form='formatted')
735	! DO ilay=1,nlay
736	! write(12,*) zkh(1,ilay), pplay(1,ilay)
737	! ENDDO
738	! CLOSE(12)
739	! endif
740	! endif
741
742	end

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

Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/turbdiff.F90 @ 263

Download in other formats: