Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

newcondens.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 26.6 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
163	c Initialisation
164	c --------------
165	c
166	IF (firstcall) THEN
167	bcond=1./tcond1mb
168	ccond=cpp/(g*latcond)
169	acond=r/latcond
170
171	firstcall=.false.
172	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
173	PRINT*,'acond,bcond,ccond',acond,bcond,acond
174
175	ico2=0
176	if (tracer) then
177	c Prepare Special treatment if one of the tracer is CO2 gas
178	do iq=1,nqmx
179	if (noms(iq).eq."co2") ico2=iq
180	enddo
181	c minimum CO2 mix. ratio below which mixing occur with layer above:
182	qco2min =0.75
183	end if
184	ENDIF
185	zcpi=1./cpp
186
187	c
188	c======================================================================
189	c Calcul of CO2 condensation sublimation
190	c ============================================================
191	c
192	c Used variable :
193	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
194	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
195	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
196	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
197	c
198	c pdtc(ngrid,nlayermx) : dT/dt due to cond/sub
199	c
200	c
201	c Tendencies set to 0 (except pdtc)
202	c -------------------------------------
203	DO l=1,nlayer
204	DO ig=1,ngrid
205	zcondicea(ig,l) = 0.
206	zfallice(ig,l) = 0.
207	pduc(ig,l) = 0
208	pdvc(ig,l) = 0
209	END DO
210	END DO
211
212	DO iq=1,nqmx
213	DO l=1,nlayer
214	DO ig=1,ngrid
215	pdqc(ig,l,iq) = 0
216	END DO
217	END DO
218	END DO
219
220	DO ig=1,ngrid
221	zfallice(ig,nlayer+1) = 0.
222	zcondices(ig) = 0.
223	pdtsrfc(ig) = 0.
224	pdpsrf(ig) = 0.
225	condsub(ig) = .false.
226	zdiceco2(ig) = 0.
227	ENDDO
228	zfallheat=0
229
230	c *************************
231	c ATMOSPHERIC CONDENSATION
232	c *************************
233
234	c forecast of atmospheric temperature zt and frost temperature ztcond
235	c --------------------------------------------------------------------
236
237	DO l=1,nlayer
238	DO ig=1,ngrid
239	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
240	ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
241	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
242	ENDDO
243	ENDDO
244
245	c Condensation/sublimation in the atmosphere
246	c ------------------------------------------
247	c (calcul of zcondicea , zfallice and pdtc)
248	c
249	DO l=nlayer , 1, -1
250	DO ig=1,ngrid
251	pdtc(ig,l)=0.
252	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
253	condsub(ig)=.true.
254	IF (zfallice(ig,l+1).gt.0) then
255	zfallheat=zfallice(ig,l+1)*
256	& (pphi(ig,l+1)-pphi(ig,l) +
257	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
258	ELSE
259	zfallheat=0.
260	ENDIF
261	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
262	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
263	& ccondpdtc(ig,l)- zfallheat
264	c Case when the ice from above sublimes entirely
265	c """""""""""""""""""""""""""""""""""""""""""""""
266	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
267	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
268	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
269	zcondicea(ig,l)= -zfallice(ig,l+1)
270	END IF
271
272	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
273	END IF
274	ENDDO
275	ENDDO
276
277	c *************************
278	c SURFACE CONDENSATION
279	c *************************
280
281	c forecast of ground temperature ztsrf and frost temperature ztcondsol
282	c --------------------------------------------------------------------
283	DO ig=1,ngrid
284	ztcondsol(ig)=1./(bcond-acondlog(.0095pplev(ig,1)))
285	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
286	ENDDO
287
288	c
289	c Condensation/sublimation on the ground
290	c --------------------------------------
291	c (calcul of zcondices , pdtsrfc)
292	c
293	DO ig=1,ngrid
294	IF(ig.GT.ngrid/2+1) THEN
295	icap=2
296	ELSE
297	icap=1
298	ENDIF
299
300	c
301	c Loop on where we have condensation/ sublimation
302	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
303	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
304	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
305	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
306	condsub(ig) = .true.
307
308	IF (zfallice(ig,1).gt.0) then
309	zfallheat=zfallice(ig,1)*
310	& (pphi(ig,1)- phisfi(ig) +
311	& cpice(ztcond(ig,1)-ztcondsol(ig)))/(latcondptimestep)
312	ELSE
313	zfallheat=0.
314	ENDIF
315
316	c condensation or partial sublimation of CO2 ice
317	c """""""""""""""""""""""""""""""""""""""""""""""
318	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
319	& /(latcond*ptimestep) - zfallheat
320	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
321
322	c If the entire CO_2 ice layer sublimes
323	c """"""""""""""""""""""""""""""""""""""""""""""""""""
324	c (including what has just condensed in the atmosphere)
325
326	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
327	& -zcondices(ig))THEN
328	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
329	pdtsrfc(ig)=(latcond/pcapcal(ig))*
330	& (zcondices(ig)+zfallheat)
331	END IF
332
333	c Changing CO2 ice amount and pressure :
334	c """"""""""""""""""""""""""""""""""""
335
336	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
337	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
338	pdpsrf(ig) = -zdiceco2(ig)*g
339
340	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
341	PRINT*,'STOP in condens'
342	PRINT*,'condensing more than total mass'
343	PRINT*,'Grid point ',ig
344	PRINT*,'Ps = ',pplev(ig,1)
345	PRINT*,'d Ps = ',pdpsrf(ig)
346	STOP
347	ENDIF
348	END IF
349	ENDDO
350
351	c ********************************************************************
352	c Surface albedo and emissivity of the surface below the snow (emisref)
353	c ********************************************************************
354	c Prepare the case where albedo varies with insolation:
355	c ----------------------------------------------------
356	if (fluxdependent) then
357
358	c Calcul du flux moyen (zonal mean)
359	do j=1,jjp1
360	Fluxmean(j)=0
361	do i=1,iim
362	ig=1+(j-2)*iim +i
363	if(j.eq.1) ig=1
364	if(j.eq.jjp1) ig=ngrid
365	Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
366	$ +fluxsurf_sw(ig,2)
367	enddo
368	Fluxmean(j)=Fluxmean(j)/float(iim)
369	enddo
370
371	c const A and B used to calculate the albedo which depends on solar flux
372	c albedice=constA+constB*Flux
373	c constT = time step to calculate the solar flux when flux decreases
374	constA=0.26
375	c constA=0.33
376	c constA=0.186
377	constB=0.00187
378	constT=10
379	endif
380
381	c Calcul de l'albedo
382	c ------------------
383	do ig =1,ngrid
384	IF(ig.GT.ngrid/2+1) THEN
385	icap=2
386	ELSE
387	icap=1
388	ENDIF
389	IF(firstcall2) THEN
390	albediceF(ig)=albedice(icap)
391	ENDIF
392	c test on the existence of co2ice ccccccccccccccccccccccc
393	if(.not.piceco2(ig).ge.0.) THEN
394	if(piceco2(ig).le.-5.e-8) print*,
395	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
396	piceco2(ig)=0.
397	endif
398	c if there is still co2ice ccccccccccccccccccccccc
399	if (piceco2(ig).gt.0) then
400	emisref(ig) = emisice(icap)
401
402	c if flux dependent albedo is used
403	c --------------------------------
404	if (fluxdependent) then
405	j=INT((ig-2)/iim)+2
406	if(ig.eq.1) j=1
407	if(ig.eq.ngrid) j=jjp1
408	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
409	c $ constA+constB*250)
410	albediceF_new(ig)=constA+constB*Fluxmean(j)
411	if (albediceF(ig).gt.albediceF_new(ig)) then
412	albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
413	$ constT)*(albediceF_new(ig)-albediceF(ig))
414	else
415	albediceF(ig)=albediceF_new(ig)
416	endif
417	c if part of the ice is transparent
418	c slopy = pente de la droite: alpha = y*co2ice/1620
419	c pour une valeur superieur a une epaisseur de glace donnee
420	c ici, epaisseur limite = 10cm
421	if (transparency) then
422	slopy=1/(1620*0.10)
423	alpha=MIN(slopy*piceco2(ig),1.)
424	psolaralb(ig,1) = alpha*albediceF(ig)
425	$ +(1-alpha)*albedodat(ig)
426	psolaralb(ig,2) = psolaralb(ig,1)
427	else
428	psolaralb(ig,1) = albediceF(ig)
429	psolaralb(ig,2) = psolaralb(ig,1)
430	endif
431	else
432	c transparency set to true and fluxdependent set to false
433	if (transparency) then
434	slopy=1/(1620*0.10)
435	alpha=MIN(slopy*piceco2(ig),1.)
436	psolaralb(ig,1) = alpha*albedice(icap)
437	$ +(1-alpha)*albedodat(ig)
438	psolaralb(ig,2) = psolaralb(ig,1)
439	else
440	c simplest case: transparency and flux dependent set to false
441	psolaralb(ig,1) = albedice(icap)
442	psolaralb(ig,2) = albedice(icap)
443	endif
444	endif
445	c no more co2ice, albedo = ground albedo
446	else
447	psolaralb(ig,1) = albedodat(ig)
448	psolaralb(ig,2) = albedodat(ig)
449	emisref(ig) = emissiv
450	pemisurf(ig) = emissiv
451	endif
452	end do ! end of the ig loop
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,nqmx
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	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
479	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
480	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
481	END DO
482
483	c Mass of each layer
484	c ------------------
485	DO l=1,nlayer
486	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
487	END DO
488
489
490	c Corresponding fluxes for T,U,V,Q
491	c """"""""""""""""""""""""""""""""
492
493	c averaging operator for TRANSPORT
494	c """"""""""""""""""""""""""""""""
495	c Value transfert at the surface interface when condensation
496	c sublimation:
497	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
498	zum(1) = 0
499	zvm(1) = 0
500	do iq=1,nqmx
501	zqm(1,iq)=0. ! most tracer do not condense !
502	enddo
503	c Special case if one of the tracer is CO2 gas
504	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
505
506	c Van Leer scheme:
507	DO l=1,nlayer+1
508	w(l)=-zmflux(l)*ptimestep
509	END DO
510	call vl1d(ztc,2.,masse,w,ztm)
511	call vl1d(zu ,2.,masse,w,zum)
512	call vl1d(zv ,2.,masse,w,zvm)
513	do iq=1,nqmx
514	do l=1,nlayer
515	zq1(l)=zq(l,iq)
516	enddo
517	zqm1(1)=zqm(1,iq)
518	call vl1d(zq1,2.,masse,w,zqm1)
519	do l=2,nlayer
520	zq( l,iq)=zq1(l)
521	zqm(l,iq)=zqm1(l)
522	enddo
523	enddo
524
525	c Surface condensation affects low winds
526	if (zmflux(1).lt.0) then
527	zum(1)= zu(1) * (w(1)/masse(1))
528	zvm(1)= zv(1) * (w(1)/masse(1))
529	if (w(1).gt.masse(1)) then ! ensure numerical stability
530	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
531	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
532	end if
533	end if
534
535	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
536	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
537	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
538	do iq=1,nqmx
539	zqm(nlayer+1,iq)= zq(nlayer,iq)
540	enddo
541
542	cc Tendencies on T, U, V, Q
543	cc """"""""""""""""""""""""
544	c DO l=1,nlayer
545	c
546	cc Tendencies on T
547	c zdtsig(ig,l) = (1/masse(l)) *
548	c & ( zmflux(l)*(ztm(l) - ztc(l))
549	c & - zmflux(l+1)*(ztm(l+1) - ztc(l))
550	c & + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
551	c pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
552	c
553	cc Tendencies on U
554	c pduc(ig,l) = (1/masse(l)) *
555	c & ( zmflux(l)*(zum(l) - zu(l))
556	c & - zmflux(l+1)*(zum(l+1) - zu(l)) )
557	c
558	c
559	cc Tendencies on V
560	c pdvc(ig,l) = (1/masse(l)) *
561	c & ( zmflux(l)*(zvm(l) - zv(l))
562	c & - zmflux(l+1)*(zvm(l+1) - zv(l)) )
563	c
564	c END DO
565
566	c Tendencies on Q
567	do iq=1,nqmx
568	if (noms(iq).eq.'co2') then
569	c SPECIAL Case when the tracer IS CO2 :
570	DO l=1,nlayer
571	pdqc(ig,l,iq)= (1/masse(l)) *
572	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
573	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
574	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
575	END DO
576	else
577	DO l=1,nlayer
578	pdqc(ig,l,iq)= (1/masse(l)) *
579	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
580	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
581	& + zcondicea(ig,l)*zq(l,iq) )
582	END DO
583	end if
584	enddo
585
586	c --------------------------------------------------------
587	c Roughly Simulate Molecular mixing when CO2 is too depleted by
588	c Surface condensation (mixing starts if qco2 < qco2min )
589	c FF 06/2004
590	c WARNING : this is now done in convadj, better (FF 02/2005)
591	c --------------------------------------------------------
592	flag=0 ! now done in convadj : must be =0
593	if (flag.eq.1) then
594	if(ico2.gt.0) then ! relevant only if one tracer is CO2
595	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
596	& .lt.qco2min) then
597	do iq=1,nqmx
598	zq(1,iq)=pq(ig,1,iq)
599	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
600	Smq(1,iq) = masse(1)*zq(1,iq)
601	end do
602	Sm(1) = masse(1)
603	do l =2,nlayermx
604	do iq=1,nqmx
605	zq(l,iq)=pq(ig,l,iq)
606	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
607	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
608	end do
609	sm(l) = sm(l-1) + masse(l)
610	if(zq(l,ico2).gt.qco2min) then
611	c mixmas: mass of atmosphere that must be mixed to reach qco2min
612	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
613	& /(zq(l,ico2)-qco2min)
614	if((mixmas.le.sm(l)))then
615	c OK if mixed mass less than mass of layers affected
616	nmix=l ! number of layer affected by mixing
617	goto 99
618	end if
619	end if
620	end do
621	99 continue
622	do iq=1,nqmx
623	qmix=zq(nmix,iq)
624	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
625	do l=1,nmix-1
626	pdqc(ig,l,iq)=
627	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
628	end do
629	c layer only partly mixed :
630	pdqc(ig,nmix,iq)=(
631	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
632	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
633
634	end do
635	end if
636	end if
637
638	endif ! (flag.eq.1)
639	end if ! if (condsub)
640	END DO ! loop on ig
641
642	c ***************************************************************
643	c CO2 snow / clouds scheme
644	c ***************************************************************
645
646	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
647	& zcondicea,zcondices,zfallice,pemisurf)
648
649	c ***************************************************************
650	c Ecriture des diagnostiques
651	c ***************************************************************
652
653	c DO l=1,nlayer
654	c DO ig=1,ngrid
655	c Taux de cond en kg.m-2.pa-1.s-1
656	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
657	c Taux de cond en kg.m-3.s-1
658	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
659	c END DO
660	c END DO
661	c call WRITEDIAGFI(ngridmx,'tconda1',
662	c &'Taux de condensation CO2 atmospherique /Pa',
663	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
664	c call WRITEDIAGFI(ngridmx,'tconda2',
665	c &'Taux de condensation CO2 atmospherique /m',
666	c & 'kg.m-3.s-1',3,tconda2)
667
668
669	!print*,'fin newcondens piceco2',piceco2(ig)
670	!print*,'ptsrf ztsrf', ptsrf(ig),ztsrf(ig)
671	!print*,'ztcondsol(ig)', ztcondsol(ig)
672	!print*,'ztcondsol(ig)', ztcondsol(ig)
673	!print,'tend. press', pdpsrf(ig)ptimestep
674
675
676	return
677	end
678
679
680
681	c *****************************************************************
682	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
683	c
684	c
685	c Operateur de moyenne inter-couche pour calcul de transport type
686	c Van-Leer " pseudo amont " dans la verticale
687	c q,w sont des arguments d'entree pour le s-pg ....
688	c masse : masse de la couche Dp/g
689	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
690	c pente_max = 2 conseillee
691	c
692	c
693	c --------------------------------------------------------------------
694	IMPLICIT NONE
695
696	#include "dimensions.h"
697
698	c
699	c
700	c
701	c Arguments:
702	c ----------
703	real masse(llm),pente_max
704	REAL q(llm),qm(llm+1)
705	REAL w(llm+1)
706	c
707	c Local
708	c ---------
709	c
710	INTEGER l
711	c
712	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
713	real sigw, Mtot, MQtot
714	integer m
715	c integer ismax,ismin
716
717
718	c On oriente tout dans le sens de la pression
719	c W > 0 WHEN DOWN !!!!!!!!!!!!!
720
721	do l=2,llm
722	dzqw(l)=q(l-1)-q(l)
723	adzqw(l)=abs(dzqw(l))
724	enddo
725
726	do l=2,llm-1
727	if(dzqw(l)*dzqw(l+1).gt.0.) then
728	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
729	else
730	dzq(l)=0.
731	endif
732	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
733	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
734	enddo
735
736	dzq(1)=0.
737	dzq(llm)=0.
738
739	do l = 1,llm-1
740
741	c Regular scheme (transfered mass < layer mass)
742	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
743	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
744	sigw=w(l+1)/masse(l+1)
745	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
746	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
747	sigw=w(l+1)/masse(l)
748	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
749
750	c Extended scheme (transfered mass > layer mass)
751	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
752	else if(w(l+1).gt.0.) then
753	m=l+1
754	Mtot = masse(m)
755	MQtot = masse(m)*q(m)
756	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
757	m=m+1
758	Mtot = Mtot + masse(m)
759	MQtot = MQtot + masse(m)*q(m)
760	end do
761	if (m.lt.llm) then
762	sigw=(w(l+1)-Mtot)/masse(m+1)
763	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
764	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
765	else
766	w(l+1) = Mtot
767	qm(l+1) = Mqtot / Mtot
768	write(,) 'top layer is disapearing !'
769	stop
770	end if
771	else ! if(w(l+1).lt.0)
772	m = l-1
773	Mtot = masse(m+1)
774	MQtot = masse(m+1)*q(m+1)
775	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
776	m=m-1
777	Mtot = Mtot + masse(m+1)
778	MQtot = MQtot + masse(m+1)*q(m+1)
779	end do
780	if (m.gt.0) then
781	sigw=(w(l+1)+Mtot)/masse(m)
782	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
783	& (q(m)-0.5(1.+sigw)dzq(m)) )
784	else
785	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
786	end if
787	end if
788	enddo
789
790	c boundary conditions (not used in newcondens !!)
791	c qm(llm+1)=0.
792	c if(w(1).gt.0.) then
793	c qm(1)=q(1)
794	c else
795	c qm(1)=0.
796	c end if
797
798	return
799	end

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

Download in other formats:

Original Format