Context Navigation

vdifc.F @ 264

Last change on this file since 264 was 227, checked in by milmd, 10 years ago
Last LMDZ version (1315) with OpenMP directives and other stuff
File size: 24.3 KB

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

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

Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/vdifc.F @ 264

Download in other formats: