Context Navigation

co2condens_mod.F @ 3126

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

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

Download in other formats:

Original Format