Context Navigation

turbdiff.F90 @ 832

Last change on this file since 832 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

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

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

Download in other formats:

Original Format