Context Navigation

vdifc.F @ 2176

Last change on this file since 2176 was 1993, checked in by jleconte, 6 years ago

29/08/2018 == JL

-watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...) which used Psat_water from watercommon.
This is now harmonized. ALl routines use Psat_water. Watersat.F has been removed, but the routine is now in watercommon for archival purpose. It is not used anymore.
-also changed the number of chars for tname in the dyn3D/infotrac.F90 to be able to run rcm1d.

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

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

Download in other formats:

Original Format