Context Navigation

vdifc.F @ 586

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

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

Download in other formats:

Original Format