Context Navigation

co2condens_mod.F @ 2903

Last change on this file since 2903 was 2903, checked in by llange, 21 months ago
PCM Fix a bug in CO2condens. Surface temperature used at the surface as a boundary condition for the VL scheme must be ztcondsol. I've directly put ztcondsol rather that ptsurf+pdtsurfc*timestep. Remove stuff about the solar albedo dependant flux to r2901, sorry about the mistake LL & FF
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
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 WRITEdiagfi(ngrid,"zfallice",
335	& "",
336	& " ",2,zfallice(ngrid,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 WRITEdiagfi(ngrid,"zfallice",
352	& "",
353	& " ",2,zdqssed_co2) ! otherwise we have not 1 day(1proc) = 1 day (24procs) test
354	ENDIF ! end of if co2clouds
355
356	call WRITEdiagfi(ngrid,"pdtc_atm",
357	& "temperature tendency due to CO2 condensation",
358	& " ",3,pdtc)
359
360	call WRITEdiagfi(ngrid,"condens_layer",
361	& "",
362	& " ",3,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 WRITEDIAGFI(ngrid,'tconda1',
719	c &'Taux de condensation CO2 atmospherique /Pa',
720	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
721	c call WRITEDIAGFI(ngrid,'tconda2',
722	c &'Taux de condensation CO2 atmospherique /m',
723	c & 'kg.m-3.s-1',3,tconda2)
724
725	! output falling co2 ice in 1st layer:
726	! call WRITEDIAGFI(ngrid,'fallice',
727	! &'Precipitation of co2 ice',
728	! & 'kg.m-2.s-1',2,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