Context Navigation

co2condens_mod.F @ 3026

Last change on this file since 3026 was 3013, checked in by llange, 18 months ago

MARS PCM
CO2 condens: Fixing bug introduced in r2977: as the cooling of the air before condensing is now done in the surface condensation loop, the boundary conditions for the VL scheme must be changed in order to avoid duplicating the cooling.

Also fixed a bug in surfini.F (2892): qsurf is now ngrid x nq x nslope: I added the dimension nslope.

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

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

Download in other formats:

Original Format