Context Navigation

co2condens_mod.F @ 2596

Last change on this file since 2596 was 2566, checked in by cmathe, 4 years ago
GCM MARS: merge co2condens subroutines
File size: 42.5 KB

Line
1	MODULE co2condens_mod
2
3	IMPLICIT NONE
4
5	logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2
6
7	CONTAINS
8
9	SUBROUTINE co2condens(ngrid,nlayer,nq,ptimestep,
10	$ pcapcal,pplay,pplev,ptsrf,pt,
11	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
12	$ piceco2,psolaralb,pemisurf,rdust,
13	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
14	$ fluxsurf_sw,zls,
15	$ zdqssed_co2,pcondicea_co2microp,
16	$ pdqsc)
17
18	use tracer_mod, only: noms, igcm_h2o_ice, igcm_h2o_vap,
19	& igcm_dust_mass, igcm_dust_number,
20	& igcm_ccn_mass, igcm_ccn_number,
21	& igcm_hdo_ice, igcm_hdo_vap,
22	& nqperes,nqfils, ! MVals: variables isotopes
23	& qperemin,masseqmin,
24	& igcm_co2
25	use surfdat_h, only: emissiv, phisfi
26	use geometry_mod, only: latitude, ! grid point latitudes (rad)
27	& longitude_deg, latitude_deg
28	use planete_h, only: obliquit
29	use comcstfi_h, only: cpp, g, r, pi
30	use dust_param_mod, only: freedust
31
32	#ifndef MESOSCALE
33	USE vertical_layers_mod, ONLY: ap, bp
34	#endif
35	IMPLICIT NONE
36	c=======================================================================
37	c subject:
38	c --------
39	c Condensation/sublimation of CO2 ice on the ground and in the
40	c atmosphere
41	c (Scheme described in Forget et al., Icarus, 1998)
42	c
43	c author: Francois Forget 1994-1996 ; updated 1996 -- 2018
44	c ------
45	c adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
46	c
47	c=======================================================================
48	c
49	c 0. Declarations :
50	c ------------------
51	c
52	include "callkeys.h"
53
54	c-----------------------------------------------------------------------
55	c Arguments :
56	c ---------
57	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
58	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
59	INTEGER,INTENT(IN) :: nq ! number of tracers
60
61	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
62	REAL,INTENT(IN) :: pcapcal(ngrid)
63	REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
64	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
65	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
66	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
67	REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
68	REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
69	! previous physical processes (K/s)
70	REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
71	! from previous physical processes
72	REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
73	! from previous physical processes
74	REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
75	! previous physical processes (K/s)
76	REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
77	REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
78	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
79	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
80	! previous physical processes
81
82	REAL,INTENT(IN) :: zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1)
83	REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1)
84
85	REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
86	REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
87	REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
88	REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius
89
90	! tendencies due to CO2 condensation/sublimation:
91	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
92	REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
93	REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
94	REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
95	REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
96	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
97	REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers
98
99	! added to calculate flux dependent albedo:
100	REAL,intent(in) :: fluxsurf_sw(ngrid,2)
101	real,intent(in) :: zls ! solar longitude (rad)
102
103	c
104	c Local variables :
105	c -----------------
106
107	INTEGER i,j
108	INTEGER l,ig,iq,icap
109	REAL zt(ngrid,nlayer)
110	REAL zcpi
111	REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
112	REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
113	REAL zdiceco2(ngrid)
114	REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
115	REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s)
116	REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
117	REAL condens_layer(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s)
118	REAL condens_column(ngrid) ! co2clouds: sum(condens_layer(ig,:)) (kg/m2/s)
119	REAL zfallheat
120	REAL zmflux(nlayer+1)
121	REAL zu(nlayer),zv(nlayer)
122	REAL zqc(nlayer,nq),zq1(nlayer)
123	REAL ztsrf(ngrid)
124	REAL ztc(nlayer), ztm(nlayer+1)
125	REAL zum(nlayer+1) , zvm(nlayer+1)
126	REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
127	REAL masse(nlayer),w(nlayer+1)
128	REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
129	REAL availco2
130	LOGICAL condsub(ngrid)
131
132	real :: emisref(ngrid)
133
134	REAL zdq_scav(ngrid,nlayer,nq) ! tendancy due to scavenging by co2
135	REAL zq(ngrid,nlayer,nq)
136
137	c variable speciale diagnostique
138	real tconda1(ngrid,nlayer)
139	real tconda2(ngrid,nlayer)
140	c REAL zdiceco2a(ngrid) ! for diagnostic only
141	real zdtsig (ngrid,nlayer)
142	real zdt (ngrid,nlayer)
143	real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
144	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
145	! then condensation temperature is computed using partial pressure of CO2
146	logical,parameter :: improved_ztcond=.true.
147
148	c local saved variables
149	integer,save :: ico2 ! index of CO2 tracer
150	real,save :: qco2,mmean
151	real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
152	real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
153	real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
154	REAL,SAVE :: acond,bcond,ccond
155	real,save :: m_co2, m_noco2, A , B
156
157	LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
158
159	c D.BARDET: for debug
160	real ztc3D(ngrid,nlayer)
161	REAL ztm3D(ngrid,nlayer)
162	REAL zmflux3D(ngrid,nlayer)
163
164	c MVals: variables isotopes
165	REAL Ratio1(nlayer),Ratiom1(nlayer+1)
166	REAL masseq(nlayer),wq(nlayer+1)
167	INTEGER ifils,iq2
168	c----------------------------------------------------------------------
169
170	c Initialisation
171	c --------------
172	c
173	! AS: firstcall OK absolute
174	IF (firstcall) THEN
175
176	bcond=1./tcond1mb
177	ccond=cpp/(g*latcond)
178	acond=r/latcond
179
180	firstcall=.false.
181	write(,) 'CO2condens: improved_ztcond=',improved_ztcond
182	PRINT*,'In co2condens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
183	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
184
185	ico2=0
186
187	if (tracer) then
188	c Prepare Special treatment if one of the tracer is CO2 gas
189	do iq=1,nq
190	if (noms(iq).eq."co2") then
191	ico2=iq
192	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
193	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
194	c Compute A and B coefficient use to compute
195	c mean molecular mass Mair defined by
196	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
197	c 1/Mair = A*q(ico2) + B
198	A =(1/m_co2 - 1/m_noco2)
199	B=1/m_noco2
200	endif
201	enddo
202	end if
203	ENDIF ! of IF (firstcall)
204	zcpi=1./cpp
205
206	c
207	c======================================================================
208	c Calcul of CO2 condensation sublimation
209	c ============================================================
210	c
211	c Used variable :
212	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
213	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
214	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
215	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
216	c
217	c pdtc(ngrid,nlayer) : dT/dt due to cond/sub
218	c
219	c
220	c Tendencies set to 0
221	c -------------------------------------
222	zcondicea(1:ngrid,1:nlayer) = 0.
223	zfallice(1:ngrid,1:nlayer+1) = 0.
224	pduc(1:ngrid,1:nlayer) = 0
225	pdvc(1:ngrid,1:nlayer) = 0
226	pdtc(1:ngrid,1:nlayer) = 0.
227	pdqsc(1:ngrid,1:nq) = 0
228
229	pdqc(1:ngrid,1:nlayer,1:nq) = 0
230
231	zcondices(1:ngrid) = 0.
232	pdtsrfc(1:ngrid) = 0.
233	pdpsrf(1:ngrid) = 0.
234	condsub(1:ngrid) = .false.
235	zdiceco2(1:ngrid) = 0.
236
237	zfallheat=0
238
239	zdq_scav(:,:,:)=0.
240
241	c Update tendencies from previous processes
242	c -------------------------------------
243	DO l=1,nlayer
244	DO ig=1,ngrid
245	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
246	do iq=1,nq
247	zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
248	enddo
249	ENDDO
250	ENDDO
251
252	c *************************
253	c ATMOSPHERIC CONDENSATION
254	c *************************
255
256	c Compute CO2 Volume mixing ratio
257	c -------------------------------
258	if (improved_ztcond.and.(ico2.ne.0)) then
259	DO l=1,nlayer
260	DO ig=1,ngrid
261	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
262	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
263	mmean=1/(A*qco2 +B)
264	vmr_co2(ig,l) = qco2*mmean/m_co2
265	ENDDO
266	ENDDO
267	else
268	DO l=1,nlayer
269	DO ig=1,ngrid
270	vmr_co2(ig,l)=0.95
271	ENDDO
272	ENDDO
273	endif
274
275	IF (.NOT. co2clouds) then
276	c forecast of atmospheric temperature zt and frost temperature ztcond
277	c --------------------------------------------------------------------
278
279	DO l=1,nlayer
280	DO ig=1,ngrid
281	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
282	if (pplay(ig,l).ge.1e-4) then
283	ztcond(ig,l)=
284	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
285	else
286	ztcond(ig,l)=0.0 !mars Monica
287	endif
288	ENDDO
289	ENDDO
290
291	ztcond(:,nlayer+1)=ztcond(:,nlayer)
292
293	c Condensation/sublimation in the atmosphere
294	c ------------------------------------------
295	c (calcul of zcondicea , zfallice and pdtc)
296	c
297	DO l=nlayer , 1, -1
298	DO ig=1,ngrid
299	pdtc(ig,l)=0.
300	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
301	condsub(ig)=.true.
302	IF (zfallice(ig,l+1).gt.0) then
303	zfallheat=zfallice(ig,l+1)*
304	& (pphi(ig,l+1)-pphi(ig,l) +
305	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
306	ELSE
307	zfallheat=0.
308	ENDIF
309	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
310	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
311	& ccondpdtc(ig,l)- zfallheat
312	c Case when the ice from above sublimes entirely
313	c """""""""""""""""""""""""""""""""""""""""""""""
314	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
315	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
316	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
317	zcondicea(ig,l)= -zfallice(ig,l+1)
318	END IF
319
320	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
321	END IF
322	ENDDO
323	ENDDO
324
325	if (scavco2cond) then
326	call scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,zq,
327	& rdust,zcondicea,zfallice,zdq_scav,pdqsc)
328	endif
329
330	ELSE ! if co2 clouds
331	condens_layer(:,:) = 0.
332	condens_column(:) = 0.
333	DO l=nlayer , 1, -1
334	DO ig=1,ngrid
335	condens_layer(ig,l) = pcondicea_co2microp(ig,l)*
336	& (pplev(ig,l) - pplev(ig,l+1))/g
337	ENDDO
338	ENDDO
339	DO ig=1,ngrid
340	condens_column(ig) = sum(condens_layer(ig,:))
341	piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep
342	ENDDO
343	ENDIF ! end of if co2clouds
344
345	call WRITEdiagfi(ngrid,"pdtc_atm",
346	& "temperature tendency due to CO2 condensation",
347	& " ",3,pdtc)
348
349	call WRITEdiagfi(ngrid,"zcondicea",
350	& "",
351	& " ",3,zcondicea)
352
353	call WRITEdiagfi(ngrid,"zfallice",
354	& "",
355	& " ",2,zfallice(ngrid,1))
356
357
358	c *************************
359	c SURFACE CONDENSATION
360	c *************************
361
362	c forecast of ground temperature ztsrf and frost temperature ztcondsol
363	c --------------------------------------------------------------------
364	DO ig=1,ngrid
365	ztcondsol(ig)=
366	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
367	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
368	ENDDO
369
370	c
371	c Condensation/sublimation on the ground
372	c --------------------------------------
373	c (compute zcondices and pdtsrfc)
374	c
375	c No microphysics of CO2 clouds
376	if (.not.co2clouds)then
377	DO ig=1,ngrid
378	IF(latitude(ig).lt.0) THEN
379	! Southern hemisphere
380	icap=2
381	ELSE
382	! Northern hemisphere
383	icap=1
384	ENDIF
385
386	c
387	c Loop on where we have condensation/ sublimation
388	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
389	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
390	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
391	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
392	condsub(ig) = .true.
393
394	IF (zfallice(ig,1).gt.0) then
395	zfallheat=zfallice(ig,1)*
396	& (pphi(ig,1)- phisfi(ig) +
397	& cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
398	ELSE
399	zfallheat=0.
400	ENDIF
401	c condensation or partial sublimation of CO2 ice
402	c """""""""""""""""""""""""""""""""""""""""""""""
403	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
404	& /(latcond*ptimestep) - zfallheat
405	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
406	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
407	#ifdef MESOSCALE
408	print*, "not enough CO2 tracer in 1st layer to condense"
409	print*, ">>> to be implemented in the mesoscale case"
410	print*, "because this uses ap levels..."
411	#else
412	c If there is not enough CO2 tracer in 1st layer to condense
413	c """"""""""""""""""""""""""""""""""""""""""""""""""""""
414	IF(ico2.ne.0) then
415	c Available CO2 tracer in layer 1 at end of timestep (kg/m2)
416	#ifndef MESOSCALE
417	availco2= pq(ig,1,ico2)*((ap(1)-ap(2))+
418	& (bp(1)-bp(2))(pplev(ig,1)/g-zdiceco2(ig)ptimestep))
419	#else
420	availco2 = pq(ig,1,igcm_co2)
421	PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT'
422	& ' CORRECTED FROM SIGMA LEVELS"
423	#endif
424	IF ((zcondices(ig) + zcondicea(ig,1))*ptimestep
425	& .gt.availco2) then
426	zcondices(ig) = availco2/ptimestep -zcondicea(ig,1)
427	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
428	pdtsrfc(ig)=(latcond/pcapcal(ig))*
429	& (zcondices(ig)+zfallheat)
430	ENDIF
431	ENDIF
432	#endif
433
434	c If the entire CO2 ice layer sublimes
435	c """"""""""""""""""""""""""""""""""""""""""""""""""""
436	c (including what has just condensed in the atmosphere)
437	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
438	& -zcondices(ig))THEN
439	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
440	pdtsrfc(ig)=(latcond/pcapcal(ig))*
441	& (zcondices(ig)+zfallheat)
442	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
443	END IF
444
445	c Changing CO2 ice amount and pressure :
446	c """"""""""""""""""""""""""""""""""""
447	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
448	pdpsrf(ig) = -zdiceco2(ig)*g
449
450	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
451	PRINT*,'STOP in condens'
452	PRINT*,'condensing more than total mass'
453	PRINT*,'Grid point ',ig
454	PRINT*,'Longitude(degrees): ',longitude_deg(ig)
455	PRINT*,'Latitude(degrees): ',latitude_deg(ig)
456	PRINT*,'Ps = ',pplev(ig,1)
457	PRINT*,'d Ps = ',pdpsrf(ig)
458	CALL abort_physic('co2condens',
459	& 'condensing more than total mass', 1)
460	ENDIF
461	END IF ! if there is condensation/sublimation
462	ENDDO ! of DO ig=1,ngrid
463	else
464	do ig=1,ngrid
465	if ((ztsrf(ig)<ztcondsol(ig)) .OR. (zdqssed_co2(ig)/=0.) .OR.
466	& ((ztsrf(ig)>ztcondsol(ig)) .AND. (piceco2(ig)/=0.))) then
467	condsub(ig) = .true.
468	zfallheat = 0.
469	!----------------------------------------------------------------------------------------------------------------------!
470	! Compute direct condensation/sublimation of CO2 ice
471	!----------------------------------------------------------------------------------------------------------------------!
472	zcondices(ig) = pcapcal(ig) * (ztcondsol(ig)-ztsrf(ig)) /
473	& (latcond*ptimestep) - zfallheat
474	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig)) / ptimestep
475	!----------------------------------------------------------------------------------------------------------------------!
476	! If there is not enough CO2 tracer in 1st layer to condense
477	!----------------------------------------------------------------------------------------------------------------------!
478	! Available CO2 tracer in layer 1 at end of timestep (kg/m2)
479	#ifndef MESOSCALE
480	availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) +
481	& (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)
482	& *ptimestep) )
483	#else
484	availco2 = pq(ig,1,igcm_co2)
485	PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT CORRECTED FROM SIGMA LEVELS"
486	#endif
487	if ( zcondices(ig) * ptimestep>availco2 ) then
488	zcondices(ig) = availco2/ptimestep
489	zdiceco2(ig) = zcondices(ig) + zdqssed_co2(ig)
490	pdtsrfc(ig) = (latcond/pcapcal(ig)) *
491	& (zcondices(ig)+zfallheat)
492	end if
493	!----------------------------------------------------------------------------------------------------------------------!
494	! If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
495	!----------------------------------------------------------------------------------------------------------------------!
496	if ( (piceco2(ig)/ptimestep) <= -zcondices(ig) ) then
497	zcondices(ig) = -piceco2(ig)/ptimestep
498	pdtsrfc(ig) = (latcond/pcapcal(ig)) *
499	& (zcondices(ig)+zfallheat)
500	end if
501	!----------------------------------------------------------------------------------------------------------------------!
502	! Changing CO2 ice amount and pressure
503	!----------------------------------------------------------------------------------------------------------------------!
504	zdiceco2(ig) = zcondices(ig) + zdqssed_co2(ig) +
505	& condens_column(ig)
506	piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
507	pdpsrf(ig) = -zdiceco2(ig) * g
508
509	if (abs(pdpsrf(ig)*ptimestep)>pplev(ig,1)) then
510	print *, 'STOP in condens'
511	print *, 'condensing more than total mass'
512	print *, 'Grid point ', ig
513	print *, 'Longitude(degrees): ', longitude_deg(ig)
514	print *, 'Latitude(degrees): ', latitude_deg(ig)
515	print *, 'Ps = ', pplev(ig,1)
516	print *, 'd Ps = ', pdpsrf(ig)
517	call abort_physic('co2condens',
518	& 'condensing more than total mass', 1)
519	end if
520	end if
521	end do
522	end if
523	c ********************************************************************
524	c Surface albedo and emissivity of the surface below the snow (emisref)
525	c ********************************************************************
526
527	! Check that amont of CO2 ice is not problematic
528	DO ig=1,ngrid
529	if(.not.piceco2(ig).ge.0.) THEN
530	if(piceco2(ig).le.-5.e-8) print*,
531	$ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig)
532	piceco2(ig)=0.
533	endif
534	ENDDO
535
536	! Set albedo and emissivity of the surface
537	! ----------------------------------------
538	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
539
540	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
541	DO ig=1,ngrid
542	if (piceco2(ig).eq.0) then
543	pemisurf(ig)=emissiv
544	endif
545	ENDDO
546
547	! firstcall2=.false.
548	c ***************************************************************
549	c Correction to account for redistribution between sigma or hybrid
550	c layers when changing surface pressure (and warming/cooling
551	c of the CO2 currently changing phase).
552	c *************************************************************
553
554	DO ig=1,ngrid
555	if (condsub(ig)) then
556	do l=1,nlayer
557	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
558	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
559	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
560	do iq=1,nq
561	zqc(l,iq)=zq(ig,l,iq)+zdq_scav(ig,l,iq)*ptimestep ! zdq_scav=0 if watercloud=false
562	enddo
563	enddo
564
565	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
566	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
567	if (.not. co2clouds) then
568	zmflux(1) = -zcondices(ig)
569	DO l=1,nlayer
570	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
571	#ifndef MESOSCALE
572	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
573	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
574	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) then
575	zmflux(l+1)=0.
576	end if
577	#else
578	zmflux(l+1) = zmflux(l) - zcondicea(ig,l)
579	if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
580	PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
581	& "dPS HAVE TO BE IMPLEMENTED"
582	#endif
583	END DO
584	else
585	zmflux(1) = - zcondices(ig) - zdqssed_co2(ig)
586	do l = 1, nlayer
587	#ifndef MESOSCALE
588	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
589	& + (bp(l)-bp(l+1)) * (-pdpsrf(ig)/g)
590	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
591	if (abs(zmflux(l+1))<1E-13.OR.bp(l+1)==0.) then
592	zmflux(l+1) = 0.
593	end if
594	#else
595	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
596	if (abs(zmflux(l+1))<1E-13) then
597	zmflux(l+1) = 0.
598	end if
599	PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
600	& "dPS HAVE TO BE IMPLEMENTED"
601	#endif
602	end do
603	end if
604	#ifdef MESOSCALE
605	print*,"absurd mass set because bp not available"
606	print*,"TO BE FIXED"
607	#else
608	c Mass of each layer at the end of timestep
609	c -----------------------------------------
610	DO l=1,nlayer
611	masse(l)=( pplev(ig,l) - pplev(ig,l+1) +
612	& (bp(l)-bp(l+1))pdpsrf(ig)ptimestep)/g
613	END DO
614	#endif
615
616	c Corresponding fluxes for T,U,V,Q
617	c """"""""""""""""""""""""""""""""
618
619	c averaging operator for TRANSPORT
620	c """"""""""""""""""""""""""""""""
621	c Value transfert at the surface interface when condensation
622	c sublimation:
623	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
624	zum(1) = 0
625	zvm(1) = 0
626	do iq=1,nq
627	zqm(1,iq)=0. ! most tracer do not condense !
628	enddo
629	c Special case if one of the tracer is CO2 gas
630	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
631
632	c Van Leer scheme:
633	DO l=1,nlayer+1
634	w(l)=-zmflux(l)*ptimestep
635	END DO
636	call vl1d(nlayer,ztc,2.,masse,w,ztm)
637	call vl1d(nlayer,zu ,2.,masse,w,zum)
638	call vl1d(nlayer,zv ,2.,masse,w,zvm)
639	! MVals: loop over the fathers ("peres")
640	do iq=1,nqperes
641	do l=1,nlayer
642	zq1(l)=zqc(l,iq)
643	enddo
644	zqm1(1)=zqm(1,iq)
645	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
646	do l=2,nlayer
647	zqc(l,iq)=zq1(l)
648	zqm(l,iq)=zqm1(l)
649	enddo
650	! MVals: loop over the sons ("fils")
651	if (nqfils(iq).gt.0) then
652	if (iq.eq.igcm_h2o_ice) then
653	iq2=igcm_hdo_ice
654	else if (iq.eq.igcm_h2o_vap) then
655	iq2=igcm_hdo_vap
656	else
657	call abort_physic("co2condens_mod","invalid isotope",1)
658	endif
659	do l=1,nlayer
660	if (zqc(l,iq).gt.qperemin) then
661	Ratio1(l)=zqc(l,iq2)/zqc(l,iq)
662	else
663	Ratio1(l)=0.
664	endif
665	masseq(l)=max(masse(l)*zqc(l,iq),masseqmin)
666	wq(l)=w(l)*zqm(l,iq)
667	enddo
668	Ratiom1(1)=zqm(1,iq2)
669	call vl1d(nlayer,Ratio1,2.,masseq,wq,Ratiom1)
670	zqm(1,iq2)=Ratiom1(1)*zqc(1,iq)
671	do l=2,nlayer
672	zqm(l,iq2)=Ratiom1(l)*zqm(l,iq)
673	enddo
674	endif !if (nqfils(iq).gt.0) then
675	enddo !iq=1,nqperes
676
677	c Surface condensation affects low winds
678	if (zmflux(1).lt.0) then
679	zum(1)= zu(1) * (w(1)/masse(1))
680	zvm(1)= zv(1) * (w(1)/masse(1))
681	if (w(1).gt.masse(1)) then ! ensure numerical stability
682	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
683	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
684	end if
685	end if
686
687	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
688	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
689	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
690	do iq=1,nq
691	zqm(nlayer+1,iq)= zqc(nlayer,iq)
692	enddo
693
694	#ifdef MESOSCALE
695	!!!! AS: This part must be commented in the mesoscale model
696	!!!! AS: ... to avoid instabilities.
697	!!!! AS: you have to compile with -DMESOSCALE to do so
698	#else
699	c Tendencies on T, U, V, Q
700	c """"""""""""""""""""""""
701	DO l=1,nlayer
702	IF(.not. co2clouds) THEN
703	c Tendencies on T
704	zdtsig(ig,l) = (1/masse(l)) *
705	& ( zmflux(l)*(ztm(l) - ztc(l))
706	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
707	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
708	ELSE
709	zdtsig(ig,l) = (1/masse(l)) *
710	& ( zmflux(l)*(ztm(l) - ztc(l))
711	& - zmflux(l+1)*(ztm(l+1) - ztc(l)))
712	ENDIF
713	c D.BARDET: for diagnotics
714	zmflux3D(ig,l)=zmflux(l)
715	ztm3D(ig,l)=ztm(l)
716	ztc3D(ig,l)=ztc(l)
717
718	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
719
720	c Tendencies on U
721	pduc(ig,l) = (1/masse(l)) *
722	& ( zmflux(l)*(zum(l) - zu(l))
723	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
724
725
726	c Tendencies on V
727	pdvc(ig,l) = (1/masse(l)) *
728	& ( zmflux(l)*(zvm(l) - zv(l))
729	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
730
731	END DO
732
733	#endif
734
735	if (.not. co2clouds) then
736	do iq=1,nq
737	! if (noms(iq).eq.'co2') then
738	if (iq.eq.ico2) then
739	c SPECIAL Case when the tracer IS CO2 :
740	DO l=1,nlayer
741	pdqc(ig,l,iq)= (1/masse(l)) *
742	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
743	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
744	& + zcondicea(ig,l)*(zqc(l,iq)-1.) )
745	END DO
746	else
747	DO l=1,nlayer
748	pdqc(ig,l,iq)= (1/masse(l)) *
749	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
750	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
751	& + zcondicea(ig,l)*zqc(l,iq) )
752
753	pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if watercloud=false
754	END DO
755	end if
756	enddo
757	else
758	c Tendencies on Q
759	do iq=1,nq
760	! if (noms(iq).eq.'co2') then
761	if (iq.eq.ico2) then
762	c SPECIAL Case when the tracer IS CO2 :
763	DO l=1,nlayer
764	pdqc(ig,l,iq)= (1/masse(l)) *
765	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
766	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
767	& + condens_layer(ig,l)*(zqc(l,iq)-1.) )
768	END DO
769	else
770	DO l=1,nlayer
771	pdqc(ig,l,iq)= (1/masse(l)) *
772	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
773	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
774	& + condens_layer(ig,l)*zqc(l,iq) )
775	END DO
776	end if
777	enddo
778	end if
779
780
781	end if ! if (condsub)
782	END DO ! loop on ig
783
784	c ***************************************************************
785	c CO2 snow / clouds scheme
786	c ***************************************************************
787	if (.not. co2clouds) then
788	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
789	& zcondicea,zcondices,zfallice,pemisurf)
790	else
791	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
792	& zcondicea,zcondices,zdqssed_co2,pemisurf)
793	end if
794	c ***************************************************************
795	c Ecriture des diagnostiques
796	c ***************************************************************
797
798	c DO l=1,nlayer
799	c DO ig=1,ngrid
800	c Taux de cond en kg.m-2.pa-1.s-1
801	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
802	c Taux de cond en kg.m-3.s-1
803	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
804	c END DO
805	c END DO
806	c call WRITEDIAGFI(ngrid,'tconda1',
807	c &'Taux de condensation CO2 atmospherique /Pa',
808	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
809	c call WRITEDIAGFI(ngrid,'tconda2',
810	c &'Taux de condensation CO2 atmospherique /m',
811	c & 'kg.m-3.s-1',3,tconda2)
812
813	! output falling co2 ice in 1st layer:
814	! call WRITEDIAGFI(ngrid,'fallice',
815	! &'Precipitation of co2 ice',
816	! & 'kg.m-2.s-1',2,zfallice(1,1))
817
818	#ifndef MESOSCALE
819	! Extra special case for surface temperature tendency pdtsrfc:
820	! we want to fix the south pole temperature to CO2 condensation temperature
821	if (caps.and.(obliquit.lt.27.)) then
822	! check if last grid point is the south pole
823	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
824	! NB: Updated surface pressure, at grid point 'ngrid', is
825	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
826	! write(,) "co2condens: South pole: latitude(ngrid)=",
827	! & latitude(ngrid)
828	ztcondsol(ngrid)=
829	& 1./(bcond-acondlog(.01vmr_co2(ngrid,1)*
830	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
831	pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
832	endif
833	endif
834	#endif
835
836	END SUBROUTINE co2condens
837
838
839
840	c *****************************************************************
841	SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
842	c
843	c
844	c Operateur de moyenne inter-couche pour calcul de transport type
845	c Van-Leer " pseudo amont " dans la verticale
846	c q,w sont des arguments d'entree pour le s-pg ....
847	c masse : masse de la couche Dp/g
848	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
849	c pente_max = 2 conseillee
850	c
851	c
852	c --------------------------------------------------------------------
853	IMPLICIT NONE
854
855	c
856	c
857	c
858	c Arguments:
859	c ----------
860	integer nlayer
861	real masse(nlayer),pente_max
862	REAL q(nlayer),qm(nlayer+1)
863	REAL w(nlayer+1)
864	c
865	c Local
866	c ---------
867	c
868	INTEGER l
869	c
870	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
871	real sigw, Mtot, MQtot
872	integer m
873	c integer ismax,ismin
874
875
876	c On oriente tout dans le sens de la pression
877	c W > 0 WHEN DOWN !!!!!!!!!!!!!
878
879	do l=2,nlayer
880	dzqw(l)=q(l-1)-q(l)
881	adzqw(l)=abs(dzqw(l))
882	enddo
883
884	do l=2,nlayer-1
885	if(dzqw(l)*dzqw(l+1).gt.0.) then
886	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
887	else
888	dzq(l)=0.
889	endif
890	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
891	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
892	enddo
893
894	dzq(1)=0.
895	dzq(nlayer)=0.
896
897	do l = 1,nlayer-1
898
899	c Regular scheme (transfered mass < layer mass)
900	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
901	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
902	sigw=w(l+1)/masse(l+1)
903	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
904	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
905	sigw=w(l+1)/masse(l)
906	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
907
908	c Extended scheme (transfered mass > layer mass)
909	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
910	else if(w(l+1).gt.0.) then
911	m=l+1
912	Mtot = masse(m)
913	MQtot = masse(m)*q(m)
914	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
915	m=m+1
916	Mtot = Mtot + masse(m)
917	MQtot = MQtot + masse(m)*q(m)
918	end do
919	if (m.lt.nlayer) then
920	sigw=(w(l+1)-Mtot)/masse(m+1)
921	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
922	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
923	else
924	w(l+1) = Mtot
925	qm(l+1) = Mqtot / Mtot
926	CALL abort_physic('co2condens',
927	& 'top layer is disapearing !', 1)
928	end if
929	else ! if(w(l+1).lt.0)
930	m = l-1
931	Mtot = masse(m+1)
932	MQtot = masse(m+1)*q(m+1)
933	if (m.gt.0) then ! because some compilers will have problems
934	! evaluating masse(0)
935	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
936	m=m-1
937	Mtot = Mtot + masse(m+1)
938	MQtot = MQtot + masse(m+1)*q(m+1)
939	if (m.eq.0) exit
940	end do
941	endif
942	if (m.gt.0) then
943	sigw=(w(l+1)+Mtot)/masse(m)
944	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
945	& (q(m)-0.5(1.+sigw)dzq(m)) )
946	else
947	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
948	end if
949	end if
950	enddo
951
952	c boundary conditions (not used in co2condens !!)
953	c qm(nlayer+1)=0.
954	c if(w(1).gt.0.) then
955	c qm(1)=q(1)
956	c else
957	c qm(1)=0.
958	c end if
959
960	END SUBROUTINE vl1d
961
962	c *****************************************************************
963	SUBROUTINE scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,pq,
964	& rdust,pcondicea,pfallice,pdq_scav,pdqsc)
965
966	c
967	c
968	c Calcul de la quantite de poussiere lessivee par les nuages de CO2
969	c
970	c --------------------------------------------------------------------
971	use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
972	& igcm_dust_mass, igcm_dust_number,
973	& igcm_ccn_mass, igcm_ccn_number,
974	& rho_dust, nuice_sed, nuice_ref,r3n_q
975	use comcstfi_h, only: g
976	use dust_param_mod, only: freedust
977	IMPLICIT NONE
978	include "callkeys.h" ! for the flags water and microphys
979	c
980	c
981	c Arguments:
982	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
983	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
984	INTEGER,INTENT(IN) :: nq ! number of tracers
985	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
986	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
987	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
988	REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius
989	REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
990	REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
991
992	REAL,INTENT(OUT) :: pdq_scav(ngrid,nlayer,nq) ! tendancy due to scavenging by co2
993	REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers
994
995	c Locals:
996	INTEGER l,ig
997	REAL scav_ratio_dust, scav_ratio_wice ! ratio of the dust/water ice mass mixing ratios in condensing CO2 ice and in air
998	REAL scav_dust_mass(nlayer+1) ! dust flux (mass) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
999	REAL scav_dust_number(nlayer+1) ! dust flux (number) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
1000	REAL scav_ccn_mass(nlayer+1) ! ccn flux (mass) scavenged towards the lower layer
1001	REAL scav_ccn_number(nlayer+1) ! ccn flux (number) scavenged towards the lower layer
1002	REAL scav_h2o_ice(nlayer+1) ! water ice flux (mass) scavenged towards the lower layer
1003	REAL massl ! mass of the layer l at point ig (kg/m2)
1004
1005	c Initialization:
1006	scav_ratio_dust = 100 !1 !10 !100 !1000
1007	scav_ratio_wice = scav_ratio_dust
1008	pdq_scav(:,:,:)=0.
1009
1010	DO ig=1,ngrid
1011	scav_dust_mass(nlayer+1)=0.
1012	scav_dust_number(nlayer+1)=0.
1013	scav_ccn_mass(nlayer+1)=0.
1014	scav_ccn_number(nlayer+1)=0.
1015	scav_h2o_ice(nlayer+1)=0.
1016
1017	DO l=nlayer , 1, -1
1018	massl=(pplev(ig,l)-pplev(ig,l+1))/g
1019	IF(pcondicea(ig,l).GT.0.)THEN ! if CO2 condenses and traps dust/water ice
1020	! Calculation of the tendencies
1021	if (freedust) then
1022	pdq_scav(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)
1023	& /ptimestep*(1-exp(
1024	& -scav_ratio_dustpcondicea(ig,l)ptimestep/massl))
1025
1026	pdq_scav(ig,l,igcm_dust_number)=pdq_scav(ig,l,igcm_dust_mass)
1027	& *r3n_q/rdust(ig,l)
1028	endif
1029	if (freedust.AND.microphys) then
1030	pdq_scav(ig,l,igcm_ccn_mass)=-pq(ig,l,igcm_ccn_mass)
1031	& /ptimestep*(1-exp(
1032	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
1033	pdq_scav(ig,l,igcm_ccn_number)=pdq_scav(ig,l,igcm_ccn_mass)
1034	& *r3n_q/rdust(ig,l)
1035	endif
1036	if (water) then
1037	pdq_scav(ig,l,igcm_h2o_ice)=-pq(ig,l,igcm_h2o_ice)
1038	& /ptimestep*(1-exp(
1039	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
1040	endif
1041
1042	ELSE IF(pcondicea(ig,l).LT.0.)THEN ! if CO2 sublimates and releases dust/water ice
1043	! Calculation of the tendencies
1044	if (freedust) then
1045	pdq_scav(ig,l,igcm_dust_mass)=-pcondicea(ig,l)/massl*
1046	& scav_dust_mass(l+1)/pfallice(ig,l+1)
1047
1048	pdq_scav(ig,l,igcm_dust_number)=-pcondicea(ig,l)/massl*
1049	& scav_dust_number(l+1)/pfallice(ig,l+1)
1050	endif
1051	if (freedust.AND.microphys) then
1052	pdq_scav(ig,l,igcm_ccn_mass)=-pcondicea(ig,l)/massl*
1053	& scav_ccn_mass(l+1)/pfallice(ig,l+1)
1054
1055	pdq_scav(ig,l,igcm_ccn_number)=-pcondicea(ig,l)/massl*
1056	& scav_ccn_number(l+1)/pfallice(ig,l+1)
1057	endif
1058	if (water) then
1059	pdq_scav(ig,l,igcm_h2o_ice)=-pcondicea(ig,l)/massl*
1060	& scav_h2o_ice(l+1)/pfallice(ig,l+1)
1061	endif
1062
1063	END IF
1064	! Calculation of the scavenged dust/wice flux towards the lower layers
1065	if (freedust) then
1066	scav_dust_mass(l)=-pdq_scav(ig,l,igcm_dust_mass)*massl
1067	& +scav_dust_mass(l+1)
1068
1069	scav_dust_number(l)=-pdq_scav(ig,l,igcm_dust_number)*massl
1070	& +scav_dust_number(l+1)
1071	endif
1072	if (freedust.AND.microphys) then
1073	scav_ccn_mass(l)=-pdq_scav(ig,l,igcm_ccn_mass)*massl
1074	& +scav_ccn_mass(l+1)
1075
1076	scav_ccn_number(l)=-pdq_scav(ig,l,igcm_ccn_number)*massl
1077	& +scav_dust_number(l+1)
1078	endif
1079	if (water) then
1080	scav_h2o_ice(l)=-pdq_scav(ig,l,igcm_h2o_ice)*massl
1081	& +scav_h2o_ice(l+1)
1082	endif
1083
1084	ENDDO
1085	! Calculation of the surface tendencies
1086	pdqsc(ig,igcm_dust_mass)=0.
1087	pdqsc(ig,igcm_dust_number)=0.
1088
1089	if (freedust) then
1090	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
1091	& +scav_dust_mass(1)
1092	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
1093	& +scav_dust_number(1)
1094	endif
1095	if (freedust.AND.microphys) then
1096	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
1097	& +scav_ccn_mass(1)
1098	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
1099	& +scav_ccn_number(1)
1100	endif
1101	if (water) then
1102	pdqsc(ig,igcm_h2o_ice)=scav_h2o_ice(1)
1103	endif
1104	ENDDO
1105
1106	END SUBROUTINE scavenging_by_co2
1107
1108	END MODULE co2condens_mod

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

Download in other formats:

Original Format