Context Navigation

co2condens_mod.F @ 2393

Last change on this file since 2393 was 2332, checked in by mvals, 5 years ago

Mars GCM:
Follow-up of the last commit for the transport of the isotopic ratio: simplification of the transmission of variables from the dynamics to the
physics:

libf/dynphy_lonlat/phymars/: iniphysiq_mod.F90: transmission of the content of 2 variables describing the isotopes instead of 4 (nqperes: number of tracers "peres", nqfils:

number of tracers "fils")

libf/phymars/: phys_state_var_init_mod.F90, tracer_mod.F: idem callsedim_mod.F: idem co2condens_mod.F: idem
libf/phymars/dyn1d: testphys1d.F: idem (the reading interface for traceur.def has been completed to fill the variables nqperes and nqfils).

File size: 36.3 KB

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

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

Download in other formats:

Original Format