Context Navigation

newcondens.F @ 1112

Last change on this file since 1112 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 31.9 KB

Line
1	SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep,
2	$ pcapcal,pplay,pplev,ptsrf,pt,
3	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
4	$ piceco2,psolaralb,pemisurf,
5	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
6	$ fluxsurf_sw,zls)
7
8	use tracer_mod, only: noms
9	use surfdat_h, only: emissiv, phisfi
10	IMPLICIT NONE
11	c=======================================================================
12	c subject:
13	c --------
14	c Condensation/sublimation of CO2 ice on the ground and in the
15	c atmosphere
16	c (Scheme described in Forget et al., Icarus, 1998)
17	c
18	c author: Francois Forget 1994-1996
19	c ------
20	c
21	c input:
22	c -----
23	c ngrid nombre de points de verticales
24	c (toutes les boucles de la physique sont au
25	c moins vectorisees sur ngrid)
26	c nlayer nombre de couches
27	c pplay(ngrid,nlayer) Pressure levels
28	c pplev(ngrid,nlayer+1) Pressure levels
29	c pt(ngrid,nlayer) temperature (en K)
30	c ptsrf(ngrid) temperature de surface
31	c
32	c \
33	c pdt(ngrid,nlayer)\ derivee temporelle physique avant condensation
34	c / ou sublimation pour pt,ptsrf
35	c pdtsrf(ngrid) /
36	c
37	c output:
38	c -------
39	c
40	c pdpsrf(ngrid) \ derivee temporelle physique (contribution de
41	c pdtc(ngrid,nlayer) / la condensation ou sublimation) pour Ps,pt,ptsrf
42	c pdtsrfc(ngrid) /
43	c
44	c Entree/sortie
45	c -------------
46	c
47	c piceco2(ngrid) : quantite de glace co2 au sol (kg/m2)
48	c psolaralb(ngrid,2) : albedo au sol
49	c pemisurf(ngrid) : emissivite du sol
50
51	c
52	c=======================================================================
53	c
54	c 0. Declarations :
55	c ------------------
56	c
57	#include "dimensions.h"
58	!#include "dimphys.h"
59	#include "comcstfi.h"
60	!#include "surfdat.h"
61	!#include "comgeomfi.h"
62	#include "comvert.h"
63	!#include "paramet.h"
64	#include "callkeys.h"
65	!#include "tracer.h"
66
67	c-----------------------------------------------------------------------
68	c Arguments :
69	c ---------
70	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
71	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
72	INTEGER,INTENT(IN) :: nq ! number of tracers
73
74	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
75	REAL,INTENT(IN) :: pcapcal(ngrid)
76	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
77	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
78	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
79	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
80	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
81	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
82	! previous physical processes (K/s)
83	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
84	! from previous physical processes
85	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
86	! from previous physical processes
87	REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
88	! previous physical processes (K/s)
89	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
90	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
91	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
92	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
93	! previous physical processes
94	REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
95	REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
96	REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
97
98	! tendencies due to CO2 condensation/sublimation:
99	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
100	REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
101	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
102	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
103	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
104	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
105
106	! added to calculate flux dependent albedo:
107	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
108	real,intent(in) :: zls ! solar longitude (rad)
109
110	c
111	c Local variables :
112	c -----------------
113
114	c variables used for albedo parametrization
115	c --------------------------------------------
116	INTEGER i,j
117	c REAL Fluxmean(jjp1)
118	INTEGER l,ig,iq,icap,nmix
119	LOGICAL transparency, fluxdependent
120	c flag transparency if you want to make the co2ice semi-transparent
121	PARAMETER(transparency=.true.)
122	c flag fluxdependent if you want the co2ice albedo to be dependent on
123	c the incident solar flux
124	PARAMETER(fluxdependent=.false.)
125	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid)
126	REAL zt(ngrid,nlayer)
127	REAL zcpi
128	REAL ztcond (ngrid,nlayer+1)
129	REAL ztcondsol(ngrid)
130	REAL zdiceco2(ngrid)
131	REAL zcondicea(ngrid,nlayer)
132	REAL zcondices(ngrid)
133	REAL zfallice(ngrid,nlayer+1) , zfallheat
134	REAL zmflux(nlayer+1)
135	REAL zu(nlayer),zv(nlayer)
136	REAL zq(nlayer,nq),zq1(nlayer)
137	REAL ztsrf(ngrid)
138	REAL ztc(nlayer), ztm(nlayer+1)
139	REAL zum(nlayer+1) , zvm(nlayer+1)
140	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
141	REAL masse(nlayer),w(nlayer+1)
142	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
143	LOGICAL condsub(ngrid)
144
145	c variable speciale diagnostique
146	real tconda1(ngrid,nlayer)
147	real tconda2(ngrid,nlayer)
148	c REAL zdiceco2a(ngrid) ! for diagnostic only
149	real zdtsig (ngrid,nlayer)
150	real zdt (ngrid,nlayer)
151	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
152	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
153	! then condensation temperature is computed using partial pressure of CO2
154	logical,parameter :: improved_ztcond=.true.
155	! Bound co2 (tracer) values...
156	logical,parameter :: bound_qco2=.false.
157	real,parameter :: qco2max=1.1
158	real,parameter :: qco2mini=0.1
159	real :: zqco2
160
161	c local saved variables
162	integer,save :: ico2 ! index of CO2 tracer
163	real,save :: qco2min,qco2,mmean
164	real,allocatable,save :: emisref(:)
165	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
166	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
167	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
168	REAL,SAVE :: acond,bcond,ccond
169	! REAL,SAVE :: albediceF(ngrid)
170	real,save :: m_co2, m_noco2, A , B
171
172	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
173
174	integer flag
175
176	c----------------------------------------------------------------------
177
178	c Initialisation
179	c --------------
180	c
181	IF (firstcall) THEN
182
183	allocate(emisref(ngrid))
184
185	bcond=1./tcond1mb
186	ccond=cpp/(g*latcond)
187	acond=r/latcond
188
189	firstcall=.false.
190	write(,) 'Newcondens: improved_ztcond=',improved_ztcond
191	write(,) 'Newcondens: bound_qco2=',bound_qco2
192	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
193	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
194
195	ico2=0
196
197	if (tracer) then
198	c Prepare Special treatment if one of the tracer is CO2 gas
199	do iq=1,nq
200	if (noms(iq).eq."co2") then
201	ico2=iq
202	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
203	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
204	c Compute A and B coefficient use to compute
205	c mean molecular mass Mair defined by
206	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
207	c 1/Mair = A*q(ico2) + B
208	A =(1/m_co2 - 1/m_noco2)
209	B=1/m_noco2
210	endif
211	enddo
212	c minimum CO2 mix. ratio below which mixing occur with layer above:
213	qco2min =0.75
214	end if
215	ENDIF ! of IF (firstcall)
216	zcpi=1./cpp
217	c
218	c======================================================================
219	c Calcul of CO2 condensation sublimation
220	c ============================================================
221	c
222	c Used variable :
223	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
224	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
225	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
226	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
227	c
228	c pdtc(ngrid,nlayer) : dT/dt due to cond/sub
229	c
230	c
231	c Tendencies set to 0 (except pdtc)
232	c -------------------------------------
233	DO l=1,nlayer
234	DO ig=1,ngrid
235	zcondicea(ig,l) = 0.
236	zfallice(ig,l) = 0.
237	pduc(ig,l) = 0
238	pdvc(ig,l) = 0
239	END DO
240	END DO
241
242	DO iq=1,nq
243	DO l=1,nlayer
244	DO ig=1,ngrid
245	pdqc(ig,l,iq) = 0
246	END DO
247	END DO
248	END DO
249
250	DO ig=1,ngrid
251	zfallice(ig,nlayer+1) = 0.
252	zcondices(ig) = 0.
253	pdtsrfc(ig) = 0.
254	pdpsrf(ig) = 0.
255	condsub(ig) = .false.
256	zdiceco2(ig) = 0.
257	ENDDO
258	zfallheat=0
259
260	c *************************
261	c ATMOSPHERIC CONDENSATION
262	c *************************
263
264	c Compute CO2 Volume mixing ratio
265	c -------------------------------
266	if (improved_ztcond.and.(ico2.ne.0)) then
267	DO l=1,nlayer
268	DO ig=1,ngrid
269	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
270	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
271	mmean=1/(A*qco2 +B)
272	vmr_co2(ig,l) = qco2*mmean/m_co2
273	ENDDO
274	ENDDO
275	else
276	DO l=1,nlayer
277	DO ig=1,ngrid
278	vmr_co2(ig,l)=0.95
279	ENDDO
280	ENDDO
281	end if
282
283	c forecast of atmospheric temperature zt and frost temperature ztcond
284	c --------------------------------------------------------------------
285
286	DO l=1,nlayer
287	DO ig=1,ngrid
288	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
289	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
290	ztcond(ig,l)=
291	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
292	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
293	ENDDO
294	ENDDO
295
296	ztcond(:,nlayer+1)=ztcond(:,nlayer)
297
298	c Condensation/sublimation in the atmosphere
299	c ------------------------------------------
300	c (calcul of zcondicea , zfallice and pdtc)
301	c
302	DO l=nlayer , 1, -1
303	DO ig=1,ngrid
304	pdtc(ig,l)=0.
305	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
306	condsub(ig)=.true.
307	IF (zfallice(ig,l+1).gt.0) then
308	zfallheat=zfallice(ig,l+1)*
309	& (pphi(ig,l+1)-pphi(ig,l) +
310	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
311	ELSE
312	zfallheat=0.
313	ENDIF
314	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
315	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
316	& ccondpdtc(ig,l)- zfallheat
317	c Case when the ice from above sublimes entirely
318	c """""""""""""""""""""""""""""""""""""""""""""""
319	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
320	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
321	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
322	zcondicea(ig,l)= -zfallice(ig,l+1)
323	END IF
324
325	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
326	END IF
327	ENDDO
328	ENDDO
329
330	c *************************
331	c SURFACE CONDENSATION
332	c *************************
333
334	c forecast of ground temperature ztsrf and frost temperature ztcondsol
335	c --------------------------------------------------------------------
336	DO ig=1,ngrid
337	ztcondsol(ig)=
338	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
339	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
340	ENDDO
341
342	c
343	c Condensation/sublimation on the ground
344	c --------------------------------------
345	c (calcul of zcondices , pdtsrfc)
346	c
347	DO ig=1,ngrid
348	IF(ig.GT.ngrid/2+1) THEN
349	icap=2
350	ELSE
351	icap=1
352	ENDIF
353
354	c
355	c Loop on where we have condensation/ sublimation
356	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
357	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
358	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
359	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
360	condsub(ig) = .true.
361
362	IF (zfallice(ig,1).gt.0) then
363	zfallheat=zfallice(ig,1)*
364	& (pphi(ig,1)- phisfi(ig) +
365	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
366	ELSE
367	zfallheat=0.
368	ENDIF
369
370	c condensation or partial sublimation of CO2 ice
371	c """""""""""""""""""""""""""""""""""""""""""""""
372	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
373	& /(latcond*ptimestep) - zfallheat
374	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
375
376	c If the entire CO_2 ice layer sublimes
377	c """"""""""""""""""""""""""""""""""""""""""""""""""""
378	c (including what has just condensed in the atmosphere)
379
380	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
381	& -zcondices(ig))THEN
382	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
383	pdtsrfc(ig)=(latcond/pcapcal(ig))*
384	& (zcondices(ig)+zfallheat)
385	END IF
386
387	c Changing CO2 ice amount and pressure :
388	c """"""""""""""""""""""""""""""""""""
389
390	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
391	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
392	pdpsrf(ig) = -zdiceco2(ig)*g
393
394	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
395	PRINT*,'STOP in condens'
396	PRINT*,'condensing more than total mass'
397	PRINT*,'Grid point ',ig
398	PRINT*,'Ps = ',pplev(ig,1)
399	PRINT*,'d Ps = ',pdpsrf(ig)
400	STOP
401	ENDIF
402	END IF ! if there is condensation/sublimmation
403	ENDDO ! of DO ig=1,ngrid
404
405	c ********************************************************************
406	c Surface albedo and emissivity of the surface below the snow (emisref)
407	c ********************************************************************
408	c Prepare the case where albedo varies with insolation:
409	c ----------------------------------------------------
410	! if (fluxdependent) then
411	!
412	c Calcul du flux moyen (zonal mean)
413	! do j=1,jjp1
414	! Fluxmean(j)=0
415	! do i=1,iim
416	! ig=1+(j-2)*iim +i
417	! if(j.eq.1) ig=1
418	! if(j.eq.jjp1) ig=ngrid
419	! Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
420	! $ +fluxsurf_sw(ig,2)
421	! enddo
422	! Fluxmean(j)=Fluxmean(j)/float(iim)
423	! enddo
424	!
425	c const A and B used to calculate the albedo which depends on solar flux
426	c albedice=constA+constB*Flux
427	c constT = time step to calculate the solar flux when flux decreases
428	! constA=0.26
429	c constA=0.33
430	c constA=0.186
431	! constB=0.00187
432	! constT=10
433	! endif ! of if (fluxdependent)
434
435	! Check that amont of CO2 ice is not problematic
436	DO ig=1,ngrid
437	if(.not.piceco2(ig).ge.0.) THEN
438	if(piceco2(ig).le.-5.e-8) print*,
439	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
440	piceco2(ig)=0.
441	endif
442	ENDDO
443
444	! Set albedo and emissivity of the surface
445	! ----------------------------------------
446	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
447
448	c Calcul de l'albedo
449	c ------------------
450	! do ig =1,ngrid
451	! IF(ig.GT.ngrid/2+1) THEN
452	! icap=2
453	! ELSE
454	! icap=1
455	! ENDIF
456	! IF(firstcall2) THEN
457	! albediceF(ig)=albedice(icap)
458	! ENDIF
459	c if there is still co2ice ccccccccccccccccccccccc
460	! if (piceco2(ig).gt.0) then
461	! emisref(ig) = emisice(icap)
462
463	c if flux dependent albedo is used
464	c --------------------------------
465	! if (fluxdependent) then
466	! j=INT((ig-2)/iim)+2
467	! if(ig.eq.1) j=1
468	! if(ig.eq.ngrid) j=jjp1
469	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
470	c $ constA+constB*250)
471	! albediceF_new(ig)=constA+constB*Fluxmean(j)
472	! if (albediceF(ig).gt.albediceF_new(ig)) then
473	! albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
474	! $ constT)*(albediceF_new(ig)-albediceF(ig))
475	! else
476	! albediceF(ig)=albediceF_new(ig)
477	! endif
478	c if part of the ice is transparent
479	c slopy = pente de la droite: alpha = y*co2ice/1620
480	c pour une valeur superieur a une epaisseur de glace donnee
481	c ici, epaisseur limite = 10cm
482	! if (transparency) then
483	! slopy=1/(1620*0.10)
484	! alpha=MIN(slopy*piceco2(ig),1.)
485	! psolaralb(ig,1) = alpha*albediceF(ig)
486	! $ +(1-alpha)*albedodat(ig)
487	! psolaralb(ig,2) = psolaralb(ig,1)
488	! else
489	! psolaralb(ig,1) = albediceF(ig)
490	! psolaralb(ig,2) = psolaralb(ig,1)
491	! endif
492	! else
493	c transparency set to true and fluxdependent set to false
494	! if (transparency) then
495	! slopy=1/(1620*0.10)
496	! alpha=MIN(slopy*piceco2(ig),1.)
497	! psolaralb(ig,1) = alpha*albedice(icap)
498	! $ +(1-alpha)*albedodat(ig)
499	! psolaralb(ig,2) = psolaralb(ig,1)
500	! else
501	c simplest case: transparency and flux dependent set to false
502	! psolaralb(ig,1) = albedice(icap)
503	! psolaralb(ig,2) = albedice(icap)
504	! endif
505	! endif
506	c no more co2ice, albedo = ground albedo
507	! else
508	! psolaralb(ig,1) = albedodat(ig)
509	! psolaralb(ig,2) = albedodat(ig)
510	! emisref(ig) = emissiv
511	! pemisurf(ig) = emissiv
512	! endif
513	! end do ! end of the ig loop
514
515	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
516	DO ig=1,ngrid
517	if (piceco2(ig).eq.0) then
518	pemisurf(ig)=emissiv
519	endif
520	ENDDO
521
522	! firstcall2=.false.
523	c ***************************************************************
524	c Correction to account for redistribution between sigma or hybrid
525	c layers when changing surface pressure (and warming/cooling
526	c of the CO2 currently changing phase).
527	c *************************************************************
528
529	DO ig=1,ngrid
530	if (condsub(ig)) then
531	do l=1,nlayer
532	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
533	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
534	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
535	do iq=1,nq
536	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
537	enddo
538	end do
539
540	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
541	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
542
543	zmflux(1) = -zcondices(ig)
544	DO l=1,nlayer
545	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
546	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
547	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
548	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
549	END DO
550
551	c Mass of each layer
552	c ------------------
553	DO l=1,nlayer
554	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
555	END DO
556
557
558	c Corresponding fluxes for T,U,V,Q
559	c """"""""""""""""""""""""""""""""
560
561	c averaging operator for TRANSPORT
562	c """"""""""""""""""""""""""""""""
563	c Value transfert at the surface interface when condensation
564	c sublimation:
565	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
566	zum(1) = 0
567	zvm(1) = 0
568	do iq=1,nq
569	zqm(1,iq)=0. ! most tracer do not condense !
570	enddo
571	c Special case if one of the tracer is CO2 gas
572	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
573
574	c Van Leer scheme:
575	DO l=1,nlayer+1
576	w(l)=-zmflux(l)*ptimestep
577	END DO
578	call vl1d(ztc,2.,masse,w,ztm)
579	call vl1d(zu ,2.,masse,w,zum)
580	call vl1d(zv ,2.,masse,w,zvm)
581	do iq=1,nq
582	do l=1,nlayer
583	zq1(l)=zq(l,iq)
584	enddo
585	zqm1(1)=zqm(1,iq)
586	call vl1d(zq1,2.,masse,w,zqm1)
587	do l=2,nlayer
588	zq( l,iq)=zq1(l)
589	zqm(l,iq)=zqm1(l)
590	enddo
591	enddo
592
593	c Surface condensation affects low winds
594	if (zmflux(1).lt.0) then
595	zum(1)= zu(1) * (w(1)/masse(1))
596	zvm(1)= zv(1) * (w(1)/masse(1))
597	if (w(1).gt.masse(1)) then ! ensure numerical stability
598	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
599	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
600	end if
601	end if
602
603	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
604	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
605	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
606	do iq=1,nq
607	zqm(nlayer+1,iq)= zq(nlayer,iq)
608	enddo
609
610	#ifdef MESOSCALE
611	!!!! AS: This part must be commented in the mesoscale model
612	!!!! AS: ... to avoid instabilities.
613	!!!! AS: you have to compile with -DMESOSCALE to do so
614	#else
615	c Tendencies on T, U, V, Q
616	c """"""""""""""""""""""""
617	DO l=1,nlayer
618
619	c Tendencies on T
620	zdtsig(ig,l) = (1/masse(l)) *
621	& ( zmflux(l)*(ztm(l) - ztc(l))
622	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
623	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
624	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
625
626	c Tendencies on U
627	pduc(ig,l) = (1/masse(l)) *
628	& ( zmflux(l)*(zum(l) - zu(l))
629	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
630
631
632	c Tendencies on V
633	pdvc(ig,l) = (1/masse(l)) *
634	& ( zmflux(l)*(zvm(l) - zv(l))
635	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
636
637	END DO
638	#endif
639
640	c Tendencies on Q
641	do iq=1,nq
642	! if (noms(iq).eq.'co2') then
643	if (iq.eq.ico2) then
644	c SPECIAL Case when the tracer IS CO2 :
645	DO l=1,nlayer
646	pdqc(ig,l,iq)= (1/masse(l)) *
647	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
648	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
649	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
650	END DO
651	else
652	DO l=1,nlayer
653	pdqc(ig,l,iq)= (1/masse(l)) *
654	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
655	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
656	& + zcondicea(ig,l)*zq(l,iq) )
657	END DO
658	end if
659	enddo
660
661	c --------------------------------------------------------
662	c Roughly Simulate Molecular mixing when CO2 is too depleted by
663	c Surface condensation (mixing starts if qco2 < qco2min )
664	c FF 06/2004
665	c WARNING : this is now done in convadj, better (FF 02/2005)
666	c --------------------------------------------------------
667	flag=0 ! now done in convadj : must be =0
668	if (flag.eq.1) then
669	if(ico2.gt.0) then ! relevant only if one tracer is CO2
670	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
671	& .lt.qco2min) then
672	do iq=1,nq
673	zq(1,iq)=pq(ig,1,iq)
674	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
675	Smq(1,iq) = masse(1)*zq(1,iq)
676	end do
677	Sm(1) = masse(1)
678	do l =2,nlayer
679	do iq=1,nq
680	zq(l,iq)=pq(ig,l,iq)
681	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
682	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
683	end do
684	sm(l) = sm(l-1) + masse(l)
685	if(zq(l,ico2).gt.qco2min) then
686	c mixmas: mass of atmosphere that must be mixed to reach qco2min
687	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
688	& /(zq(l,ico2)-qco2min)
689	if((mixmas.le.sm(l)))then
690	c OK if mixed mass less than mass of layers affected
691	nmix=l ! number of layer affected by mixing
692	goto 99
693	end if
694	end if
695	end do
696	99 continue
697	do iq=1,nq
698	qmix=zq(nmix,iq)
699	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
700	do l=1,nmix-1
701	pdqc(ig,l,iq)=
702	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
703	end do
704	c layer only partly mixed :
705	pdqc(ig,nmix,iq)=(
706	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
707	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
708
709	end do
710	end if
711	end if
712
713	endif ! (flag.eq.1)
714	end if ! if (condsub)
715	END DO ! loop on ig
716
717	c ***************************************************************
718	c CO2 snow / clouds scheme
719	c ***************************************************************
720
721	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
722	& zcondicea,zcondices,zfallice,pemisurf)
723
724	c ***************************************************************
725	c Ecriture des diagnostiques
726	c ***************************************************************
727
728	c DO l=1,nlayer
729	c DO ig=1,ngrid
730	c Taux de cond en kg.m-2.pa-1.s-1
731	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
732	c Taux de cond en kg.m-3.s-1
733	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
734	c END DO
735	c END DO
736	c call WRITEDIAGFI(ngrid,'tconda1',
737	c &'Taux de condensation CO2 atmospherique /Pa',
738	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
739	c call WRITEDIAGFI(ngrid,'tconda2',
740	c &'Taux de condensation CO2 atmospherique /m',
741	c & 'kg.m-3.s-1',3,tconda2)
742
743	! output falling co2 ice in 1st layer:
744	! call WRITEDIAGFI(ngrid,'fallice',
745	! &'Precipitation of co2 ice',
746	! & 'kg.m-2.s-1',2,zfallice(1,1))
747
748	!! Specific stuff to bound co2 tracer values ....
749	if (bound_qco2.and.(ico2.ne.0)) then
750	do ig=1,ngrid
751	do l=1,nlayer
752	zqco2=pq(ig,l,ico2)
753	& +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
754	if (zqco2.gt.qco2max) then
755	! correct pdqc:
756	pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
757	& -pdq(ig,l,ico2)
758	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
759	& " so that co2 conc. does not exceed",qco2max
760	write(,) " ig:",ig," l:",l
761	endif ! of if (zqco2.gt.qco2max)
762	if (zqco2.lt.qco2mini) then
763	! correct pdqc:
764	pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
765	& -pdq(ig,l,ico2)
766	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
767	& " so that co2 conc. is not less than",qco2mini
768	write(,) " ig:",ig," l:",l
769	endif ! of if (zqco2.lt.qco2mini)
770	end do
771	enddo
772	endif ! of if (bound_qco2.and.(ico2.ne.0)) then
773
774	return
775	end
776
777
778
779	c *****************************************************************
780	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
781	c
782	c
783	c Operateur de moyenne inter-couche pour calcul de transport type
784	c Van-Leer " pseudo amont " dans la verticale
785	c q,w sont des arguments d'entree pour le s-pg ....
786	c masse : masse de la couche Dp/g
787	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
788	c pente_max = 2 conseillee
789	c
790	c
791	c --------------------------------------------------------------------
792	IMPLICIT NONE
793
794	#include "dimensions.h"
795
796	c
797	c
798	c
799	c Arguments:
800	c ----------
801	real masse(llm),pente_max
802	REAL q(llm),qm(llm+1)
803	REAL w(llm+1)
804	c
805	c Local
806	c ---------
807	c
808	INTEGER l
809	c
810	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
811	real sigw, Mtot, MQtot
812	integer m
813	c integer ismax,ismin
814
815
816	c On oriente tout dans le sens de la pression
817	c W > 0 WHEN DOWN !!!!!!!!!!!!!
818
819	do l=2,llm
820	dzqw(l)=q(l-1)-q(l)
821	adzqw(l)=abs(dzqw(l))
822	enddo
823
824	do l=2,llm-1
825	if(dzqw(l)*dzqw(l+1).gt.0.) then
826	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
827	else
828	dzq(l)=0.
829	endif
830	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
831	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
832	enddo
833
834	dzq(1)=0.
835	dzq(llm)=0.
836
837	do l = 1,llm-1
838
839	c Regular scheme (transfered mass < layer mass)
840	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
841	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
842	sigw=w(l+1)/masse(l+1)
843	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
844	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
845	sigw=w(l+1)/masse(l)
846	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
847
848	c Extended scheme (transfered mass > layer mass)
849	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
850	else if(w(l+1).gt.0.) then
851	m=l+1
852	Mtot = masse(m)
853	MQtot = masse(m)*q(m)
854	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
855	m=m+1
856	Mtot = Mtot + masse(m)
857	MQtot = MQtot + masse(m)*q(m)
858	end do
859	if (m.lt.llm) then
860	sigw=(w(l+1)-Mtot)/masse(m+1)
861	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
862	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
863	else
864	w(l+1) = Mtot
865	qm(l+1) = Mqtot / Mtot
866	write(,) 'top layer is disapearing !'
867	stop
868	end if
869	else ! if(w(l+1).lt.0)
870	m = l-1
871	Mtot = masse(m+1)
872	MQtot = masse(m+1)*q(m+1)
873	if (m.gt.0) then ! because some compilers will have problems
874	! evaluating masse(0)
875	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
876	m=m-1
877	Mtot = Mtot + masse(m+1)
878	MQtot = MQtot + masse(m+1)*q(m+1)
879	if (m.eq.0) exit
880	end do
881	endif
882	if (m.gt.0) then
883	sigw=(w(l+1)+Mtot)/masse(m)
884	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
885	& (q(m)-0.5(1.+sigw)dzq(m)) )
886	else
887	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
888	end if
889	end if
890	enddo
891
892	c boundary conditions (not used in newcondens !!)
893	c qm(llm+1)=0.
894	c if(w(1).gt.0.) then
895	c qm(1)=q(1)
896	c else
897	c qm(1)=0.
898	c end if
899
900	return
901	end

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

Download in other formats:

Original Format