Context Navigation

vdifc.F @ 601

Last change on this file since 601 was 600, checked in by jleconte, 13 years ago

Corrects the computation of planck function at the surface in sfluxi

so that its integral is equal to sigma Tsurf^4.

This ensure that no flux is lost due to:

-truncation of the planck function at high/low wavenumber
-numerical error during first spectral computation of the planck function
-discrepancy between Tsurf and NTS/NTfac in sfluxi

OLR now equal to LW net heating/cooling at equilibrium!
As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the

corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.

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

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

Download in other formats:

Original Format