Context Navigation

co2condens_mod.F @ 2141

Last change on this file since 2141 was 2124, checked in by emillour, 6 years ago
Mars GCM: Updated co2condens to correctly conserve tracer mass EM+FF
File size: 27.5 KB

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

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

Download in other formats:

Original Format