Context Navigation

vdifc.F @ 1834

Last change on this file since 1834 was 1542, checked in by emillour, 9 years ago

Generic GCM:

fix for 1D in writediagfi to enable writing at "ecritphy" rate.
move iniprint.h to "misc"
Some code cleanup in anticipation of future updates:
- changed variable names in comgeomphy.F90: give them more explicit names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
- removed long(), lati() and area() from comgeomfi_h.F90, use longitude(), latitude() and cell_are() from comgeomphy.F90 instead

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

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

Download in other formats:

Original Format