Context Navigation

co2condens_mod.F @ 2685

Last change on this file since 2685 was 2669, checked in by abierjon, 3 years ago

Mars GCM:
Put the CO2 ice scavenging ratio for dust and water ice to 20 instead of 100, to prevent drying of the water cycle.

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

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

Download in other formats:

Original Format