Context Navigation

vdifc.F @ 269

Last change on this file since 269 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

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

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

Download in other formats:

Original Format