Context Navigation

turbdiff.F90 @ 726

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

17/07/2012 == JL for LK

Generalization of aerosol scheme:
- any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
- addition of a module with the id numbers for aerosols (aerosol_mod.F90).
- initialization of aerosols id numbers in iniaerosol.F90
- compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

If starting from an old start file, recreate start file with the q=0 option in newstart.e.
update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

To add other aerosols, see Laura Kerber.

File size: 24.6 KB

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

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

Download in other formats:

Original Format