Context Navigation

newcondens.F @ 317

Last change on this file since 317 was 120, checked in by emillour, 14 years ago

-Mars GCM:

set internal computations using double precision in growthrate.F and

watercloud.F (otherwise we sometimes end up with Nans).

add extra checks in newcondens.F to avoid possibility of out of bounds

evaluation of array masse()

File size: 30.2 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,zls)
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	c-----------------------------------------------------------------------
66	c Arguments :
67	c ---------
68	INTEGER ngrid, nlayer, nq
69
70	REAL ptimestep
71	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
72	REAL pcapcal(ngrid)
73	REAL pt(ngrid,nlayer)
74	REAL ptsrf(ngrid)
75	REAL pphi(ngrid,nlayer)
76	REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
77	REAL pdtsrfc(ngrid),pdpsrf(ngrid)
78	REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
79
80	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
81	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
82	REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
83	REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
84	REAL pdqc(ngrid,nlayer,nq)
85	REAL fluxsurf_sw(ngrid,2) ! added to calculate flux dependent albedo
86	real zls ! solar longitude (rad)
87	c
88	c Local variables :
89	c -----------------
90
91	c variables used for albedo parametrization
92	c --------------------------------------------
93	INTEGER i,j
94	REAL Fluxmean(jjp1)
95	INTEGER l,ig,iq,icap,nmix
96	LOGICAL transparency, fluxdependent
97	c flag transparency if you want to make the co2ice semi-transparent
98	PARAMETER(transparency=.true.)
99	c flag fluxdependent if you want the co2ice albedo to be dependent on
100	c the incident solar flux
101	PARAMETER(fluxdependent=.false.)
102	REAL slopy,alpha,constA,constB,constT,albediceF_new(ngridmx)
103	REAL zt(ngridmx,nlayermx)
104	REAL zcpi
105	REAL ztcond (ngridmx,nlayermx)
106	REAL ztcondsol(ngridmx)
107	REAL zdiceco2(ngridmx)
108	REAL zcondicea(ngridmx,nlayermx)
109	REAL zcondices(ngridmx)
110	REAL zfallice(ngridmx,nlayermx+1) , zfallheat
111	REAL zmflux(nlayermx+1)
112	REAL zu(nlayermx),zv(nlayermx)
113	REAL zq(nlayermx,nqmx),zq1(nlayermx)
114	REAL ztsrf(ngridmx)
115	REAL ztc(nlayermx), ztm(nlayermx+1)
116	REAL zum(nlayermx+1) , zvm(nlayermx+1)
117	REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
118	REAL masse(nlayermx),w(nlayermx+1)
119	REAL Sm(nlayermx),Smq(nlayermx,nqmx),mixmas,qmix
120	LOGICAL condsub(ngridmx)
121
122	c variable speciale diagnostique
123	real tconda1(ngridmx,nlayermx)
124	real tconda2(ngridmx,nlayermx)
125	c REAL zdiceco2a(ngridmx) ! for diagnostic only
126	real zdtsig (ngridmx,nlayermx)
127	real zdt (ngridmx,nlayermx)
128	real vmr_co2(ngridmx,nlayermx) ! co2 volume mixing ratio
129	! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
130	! then condensation temperature is computed using partial pressure of CO2
131	logical,parameter :: improved_ztcond=.true.
132	! Bound co2 (tracer) values...
133	logical,parameter :: bound_qco2=.false.
134	real,parameter :: qco2max=1.1
135	real,parameter :: qco2mini=0.1
136	real :: zqco2
137
138	c local saved variables
139	integer ico2
140	real qco2min,qco2,mmean
141	save ico2,qco2min
142	REAL emisref(ngridmx)
143	REAL latcond,tcond1mb
144	REAL acond,bcond,ccond,cpice
145	! REAL albediceF(ngridmx)
146	SAVE emisref
147	SAVE latcond,acond,bcond,ccond,cpice
148	! SAVE albediceF
149	real m_co2, m_noco2, A , B
150	save A, B, m_co2, m_noco2
151
152	LOGICAL firstcall !,firstcall2
153	SAVE firstcall !,firstcall2
154	! REAL SSUM
155	! EXTERNAL SSUM
156
157	DATA latcond,tcond1mb/5.9e5,136.27/
158	DATA cpice /1000./
159	DATA firstcall/.true./
160	! DATA firstcall2/.true./
161
162	integer flag
163
164	c----------------------------------------------------------------------
165
166	c Initialisation
167	c --------------
168	c
169	IF (firstcall) THEN
170	bcond=1./tcond1mb
171	ccond=cpp/(g*latcond)
172	acond=r/latcond
173
174	firstcall=.false.
175	write(,) 'Newcondens: improved_ztcond=',improved_ztcond
176	write(,) 'Newcondens: bound_qco2=',bound_qco2
177	PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
178	PRINT*,'acond,bcond,ccond',acond,bcond,ccond
179
180	ico2=0
181
182	if (tracer) then
183	c Prepare Special treatment if one of the tracer is CO2 gas
184	do iq=1,nqmx
185	if (noms(iq).eq."co2") then
186	ico2=iq
187	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
188	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
189	c Compute A and B coefficient use to compute
190	c mean molecular mass Mair defined by
191	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
192	c 1/Mair = A*q(ico2) + B
193	A =(1/m_co2 - 1/m_noco2)
194	B=1/m_noco2
195	endif
196	enddo
197	c minimum CO2 mix. ratio below which mixing occur with layer above:
198	qco2min =0.75
199	end if
200	ENDIF
201	zcpi=1./cpp
202	c
203	c======================================================================
204	c Calcul of CO2 condensation sublimation
205	c ============================================================
206	c
207	c Used variable :
208	c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2)
209	c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s)
210	c zcondices(ngrid): condensation rate on the ground (kg/m2/s)
211	c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
212	c
213	c pdtc(ngrid,nlayermx) : dT/dt due to cond/sub
214	c
215	c
216	c Tendencies set to 0 (except pdtc)
217	c -------------------------------------
218	DO l=1,nlayer
219	DO ig=1,ngrid
220	zcondicea(ig,l) = 0.
221	zfallice(ig,l) = 0.
222	pduc(ig,l) = 0
223	pdvc(ig,l) = 0
224	END DO
225	END DO
226
227	DO iq=1,nqmx
228	DO l=1,nlayer
229	DO ig=1,ngrid
230	pdqc(ig,l,iq) = 0
231	END DO
232	END DO
233	END DO
234
235	DO ig=1,ngrid
236	zfallice(ig,nlayer+1) = 0.
237	zcondices(ig) = 0.
238	pdtsrfc(ig) = 0.
239	pdpsrf(ig) = 0.
240	condsub(ig) = .false.
241	zdiceco2(ig) = 0.
242	ENDDO
243	zfallheat=0
244
245	c *************************
246	c ATMOSPHERIC CONDENSATION
247	c *************************
248
249	c Compute CO2 Volume mixing ratio
250	c -------------------------------
251	if (improved_ztcond.and.(ico2.ne.0)) then
252	DO l=1,nlayer
253	DO ig=1,ngrid
254	qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
255	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
256	mmean=1/(A*qco2 +B)
257	vmr_co2(ig,l) = qco2*mmean/m_co2
258	ENDDO
259	ENDDO
260	else
261	DO l=1,nlayer
262	DO ig=1,ngrid
263	vmr_co2(ig,l)=0.95
264	ENDDO
265	ENDDO
266	end if
267
268	c forecast of atmospheric temperature zt and frost temperature ztcond
269	c --------------------------------------------------------------------
270
271	DO l=1,nlayer
272	DO ig=1,ngrid
273	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
274	! ztcond(ig,l)=1./(bcond-acondlog(.0095pplay(ig,l)))
275	ztcond(ig,l)=
276	& 1./(bcond-acondlog(.01vmr_co2(ig,l)*pplay(ig,l)))
277	if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
278	ENDDO
279	ENDDO
280
281	c Condensation/sublimation in the atmosphere
282	c ------------------------------------------
283	c (calcul of zcondicea , zfallice and pdtc)
284	c
285	DO l=nlayer , 1, -1
286	DO ig=1,ngrid
287	pdtc(ig,l)=0.
288	IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
289	condsub(ig)=.true.
290	IF (zfallice(ig,l+1).gt.0) then
291	zfallheat=zfallice(ig,l+1)*
292	& (pphi(ig,l+1)-pphi(ig,l) +
293	& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
294	ELSE
295	zfallheat=0.
296	ENDIF
297	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
298	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
299	& ccondpdtc(ig,l)- zfallheat
300	c Case when the ice from above sublimes entirely
301	c """""""""""""""""""""""""""""""""""""""""""""""
302	IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
303	pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
304	& (ccond*(pplev(ig,l)-pplev(ig,l+1)))
305	zcondicea(ig,l)= -zfallice(ig,l+1)
306	END IF
307
308	zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
309	END IF
310	ENDDO
311	ENDDO
312
313	c *************************
314	c SURFACE CONDENSATION
315	c *************************
316
317	c forecast of ground temperature ztsrf and frost temperature ztcondsol
318	c --------------------------------------------------------------------
319	DO ig=1,ngrid
320	ztcondsol(ig)=
321	& 1./(bcond-acondlog(.01vmr_co2(ig,1)*pplev(ig,1)))
322	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
323	ENDDO
324
325	c
326	c Condensation/sublimation on the ground
327	c --------------------------------------
328	c (calcul of zcondices , pdtsrfc)
329	c
330	DO ig=1,ngrid
331	IF(ig.GT.ngrid/2+1) THEN
332	icap=2
333	ELSE
334	icap=1
335	ENDIF
336
337	c
338	c Loop on where we have condensation/ sublimation
339	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond
340	$ (zfallice(ig,1).NE.0.) .OR. ! falling snow
341	$ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim.
342	$ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
343	condsub(ig) = .true.
344
345	IF (zfallice(ig,1).gt.0) then
346	zfallheat=zfallice(ig,1)*
347	& (pphi(ig,1)- phisfi(ig) +
348	& cpice(ztcond(ig,1)-ztcondsol(ig)))/(latcondptimestep)
349	ELSE
350	zfallheat=0.
351	ENDIF
352
353	c condensation or partial sublimation of CO2 ice
354	c """""""""""""""""""""""""""""""""""""""""""""""
355	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
356	& /(latcond*ptimestep) - zfallheat
357	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
358
359	c If the entire CO_2 ice layer sublimes
360	c """"""""""""""""""""""""""""""""""""""""""""""""""""
361	c (including what has just condensed in the atmosphere)
362
363	IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
364	& -zcondices(ig))THEN
365	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
366	pdtsrfc(ig)=(latcond/pcapcal(ig))*
367	& (zcondices(ig)+zfallheat)
368	END IF
369
370	c Changing CO2 ice amount and pressure :
371	c """"""""""""""""""""""""""""""""""""
372
373	zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
374	piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
375	pdpsrf(ig) = -zdiceco2(ig)*g
376
377	IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
378	PRINT*,'STOP in condens'
379	PRINT*,'condensing more than total mass'
380	PRINT*,'Grid point ',ig
381	PRINT*,'Ps = ',pplev(ig,1)
382	PRINT*,'d Ps = ',pdpsrf(ig)
383	STOP
384	ENDIF
385	END IF ! if there is condensation/sublimmation
386	ENDDO ! of DO ig=1,ngrid
387
388	c ********************************************************************
389	c Surface albedo and emissivity of the surface below the snow (emisref)
390	c ********************************************************************
391	c Prepare the case where albedo varies with insolation:
392	c ----------------------------------------------------
393	! if (fluxdependent) then
394	!
395	c Calcul du flux moyen (zonal mean)
396	! do j=1,jjp1
397	! Fluxmean(j)=0
398	! do i=1,iim
399	! ig=1+(j-2)*iim +i
400	! if(j.eq.1) ig=1
401	! if(j.eq.jjp1) ig=ngrid
402	! Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
403	! $ +fluxsurf_sw(ig,2)
404	! enddo
405	! Fluxmean(j)=Fluxmean(j)/float(iim)
406	! enddo
407	!
408	c const A and B used to calculate the albedo which depends on solar flux
409	c albedice=constA+constB*Flux
410	c constT = time step to calculate the solar flux when flux decreases
411	! constA=0.26
412	c constA=0.33
413	c constA=0.186
414	! constB=0.00187
415	! constT=10
416	! endif ! of if (fluxdependent)
417
418	! Check that amont of CO2 ice is not problematic
419	DO ig=1,ngrid
420	if(.not.piceco2(ig).ge.0.) THEN
421	if(piceco2(ig).le.-5.e-8) print*,
422	$ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
423	piceco2(ig)=0.
424	endif
425	ENDDO
426
427	! Set albedo and emissivity of the surface
428	! ----------------------------------------
429	CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
430
431	c Calcul de l'albedo
432	c ------------------
433	! do ig =1,ngrid
434	! IF(ig.GT.ngrid/2+1) THEN
435	! icap=2
436	! ELSE
437	! icap=1
438	! ENDIF
439	! IF(firstcall2) THEN
440	! albediceF(ig)=albedice(icap)
441	! ENDIF
442	c if there is still co2ice ccccccccccccccccccccccc
443	! if (piceco2(ig).gt.0) then
444	! emisref(ig) = emisice(icap)
445
446	c if flux dependent albedo is used
447	c --------------------------------
448	! if (fluxdependent) then
449	! j=INT((ig-2)/iim)+2
450	! if(ig.eq.1) j=1
451	! if(ig.eq.ngrid) j=jjp1
452	c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
453	c $ constA+constB*250)
454	! albediceF_new(ig)=constA+constB*Fluxmean(j)
455	! if (albediceF(ig).gt.albediceF_new(ig)) then
456	! albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
457	! $ constT)*(albediceF_new(ig)-albediceF(ig))
458	! else
459	! albediceF(ig)=albediceF_new(ig)
460	! endif
461	c if part of the ice is transparent
462	c slopy = pente de la droite: alpha = y*co2ice/1620
463	c pour une valeur superieur a une epaisseur de glace donnee
464	c ici, epaisseur limite = 10cm
465	! if (transparency) then
466	! slopy=1/(1620*0.10)
467	! alpha=MIN(slopy*piceco2(ig),1.)
468	! psolaralb(ig,1) = alpha*albediceF(ig)
469	! $ +(1-alpha)*albedodat(ig)
470	! psolaralb(ig,2) = psolaralb(ig,1)
471	! else
472	! psolaralb(ig,1) = albediceF(ig)
473	! psolaralb(ig,2) = psolaralb(ig,1)
474	! endif
475	! else
476	c transparency set to true and fluxdependent set to false
477	! if (transparency) then
478	! slopy=1/(1620*0.10)
479	! alpha=MIN(slopy*piceco2(ig),1.)
480	! psolaralb(ig,1) = alpha*albedice(icap)
481	! $ +(1-alpha)*albedodat(ig)
482	! psolaralb(ig,2) = psolaralb(ig,1)
483	! else
484	c simplest case: transparency and flux dependent set to false
485	! psolaralb(ig,1) = albedice(icap)
486	! psolaralb(ig,2) = albedice(icap)
487	! endif
488	! endif
489	c no more co2ice, albedo = ground albedo
490	! else
491	! psolaralb(ig,1) = albedodat(ig)
492	! psolaralb(ig,2) = albedodat(ig)
493	! emisref(ig) = emissiv
494	! pemisurf(ig) = emissiv
495	! endif
496	! end do ! end of the ig loop
497
498	! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
499	DO ig=1,ngrid
500	if (piceco2(ig).eq.0) then
501	pemisurf(ig)=emissiv
502	endif
503	ENDDO
504
505	! firstcall2=.false.
506	c ***************************************************************
507	c Correction to account for redistribution between sigma or hybrid
508	c layers when changing surface pressure (and warming/cooling
509	c of the CO2 currently changing phase).
510	c *************************************************************
511
512	DO ig=1,ngrid
513	if (condsub(ig)) then
514	do l=1,nlayer
515	ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep
516	zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep
517	zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep
518	do iq=1,nqmx
519	zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
520	enddo
521	end do
522
523	c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up)
524	c """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
525
526	zmflux(1) = -zcondices(ig)
527	DO l=1,nlayer
528	zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
529	& + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
530	c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
531	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
532	END DO
533
534	c Mass of each layer
535	c ------------------
536	DO l=1,nlayer
537	masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
538	END DO
539
540
541	c Corresponding fluxes for T,U,V,Q
542	c """"""""""""""""""""""""""""""""
543
544	c averaging operator for TRANSPORT
545	c """"""""""""""""""""""""""""""""
546	c Value transfert at the surface interface when condensation
547	c sublimation:
548	ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
549	zum(1) = 0
550	zvm(1) = 0
551	do iq=1,nqmx
552	zqm(1,iq)=0. ! most tracer do not condense !
553	enddo
554	c Special case if one of the tracer is CO2 gas
555	if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
556
557	c Van Leer scheme:
558	DO l=1,nlayer+1
559	w(l)=-zmflux(l)*ptimestep
560	END DO
561	call vl1d(ztc,2.,masse,w,ztm)
562	call vl1d(zu ,2.,masse,w,zum)
563	call vl1d(zv ,2.,masse,w,zvm)
564	do iq=1,nqmx
565	do l=1,nlayer
566	zq1(l)=zq(l,iq)
567	enddo
568	zqm1(1)=zqm(1,iq)
569	call vl1d(zq1,2.,masse,w,zqm1)
570	do l=2,nlayer
571	zq( l,iq)=zq1(l)
572	zqm(l,iq)=zqm1(l)
573	enddo
574	enddo
575
576	c Surface condensation affects low winds
577	if (zmflux(1).lt.0) then
578	zum(1)= zu(1) * (w(1)/masse(1))
579	zvm(1)= zv(1) * (w(1)/masse(1))
580	if (w(1).gt.masse(1)) then ! ensure numerical stability
581	zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
582	zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
583	end if
584	end if
585
586	ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
587	zum(nlayer+1)= zu(nlayer) ! should not be used, but...
588	zvm(nlayer+1)= zv(nlayer) ! should not be used, but...
589	do iq=1,nqmx
590	zqm(nlayer+1,iq)= zq(nlayer,iq)
591	enddo
592
593	#ifdef MESOSCALE
594	!!!! AS: This part must be commented in the mesoscale model
595	!!!! AS: ... to avoid instabilities.
596	!!!! AS: you have to compile with -DMESOSCALE to do so
597	#else
598	c Tendencies on T, U, V, Q
599	c """"""""""""""""""""""""
600	DO l=1,nlayer
601
602	c Tendencies on T
603	zdtsig(ig,l) = (1/masse(l)) *
604	& ( zmflux(l)*(ztm(l) - ztc(l))
605	& - zmflux(l+1)*(ztm(l+1) - ztc(l))
606	& + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) )
607	pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l)
608
609	c Tendencies on U
610	pduc(ig,l) = (1/masse(l)) *
611	& ( zmflux(l)*(zum(l) - zu(l))
612	& - zmflux(l+1)*(zum(l+1) - zu(l)) )
613
614
615	c Tendencies on V
616	pdvc(ig,l) = (1/masse(l)) *
617	& ( zmflux(l)*(zvm(l) - zv(l))
618	& - zmflux(l+1)*(zvm(l+1) - zv(l)) )
619
620	END DO
621	#endif
622
623	c Tendencies on Q
624	do iq=1,nqmx
625	! if (noms(iq).eq.'co2') then
626	if (iq.eq.ico2) then
627	c SPECIAL Case when the tracer IS CO2 :
628	DO l=1,nlayer
629	pdqc(ig,l,iq)= (1/masse(l)) *
630	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
631	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
632	& + zcondicea(ig,l)*(zq(l,iq)-1.) )
633	END DO
634	else
635	DO l=1,nlayer
636	pdqc(ig,l,iq)= (1/masse(l)) *
637	& ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
638	& - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))
639	& + zcondicea(ig,l)*zq(l,iq) )
640	END DO
641	end if
642	enddo
643
644	c --------------------------------------------------------
645	c Roughly Simulate Molecular mixing when CO2 is too depleted by
646	c Surface condensation (mixing starts if qco2 < qco2min )
647	c FF 06/2004
648	c WARNING : this is now done in convadj, better (FF 02/2005)
649	c --------------------------------------------------------
650	flag=0 ! now done in convadj : must be =0
651	if (flag.eq.1) then
652	if(ico2.gt.0) then ! relevant only if one tracer is CO2
653	if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
654	& .lt.qco2min) then
655	do iq=1,nqmx
656	zq(1,iq)=pq(ig,1,iq)
657	& +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
658	Smq(1,iq) = masse(1)*zq(1,iq)
659	end do
660	Sm(1) = masse(1)
661	do l =2,nlayermx
662	do iq=1,nqmx
663	zq(l,iq)=pq(ig,l,iq)
664	& +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
665	smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
666	end do
667	sm(l) = sm(l-1) + masse(l)
668	if(zq(l,ico2).gt.qco2min) then
669	c mixmas: mass of atmosphere that must be mixed to reach qco2min
670	mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
671	& /(zq(l,ico2)-qco2min)
672	if((mixmas.le.sm(l)))then
673	c OK if mixed mass less than mass of layers affected
674	nmix=l ! number of layer affected by mixing
675	goto 99
676	end if
677	end if
678	end do
679	99 continue
680	do iq=1,nqmx
681	qmix=zq(nmix,iq)
682	& +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
683	do l=1,nmix-1
684	pdqc(ig,l,iq)=
685	& (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
686	end do
687	c layer only partly mixed :
688	pdqc(ig,nmix,iq)=(
689	& qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
690	& -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
691
692	end do
693	end if
694	end if
695
696	endif ! (flag.eq.1)
697	end if ! if (condsub)
698	END DO ! loop on ig
699
700	c ***************************************************************
701	c CO2 snow / clouds scheme
702	c ***************************************************************
703
704	call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
705	& zcondicea,zcondices,zfallice,pemisurf)
706
707	c ***************************************************************
708	c Ecriture des diagnostiques
709	c ***************************************************************
710
711	c DO l=1,nlayer
712	c DO ig=1,ngrid
713	c Taux de cond en kg.m-2.pa-1.s-1
714	c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
715	c Taux de cond en kg.m-3.s-1
716	c tconda2(ig,l)=tconda1(ig,l)pplay(ig,l)g/(r*pt(ig,l))
717	c END DO
718	c END DO
719	c call WRITEDIAGFI(ngridmx,'tconda1',
720	c &'Taux de condensation CO2 atmospherique /Pa',
721	c & 'kg.m-2.Pa-1.s-1',3,tconda1)
722	c call WRITEDIAGFI(ngridmx,'tconda2',
723	c &'Taux de condensation CO2 atmospherique /m',
724	c & 'kg.m-3.s-1',3,tconda2)
725
726	! output falling co2 ice in 1st layer:
727	! call WRITEDIAGFI(ngridmx,'fallice',
728	! &'Precipitation of co2 ice',
729	! & 'kg.m-2.s-1',2,zfallice(1,1))
730
731	!! Specific stuff to bound co2 tracer values ....
732	if (bound_qco2.and.(ico2.ne.0)) then
733	do ig=1,ngridmx
734	do l=1,nlayermx
735	zqco2=pq(ig,l,ico2)
736	& +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
737	if (zqco2.gt.qco2max) then
738	! correct pdqc:
739	pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
740	& -pdq(ig,l,ico2)
741	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
742	& " so that co2 conc. does not exceed",qco2max
743	write(,) " ig:",ig," l:",l
744	endif ! of if (zqco2.gt.qco2max)
745	if (zqco2.lt.qco2mini) then
746	! correct pdqc:
747	pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
748	& -pdq(ig,l,ico2)
749	write(,) "newcondens: adapting pdqc(ig,l,ico2)",
750	& " so that co2 conc. is not less than",qco2mini
751	write(,) " ig:",ig," l:",l
752	endif ! of if (zqco2.lt.qco2mini)
753	end do
754	enddo
755	endif ! of if (bound_qco2.and.(ico2.ne.0)) then
756
757	return
758	end
759
760
761
762	c *****************************************************************
763	SUBROUTINE vl1d(q,pente_max,masse,w,qm)
764	c
765	c
766	c Operateur de moyenne inter-couche pour calcul de transport type
767	c Van-Leer " pseudo amont " dans la verticale
768	c q,w sont des arguments d'entree pour le s-pg ....
769	c masse : masse de la couche Dp/g
770	c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
771	c pente_max = 2 conseillee
772	c
773	c
774	c --------------------------------------------------------------------
775	IMPLICIT NONE
776
777	#include "dimensions.h"
778
779	c
780	c
781	c
782	c Arguments:
783	c ----------
784	real masse(llm),pente_max
785	REAL q(llm),qm(llm+1)
786	REAL w(llm+1)
787	c
788	c Local
789	c ---------
790	c
791	INTEGER l
792	c
793	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
794	real sigw, Mtot, MQtot
795	integer m
796	c integer ismax,ismin
797
798
799	c On oriente tout dans le sens de la pression
800	c W > 0 WHEN DOWN !!!!!!!!!!!!!
801
802	do l=2,llm
803	dzqw(l)=q(l-1)-q(l)
804	adzqw(l)=abs(dzqw(l))
805	enddo
806
807	do l=2,llm-1
808	if(dzqw(l)*dzqw(l+1).gt.0.) then
809	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
810	else
811	dzq(l)=0.
812	endif
813	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
814	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
815	enddo
816
817	dzq(1)=0.
818	dzq(llm)=0.
819
820	do l = 1,llm-1
821
822	c Regular scheme (transfered mass < layer mass)
823	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
824	if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
825	sigw=w(l+1)/masse(l+1)
826	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
827	else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
828	sigw=w(l+1)/masse(l)
829	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
830
831	c Extended scheme (transfered mass > layer mass)
832	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
833	else if(w(l+1).gt.0.) then
834	m=l+1
835	Mtot = masse(m)
836	MQtot = masse(m)*q(m)
837	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
838	m=m+1
839	Mtot = Mtot + masse(m)
840	MQtot = MQtot + masse(m)*q(m)
841	end do
842	if (m.lt.llm) then
843	sigw=(w(l+1)-Mtot)/masse(m+1)
844	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)
845	& (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
846	else
847	w(l+1) = Mtot
848	qm(l+1) = Mqtot / Mtot
849	write(,) 'top layer is disapearing !'
850	stop
851	end if
852	else ! if(w(l+1).lt.0)
853	m = l-1
854	Mtot = masse(m+1)
855	MQtot = masse(m+1)*q(m+1)
856	if (m.gt.0) then ! because some compilers will have problems
857	! evaluating masse(0)
858	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
859	m=m-1
860	Mtot = Mtot + masse(m+1)
861	MQtot = MQtot + masse(m+1)*q(m+1)
862	if (m.eq.0) exit
863	end do
864	endif
865	if (m.gt.0) then
866	sigw=(w(l+1)+Mtot)/masse(m)
867	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)
868	& (q(m)-0.5(1.+sigw)dzq(m)) )
869	else
870	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
871	end if
872	end if
873	enddo
874
875	c boundary conditions (not used in newcondens !!)
876	c qm(llm+1)=0.
877	c if(w(1).gt.0.) then
878	c qm(1)=q(1)
879	c else
880	c qm(1)=0.
881	c end if
882
883	return
884	end

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

Download in other formats:

Original Format