Context Navigation

turbdiff.F90 @ 730

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

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

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

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

Download in other formats:

Original Format