Context Navigation

vdifc.F @ 1145

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

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

Download in other formats:

Original Format