Context Navigation

co2condens_mod.F @ 2932

Last change on this file since 2932 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

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

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

Download in other formats:

Original Format