Context Navigation

vdifc.F @ 1448

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