Context Navigation

vdifc.F @ 1351

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

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

Download in other formats:

Original Format