Context Navigation

vdifc_mod.F @ 3576

Last change on this file since 3576 was 3236, checked in by gmilcareck, 21 months ago
Add virtual correction for convective adjustment and turbulent diffusion (turbdiff and vdifc) + correction of allocated tables in physiq_mod for varspec.
File size: 23.8 KB

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

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

Download in other formats:

Original Format