Context Navigation

co2condens_mod.F @ 2953

Last change on this file since 2953 was 2953, checked in by romain.vande, 19 months ago
Mars PCM : Adapt vdifc, co2condens, rocketduststorm and topmons routines to the subslope parametrisation. RV
File size: 42.0 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	c & /((latcond+cpp(zt(ig,1)-ztcondsol(ig)))ptimestep)
446	& /(latcond*ptimestep)
447	& - zfallheat
448	else
449	c sublimation
450	zcondices(ig,islope)=pcapcal(ig,islope)
451	& *(ztcondsol(ig)-ztsrf(ig,islope))
452	& /(latcond*ptimestep)
453	& - zfallheat
454	endif
455	pdtsrfc(ig,islope) = (ztcondsol(ig) - ztsrf(ig,islope))
456	& /ptimestep
457	#ifdef MESOSCALE
458	print*, "not enough CO2 tracer in 1st layer to condense"
459	print*, ">>> to be implemented in the mesoscale case"
460	print*, "because this uses ap levels..."
461	#else
462	c If there is not enough CO2 tracer in 1st layer to condense
463	c """"""""""""""""""""""""""""""""""""""""""""""""""""""
464	IF(ico2.ne.0) then
465	c Available CO2 tracer in layer 1 at end of timestep (kg/m2)
466	#ifndef MESOSCALE
467	availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+
468	& (bp(1)-bp(2))*(pplev(ig,1)/g -
469	& (zcondices(ig,islope) + zfallice(ig,1))
470	& *ptimestep))
471	if ((zcondices(ig,islope) + condens_layer(ig,1))
472	& *ptimestep
473	& .gt.availco2) then
474	zcondices(ig,islope) = availco2/ptimestep -
475	& condens_layer(ig,1)
476	pdtsrfc(ig,islope) = (latcond/pcapcal(ig,islope))*
477	& (zcondices(ig,islope)+zfallheat)
478	end if
479	#else
480	availco2 = pq(ig,1,igcm_co2)
481	PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT'
482	& ' CORRECTED FROM SIGMA LEVELS"
483	#endif
484	ENDIF
485	#endif
486
487	c If the entire CO2 ice layer sublimes on the slope
488	c """"""""""""""""""""""""""""""""""""""""""""""""""""
489	c (including what has just condensed in the atmosphere)
490	if (co2clouds) then
491	IF((piceco2(ig,islope)/ptimestep).LE.
492	& -zcondices(ig,islope))THEN
493	zcondices(ig,islope) = -piceco2(ig,islope)/ptimestep
494	pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) *
495	& (zcondices(ig,islope)+zfallheat)
496	END IF
497	else
498	IF((piceco2(ig,islope)/ptimestep+zfallice(ig,1)).LE.
499	& -zcondices(ig,islope))THEN
500	zcondices(ig,islope) = -piceco2(ig,islope)
501	& /ptimestep - zfallice(ig,1)
502	pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) *
503	& (zcondices(ig,islope)+zfallheat)
504	END IF
505	end if
506
507	c Changing CO2 ice amount and pressure per slope:
508	c """"""""""""""""""""""""""""""""""""
509	zdiceco2(ig,islope) = zcondices(ig,islope)+zfallice(ig,1)
510	& + condens_column(ig)
511	if (co2clouds) then
512	! add here only direct condensation/sublimation
513	piceco2(ig,islope) = piceco2(ig,islope) +
514	& zcondices(ig,islope)*ptimestep
515	else
516	! add here also CO2 ice in the atmosphere
517	piceco2(ig,islope) = piceco2(ig,islope) +
518	& zdiceco2(ig,islope)*ptimestep
519	end if
520
521	zcondices_mesh_avg(ig) = zcondices_mesh_avg(ig) +
522	& zcondices(ig,islope)* subslope_dist(ig,islope)
523
524	zdiceco2_mesh_avg(ig) = zdiceco2_mesh_avg(ig) +
525	& zdiceco2(ig,islope)* subslope_dist(ig,islope)
526
527	END IF ! if there is condensation/sublimation
528
529	piceco2(ig,islope) = piceco2(ig,islope)
530	& * cos(pi*def_slope_mean(islope)/180.)
531
532	ENDDO !islope
533
534	pdpsrf(ig) = -zdiceco2_mesh_avg(ig)*g
535
536	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
537	PRINT*,'STOP in condens'
538	PRINT*,'condensing more than total mass'
539	PRINT*,'Grid point ',ig
540	PRINT*,'Longitude(degrees): ',longitude_deg(ig)
541	PRINT*,'Latitude(degrees): ',latitude_deg(ig)
542	PRINT*,'Ps = ',pplev(ig,1)
543	PRINT*,'d Ps = ',pdpsrf(ig)
544	CALL abort_physic('co2condens',
545	& 'condensing more than total mass', 1)
546	ENDIF
547
548	ENDDO ! of DO ig=1,ngrid
549
550
551	c ********************************************************************
552	c Surface albedo and emissivity of the surface below the snow (emisref)
553	c ********************************************************************
554
555	! Check that amont of CO2 ice is not problematic
556	DO ig=1,ngrid
557	DO islope = 1,nslope
558	if(.not.piceco2(ig,islope).ge.0.) THEN
559	if(piceco2(ig,islope).le.-5.e-8) print*,
560	$ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig,islope)
561	piceco2(ig,islope)=0.
562	endif
563	ENDDO
564	ENDDO
565
566	! Set albedo and emissivity of the surface
567	! ----------------------------------------
568	DO islope = 1,nslope
569	piceco2_tmp(:) = piceco2(:,islope)
570	alb_tmp(:,:) = psolaralb(:,:,islope)
571	emisref_tmp(:) = 0.
572	CALL albedocaps(zls,ngrid,piceco2_tmp,alb_tmp,emisref_tmp)
573	psolaralb(:,1,islope) = alb_tmp(:,1)
574	psolaralb(:,2,islope) = alb_tmp(:,2)
575	emisref(:,islope) = emisref_tmp(:)
576	ENDDO
577
578	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
579	DO ig=1,ngrid
580	DO islope = 1,nslope
581	if (piceco2(ig,islope).eq.0) then
582	pemisurf(ig,islope)=emissiv
583	endif
584	ENDDO
585	ENDDO
586
587	! firstcall2=.false.
588	c ***************************************************************
589	c Correction to account for redistribution between sigma or hybrid
590	c layers when changing surface pressure (and warming/cooling
591	c of the CO2 currently changing phase).
592	c *************************************************************
593
594	DO ig=1,ngrid
595	if (any(condsub(ig,:))) then
596	do l=1,nlayer
597	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
598	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
599	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
600	do iq=1,nq
601	zqc(l,iq)=zq(ig,l,iq)+zdq_scav(ig,l,iq)*ptimestep ! zdq_scav=0 if co2clouds=true
602	enddo
603	enddo
604
605	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
606	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
607	zmflux(1) = -zcondices_mesh_avg(ig) - zdqssed_co2(ig)
608	DO l=1,nlayer
609	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
610	#ifndef MESOSCALE
611	& + (bp(l)-bp(l+1))*(-pdpsrf(ig)/g)
612	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
613	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) then
614	zmflux(l+1)=0.
615	end if
616	#else
617	zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
618	if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
619	PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
620	& "dPS HAVE TO BE IMPLEMENTED"
621	#endif
622	END DO
623	#ifdef MESOSCALE
624	print*,"absurd mass set because bp not available"
625	print*,"TO BE FIXED"
626	#else
627	c Mass of each layer at the end of timestep
628	c -----------------------------------------
629	DO l=1,nlayer
630	masse(l)=( pplev(ig,l) - pplev(ig,l+1) +
631	& (bp(l)-bp(l+1))pdpsrf(ig)ptimestep)/g
632	END DO
633	#endif
634
635	c Corresponding fluxes for T,U,V,Q
636	c """"""""""""""""""""""""""""""""
637
638	c averaging operator for TRANSPORT
639	c """"""""""""""""""""""""""""""""
640	c Value transfert at the surface interface when condensation
641	c sublimation:
642	ztm(1) = ztcondsol(ig)
643	zum(1) = 0
644	zvm(1) = 0
645	do iq=1,nq
646	zqm(1,iq)=0. ! most tracer do not condense !
647	enddo
648	c Special case if one of the tracer is CO2 gas
649	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
650
651	c Van Leer scheme:
652	DO l=1,nlayer+1
653	w(l)=-zmflux(l)*ptimestep
654	END DO
655	call vl1d(nlayer,ztc,2.,masse,w,ztm)
656	call vl1d(nlayer,zu ,2.,masse,w,zum)
657	call vl1d(nlayer,zv ,2.,masse,w,zvm)
658	! MVals: loop over the fathers ("peres")
659	do iq=1,nqperes
660	do l=1,nlayer
661	zq1(l)=zqc(l,iq)
662	enddo
663	zqm1(1)=zqm(1,iq)
664	call vl1d(nlayer,zq1,2.,masse,w,zqm1)
665	do l=2,nlayer
666	zqc(l,iq)=zq1(l)
667	zqm(l,iq)=zqm1(l)
668	enddo
669	! MVals: loop over the sons ("fils")
670	if (nqfils(iq).gt.0) then
671	if (iq.eq.igcm_h2o_ice) then
672	iq2=igcm_hdo_ice
673	else if (iq.eq.igcm_h2o_vap) then
674	iq2=igcm_hdo_vap
675	else
676	call abort_physic("co2condens_mod","invalid isotope",1)
677	endif
678	do l=1,nlayer
679	if (zqc(l,iq).gt.qperemin) then
680	Ratio1(l)=zqc(l,iq2)/zqc(l,iq)
681	else
682	Ratio1(l)=0.
683	endif
684	masseq(l)=max(masse(l)*zqc(l,iq),masseqmin)
685	wq(l)=w(l)*zqm(l,iq)
686	enddo
687	Ratiom1(1)=zqm(1,iq2)
688	call vl1d(nlayer,Ratio1,2.,masseq,wq,Ratiom1)
689	zqm(1,iq2)=Ratiom1(1)*zqc(1,iq)
690	do l=2,nlayer
691	zqm(l,iq2)=Ratiom1(l)*zqm(l,iq)
692	enddo
693	endif !if (nqfils(iq).gt.0) then
694	enddo !iq=1,nqperes
695
696	c Surface condensation affects low winds
697	if (zmflux(1).lt.0) then
698	zum(1)= zu(1) * (w(1)/masse(1))
699	zvm(1)= zv(1) * (w(1)/masse(1))
700	if (w(1).gt.masse(1)) then ! ensure numerical stability
701	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
702	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
703	end if
704	end if
705
706	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
707	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
708	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
709	do iq=1,nq
710	zqm(nlayer+1,iq)= zqc(nlayer,iq)
711	enddo
712
713	#ifdef MESOSCALE
714	!!!! AS: This part must be commented in the mesoscale model
715	!!!! AS: ... to avoid instabilities.
716	!!!! AS: you have to compile with -DMESOSCALE to do so
717	#else
718	c Tendencies on T, U, V, Q
719	c """"""""""""""""""""""""
720	DO l=1,nlayer
721	IF(.not. co2clouds) THEN
722	c Tendencies on T
723	zdtsig(ig,l) = (1/masse(l)) *
724	& ( zmflux(l)*(ztm(l) - ztc(l))
725	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
726	& + condens_layer(ig,l)*(ztcond(ig,l)-ztc(l)) )
727	ELSE
728	zdtsig(ig,l) = (1/masse(l)) *
729	& ( zmflux(l)*(ztm(l) - ztc(l))
730	& - zmflux(l+1)*(ztm(l+1) - ztc(l)))
731	ENDIF
732	c D.BARDET: for diagnotics
733	zmflux3D(ig,l)=zmflux(l)
734	ztm3D(ig,l)=ztm(l)
735	ztc3D(ig,l)=ztc(l)
736
737	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
738
739	c Tendencies on U
740	pduc(ig,l) = (1/masse(l)) *
741	& ( zmflux(l)*(zum(l) - zu(l))
742	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
743
744
745	c Tendencies on V
746	pdvc(ig,l) = (1/masse(l)) *
747	& ( zmflux(l)*(zvm(l) - zv(l))
748	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
749
750	END DO
751
752	#endif
753
754	do iq=1,nq
755	! if (noms(iq).eq.'co2') then
756	if (iq.eq.ico2) then
757	c SPECIAL Case when the tracer IS CO2 :
758	DO l=1,nlayer
759	pdqc(ig,l,iq)= (1/masse(l)) *
760	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
761	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
762	& + condens_layer(ig,l)*(zqc(l,iq)-1.) )
763	END DO
764	else
765	DO l=1,nlayer
766	pdqc(ig,l,iq)= (1/masse(l)) *
767	& ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
768	& - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
769	& + condens_layer(ig,l)*zqc(l,iq) )
770
771	pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if co2clouds=true
772	END DO
773	end if
774	enddo
775
776	end if ! if (condsub)
777	END DO ! loop on ig
778
779	c ***************************************************************
780	c CO2 snow / clouds scheme
781	c ***************************************************************
782	DO islope = 1,nslope
783	emisref_tmp(:) = emisref(:,islope)
784	condsub_tmp(:) = condsub(:,islope)
785	condens_layer_tmp(:,:) = condens_layer(:,:)*
786	& cos(pi*def_slope_mean(islope)/180.)
787	zcondices_tmp(:) = zcondices(:,islope)*
788	& cos(pi*def_slope_mean(islope)/180.)
789	zfallice_tmp(:,:) = zfallice(:,:)*
790	& cos(pi*def_slope_mean(islope)/180.)
791	pemisurf_tmp(:) = pemisurf(:,islope)
792
793	call co2snow(ngrid,nlayer,ptimestep,emisref_tmp,condsub_tmp,
794	& pplev,condens_layer_tmp,zcondices_tmp,zfallice_tmp,
795	& pemisurf_tmp)
796	pemisurf(:,islope) = pemisurf_tmp(:)
797
798	ENDDO
799	c ***************************************************************
800	c Ecriture des diagnostiques
801	c ***************************************************************
802
803	c DO l=1,nlayer
804	c DO ig=1,ngrid
805	c Taux de cond en kg.m-2.pa-1.s-1
806	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
807	c Taux de cond en kg.m-3.s-1
808	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
809	c END DO
810	c END DO
811	c call write_output('tconda1',
812	c &'Taux de condensation CO2 atmospherique /Pa',
813	c & 'kg.m-2.Pa-1.s-1',tconda1)
814	c call write_output('tconda2',
815	c &'Taux de condensation CO2 atmospherique /m',
816	c & 'kg.m-3.s-1',tconda2)
817
818	! output falling co2 ice in 1st layer:
819	! call write_output('fallice',
820	! &'Precipitation of co2 ice',
821	! & 'kg.m-2.s-1',zfallice(1,1))
822
823	#ifndef MESOSCALE
824	! Extra special case for surface temperature tendency pdtsrfc:
825	! we want to fix the south pole temperature to CO2 condensation temperature
826	if (caps.and.(obliquit.lt.27.)) then
827	! check if last grid point is the south pole
828	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
829	! NB: Updated surface pressure, at grid point 'ngrid', is
830	! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
831	ztcondsol(ngrid)=
832	& 1./(bcond-acondlog(.01vmr_co2(ngrid,1)*
833	& (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
834	DO islope = 1,nslope
835	pdtsrfc(ngrid,islope)=(ztcondsol(ngrid)-
836	& ztsrf(ngrid,islope))/ptimestep
837	ENDDO ! islope = 1,nslope
838	endif
839	endif
840	#endif
841
842	END SUBROUTINE co2condens
843
844
845
846	c *****************************************************************
847	SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
848	c
849	c
850	c Operateur de moyenne inter-couche pour calcul de transport type
851	c Van-Leer " pseudo amont " dans la verticale
852	c q,w sont des arguments d'entree pour le s-pg ....
853	c masse : masse de la couche Dp/g
854	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
855	c pente_max = 2 conseillee
856	c
857	c
858	c --------------------------------------------------------------------
859	IMPLICIT NONE
860
861	c
862	c
863	c
864	c Arguments:
865	c ----------
866	integer nlayer
867	real masse(nlayer),pente_max
868	REAL q(nlayer),qm(nlayer+1)
869	REAL w(nlayer+1)
870	c
871	c Local
872	c ---------
873	c
874	INTEGER l
875	c
876	real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
877	real sigw, Mtot, MQtot
878	integer m
879	c integer ismax,ismin
880
881
882	c On oriente tout dans le sens de la pression
883	c W > 0 WHEN DOWN !!!!!!!!!!!!!
884
885	do l=2,nlayer
886	dzqw(l)=q(l-1)-q(l)
887	adzqw(l)=abs(dzqw(l))
888	enddo
889
890	do l=2,nlayer-1
891	if(dzqw(l)*dzqw(l+1).gt.0.) then
892	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
893	else
894	dzq(l)=0.
895	endif
896	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
897	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
898	enddo
899
900	dzq(1)=0.
901	dzq(nlayer)=0.
902
903	do l = 1,nlayer-1
904
905	c Regular scheme (transfered mass < layer mass)
906	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
907	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
908	sigw=w(l+1)/masse(l+1)
909	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
910	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
911	sigw=w(l+1)/masse(l)
912	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
913
914	c Extended scheme (transfered mass > layer mass)
915	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
916	else if(w(l+1).gt.0.) then
917	m=l+1
918	Mtot = masse(m)
919	MQtot = masse(m)*q(m)
920	do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
921	m=m+1
922	Mtot = Mtot + masse(m)
923	MQtot = MQtot + masse(m)*q(m)
924	end do
925	if (m.lt.nlayer) then
926	sigw=(w(l+1)-Mtot)/masse(m+1)
927	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
928	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
929	else
930	w(l+1) = Mtot
931	qm(l+1) = Mqtot / Mtot
932	CALL abort_physic('co2condens',
933	& 'top layer is disapearing !', 1)
934	end if
935	else ! if(w(l+1).lt.0)
936	m = l-1
937	Mtot = masse(m+1)
938	MQtot = masse(m+1)*q(m+1)
939	if (m.gt.0) then ! because some compilers will have problems
940	! evaluating masse(0)
941	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
942	m=m-1
943	Mtot = Mtot + masse(m+1)
944	MQtot = MQtot + masse(m+1)*q(m+1)
945	if (m.eq.0) exit
946	end do
947	endif
948	if (m.gt.0) then
949	sigw=(w(l+1)+Mtot)/masse(m)
950	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
951	& (q(m)-0.5(1.+sigw)dzq(m)) )
952	else
953	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
954	end if
955	end if
956	enddo
957
958	c boundary conditions (not used in co2condens !!)
959	c qm(nlayer+1)=0.
960	c if(w(1).gt.0.) then
961	c qm(1)=q(1)
962	c else
963	c qm(1)=0.
964	c end if
965
966	END SUBROUTINE vl1d
967
968	c *****************************************************************
969	SUBROUTINE scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,pq,
970	& rdust,pcondicea,pfallice,pdq_scav,pdqsc)
971
972	c
973	c
974	c Calcul de la quantite de poussiere lessivee par les nuages de CO2
975	c
976	c --------------------------------------------------------------------
977	use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
978	& igcm_dust_mass, igcm_dust_number,
979	& igcm_ccn_mass, igcm_ccn_number,
980	& rho_dust, nuice_sed, nuice_ref,r3n_q
981	use comcstfi_h, only: g
982	use dust_param_mod, only: freedust
983	IMPLICIT NONE
984	include "callkeys.h" ! for the flags water and microphys
985	c
986	c
987	c Arguments:
988	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
989	INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
990	INTEGER,INTENT(IN) :: nq ! number of tracers
991	REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
992	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
993	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
994	REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius
995	REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s)
996	REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
997
998	REAL,INTENT(OUT) :: pdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2
999	REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers
1000
1001	c Locals:
1002	INTEGER l,ig
1003	REAL scav_ratio_dust, scav_ratio_wice ! ratio of the dust/water ice mass mixing ratios in condensing CO2 ice and in air
1004	REAL scav_dust_mass(nlayer+1) ! dust flux (mass) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
1005	REAL scav_dust_number(nlayer+1) ! dust flux (number) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD)
1006	REAL scav_ccn_mass(nlayer+1) ! ccn flux (mass) scavenged towards the lower layer
1007	REAL scav_ccn_number(nlayer+1) ! ccn flux (number) scavenged towards the lower layer
1008	REAL scav_h2o_ice(nlayer+1) ! water ice flux (mass) scavenged towards the lower layer
1009	REAL massl ! mass of the layer l at point ig (kg/m2)
1010
1011	c Initialization:
1012	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
1013	scav_ratio_wice = scav_ratio_dust ! while not drying up the water cycle (which occurs at scav_ratio_wice values above 50 at least)
1014	pdq_scav(:,:,:)=0.
1015	pdqsc(:,:)=0.
1016
1017	DO ig=1,ngrid
1018	scav_dust_mass(nlayer+1)=0.
1019	scav_dust_number(nlayer+1)=0.
1020	scav_ccn_mass(nlayer+1)=0.
1021	scav_ccn_number(nlayer+1)=0.
1022	scav_h2o_ice(nlayer+1)=0.
1023
1024	DO l=nlayer , 1, -1
1025	massl=(pplev(ig,l)-pplev(ig,l+1))/g
1026	IF(pcondicea(ig,l).GT.0.)THEN ! if CO2 condenses and traps dust/water ice
1027	! Calculation of the tendencies
1028	if (freedust) then
1029	pdq_scav(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)
1030	& /ptimestep*(1-exp(
1031	& -scav_ratio_dustpcondicea(ig,l)ptimestep/massl))
1032
1033	pdq_scav(ig,l,igcm_dust_number)=pdq_scav(ig,l,igcm_dust_mass)
1034	& *r3n_q/rdust(ig,l)
1035	endif
1036	if (freedust.AND.microphys) then
1037	pdq_scav(ig,l,igcm_ccn_mass)=-pq(ig,l,igcm_ccn_mass)
1038	& /ptimestep*(1-exp(
1039	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
1040	pdq_scav(ig,l,igcm_ccn_number)=pdq_scav(ig,l,igcm_ccn_mass)
1041	& *r3n_q/rdust(ig,l)
1042	endif
1043	if (water) then
1044	pdq_scav(ig,l,igcm_h2o_ice)=-pq(ig,l,igcm_h2o_ice)
1045	& /ptimestep*(1-exp(
1046	& -scav_ratio_wicepcondicea(ig,l)ptimestep/massl))
1047	endif
1048
1049	ELSE IF(pcondicea(ig,l).LT.0.)THEN ! if CO2 sublimates and releases dust/water ice
1050	! Calculation of the tendencies
1051	if (freedust) then
1052	pdq_scav(ig,l,igcm_dust_mass)=-pcondicea(ig,l)/massl*
1053	& scav_dust_mass(l+1)/pfallice(ig,l+1)
1054
1055	pdq_scav(ig,l,igcm_dust_number)=-pcondicea(ig,l)/massl*
1056	& scav_dust_number(l+1)/pfallice(ig,l+1)
1057	endif
1058	if (freedust.AND.microphys) then
1059	pdq_scav(ig,l,igcm_ccn_mass)=-pcondicea(ig,l)/massl*
1060	& scav_ccn_mass(l+1)/pfallice(ig,l+1)
1061
1062	pdq_scav(ig,l,igcm_ccn_number)=-pcondicea(ig,l)/massl*
1063	& scav_ccn_number(l+1)/pfallice(ig,l+1)
1064	endif
1065	if (water) then
1066	pdq_scav(ig,l,igcm_h2o_ice)=-pcondicea(ig,l)/massl*
1067	& scav_h2o_ice(l+1)/pfallice(ig,l+1)
1068	endif
1069
1070	END IF
1071	! Calculation of the scavenged dust/wice flux towards the lower layers
1072	if (freedust) then
1073	scav_dust_mass(l)=-pdq_scav(ig,l,igcm_dust_mass)*massl
1074	& +scav_dust_mass(l+1)
1075
1076	scav_dust_number(l)=-pdq_scav(ig,l,igcm_dust_number)*massl
1077	& +scav_dust_number(l+1)
1078	endif
1079	if (freedust.AND.microphys) then
1080	scav_ccn_mass(l)=-pdq_scav(ig,l,igcm_ccn_mass)*massl
1081	& +scav_ccn_mass(l+1)
1082
1083	scav_ccn_number(l)=-pdq_scav(ig,l,igcm_ccn_number)*massl
1084	& +scav_dust_number(l+1)
1085	endif
1086	if (water) then
1087	scav_h2o_ice(l)=-pdq_scav(ig,l,igcm_h2o_ice)*massl
1088	& +scav_h2o_ice(l+1)
1089	endif
1090
1091	ENDDO
1092	! Calculation of the surface tendencies
1093	if (freedust) then
1094	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
1095	& +scav_dust_mass(1)
1096	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
1097	& +scav_dust_number(1)
1098	endif
1099	if (freedust.AND.microphys) then
1100	pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass)
1101	& +scav_ccn_mass(1)
1102	pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number)
1103	& +scav_ccn_number(1)
1104	endif
1105	if (water) then
1106	pdqsc(ig,igcm_h2o_ice)=scav_h2o_ice(1)
1107	endif
1108
1109	ENDDO ! loop on ngrid
1110
1111	END SUBROUTINE scavenging_by_co2
1112
1113	END MODULE co2condens_mod

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

Download in other formats:

Original Format