Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/vdifc.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 30.1 KB

Line
1	subroutine vdifc(ngrid,nlay,nq,rnat,ppopsk,
2	& ptimestep,pcapcal,lecrit,
3	& pplay,pplev,pzlay,pzlev,pz0,
4	& pu,pv,ph,pq,ptsrf,pemis,pqsurf,
5	& pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
6	& pdudif,pdvdif,pdhdif,pdtsrf,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	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	!
34	!==================================================================
35
36	!-----------------------------------------------------------------------
37	! declarations
38	! ------------
39
40	!-----------------------------------------------------------------------
41	! INCLUDE 'dimensions.h'
42	!
43	! dimensions.h contient les dimensions du modele
44	! ndm est tel que iim=2**ndm
45	!-----------------------------------------------------------------------
46
47	INTEGER iim,jjm,llm,ndm
48
49	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
50
51	!-----------------------------------------------------------------------
52	!-----------------------------------------------------------------------
53	! INCLUDE 'dimphys.h'
54
55	! ngridmx : number of horizontal grid points
56	! note: the -1/jjm term will be 0; unless jj=1
57	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
58	! nlayermx : number of atmospheric layers
59	integer, parameter :: nlayermx = llm
60	! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
61	!integer, parameter :: nsoilmx = 4 ! for a test
62	!integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
63	!integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
64	!-----------------------------------------------------------------------
65
66	!-----------------------------------------------------------------------
67	! INCLUDE "comcstfi.h"
68
69	common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
70	common/comcstfi/avocado!,molrad,visc
71
72	real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
73	real avocado!,molrad,visc
74
75	!
76	! For Fortran 77/Fortran 90 compliance always use line continuation
77	! symbols '&' in columns 73 and 6
78	!
79	! Group commons according to their type for minimal performance impact
80
81	COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj &
82	& , co2cond,callsoil &
83	& , season,diurnal,tlocked,rings_shadow,lwrite &
84	& , callstats,calleofdump &
85	& , enertest &
86	& , callgasvis,continuum,H2Ocont_simple,graybody &
87	& , radfixed &
88	& , meanOLR, specOLR &
89	& , kastprof &
90	& , nosurf, oblate &
91	& , newtonian, testradtimes &
92	& , check_cpp_match, force_cpp &
93	& , rayleigh &
94	& , stelbbody &
95	& , nearco2cond &
96	& , tracer, mass_redistrib, varactive, varfixed &
97	& , sedimentation,water,watercond,waterrain &
98	& , aeroco2,aeroh2o,aeroh2so4,aeroback2lay &
99	& , aerofixco2,aerofixh2o &
100	& , hydrology, sourceevol &
101	& , CLFvarying &
102	& , strictboundcorrk &
103	& , ok_slab_ocean &
104	& , ok_slab_sic &
105	& , ok_slab_heat_transp
106
107
108	COMMON/callkeys_i/iaervar,iddist,iradia,startype
109
110	COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb, &
111	& Tstrat,tplanet,obs_tau_col_tropo, &
112	& obs_tau_col_strato,pres_bottom_tropo, &
113	& pres_top_tropo,pres_bottom_strato, &
114	& pres_top_strato,size_tropo,size_strato,satval, &
115	& CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
116	& maxicethick,Tsaldiff,tau_relax,cloudlvl, &
117	& icetstep,intheat,flatten,Rmean,J2,MassPlanet
118
119	logical callrad,corrk,calldifv,UseTurbDiff &
120	& , calladj,co2cond,callsoil &
121	& , season,diurnal,tlocked,rings_shadow,lwrite &
122	& , callstats,calleofdump &
123	& , callgasvis,continuum,H2Ocont_simple,graybody &
124	& , strictboundcorrk
125
126	logical enertest
127	logical nonideal
128	logical meanOLR
129	logical specOLR
130	logical kastprof
131	logical newtonian
132	logical check_cpp_match
133	logical force_cpp
134	logical testradtimes
135	logical rayleigh
136	logical stelbbody
137	logical ozone
138	logical nearco2cond
139	logical tracer
140	logical mass_redistrib
141	logical varactive
142	logical varfixed
143	logical radfixed
144	logical sedimentation
145	logical water,watercond,waterrain
146	logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
147	logical aerofixco2,aerofixh2o
148	logical hydrology
149	logical sourceevol
150	logical CLFvarying
151	logical nosurf
152	logical oblate
153	logical ok_slab_ocean
154	logical ok_slab_sic
155	logical ok_slab_heat_transp
156
157	integer iddist
158	integer iaervar
159	integer iradia
160	integer startype
161
162	real topdustref
163	real Nmix_co2
164	real dusttau
165	real Fat1AU
166	real stelTbb
167	real Tstrat
168	real tplanet
169	real obs_tau_col_tropo
170	real obs_tau_col_strato
171	real pres_bottom_tropo
172	real pres_top_tropo
173	real pres_bottom_strato
174	real pres_top_strato
175	real size_tropo
176	real size_strato
177	real satval
178	real CLFfixval
179	real n2mixratio
180	real co2supsat
181	real pceil
182	real albedosnow
183	real maxicethick
184	real Tsaldiff
185	real tau_relax
186	real cloudlvl
187	real icetstep
188	real intheat
189	real flatten
190	real Rmean
191	real J2
192	real MassPlanet
193
194	! arguments
195	! ---------
196	INTEGER ngrid,nlay
197	REAL ptimestep
198	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
199	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
200	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
201	REAL ptsrf(ngrid),pemis(ngrid)
202	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
203	REAL pfluxsrf(ngrid)
204	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
205	REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
206	REAL pq2(ngrid,nlay+1)
207
208	real rnat(ngrid)
209
210	! Arguments added for condensation
211	REAL ppopsk(ngrid,nlay)
212	logical lecrit
213	REAL pz0
214
215	! Tracers
216	! --------
217	integer nq
218	real pqsurf(ngrid,nq)
219	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
220	real pdqdif(ngrid,nlay,nq)
221	real pdqsdif(ngrid,nq)
222
223	! local
224	! -----
225	integer ilev,ig,ilay,nlev
226
227	REAL z4st,zdplanck(ngrid)
228	REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
229	REAL zcdv(ngrid),zcdh(ngrid)
230	REAL zcdv_true(ngrid),zcdh_true(ngrid)
231	REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
232	REAL zh(ngrid,nlayermx)
233	REAL ztsrf2(ngrid)
234	REAL z1(ngrid),z2(ngrid)
235	REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
236	REAL zb0(ngrid,nlayermx)
237	REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
238	REAL zcst1
239	REAL zu2!, a
240	REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx)
241	REAL evap(ngrid)
242	REAL zcq0(ngrid),zdq0(ngrid)
243	REAL zx_alf1(ngrid),zx_alf2(ngrid)
244
245	LOGICAL firstcall
246	SAVE firstcall
247
248	LOGICAL lastcall
249
250	! variables added for CO2 condensation
251	! ------------------------------------
252	REAL hh !, zhcond(ngrid,nlayermx)
253	! REAL latcond,tcond1mb
254	! REAL acond,bcond
255	! SAVE acond,bcond
256	! DATA latcond,tcond1mb/5.9e5,136.27/
257
258	! Tracers
259	! -------
260	INTEGER iq
261	REAL zq(ngrid,nlayermx,nq)
262	REAL zq1temp(ngrid)
263	REAL rho(ngrid) ! near-surface air density
264	REAL qsat(ngrid)
265	DATA firstcall/.true./
266	REAL kmixmin
267
268	! Variables added for implicit latent heat inclusion
269	! --------------------------------------------------
270	real latconst, dqsat(ngrid), qsat_temp1, qsat_temp2
271	real z1_Tdry(ngrid), z2_Tdry(ngrid)
272	real z1_T(ngrid), z2_T(ngrid)
273	real zb_T(ngrid)
274	real zc_T(ngrid,nlayermx)
275	real zd_T(ngrid,nlayermx)
276	real lat1(ngrid), lat2(ngrid)
277	real surfh2otot
278	logical surffluxdiag
279	integer isub ! sub-loop for precision
280
281	integer ivap, iice ! also make liq for clarity on surface...
282	save ivap, iice
283
284	real, parameter :: karman=0.4
285	real cd0, roughratio
286
287	logical forceWC
288	real masse, Wtot, Wdiff
289
290	real dqsdif_total(ngrid)
291	real zq0(ngrid)
292
293	forceWC=.true.
294	! forceWC=.false.
295
296
297	! Coherence test
298	! --------------
299
300	IF (firstcall) THEN
301	! To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
302	! bcond=1./tcond1mb
303	! acond=r/latcond
304	! PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
305	! PRINT*,' acond,bcond',acond,bcond
306
307	if(water)then
308	! iliq=igcm_h2o_vap
309	ivap=igcm_h2o_vap
310	iice=igcm_h2o_ice ! simply to make the code legible
311	! to be generalised later
312	endif
313
314	firstcall=.false.
315	ENDIF
316
317	!-----------------------------------------------------------------------
318	! 1. Initialisation
319	! -----------------
320
321	nlev=nlay+1
322
323	! Calculate rhodz and dtrho/dz=dtrho*2 g/dp
324	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
325	! ---------------------------------------------
326
327	DO ilay=1,nlay
328	DO ig=1,ngrid
329	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
330	ENDDO
331	ENDDO
332
333	zcst1=4.gptimestep/(R*R)
334	DO ilev=2,nlev-1
335	DO ig=1,ngrid
336	zb0(ig,ilev)=pplev(ig,ilev)*
337	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
338	s (ph(ig,ilev-1)+ph(ig,ilev))
339	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
340	s (pplay(ig,ilev-1)-pplay(ig,ilev))
341	ENDDO
342	ENDDO
343	DO ig=1,ngrid
344	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
345	ENDDO
346
347	dqsdif_total(:)=0.0
348
349	!-----------------------------------------------------------------------
350	! 2. Add the physical tendencies computed so far
351	! ----------------------------------------------
352
353	DO ilev=1,nlay
354	DO ig=1,ngrid
355	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
356	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
357	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
358	ENDDO
359	ENDDO
360	if(tracer) then
361	DO iq =1, nq
362	DO ilev=1,nlay
363	DO ig=1,ngrid
364	zq(ig,ilev,iq)=pq(ig,ilev,iq) +
365	& pdqfi(ig,ilev,iq)*ptimestep
366	ENDDO
367	ENDDO
368	ENDDO
369	end if
370
371	!-----------------------------------------------------------------------
372	! 3. Turbulence scheme
373	! --------------------
374	!
375	! Source of turbulent kinetic energy at the surface
376	! -------------------------------------------------
377	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
378
379	DO ig=1,ngrid
380	roughratio = 1.E+0 + pzlay(ig,1)/pz0
381	cd0 = karman/log(roughratio)
382	cd0 = cd0*cd0
383	zcdv_true(ig) = cd0
384	zcdh_true(ig) = cd0
385	if (nosurf) then
386	zcdv_true(ig) = 0. !! disable sensible momentum flux
387	zcdh_true(ig) = 0. !! disable sensible heat flux
388	endif
389	ENDDO
390
391	DO ig=1,ngrid
392	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
393	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
394	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
395	ENDDO
396
397	! Turbulent diffusion coefficients in the boundary layer
398	! ------------------------------------------------------
399
400	call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay
401	& ,pu,pv,ph,zcdv_true
402	& ,pq2,zkv,zkh)
403
404	! Adding eddy mixing to mimic 3D general circulation in 1D
405	! R. Wordsworth & F. Forget (2010)
406	if ((ngrid.eq.1)) then
407	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
408	do ilev=1,nlay
409	do ig=1,ngrid
410	!zkh(ig,ilev) = 1.0
411	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
412	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
413	end do
414	end do
415	end if
416
417	!-----------------------------------------------------------------------
418	! 4. Implicit inversion of u
419	! --------------------------
420
421	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
422	! avec
423	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
424	! et
425	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
426	! donc les entrees sont /zcdv/ pour la condition a la limite sol
427	! et /zkv/ = Ku
428
429	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
430	CALL multipl(ngrid,zcdv,zb0,zb)
431
432	DO ig=1,ngrid
433	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
434	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
435	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
436	ENDDO
437
438	DO ilay=nlay-1,1,-1
439	DO ig=1,ngrid
440	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
441	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
442	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
443	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
444	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
445	ENDDO
446	ENDDO
447
448	DO ig=1,ngrid
449	zu(ig,1)=zc(ig,1)
450	ENDDO
451	DO ilay=2,nlay
452	DO ig=1,ngrid
453	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
454	ENDDO
455	ENDDO
456
457	!-----------------------------------------------------------------------
458	! 5. Implicit inversion of v
459	! --------------------------
460
461	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
462	! avec
463	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
464	! et
465	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
466	! donc les entrees sont /zcdv/ pour la condition a la limite sol
467	! et /zkv/ = Kv
468
469	DO ig=1,ngrid
470	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
471	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
472	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
473	ENDDO
474
475	DO ilay=nlay-1,1,-1
476	DO ig=1,ngrid
477	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
478	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
479	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
480	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
481	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
482	ENDDO
483	ENDDO
484
485	DO ig=1,ngrid
486	zv(ig,1)=zc(ig,1)
487	ENDDO
488	DO ilay=2,nlay
489	DO ig=1,ngrid
490	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
491	ENDDO
492	ENDDO
493
494	!----------------------------------------------------------------------------
495	! 6. Implicit inversion of h, not forgetting the coupling with the ground
496
497	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
498	! avec
499	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
500	! et
501	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
502	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
503	! et /zkh/ = Kh
504
505	! Using the wind modified by friction for lifting and sublimation
506	! ---------------------------------------------------------------
507	DO ig=1,ngrid
508	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
509	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
510	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
511	ENDDO
512
513	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
514	CALL multipl(ngrid,zcdh,zb0,zb)
515
516	DO ig=1,ngrid
517	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
518	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
519	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
520	ENDDO
521
522	DO ilay=nlay-1,2,-1
523	DO ig=1,ngrid
524	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
525	& zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
526	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
527	& zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
528	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
529	ENDDO
530	ENDDO
531
532	DO ig=1,ngrid
533	z1(ig)=1./(za(ig,1)+zb(ig,1)+
534	& zb(ig,2)*(1.-zd(ig,2)))
535	zc(ig,1)=(za(ig,1)*zh(ig,1)+
536	& zb(ig,2)zc(ig,2))z1(ig)
537	zd(ig,1)=zb(ig,1)*z1(ig)
538	ENDDO
539
540	! Calculate (d Planck / dT) at the interface temperature
541	! ------------------------------------------------------
542
543	z4st=4.0sigmaptimestep
544	DO ig=1,ngrid
545	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
546	ENDDO
547
548	! Calculate temperature tendency at the interface (dry case)
549	! ----------------------------------------------------------
550	! Sum of fluxes at interface at time t + \delta t gives change in T:
551	! radiative fluxes
552	! turbulent convective (sensible) heat flux
553	! flux (if any) from subsurface
554
555	if(.not.water) then
556
557	DO ig=1,ngrid
558
559	z1(ig) = pcapcal(ig)ptsrf(ig) + cppzb(ig,1)*zc(ig,1)
560	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
561	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
562	& +zdplanck(ig)
563	ztsrf2(ig) = z1(ig) / z2(ig)
564	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
565	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
566	ENDDO
567
568	! Recalculate temperature to top of atmosphere, starting from ground
569	! ------------------------------------------------------------------
570
571	DO ilay=2,nlay
572	DO ig=1,ngrid
573	hh = zh(ig,ilay-1)
574	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
575	ENDDO
576	ENDDO
577
578	endif ! not water
579
580	!-----------------------------------------------------------------------
581	! TRACERS (no vapour)
582	! -------
583
584	if(tracer) then
585
586	! Calculate vertical flux from the bottom to the first layer (dust)
587	! -----------------------------------------------------------------
588	do ig=1,ngrid
589	rho(ig) = zb0(ig,1) /ptimestep
590	end do
591
592	call zerophys(ngrid*nq,pdqsdif)
593
594	! Implicit inversion of q
595	! -----------------------
596	do iq=1,nq
597
598	if (iq.ne.igcm_h2o_vap) then
599
600	DO ig=1,ngrid
601	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
602	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
603	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
604	ENDDO
605
606	DO ilay=nlay-1,2,-1
607	DO ig=1,ngrid
608	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
609	& zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
610	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
611	& zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
612	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
613	ENDDO
614	ENDDO
615
616	if ((water).and.(iq.eq.iice)) then
617	! special case for water ice tracer: do not include
618	! h2o ice tracer from surface (which is set when handling
619	! h2o vapour case (see further down).
620	! zb(ig,1)=0 if iq ne ivap
621	DO ig=1,ngrid
622	z1(ig)=1./(za(ig,1)+
623	& zb(ig,2)*(1.-zdq(ig,2)))
624	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
625	& zb(ig,2)zcq(ig,2))z1(ig)
626	ENDDO
627	else ! general case
628	DO ig=1,ngrid
629	z1(ig)=1./(za(ig,1)+
630	& zb(ig,2)*(1.-zdq(ig,2)))
631	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
632	& zb(ig,2)*zcq(ig,2)
633	& +(-pdqsdif(ig,iq))ptimestep)z1(ig)
634	! tracer flux from surface
635	! currently pdqsdif always zero here,
636	! so last line is superfluous
637	enddo
638	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
639
640
641	! Starting upward calculations for simple tracer mixing (e.g., dust)
642	do ig=1,ngrid
643	zq(ig,1,iq)=zcq(ig,1)
644	end do
645
646	do ilay=2,nlay
647	do ig=1,ngrid
648	zq(ig,ilay,iq)=zcq(ig,ilay)+
649	$ zdq(ig,ilay)*zq(ig,ilay-1,iq)
650	end do
651	end do
652
653	endif ! if (iq.ne.igcm_h2o_vap)
654
655	! Calculate temperature tendency including latent heat term
656	! and assuming an infinite source of water on the ground
657	! ------------------------------------------------------------------
658
659	if (water.and.(iq.eq.igcm_h2o_vap)) then
660
661	! compute evaporation efficiency
662	do ig = 1, ngrid
663	if(nint(rnat(ig)).eq.1)then
664	dryness(ig)=pqsurf(ig,ivap)+pqsurf(ig,iice)
665	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
666	dryness(ig)=MAX(0.,dryness(ig))
667	endif
668	enddo
669
670	do ig=1,ngrid
671
672	! Calculate the value of qsat at the surface (water)
673	call watersat(ptsrf(ig),pplev(ig,1),qsat(ig))
674	call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1)
675	call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2)
676	dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
677	! calculate dQsat / dT by finite differences
678	! we cannot use the updated temperature value yet...
679
680	enddo
681
682	! coefficients for q
683
684	do ig=1,ngrid
685	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
686	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
687	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
688	enddo
689
690	do ilay=nlay-1,2,-1
691	do ig=1,ngrid
692	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
693	$ zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
694	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
695	$ zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
696	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
697	enddo
698	enddo
699
700	do ig=1,ngrid
701	z1(ig)=1./(za(ig,1)+zb(ig,1)*dryness(ig)+
702	$ zb(ig,2)*(1.-zdq(ig,2)))
703	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
704	$ zb(ig,2)zcq(ig,2))z1(ig)
705	zdq(ig,1)=dryness(ig)zb(ig,1)z1(ig)
706	enddo
707
708	! calculation of h0 and h1
709	do ig=1,ngrid
710	zdq0(ig) = dqsat(ig)
711	zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
712
713	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzb(ig,1)*zc(ig,1)
714	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
715	& + zb(ig,1)dryness(ig)RLVTT*
716	& ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
717
718	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
719	& +zdplanck(ig)
720	& +zb(ig,1)dryness(ig)RLVTTzdq0(ig)
721	& (1.0-zdq(ig,1))
722
723	ztsrf2(ig) = z1(ig) / z2(ig)
724	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
725	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
726	enddo
727
728	! calculation of qs and q1
729	do ig=1,ngrid
730	zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf2(ig)
731	zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
732	enddo
733
734	! calculation of evaporation
735	do ig=1,ngrid
736	evap(ig)= zb(ig,1)dryness(ig)(zq(ig,1,ivap)-zq0(ig))
737	dqsdif_total(ig)=evap(ig)
738	enddo
739
740	! recalculate temperature and q(vap) to top of atmosphere, starting from ground
741	do ilay=2,nlay
742	do ig=1,ngrid
743	zq(ig,ilay,iq)=zcq(ig,ilay)
744	& +zdq(ig,ilay)*zq(ig,ilay-1,iq)
745	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
746	end do
747	end do
748
749	do ig=1,ngrid
750
751	! --------------------------------------------------------------------------
752	! On the ocean, if T > 0 C then the vapour tendency must replace the ice one
753	! The surface vapour tracer is actually liquid. To make things difficult.
754
755	if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
756
757	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
758	pdqsdif(ig,iice)=0.0
759
760
761	elseif (nint(rnat(ig)).eq.1) then ! (continent)
762
763	! --------------------------------------------------------
764	! Now check if we've taken too much water from the surface
765	! This can only occur on the continent
766
767	! If water is evaporating / subliming, we take it from ice before liquid
768	! -- is this valid??
769	if(dqsdif_total(ig).lt.0)then
770	pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
771	pdqsdif(ig,iice)=max(-pqsurf(ig,iice)/ptimestep
772	& ,pdqsdif(ig,iice))
773	endif
774	! sublimation only greater than qsurf(ice)
775	! ----------------------------------------
776	! we just convert some liquid to vapour too
777	! if latent heats are the same, no big deal
778	if (-dqsdif_total(ig).gt.pqsurf(ig,iice))then
779	pdqsdif(ig,iice) = -pqsurf(ig,iice)/ptimestep ! removes all the ice!
780	pdqsdif(ig,ivap) = dqsdif_total(ig)/ptimestep
781	& - pdqsdif(ig,iice) ! take the remainder from the liquid instead
782	pdqsdif(ig,ivap) = max(-pqsurf(ig,ivap)/ptimestep
783	& ,pdqsdif(ig,ivap))
784	endif
785
786	endif ! if (rnat.ne.1)
787
788	! If water vapour is condensing, we must decide whether it forms ice or liquid.
789	if(dqsdif_total(ig).gt.0)then ! a bug was here!
790	if(ztsrf2(ig).gt.T_h2O_ice_liq)then
791	pdqsdif(ig,iice)=0.0
792	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
793	else
794	pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
795	pdqsdif(ig,ivap)=0.0
796	endif
797	endif
798
799	end do ! of DO ig=1,ngrid
800	endif ! if (water et iq=ivap)
801	end do ! of do iq=1,nq
802	endif ! traceur
803
804
805	!-----------------------------------------------------------------------
806	! 8. Final calculation of the vertical diffusion tendencies
807	! -----------------------------------------------------------------
808
809	do ilev = 1, nlay
810	do ig=1,ngrid
811	pdudif(ig,ilev)=(zu(ig,ilev)-
812	& (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
813	pdvdif(ig,ilev)=(zv(ig,ilev)-
814	& (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
815	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
816
817	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
818	enddo
819	enddo
820
821	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
822	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zh(ig,1)-ztsrf2(ig))
823	ENDDO
824
825	if (tracer) then
826	do iq = 1, nq
827	do ilev = 1, nlay
828	do ig=1,ngrid
829	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
830	& (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/
831	& ptimestep
832	enddo
833	enddo
834	enddo
835
836	if(water.and.forceWC)then ! force water conservation in model
837	! we calculate the difference and add it to the ground
838	! this is ugly and should be improved in the future
839	do ig=1,ngrid
840	Wtot=0.0
841	do ilay=1,nlay
842	masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
843	! Wtot=Wtot+masse*(zq(ig,ilay,iice)-
844	! & (pq(ig,ilay,iice)+pdqfi(ig,ilay,iice)*ptimestep))
845	Wtot=Wtot+masse*(zq(ig,ilay,ivap)-
846	& (pq(ig,ilay,ivap)+pdqfi(ig,ilay,ivap)*ptimestep))
847	enddo
848	Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
849
850	if(ztsrf2(ig).gt.T_h2O_ice_liq)then
851	pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
852	else
853	pdqsdif(ig,iice)=pdqsdif(ig,iice)-Wdiff
854	endif
855	enddo
856
857	endif
858
859	endif
860
861	if(water)then
862	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
863	endif
864
865	! if(lastcall)then
866	! if(ngrid.eq.1)then
867	! print*,'Saving k.out...'
868	! OPEN(12,file='k.out',form='formatted')
869	! DO ilay=1,nlay
870	! write(12,*) zkh(1,ilay), pplay(1,ilay)
871	! ENDDO
872	! CLOSE(12)
873	! endif
874	! endif
875
876
877	return
878	end

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

Download in other formats: