Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/mass_redistribution.f90 @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 22.0 KB

Line
1	SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep, &
2	rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs, &
3	pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr, &
4	pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
5
6	USE watercommon_h, Only: Tsat_water,RLVTT
7	use surfdat_h
8	use radcommon_h, only: glat
9	USE comgeomfi_h
10	USE tracer_h
11
12	IMPLICIT NONE
13	!=======================================================================
14	! subject:
15	! --------
16	! Mass and momentum fluxes through sigma levels as the surface pressure is modified are also taken into account
17	!
18	! author: Jeremy Leconte 2012 (from F.Forget 1998)
19	! ------
20	!
21	! input:
22	! -----
23	! ngrid nombre de points de verticales
24	! (toutes les boucles de la physique sont au
25	! moins vectorisees sur ngrid)
26	! nlayer nombre de couches
27	! pplay(ngrid,nlayer) Pressure levels
28	! pplev(ngrid,nlayer+1) Pressure levels
29	! nq Number of tracers
30	!
31	! pt(ngrid,nlayer) temperature (en K)
32	! pq(ngrid,nlayer,nq) tracer specific concentration (kg/kg of air)
33	! pu,pv (ngrid,nlayer) wind velocity (m/s)
34	!
35	!
36	! pdX physical tendency of X before mass redistribution
37	!
38	! pdmassmr air Mass added to the atmosphere in each layer (kg/m2/s)
39	!
40	! output:
41	! -------
42	!
43	! pdXmr(ngrid) physical tendency of X after mass redistribution
44	!
45	!
46	!
47	!=======================================================================
48	!
49	! 0. Declarations :
50	! ------------------
51	!
52	!-----------------------------------------------------------------------
53	! INCLUDE 'dimensions.h'
54	!
55	! dimensions.h contient les dimensions du modele
56	! ndm est tel que iim=2**ndm
57	!-----------------------------------------------------------------------
58
59	INTEGER iim,jjm,llm,ndm
60
61	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
62
63	!-----------------------------------------------------------------------
64	!-----------------------------------------------------------------------
65	! INCLUDE 'dimphys.h'
66
67	! ngridmx : number of horizontal grid points
68	! note: the -1/jjm term will be 0; unless jj=1
69	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
70	! nlayermx : number of atmospheric layers
71	integer, parameter :: nlayermx = llm
72	! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
73	!integer, parameter :: nsoilmx = 4 ! for a test
74	!integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
75	!integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
76	!-----------------------------------------------------------------------
77
78	!-----------------------------------------------------------------------
79	! INCLUDE "comcstfi.h"
80
81	common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
82	common/comcstfi/avocado!,molrad,visc
83
84	real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
85	real avocado!,molrad,visc
86
87	!
88	! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
89	!
90	!-----------------------------------------------------------------------
91	! INCLUDE 'comvert.h'
92
93	COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), &
94	& pa,preff,nivsigs(llm),nivsig(llm+1), &
95	& aps(llm),bps(llm),scaleheight,pseudoalt(llm)
96
97	common/comverti/disvert_type, pressure_exner
98
99	real ap ! hybrid pressure contribution at interlayers
100	real bp ! hybrid sigma contribution at interlayer
101	real presnivs ! (reference) pressure at mid-layers
102	real dpres
103	real pa ! reference pressure (Pa) at which hybrid coordinates
104	! become purely pressure
105	real preff ! reference surface pressure (Pa)
106	real nivsigs
107	real nivsig
108	real aps ! hybrid pressure contribution at mid-layers
109	real bps ! hybrid sigma contribution at mid-layers
110	real scaleheight ! atmospheric (reference) scale height (km)
111	real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
112	! preff and scaleheight
113
114	integer disvert_type ! type of vertical discretization:
115	! 1: Earth (default for planet_type==earth),
116	! automatic generation
117	! 2: Planets (default for planet_type!=earth),
118	! using 'z2sig.def' (or 'esasig.def) file
119
120	logical pressure_exner
121	! compute pressure inside layers using Exner function, else use mean
122	! of pressure values at interfaces
123
124	!-----------------------------------------------------------------------
125	!
126	! $Header$
127	!
128	!
129	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
130	! veillez n'utiliser que des ! pour les commentaires
131	! et bien positionner les & des lignes de continuation
132	! (les placer en colonne 6 et en colonne 73)
133	!
134	!
135	!-----------------------------------------------------------------------
136	! INCLUDE 'paramet.h'
137
138	INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
139	INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
140	INTEGER ijmllm,mvar
141	INTEGER jcfil,jcfllm
142
143	PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 &
144	& ,jjp1=jjm+1-1/jjm)
145	PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 )
146	PARAMETER( kftd = iim/2 -ndm )
147	PARAMETER( ip1jm = iip1jjm, ip1jmp1= iip1jjp1 )
148	PARAMETER( ip1jmi1= ip1jm - iip1 )
149	PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
150	PARAMETER( mvar= ip1jmp1( 2llm+1) + ijmllm )
151	PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
152
153	!-----------------------------------------------------------------------
154	!
155	! For Fortran 77/Fortran 90 compliance always use line continuation
156	! symbols '&' in columns 73 and 6
157	!
158	! Group commons according to their type for minimal performance impact
159
160	COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj &
161	& , co2cond,callsoil &
162	& , season,diurnal,tlocked,rings_shadow,lwrite &
163	& , callstats,calleofdump &
164	& , enertest &
165	& , callgasvis,continuum,H2Ocont_simple,graybody &
166	& , radfixed &
167	& , meanOLR, specOLR &
168	& , kastprof &
169	& , nosurf, oblate &
170	& , newtonian, testradtimes &
171	& , check_cpp_match, force_cpp &
172	& , rayleigh &
173	& , stelbbody &
174	& , nearco2cond &
175	& , tracer, mass_redistrib, varactive, varfixed &
176	& , sedimentation,water,watercond,waterrain &
177	& , aeroco2,aeroh2o,aeroh2so4,aeroback2lay &
178	& , aerofixco2,aerofixh2o &
179	& , hydrology, sourceevol &
180	& , CLFvarying &
181	& , strictboundcorrk &
182	& , ok_slab_ocean &
183	& , ok_slab_sic &
184	& , ok_slab_heat_transp
185
186
187	COMMON/callkeys_i/iaervar,iddist,iradia,startype
188
189	COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb, &
190	& Tstrat,tplanet,obs_tau_col_tropo, &
191	& obs_tau_col_strato,pres_bottom_tropo, &
192	& pres_top_tropo,pres_bottom_strato, &
193	& pres_top_strato,size_tropo,size_strato,satval, &
194	& CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
195	& maxicethick,Tsaldiff,tau_relax,cloudlvl, &
196	& icetstep,intheat,flatten,Rmean,J2,MassPlanet
197
198	logical callrad,corrk,calldifv,UseTurbDiff &
199	& , calladj,co2cond,callsoil &
200	& , season,diurnal,tlocked,rings_shadow,lwrite &
201	& , callstats,calleofdump &
202	& , callgasvis,continuum,H2Ocont_simple,graybody &
203	& , strictboundcorrk
204
205	logical enertest
206	logical nonideal
207	logical meanOLR
208	logical specOLR
209	logical kastprof
210	logical newtonian
211	logical check_cpp_match
212	logical force_cpp
213	logical testradtimes
214	logical rayleigh
215	logical stelbbody
216	logical ozone
217	logical nearco2cond
218	logical tracer
219	logical mass_redistrib
220	logical varactive
221	logical varfixed
222	logical radfixed
223	logical sedimentation
224	logical water,watercond,waterrain
225	logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
226	logical aerofixco2,aerofixh2o
227	logical hydrology
228	logical sourceevol
229	logical CLFvarying
230	logical nosurf
231	logical oblate
232	logical ok_slab_ocean
233	logical ok_slab_sic
234	logical ok_slab_heat_transp
235
236	integer iddist
237	integer iaervar
238	integer iradia
239	integer startype
240
241	real topdustref
242	real Nmix_co2
243	real dusttau
244	real Fat1AU
245	real stelTbb
246	real Tstrat
247	real tplanet
248	real obs_tau_col_tropo
249	real obs_tau_col_strato
250	real pres_bottom_tropo
251	real pres_top_tropo
252	real pres_bottom_strato
253	real pres_top_strato
254	real size_tropo
255	real size_strato
256	real satval
257	real CLFfixval
258	real n2mixratio
259	real co2supsat
260	real pceil
261	real albedosnow
262	real maxicethick
263	real Tsaldiff
264	real tau_relax
265	real cloudlvl
266	real icetstep
267	real intheat
268	real flatten
269	real Rmean
270	real J2
271	real MassPlanet
272
273	!-----------------------------------------------------------------------
274	! Arguments :
275	! ---------
276	INTEGER ngrid, nlayer, nq
277	REAL ptimestep
278	REAL pcapcal(ngrid)
279	INTEGER rnat(ngrid)
280	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
281	REAL pt(ngrid,nlayer),pdt(ngrid,nlayer)
282	REAL ptsrf(ngrid),pdtsrf(ngrid)
283	REAL pdtmr(ngrid,nlayer)
284	REAL pu(ngrid,nlayer) , pv(ngrid,nlayer)
285	REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
286	REAL pdmassmr(ngrid,nlayer)
287	REAL pdumr(ngrid,nlayer) , pdvmr(ngrid,nlayer)
288	REAL pq(ngrid,nlayer,nq),pdq(ngrid,nlayer,nq)
289	REAL pqs(ngrid,nq)
290	REAL pdqmr(ngrid,nlayer,nq),pdqsmr(ngrid,nq)
291	REAL pdpsrfmr(ngrid),pdtsrfmr(ngrid)
292	!
293	! Local variables :
294	! -----------------
295
296	! Boiling/sublimation
297	REAL Tsat(ngrid),zmassboil(ngrid)
298
299	! vertical reorganization of sigma levels
300	REAL zzu(nlayermx),zzv(nlayermx)
301	REAL zzq(nlayermx,nq),zzt(nlayermx)
302	! Dummy variables
303	INTEGER n,l,ig,iq
304	REAL zdtsig(ngrid,nlayermx)
305	REAL zmass(ngrid,nlayermx),zzmass(nlayermx),w(nlayermx+1)
306	REAL zdmass_sum(ngrid,nlayermx+1)
307	REAL zmflux(nlayermx+1)
308	REAL zq1(nlayermx)
309	REAL ztsrf(ngrid)
310	REAL ztm(nlayermx+1)
311	REAL zum(nlayermx+1) , zvm(nlayermx+1)
312	REAL zqm(nlayermx+1,nq),zqm1(nlayermx+1)
313
314	! local saved variables
315	LOGICAL, SAVE :: firstcall=.true.
316
317	!----------------------------------------------------------------------
318
319	! Initialisation
320	! --------------
321	!
322	IF (firstcall) THEN
323	firstcall=.false.
324	ENDIF
325	!
326	!======================================================================
327	! Calcul of h2o condensation
328	! ============================================================
329	!
330	! Used variable :
331	! pdmassmr : air Mass added to the atmosphere in each layer per unit time (kg/m2/s)
332	! zdmass_sum(ngrid,l) : total air mass added to the atm above layer l per unit time (kg/m2/s)
333	!
334	!
335	! Surface tracer Tendencies set to 0
336	! -------------------------------------
337	pdqsmr(1:ngrid,1:nq)=0.
338
339	ztsrf(1:ngrid) = ptsrf(1:ngrid) + pdtsrf(1:ngrid)*ptimestep
340
341
342	DO ig=1,ngrid
343	zdmass_sum(ig,nlayermx+1)=0.
344	DO l = nlayermx, 1, -1
345	zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/glat(ig)
346	zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
347	END DO
348	END DO
349
350
351	if (water) then
352	do ig=1,ngrid
353	call Tsat_water(pplev(ig,1)+zdmass_sum(ig,1)gptimestep,Tsat(ig))
354	enddo
355	call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat)
356
357	do ig=1,ngrid
358	if (ztsrf(ig).gt.Tsat(ig)) then
359	zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT/ptimestep
360	if ((zmassboil(ig)*ptimestep.gt.pqs(ig,igcm_h2o_vap)).and.(rnat(ig).eq.1)) then
361	zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep
362	endif
363	zmassboil(ig)=zmassboil(ig)*0.0 !momentary, should be 1. JL12
364	pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig)
365	pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig)
366	ztsrf(ig)=ptsrf(ig)+pdtsrfmr(ig)*ptimestep
367	else
368	zmassboil(ig)=0.
369	pdtsrfmr(ig)=0.
370	endif
371	enddo
372	endif
373
374	! *************************
375	! UPDATE SURFACE
376	! *************************
377	! Changing pressure at the surface:
378	! """"""""""""""""""""""""""""""""""""
379
380	pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g
381
382	do ig = 1, ngrid
383	IF(ABS(pdpsrfmr(ig)*ptimestep).GT.pplev(ig,1)) THEN
384	PRINT*,'STOP in condens'
385	PRINT*,'condensing more than total mass'
386	PRINT*,'Grid point ',ig
387	PRINT*,'Ps = ',pplev(ig,1)
388	PRINT,'d Ps = ',pdpsrfmr(ig)ptimestep
389	STOP
390	ENDIF
391	enddo ! of DO ig=1,ngrid
392
393
394	! ***************************************************************
395	! Correction to account for redistribution between sigma or hybrid
396	! layers when changing surface pressure
397	! zzx quantities have dimension (nlayermx) to speed up calculation
398	! *************************************************************
399
400	DO ig=1,ngrid
401	zzt(1:nlayermx) = pt(ig,1:nlayermx) + pdt(ig,1:nlayermx) * ptimestep
402	zzu(1:nlayermx) = pu(ig,1:nlayermx) + pdu(ig,1:nlayermx) * ptimestep
403	zzv(1:nlayermx) = pv(ig,1:nlayermx) + pdv(ig,1:nlayermx) * ptimestep
404	zzq(1:nlayermx,1:nq)=pq(ig,1:nlayermx,1:nq)+pdq(ig,1:nlayermx,1:nq)*ptimestep ! must add the water that has fallen???
405
406	! Mass fluxes of air through the sigma levels (kg.m-2.s-1) (>0 when up)
407	! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
408
409	zmflux(1) = zmassboil(ig)
410	DO l=1,nlayer
411	zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
412	! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
413	if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
414	END DO
415
416	! Mass of each layer
417	! ------------------
418	zzmass(1:nlayermx)=zmass(ig,1:nlayermx)(1.+pdpsrfmr(ig)ptimestep/pplev(ig,1))
419
420
421	! Corresponding fluxes for T,U,V,Q
422	! """"""""""""""""""""""""""""""""
423
424	! averaging operator for TRANSPORT
425	! """"""""""""""""""""""""""""""""
426	! Value transfert at the surface interface when condensation
427	! sublimation:
428	ztm(1) = ztsrf(ig)
429	zum(1) = 0.
430	zvm(1) = 0.
431	zqm(1,1:nq)=0. ! most tracer do not condense !
432	if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
433
434	! Van Leer scheme:
435	w(1:nlayermx+1)=-zmflux(1:nlayermx+1)*ptimestep
436	call vl1d(zzt,2.,zzmass,w,ztm)
437	call vl1d(zzu ,2.,zzmass,w,zum)
438	call vl1d(zzv ,2.,zzmass,w,zvm)
439	do iq=1,nq
440	zq1(1:nlayermx)=zzq(1:nlayermx,iq)
441	zqm1(1)=zqm(1,iq)
442	! print*,iq
443	! print*,zq1
444	call vl1d(zq1,2.,zzmass,w,zqm1)
445	do l=2,nlayer
446	zzq(l,iq)=zq1(l)
447	zqm(l,iq)=zqm1(l)
448	enddo
449	enddo
450
451	! Surface condensation affects low winds
452	if (zmflux(1).lt.0) then
453	zum(1)= zzu(1) * (w(1)/zzmass(1))
454	zvm(1)= zzv(1) * (w(1)/zzmass(1))
455	if (w(1).gt.zzmass(1)) then ! ensure numerical stability
456	zum(1)= (zzu(1)-zum(2))*zzmass(1)/w(1) + zum(2)
457	zvm(1)= (zzv(1)-zvm(2))*zzmass(1)/w(1) + zvm(2)
458	end if
459	end if
460
461	ztm(nlayer+1)= zzt(nlayer) ! should not be used, but...
462	zum(nlayer+1)= zzu(nlayer) ! should not be used, but...
463	zvm(nlayer+1)= zzv(nlayer) ! should not be used, but...
464	zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq)
465
466	! Tendencies on T, U, V, Q
467	! """"""""""""""""""""""""
468	DO l=1,nlayer
469
470	! Tendencies on T
471	pdtmr(ig,l) = (1/zzmass(l)) * &
472	(zmflux(l)(ztm(l) - zzt(l))-zmflux(l+1)(ztm(l+1)-zzt(l)))
473	!JL12 the last term in Newcondens has been set to zero because we are only dealing with redistribution here
474
475	! Tendencies on U
476	pdumr(ig,l) = (1/zzmass(l)) ( zmflux(l)(zum(l) - zzu(l)) - zmflux(l+1)*(zum(l+1) - zzu(l)) )
477
478	! Tendencies on V
479	pdvmr(ig,l) = (1/zzmass(l)) ( zmflux(l)(zvm(l) - zzv(l)) - zmflux(l+1)*(zvm(l+1) - zzv(l)) )
480
481	END DO
482
483	! Tendencies on Q
484	do iq=1,nq
485	DO l=1,nlayer
486	pdqmr(ig,l,iq)= (1/zzmass(l)) * &
487	(zmflux(l)(zqm(l,iq)-zzq(l,iq))- zmflux(l+1)(zqm(l+1,iq)-zzq(l,iq)) - pdmassmr(ig,l)*zzq(l,iq))
488	END DO
489	enddo
490
491	END DO ! loop on ig
492
493	return
494	end
495
496
497
498	! *****************************************************************
499	SUBROUTINE vl1d(q,pente_max,zzmass,w,qm)
500	!
501	!
502	! Operateur de moyenne inter-couche pour calcul de transport type
503	! Van-Leer " pseudo amont " dans la verticale
504	! q,w sont des arguments d'entree pour le s-pg ....
505	! masse : masse de la couche Dp/g
506	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
507	! pente_max = 2 conseillee
508	!
509	!
510	! --------------------------------------------------------------------
511	IMPLICIT NONE
512
513	!-----------------------------------------------------------------------
514	! INCLUDE 'dimensions.h'
515	!
516	! dimensions.h contient les dimensions du modele
517	! ndm est tel que iim=2**ndm
518	!-----------------------------------------------------------------------
519
520	INTEGER iim,jjm,llm,ndm
521
522	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
523
524	!-----------------------------------------------------------------------
525
526	! Arguments:
527	! ----------
528	real zzmass(llm),pente_max
529	REAL q(llm),qm(llm+1)
530	REAL w(llm+1)
531	!
532	! Local
533	! ---------
534	!
535	INTEGER l
536	!
537	real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
538	real sigw, Mtot, MQtot
539	integer m
540	! integer ismax,ismin
541
542
543	! On oriente tout dans le sens de la pression
544	! W > 0 WHEN DOWN !!!!!!!!!!!!!
545
546	do l=2,llm
547	dzqw(l)=q(l-1)-q(l)
548	adzqw(l)=abs(dzqw(l))
549	enddo
550
551	do l=2,llm-1
552	if(dzqw(l)*dzqw(l+1).gt.0.) then
553	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
554	else
555	dzq(l)=0.
556	endif
557	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
558	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
559	enddo
560
561	dzq(1)=0.
562	dzq(llm)=0.
563
564	do l = 1,llm-1
565
566	! Regular scheme (transfered mass < layer mass)
567	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
568	if(w(l+1).gt.0. .and. w(l+1).le.zzmass(l+1)) then
569	sigw=w(l+1)/zzmass(l+1)
570	qm(l+1)=(q(l+1)+0.5(1.-sigw)dzq(l+1))
571	else if(w(l+1).le.0. .and. -w(l+1).le.zzmass(l)) then
572	sigw=w(l+1)/zzmass(l)
573	qm(l+1)=(q(l)-0.5(1.+sigw)dzq(l))
574
575	! Extended scheme (transfered mass > layer mass)
576	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
577	else if(w(l+1).gt.0.) then
578	m=l+1
579	Mtot = zzmass(m)
580	MQtot = zzmass(m)*q(m)
581	do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+zzmass(m+1))))
582	m=m+1
583	Mtot = Mtot + zzmass(m)
584	MQtot = MQtot + zzmass(m)*q(m)
585	end do
586	if (m.lt.llm) then
587	sigw=(w(l+1)-Mtot)/zzmass(m+1)
588	qm(l+1)= (1/w(l+1))(MQtot + (w(l+1)-Mtot)(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
589	else
590	! w(l+1) = Mtot
591	! qm(l+1) = Mqtot / Mtot
592	write(,) 'top layer is disappearing !',l,Mtot,w(l+1),qm(l+1)
593	print*,zzmass
594	print*,w
595	print*,q
596	print*,qm
597	stop
598	end if
599	else ! if(w(l+1).lt.0)
600	m = l-1
601	Mtot = zzmass(m+1)
602	MQtot = zzmass(m+1)*q(m+1)
603	if (m.gt.0) then ! because some compilers will have problems
604	! evaluating zzmass(0)
605	do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+zzmass(m))))
606	m=m-1
607	Mtot = Mtot + zzmass(m+1)
608	MQtot = MQtot + zzmass(m+1)*q(m+1)
609	if (m.eq.0) exit
610	end do
611	endif
612	if (m.gt.0) then
613	sigw=(w(l+1)+Mtot)/zzmass(m)
614	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)(q(m)-0.5(1.+sigw)dzq(m)) )
615	else
616	qm(l+1)= (-1/w(l+1))(MQtot + (-w(l+1)-Mtot)qm(1))
617	end if
618	end if
619	enddo
620
621	! boundary conditions (not used in newcondens !!)
622	! qm(llm+1)=0.
623	! if(w(1).gt.0.) then
624	! qm(1)=q(1)
625	! else
626	! qm(1)=0.
627	! end if
628
629	return
630	end

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

Download in other formats: