Context Navigation

turbdiff.F90 @ 601

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

Corrects the computation of planck function at the surface in sfluxi

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

This ensure that no flux is lost due to:

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

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

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

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

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

Download in other formats:

Original Format