Context Navigation

vdifc.F @ 486

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

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

Download in other formats:

Original Format