Context Navigation

co2condens_mod.F @ 3493

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

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

Download in other formats:

Original Format