Context Navigation

co2condens_mod.F @ 2997

Last change on this file since 2997 was 2977, checked in by llange, 19 months ago
Mars PCM Change the latent heat when computing the amount of CO2 that condenses on a slope: a warm air must be cooled off before condensing on the slope. Adapted from pluton PCM (see Eq.12 from Forget et al., 2017) LL
File size: 41.9 KB

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

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

Download in other formats:

Original Format