Context Navigation

co2condens_mod.F @ 2740

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

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

Download in other formats:

Original Format