Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

newcondens.F @ 134

Last change on this file since 134 was 57, checked in by aslmd, 14 years ago
mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 26.4 KB

Line
1	SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep,
2	$ pcapcal,pplay,pplev,ptsrf,pt,
3	$ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
4	$ piceco2,psolaralb,pemisurf,
5	$ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
6	$ fluxsurf_sw)
7
8	IMPLICIT NONE
9	c=======================================================================
10	c subject:
11	c --------
12	c Condensation/sublimation of CO2 ice on the ground and in the
13	c atmosphere
14	c (Scheme described in Forget et al., Icarus, 1998)
15	c
16	c author: Francois Forget 1994-1996
17	c ------
18	c
19	c input:
20	c -----
21	c ngrid nombre de points de verticales
22	c (toutes les boucles de la physique sont au
23	c moins vectorisees sur ngrid)
24	c nlayer nombre de couches
25	c pplay(ngrid,nlayer) Pressure levels
26	c pplev(ngrid,nlayer+1) Pressure levels
27	c pt(ngrid,nlayer) temperature (en K)
28	c ptsrf(ngrid) temperature de surface
29	c
30	c \
31	c pdt(ngrid,nlayermx) \ derivee temporelle physique avant condensation
32	c / ou sublimation pour pt,ptsrf
33	c pdtsrf(ngrid) /
34	c
35	c output:
36	c -------
37	c
38	c pdpsrf(ngrid) \ derivee temporelle physique (contribution de
39	c pdtc(ngrid,nlayermx) / la condensation ou sublimation) pour Ps,pt,ptsrf
40	c pdtsrfc(ngrid) /
41	c
42	c Entree/sortie
43	c -------------
44	c
45	c piceco2(ngrid) : quantite de glace co2 au sol (kg/m2)
46	c psolaralb(ngrid,2) : albedo au sol
47	c pemisurf(ngrid) : emissivite du sol
48
49	c
50	c=======================================================================
51	c
52	c 0. Declarations :
53	c ------------------
54	c
55	#include "dimensions.h"
56	#include "dimphys.h"
57	#include "comcstfi.h"
58	#include "surfdat.h"
59	#include "comgeomfi.h"
60	#include "comvert.h"
61	#include "paramet.h"
62	#include "callkeys.h"
63	#include "tracer.h"
64
65	#include "fisice.h"
66
67	c-----------------------------------------------------------------------
68	c Arguments :
69	c ---------
70	INTEGER ngrid, nlayer, nq
71
72	REAL ptimestep
73	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
74	REAL pcapcal(ngrid)
75	REAL pt(ngrid,nlayer)
76	REAL ptsrf(ngrid)
77	REAL pphi(ngrid,nlayer)
78	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
79	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
80	REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
81
82	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
83	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
84	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
85	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
86	REAL pdqc(ngrid,nlayer,nq)
87	REAL fluxsurf_sw(ngrid,2) ! added to calculate flux dependent albedo
88	c
89	c Local variables :
90	c -----------------
91
92	c variables used for albedo parametrization
93	c --------------------------------------------
94	INTEGER i,j
95	REAL Fluxmean(jjp1)
96	INTEGER l,ig,iq,icap,nmix
97	LOGICAL transparency, fluxdependent
98	c flag transparency if you want to make the co2ice semi-transparent
99	PARAMETER(transparency=.true.)
100	c flag fluxdependent if you want the co2ice albedo to be dependent on
101	c the incident solar flux
102	PARAMETER(fluxdependent=.false.)
103	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngridmx)
104	REAL zt(ngridmx,nlayermx)
105	REAL zcpi
106	REAL ztcond (ngridmx,nlayermx)
107	REAL ztcondsol(ngridmx)
108	REAL zdiceco2(ngridmx)
109	c REAL zcondicea(ngridmx,nlayermx) ! Already defined in fisice.h
110	REAL zcondices(ngridmx)
111	REAL zfallice(ngridmx,nlayermx+1) , zfallheat
112	REAL zmflux(nlayermx+1)
113	REAL zu(nlayermx),zv(nlayermx)
114	REAL zq(nlayermx,nqmx),zq1(nlayermx)
115	REAL ztsrf(ngridmx)
116	REAL ztc(nlayermx), ztm(nlayermx+1)
117	REAL zum(nlayermx+1) , zvm(nlayermx+1)
118	REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
119	REAL masse(nlayermx),w(nlayermx+1)
120	REAL Sm(nlayermx),Smq(nlayermx,nqmx),mixmas,qmix
121	LOGICAL condsub(ngridmx)
122
123	c variable speciale diagnostique
124	real tconda1(ngridmx,nlayermx)
125	real tconda2(ngridmx,nlayermx)
126	c REAL zdiceco2a(ngridmx) ! for diagnostic only
127	real zdtsig (ngridmx,nlayermx)
128	real zdt (ngridmx,nlayermx)
129
130
131	c local saved variables
132	integer ico2
133	real qco2min
134	save ico2,qco2min
135	REAL emisref(ngridmx)
136	REAL latcond,tcond1mb
137	REAL acond,bcond,ccond,cpice
138	REAL albediceF(ngridmx)
139	SAVE emisref
140	SAVE latcond,acond,bcond,ccond,cpice
141	SAVE albediceF
142
143
144	LOGICAL firstcall,firstcall2
145	SAVE firstcall,firstcall2
146	REAL SSUM
147	EXTERNAL SSUM
148
149	c common/scratch/zt,ztcond,zcondicea,zfallice,tconda1
150	c , ,tconda2,zdtsig,zdt,zu,zv
151
152
153	DATA latcond,tcond1mb/5.9e5,136.27/
154	DATA cpice /1000./
155	DATA firstcall/.true./
156	DATA firstcall2/.true./
157
158	integer flag
159
160	c----------------------------------------------------------------------
161
162	c Initialisation
163	c --------------
164	c
165	IF (firstcall) THEN
166	bcond=1./tcond1mb
167	ccond=cpp/(g*latcond)
168	acond=r/latcond
169
170	firstcall=.false.
171	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
172	PRINT*,'acond,bcond,ccond',acond,bcond,acond
173
174	ico2=0
175	if (tracer) then
176	c Prepare Special treatment if one of the tracer is CO2 gas
177	do iq=1,nqmx
178	if (noms(iq).eq."co2") ico2=iq
179	enddo
180	c minimum CO2 mix. ratio below which mixing occur with layer above:
181	qco2min =0.75
182	end if
183	ENDIF
184	zcpi=1./cpp
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,nlayermx) : 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,nqmx
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 forecast of atmospheric temperature zt and frost temperature ztcond
233	c --------------------------------------------------------------------
234
235	DO l=1,nlayer
236	DO ig=1,ngrid
237	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
238	ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
239	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
240	ENDDO
241	ENDDO
242
243	c Condensation/sublimation in the atmosphere
244	c ------------------------------------------
245	c (calcul of zcondicea , zfallice and pdtc)
246	c
247	DO l=nlayer , 1, -1
248	DO ig=1,ngrid
249	pdtc(ig,l)=0.
250	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
251	condsub(ig)=.true.
252	IF (zfallice(ig,l+1).gt.0) then
253	zfallheat=zfallice(ig,l+1)*
254	& (pphi(ig,l+1)-pphi(ig,l) +
255	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
256	ELSE
257	zfallheat=0.
258	ENDIF
259	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
260	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
261	& ccondpdtc(ig,l)- zfallheat
262	c Case when the ice from above sublimes entirely
263	c """""""""""""""""""""""""""""""""""""""""""""""
264	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
265	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
266	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
267	zcondicea(ig,l)= -zfallice(ig,l+1)
268	END IF
269
270	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
271	END IF
272	ENDDO
273	ENDDO
274
275	c *************************
276	c SURFACE CONDENSATION
277	c *************************
278
279	c forecast of ground temperature ztsrf and frost temperature ztcondsol
280	c --------------------------------------------------------------------
281	DO ig=1,ngrid
282	ztcondsol(ig)=1./(bcond-acondlog(.0095pplev(ig,1)))
283	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
284	ENDDO
285
286	c
287	c Condensation/sublimation on the ground
288	c --------------------------------------
289	c (calcul of zcondices , pdtsrfc)
290	c
291	DO ig=1,ngrid
292	IF(ig.GT.ngrid/2+1) THEN
293	icap=2
294	ELSE
295	icap=1
296	ENDIF
297
298	c
299	c Loop on where we have condensation/ sublimation
300	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
301	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
302	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
303	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
304	condsub(ig) = .true.
305
306	IF (zfallice(ig,1).gt.0) then
307	zfallheat=zfallice(ig,1)*
308	& (pphi(ig,1)- phisfi(ig) +
309	& cpice(ztcond(ig,1)-ztcondsol(ig)))/(latcondptimestep)
310	ELSE
311	zfallheat=0.
312	ENDIF
313
314	c condensation or partial sublimation of CO2 ice
315	c """""""""""""""""""""""""""""""""""""""""""""""
316	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
317	& /(latcond*ptimestep) - zfallheat
318	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
319
320	c If the entire CO_2 ice layer sublimes
321	c """"""""""""""""""""""""""""""""""""""""""""""""""""
322	c (including what has just condensed in the atmosphere)
323
324	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
325	& -zcondices(ig))THEN
326	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
327	pdtsrfc(ig)=(latcond/pcapcal(ig))*
328	& (zcondices(ig)+zfallheat)
329	END IF
330
331	c Changing CO2 ice amount and pressure :
332	c """"""""""""""""""""""""""""""""""""
333
334	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
335	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
336	pdpsrf(ig) = -zdiceco2(ig)*g
337
338	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
339	PRINT*,'STOP in condens'
340	PRINT*,'condensing more than total mass'
341	PRINT*,'Grid point ',ig
342	PRINT*,'Ps = ',pplev(ig,1)
343	PRINT*,'d Ps = ',pdpsrf(ig)
344	STOP
345	ENDIF
346	END IF
347	ENDDO
348
349	c ********************************************************************
350	c Surface albedo and emissivity of the surface below the snow (emisref)
351	c ********************************************************************
352	c Prepare the case where albedo varies with insolation:
353	c ----------------------------------------------------
354	if (fluxdependent) then
355
356	c Calcul du flux moyen (zonal mean)
357	do j=1,jjp1
358	Fluxmean(j)=0
359	do i=1,iim
360	ig=1+(j-2)*iim +i
361	if(j.eq.1) ig=1
362	if(j.eq.jjp1) ig=ngrid
363	Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
364	$ +fluxsurf_sw(ig,2)
365	enddo
366	Fluxmean(j)=Fluxmean(j)/float(iim)
367	enddo
368
369	c const A and B used to calculate the albedo which depends on solar flux
370	c albedice=constA+constB*Flux
371	c constT = time step to calculate the solar flux when flux decreases
372	constA=0.26
373	c constA=0.33
374	c constA=0.186
375	constB=0.00187
376	constT=10
377	endif
378
379	c Calcul de l'albedo
380	c ------------------
381	do ig =1,ngrid
382	IF(ig.GT.ngrid/2+1) THEN
383	icap=2
384	ELSE
385	icap=1
386	ENDIF
387	IF(firstcall2) THEN
388	albediceF(ig)=albedice(icap)
389	ENDIF
390	c test on the existence of co2ice ccccccccccccccccccccccc
391	if(.not.piceco2(ig).ge.0.) THEN
392	if(piceco2(ig).le.-5.e-8) print*,
393	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
394	piceco2(ig)=0.
395	endif
396	c if there is still co2ice ccccccccccccccccccccccc
397	if (piceco2(ig).gt.0) then
398	emisref(ig) = emisice(icap)
399
400	c if flux dependent albedo is used
401	c --------------------------------
402	if (fluxdependent) then
403	j=INT((ig-2)/iim)+2
404	if(ig.eq.1) j=1
405	if(ig.eq.ngrid) j=jjp1
406	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
407	c $ constA+constB*250)
408	albediceF_new(ig)=constA+constB*Fluxmean(j)
409	if (albediceF(ig).gt.albediceF_new(ig)) then
410	albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
411	$ constT)*(albediceF_new(ig)-albediceF(ig))
412	else
413	albediceF(ig)=albediceF_new(ig)
414	endif
415	c if part of the ice is transparent
416	c slopy = pente de la droite: alpha = y*co2ice/1620
417	c pour une valeur superieur a une epaisseur de glace donnee
418	c ici, epaisseur limite = 10cm
419	if (transparency) then
420	slopy=1/(1620*0.10)
421	alpha=MIN(slopy*piceco2(ig),1.)
422	psolaralb(ig,1) = alpha*albediceF(ig)
423	$ +(1-alpha)*albedodat(ig)
424	psolaralb(ig,2) = psolaralb(ig,1)
425	else
426	psolaralb(ig,1) = albediceF(ig)
427	psolaralb(ig,2) = psolaralb(ig,1)
428	endif
429	else
430	c transparency set to true and fluxdependent set to false
431	if (transparency) then
432	slopy=1/(1620*0.10)
433	alpha=MIN(slopy*piceco2(ig),1.)
434	psolaralb(ig,1) = alpha*albedice(icap)
435	$ +(1-alpha)*albedodat(ig)
436	psolaralb(ig,2) = psolaralb(ig,1)
437	else
438	c simplest case: transparency and flux dependent set to false
439	psolaralb(ig,1) = albedice(icap)
440	psolaralb(ig,2) = albedice(icap)
441	endif
442	endif
443	c no more co2ice, albedo = ground albedo
444	else
445	psolaralb(ig,1) = albedodat(ig)
446	psolaralb(ig,2) = albedodat(ig)
447	emisref(ig) = emissiv
448	pemisurf(ig) = emissiv
449	endif
450	end do ! end of the ig loop
451
452	firstcall2=.false.
453	c ***************************************************************
454	c Correction to account for redistribution between sigma or hybrid
455	c layers when changing surface pressure (and warming/cooling
456	c of the CO2 currently changing phase).
457	c *************************************************************
458
459	DO ig=1,ngrid
460	if (condsub(ig)) then
461	do l=1,nlayer
462	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
463	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
464	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
465	do iq=1,nqmx
466	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
467	enddo
468	end do
469
470	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
471	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
472
473	zmflux(1) = -zcondices(ig)
474	DO l=1,nlayer
475	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
476	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
477	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
478	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
479	END DO
480
481	c Mass of each layer
482	c ------------------
483	DO l=1,nlayer
484	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
485	END DO
486
487
488	c Corresponding fluxes for T,U,V,Q
489	c """"""""""""""""""""""""""""""""
490
491	c averaging operator for TRANSPORT
492	c """"""""""""""""""""""""""""""""
493	c Value transfert at the surface interface when condensation
494	c sublimation:
495	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
496	zum(1) = 0
497	zvm(1) = 0
498	do iq=1,nqmx
499	zqm(1,iq)=0. ! most tracer do not condense !
500	enddo
501	c Special case if one of the tracer is CO2 gas
502	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
503
504	c Van Leer scheme:
505	DO l=1,nlayer+1
506	w(l)=-zmflux(l)*ptimestep
507	END DO
508	call vl1d(ztc,2.,masse,w,ztm)
509	call vl1d(zu ,2.,masse,w,zum)
510	call vl1d(zv ,2.,masse,w,zvm)
511	do iq=1,nqmx
512	do l=1,nlayer
513	zq1(l)=zq(l,iq)
514	enddo
515	zqm1(1)=zqm(1,iq)
516	call vl1d(zq1,2.,masse,w,zqm1)
517	do l=2,nlayer
518	zq( l,iq)=zq1(l)
519	zqm(l,iq)=zqm1(l)
520	enddo
521	enddo
522
523	c Surface condensation affects low winds
524	if (zmflux(1).lt.0) then
525	zum(1)= zu(1) * (w(1)/masse(1))
526	zvm(1)= zv(1) * (w(1)/masse(1))
527	if (w(1).gt.masse(1)) then ! ensure numerical stability
528	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
529	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
530	end if
531	end if
532
533	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
534	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
535	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
536	do iq=1,nqmx
537	zqm(nlayer+1,iq)= zq(nlayer,iq)
538	enddo
539
540	c Tendencies on T, U, V, Q
541	c """"""""""""""""""""""""
542	DO l=1,nlayer
543
544	c Tendencies on T
545	zdtsig(ig,l) = (1/masse(l)) *
546	& ( zmflux(l)*(ztm(l) - ztc(l))
547	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
548	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
549	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
550
551	c Tendencies on U
552	pduc(ig,l) = (1/masse(l)) *
553	& ( zmflux(l)*(zum(l) - zu(l))
554	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
555
556
557	c Tendencies on V
558	pdvc(ig,l) = (1/masse(l)) *
559	& ( zmflux(l)*(zvm(l) - zv(l))
560	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
561
562	END DO
563
564	c Tendencies on Q
565	do iq=1,nqmx
566	if (noms(iq).eq.'co2') then
567	c SPECIAL Case when the tracer IS CO2 :
568	DO l=1,nlayer
569	pdqc(ig,l,iq)= (1/masse(l)) *
570	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
571	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
572	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
573	END DO
574	else
575	DO l=1,nlayer
576	pdqc(ig,l,iq)= (1/masse(l)) *
577	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
578	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
579	& + zcondicea(ig,l)*zq(l,iq) )
580	END DO
581	end if
582	enddo
583
584	c --------------------------------------------------------
585	c Roughly Simulate Molecular mixing when CO2 is too depleted by
586	c Surface condensation (mixing starts if qco2 < qco2min )
587	c FF 06/2004
588	c WARNING : this is now done in convadj, better (FF 02/2005)
589	c --------------------------------------------------------
590	flag=0 ! now done in convadj : must be =0
591	if (flag.eq.1) then
592	if(ico2.gt.0) then ! relevant only if one tracer is CO2
593	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
594	& .lt.qco2min) then
595	do iq=1,nqmx
596	zq(1,iq)=pq(ig,1,iq)
597	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
598	Smq(1,iq) = masse(1)*zq(1,iq)
599	end do
600	Sm(1) = masse(1)
601	do l =2,nlayermx
602	do iq=1,nqmx
603	zq(l,iq)=pq(ig,l,iq)
604	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
605	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
606	end do
607	sm(l) = sm(l-1) + masse(l)
608	if(zq(l,ico2).gt.qco2min) then
609	c mixmas: mass of atmosphere that must be mixed to reach qco2min
610	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
611	& /(zq(l,ico2)-qco2min)
612	if((mixmas.le.sm(l)))then
613	c OK if mixed mass less than mass of layers affected
614	nmix=l ! number of layer affected by mixing
615	goto 99
616	end if
617	end if
618	end do
619	99 continue
620	do iq=1,nqmx
621	qmix=zq(nmix,iq)
622	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
623	do l=1,nmix-1
624	pdqc(ig,l,iq)=
625	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
626	end do
627	c layer only partly mixed :
628	pdqc(ig,nmix,iq)=(
629	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
630	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
631
632	end do
633	end if
634	end if
635
636	endif ! (flag.eq.1)
637	end if ! if (condsub)
638	END DO ! loop on ig
639
640	c ***************************************************************
641	c CO2 snow / clouds scheme
642	c ***************************************************************
643
644	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
645	& zcondicea,zcondices,zfallice,pemisurf)
646
647	c ***************************************************************
648	c Ecriture des diagnostiques
649	c ***************************************************************
650
651	c DO l=1,nlayer
652	c DO ig=1,ngrid
653	c Taux de cond en kg.m-2.pa-1.s-1
654	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
655	c Taux de cond en kg.m-3.s-1
656	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
657	c END DO
658	c END DO
659	c call WRITEDIAGFI(ngridmx,'tconda1',
660	c &'Taux de condensation CO2 atmospherique /Pa',
661	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
662	c call WRITEDIAGFI(ngridmx,'tconda2',
663	c &'Taux de condensation CO2 atmospherique /m',
664	c & 'kg.m-3.s-1',3,tconda2)
665
666	return
667	end
668
669
670
671	c *****************************************************************
672	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
673	c
674	c
675	c Operateur de moyenne inter-couche pour calcul de transport type
676	c Van-Leer " pseudo amont " dans la verticale
677	c q,w sont des arguments d'entree pour le s-pg ....
678	c masse : masse de la couche Dp/g
679	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
680	c pente_max = 2 conseillee
681	c
682	c
683	c --------------------------------------------------------------------
684	IMPLICIT NONE
685
686	#include "dimensions.h"
687
688	c
689	c
690	c
691	c Arguments:
692	c ----------
693	real masse(llm),pente_max
694	REAL q(llm),qm(llm+1)
695	REAL w(llm+1)
696	c
697	c Local
698	c ---------
699	c
700	INTEGER l
701	c
702	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
703	real sigw, Mtot, MQtot
704	integer m
705	c integer ismax,ismin
706
707
708	c On oriente tout dans le sens de la pression
709	c W > 0 WHEN DOWN !!!!!!!!!!!!!
710
711	do l=2,llm
712	dzqw(l)=q(l-1)-q(l)
713	adzqw(l)=abs(dzqw(l))
714	enddo
715
716	do l=2,llm-1
717	if(dzqw(l)*dzqw(l+1).gt.0.) then
718	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
719	else
720	dzq(l)=0.
721	endif
722	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
723	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
724	enddo
725
726	dzq(1)=0.
727	dzq(llm)=0.
728
729	do l = 1,llm-1
730
731	c Regular scheme (transfered mass < layer mass)
732	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
733	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
734	sigw=w(l+1)/masse(l+1)
735	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
736	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
737	sigw=w(l+1)/masse(l)
738	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
739
740	c Extended scheme (transfered mass > layer mass)
741	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
742	else if(w(l+1).gt.0.) then
743	m=l+1
744	Mtot = masse(m)
745	MQtot = masse(m)*q(m)
746	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
747	m=m+1
748	Mtot = Mtot + masse(m)
749	MQtot = MQtot + masse(m)*q(m)
750	end do
751	if (m.lt.llm) then
752	sigw=(w(l+1)-Mtot)/masse(m+1)
753	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
754	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
755	else
756	w(l+1) = Mtot
757	qm(l+1) = Mqtot / Mtot
758	write(,) 'top layer is disapearing !'
759	stop
760	end if
761	else ! if(w(l+1).lt.0)
762	m = l-1
763	Mtot = masse(m+1)
764	MQtot = masse(m+1)*q(m+1)
765	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
766	m=m-1
767	Mtot = Mtot + masse(m+1)
768	MQtot = MQtot + masse(m+1)*q(m+1)
769	end do
770	if (m.gt.0) then
771	sigw=(w(l+1)+Mtot)/masse(m)
772	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
773	& (q(m)-0.5(1.+sigw)dzq(m)) )
774	else
775	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
776	end if
777	end if
778	enddo
779
780	c boundary conditions (not used in newcondens !!)
781	c qm(llm+1)=0.
782	c if(w(1).gt.0.) then
783	c qm(1)=q(1)
784	c else
785	c qm(1)=0.
786	c end if
787
788	return
789	end

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

Download in other formats:

Original Format