source: trunk/LMDZ.PLUTO/libf/phypluto/mass_redistribution_mod.F90 @ 3380

Last change on this file since 3380 was 3275, checked in by afalco, 8 months ago

Pluto PCM:
Changed _vap to _gas;
Included surfprop.F90;
callcorrk includes methane
AF

File size: 12.7 KB
Line 
1MODULE mass_redistribution_mod
2
3CONTAINS
4
5SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
6                       pcapcal,pplay,pplev,pt,ptsrf,pq,pqs,     &
7                       pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
8                       pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
9
10       use radcommon_h, only: glat
11       USE planete_mod, only: bp
12       use comcstfi_mod, only: g
13
14       IMPLICIT NONE
15!=======================================================================
16!   subject:
17!   --------
18!     Mass and momentum fluxes through sigma levels as the surface pressure is modified are also taken into account
19!
20!   author:   Jeremy Leconte 2012 (from F.Forget 1998)
21!   ------
22!
23!   input:
24!   -----
25!    ngrid                 nombre de points de verticales
26!                          (toutes les boucles de la physique sont au
27!                          moins vectorisees sur ngrid)
28!    nlayer                nombre de couches
29!    pplay(ngrid,nlayer)   Pressure levels
30!    pplev(ngrid,nlayer+1) Pressure levels
31!    nq                    Number of tracers
32!
33!    pt(ngrid,nlayer)      temperature (en K)
34!    pq(ngrid,nlayer,nq)   tracer specific concentration (kg/kg of air)
35!    pu,pv (ngrid,nlayer)  wind velocity (m/s)
36!
37!
38!    pdX                   physical tendency of X before mass redistribution
39!
40!    pdmassmr                air Mass added to the atmosphere in each layer (kg/m2/s)
41!
42!   output:
43!   -------
44!
45!    pdXmr(ngrid)           physical tendency of X after mass redistribution
46!
47!
48!
49!=======================================================================
50!
51!    0.  Declarations :
52!    ------------------
53
54!-----------------------------------------------------------------------
55!    Arguments :
56!    ---------
57      INTEGER,INTENT(IN) :: ngrid, nlayer, nq
58      REAL,INTENT(IN) :: ptimestep
59      REAL,INTENT(IN) :: pcapcal(ngrid)
60      REAL,INTENT(IN) :: pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
61      REAL,INTENT(IN) :: pt(ngrid,nlayer),pdt(ngrid,nlayer)
62      REAL,INTENT(IN) :: ptsrf(ngrid),pdtsrf(ngrid)
63      REAL,INTENT(OUT) :: pdtmr(ngrid,nlayer)
64      REAL,INTENT(IN) :: pu(ngrid,nlayer) ,  pv(ngrid,nlayer)
65      REAL,INTENT(IN) :: pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
66      REAL,INTENT(IN) :: pdmassmr(ngrid,nlayer)
67      REAL,INTENT(OUT) :: pdumr(ngrid,nlayer) , pdvmr(ngrid,nlayer)
68      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq),pdq(ngrid,nlayer,nq)
69      REAL,INTENT(IN) :: pqs(ngrid,nq)
70      REAL,INTENT(OUT) :: pdqmr(ngrid,nlayer,nq),pdqsmr(ngrid,nq)
71      REAL,INTENT(OUT) :: pdpsrfmr(ngrid),pdtsrfmr(ngrid)
72!
73!    Local variables :
74!    -----------------
75
76!    Boiling/sublimation
77      REAL Tsat(ngrid),zmassboil(ngrid)
78
79!    vertical reorganization of sigma levels
80      REAL zzu(nlayer),zzv(nlayer)
81      REAL zzq(nlayer,nq),zzt(nlayer)
82!    Dummy variables
83      INTEGER n,l,ig,iq
84      REAL zdtsig(ngrid,nlayer)
85      REAL zmass(ngrid,nlayer),zzmass(nlayer),w(nlayer+1)
86      REAL zdmass_sum(ngrid,nlayer+1)
87      REAL zmflux(nlayer+1)
88      REAL zq1(nlayer)
89      REAL ztsrf(ngrid)
90      REAL ztm(nlayer+1)
91      REAL zum(nlayer+1) , zvm(nlayer+1)
92      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
93      REAL sigma(nlayer+1)
94
95!   local saved variables
96      LOGICAL, SAVE :: firstcall=.true.
97!$OMP THREADPRIVATE(firstcall)
98
99!----------------------------------------------------------------------
100
101!   Initialisation
102!   --------------
103!
104      IF (firstcall) THEN
105         firstcall=.false.
106      ENDIF
107!
108!======================================================================
109!    Calcul of h2o condensation
110!    ============================================================
111!
112!    Used variable :
113!       pdmassmr      : air Mass added to the atmosphere in each layer per unit time (kg/m2/s)
114!       zdmass_sum(ngrid,l) : total air mass added to the atm above layer l per unit time (kg/m2/s)
115!
116!
117!     Surface tracer Tendencies set to 0
118!     -------------------------------------
119      pdqsmr(1:ngrid,1:nq)=0.
120
121      ztsrf(1:ngrid) = ptsrf(1:ngrid) + pdtsrf(1:ngrid)*ptimestep
122
123
124      DO ig=1,ngrid
125         zdmass_sum(ig,nlayer+1)=0.
126         DO l = nlayer, 1, -1
127           zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/glat(ig)
128           zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
129         END DO
130      END DO
131
132!     *************************
133!           UPDATE SURFACE
134!     *************************
135!    Changing pressure at the surface:
136!    """"""""""""""""""""""""""""""""""""
137
138      zmassboil(1:ngrid) = 0 ! AF: TODO: update this for pluto
139      pdtsrfmr(1:ngrid) = 0
140
141      pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g
142
143      do ig = 1, ngrid
144        IF(ABS(pdpsrfmr(ig)*ptimestep).GT.pplev(ig,1)) THEN
145         PRINT*,'STOP in condens in mass_redistribution'
146         PRINT*,'condensing more than total mass'
147         PRINT*,'Grid point ',ig
148         PRINT*,'Ps = ',pplev(ig,1)
149         PRINT*,'d Ps = ',pdpsrfmr(ig)*ptimestep
150         STOP
151        ENDIF
152      enddo ! of DO ig=1,ngrid
153
154
155! ***************************************************************
156!  Correction to account for redistribution between sigma or hybrid
157!  layers when changing surface pressure
158!  zzx quantities have dimension (nlayer) to speed up calculation
159! *************************************************************
160
161      DO ig=1,ngrid
162         zzt(1:nlayer)  = pt(ig,1:nlayer) + pdt(ig,1:nlayer) * ptimestep
163         zzu(1:nlayer)  = pu(ig,1:nlayer) + pdu(ig,1:nlayer) * ptimestep
164         zzv(1:nlayer)  = pv(ig,1:nlayer) + pdv(ig,1:nlayer) * ptimestep
165         zzq(1:nlayer,1:nq)=pq(ig,1:nlayer,1:nq)+pdq(ig,1:nlayer,1:nq)*ptimestep ! must add the water that has fallen???
166
167!  Mass fluxes of air through the sigma levels (kg.m-2.s-1)  (>0 when up)
168!  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
169         zmflux(1) = zmassboil(ig)
170         sigma(1)=1
171         DO l=1,nlayer
172           ! Ehouarn: shouldn't we rather compute sigma levels than use bp()?
173!           sigma(l+1)=pplev(ig,l+1)/pplev(ig,1)
174!           zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - &
175!                        (sigma(l)-sigma(l+1))*(zdmass_sum(ig,1)+zmflux(1))
176!            if (abs(zmflux(l+1)).lt.1E-13.OR.sigma(l+1).eq.0.) zmflux(l+1)=0.
177           ! Ehouarn: but for now leave things as before
178            zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
179! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
180            if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
181         END DO
182
183! Mass of each layer
184! ------------------
185         zzmass(1:nlayer)=zmass(ig,1:nlayer)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1))
186
187
188!  Corresponding fluxes for T,U,V,Q
189!  """"""""""""""""""""""""""""""""
190
191!        averaging operator for TRANSPORT
192!        """"""""""""""""""""""""""""""""
193!        Value transfert at the surface interface when condensation
194!        sublimation:
195         ztm(1) = ztsrf(ig)
196         zum(1) = 0.
197         zvm(1) = 0.
198         zqm(1,1:nq)=0. ! most tracer do not condense !
199         !! if (water) zqm(1,igcm_h2o_gas)=1. ! flux is 100% h2o at surface
200
201!        Van Leer scheme:
202         w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
203         call vl1d(nlayer,zzt,2.,zzmass,w,ztm)
204         call vl1d(nlayer,zzu,2.,zzmass,w,zum)
205         call vl1d(nlayer,zzv,2.,zzmass,w,zvm)
206         do iq=1,nq
207           zq1(1:nlayer)=zzq(1:nlayer,iq)
208           zqm1(1)=zqm(1,iq)
209!               print*,iq
210!               print*,zq1
211           call vl1d(nlayer,zq1,2.,zzmass,w,zqm1)
212           do l=2,nlayer
213              zzq(l,iq)=zq1(l)
214              zqm(l,iq)=zqm1(l)
215           enddo
216         enddo
217
218!        Surface condensation affects low winds
219         if (zmflux(1).lt.0) then
220            zum(1)= zzu(1) *  (w(1)/zzmass(1))
221            zvm(1)= zzv(1) *  (w(1)/zzmass(1))
222            if (w(1).gt.zzmass(1)) then ! ensure numerical stability
223               zum(1)= (zzu(1)-zum(2))*zzmass(1)/w(1) + zum(2)
224               zvm(1)= (zzv(1)-zvm(2))*zzmass(1)/w(1) + zvm(2)
225            end if
226         end if
227
228         ztm(nlayer+1)= zzt(nlayer) ! should not be used, but...
229         zum(nlayer+1)= zzu(nlayer)  ! should not be used, but...
230         zvm(nlayer+1)= zzv(nlayer)  ! should not be used, but...
231         zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq)
232
233!        Tendencies on T, U, V, Q
234!        """"""""""""""""""""""""
235         DO l=1,nlayer
236
237!           Tendencies on T
238            pdtmr(ig,l) = (1/zzmass(l)) *   &
239                (zmflux(l)*(ztm(l) - zzt(l))-zmflux(l+1)*(ztm(l+1)-zzt(l)))
240                  !JL12 the last term in Newcondens has been set to zero because we are only dealing with redistribution here
241
242!           Tendencies on U
243            pdumr(ig,l)   = (1/zzmass(l)) *( zmflux(l)*(zum(l) - zzu(l)) - zmflux(l+1)*(zum(l+1) - zzu(l)) )
244
245!           Tendencies on V
246            pdvmr(ig,l)   = (1/zzmass(l)) *( zmflux(l)*(zvm(l) - zzv(l)) - zmflux(l+1)*(zvm(l+1) - zzv(l)) )
247
248         END DO
249
250!        Tendencies on Q
251         do iq=1,nq
252            DO l=1,nlayer
253               pdqmr(ig,l,iq)= (1/zzmass(l)) *   &
254                   (zmflux(l)*(zqm(l,iq)-zzq(l,iq))- zmflux(l+1)*(zqm(l+1,iq)-zzq(l,iq)) - pdmassmr(ig,l)*zzq(l,iq))
255            END DO
256         enddo
257
258      END DO  ! loop on ig
259
260CONTAINS
261
262! *****************************************************************
263      SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm)
264!
265!
266!     Operateur de moyenne inter-couche pour calcul de transport type
267!     Van-Leer " pseudo amont " dans la verticale
268!    q,w sont des arguments d'entree  pour le s-pg ....
269!    masse : masse de la couche Dp/g
270!    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
271!    pente_max = 2 conseillee
272!
273!
274!   --------------------------------------------------------------------
275
276      IMPLICIT NONE
277
278!   Arguments:
279!   ----------
280      integer,intent(in) :: llm
281      real zzmass(llm),pente_max
282      REAL q(llm),qm(llm+1)
283      REAL w(llm+1)
284!
285!      Local
286!   ---------
287!
288      INTEGER l
289!
290      real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
291      real sigw, Mtot, MQtot
292      integer m
293!     integer ismax,ismin
294
295
296!    On oriente tout dans le sens de la pression
297!     W > 0 WHEN DOWN !!!!!!!!!!!!!
298
299      do l=2,llm
300            dzqw(l)=q(l-1)-q(l)
301            adzqw(l)=abs(dzqw(l))
302      enddo
303
304      do l=2,llm-1
305            if(dzqw(l)*dzqw(l+1).gt.0.) then
306                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
307            else
308                dzq(l)=0.
309            endif
310            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
311            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
312      enddo
313
314         dzq(1)=0.
315         dzq(llm)=0.
316
317       do l = 1,llm-1
318
319!         Regular scheme (transfered mass < layer mass)
320!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
321          if(w(l+1).gt.0. .and. w(l+1).le.zzmass(l+1)) then
322             sigw=w(l+1)/zzmass(l+1)
323             qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
324          else if(w(l+1).le.0. .and. -w(l+1).le.zzmass(l)) then
325             sigw=w(l+1)/zzmass(l)
326             qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
327
328!         Extended scheme (transfered mass > layer mass)
329!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
330          else if(w(l+1).gt.0.) then
331             m=l+1
332             Mtot = zzmass(m)
333             MQtot = zzmass(m)*q(m)
334             do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+zzmass(m+1))))
335                m=m+1
336                Mtot = Mtot + zzmass(m)
337                MQtot = MQtot + zzmass(m)*q(m)
338             end do
339             if (m.lt.llm) then
340                sigw=(w(l+1)-Mtot)/zzmass(m+1)
341                qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*(q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
342             else
343!                w(l+1) = Mtot
344!                qm(l+1) = Mqtot / Mtot
345                write(*,*) 'top layer is disappearing !',l,Mtot,w(l+1),qm(l+1)
346                print*,zzmass
347                print*,w
348                print*,q
349                print*,qm
350               stop
351             end if
352          else      ! if(w(l+1).lt.0)
353             m = l-1
354             Mtot = zzmass(m+1)
355             MQtot = zzmass(m+1)*q(m+1)
356             if (m.gt.0) then ! because some compilers will have problems
357                              ! evaluating zzmass(0)
358              do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+zzmass(m))))
359                m=m-1
360                Mtot = Mtot + zzmass(m+1)
361                MQtot = MQtot + zzmass(m+1)*q(m+1)
362                if (m.eq.0) exit
363              end do
364             endif
365             if (m.gt.0) then
366                sigw=(w(l+1)+Mtot)/zzmass(m)
367                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*(q(m)-0.5*(1.+sigw)*dzq(m))  )
368             else
369                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
370             end if
371          end if
372       enddo
373
374!     boundary conditions (not used in newcondens !!)
375!         qm(llm+1)=0.
376!         if(w(1).gt.0.) then
377!            qm(1)=q(1)
378!         else
379!           qm(1)=0.
380!         end if
381
382       END SUBROUTINE vl1d
383
384END SUBROUTINE mass_redistribution
385
386END MODULE mass_redistribution_mod
Note: See TracBrowser for help on using the repository browser.