Context Navigation

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

Last change on this file since 222 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

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

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

Download in other formats: