Sources/dyn3d/vlyqs.f

      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
!
!     Auteurs:   P.Le Van, F.Hourdin, F.Forget
!
!    ********************************************************************
!     Shema  d'advection " pseudo amont " .
!    ********************************************************************
!     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
!     qsat            est   un argument de sortie pour le s-pg ....
!
!
!   --------------------------------------------------------------------
!
      use dimens_m
      use paramet_m
      use comconst
      use disvert_m
      use conf_gcm_m
      use comgeom
      USE nr_util, ONLY : pi
      IMPLICIT NONE
!
!
!   Arguments:
!   ----------
      REAL masse(ip1jmp1,llm),pente_max
      REAL masse_adv_v( ip1jm,llm)
      REAL q(ip1jmp1,llm)
      REAL qsat(ip1jmp1,llm)
!
!      Local
!   ---------
!
      INTEGER i,ij,l
!
      REAL airej2,airejjm,airescb(iim),airesch(iim)
      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
      REAL qbyv(ip1jm,llm)

      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
!     REAL newq,oldmasse
      Logical first,testcpu
      REAL temps0,temps1,temps2,temps3,temps4,temps5
      SAVE temps0,temps1,temps2,temps3,temps4,temps5
      SAVE first,testcpu

      REAL convpn,convps,convmpn,convmps
      REAL sinlon(iip1),sinlondlon(iip1)
      REAL coslon(iip1),coslondlon(iip1)
      SAVE sinlon,coslon,sinlondlon,coslondlon
      SAVE airej2,airejjm
!
!
      REAL      SSUM

      DATA first,testcpu/.true.,.false./
      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./

      IF(first) THEN
         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
         first=.false.
         do i=2,iip1
            coslon(i)=cos(rlonv(i))
            sinlon(i)=sin(rlonv(i))
            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
         ENDDO
         coslon(1)=coslon(iip1)
         coslondlon(1)=coslondlon(iip1)
         sinlon(1)=sinlon(iip1)
         sinlondlon(1)=sinlondlon(iip1)
         airej2 = SSUM( iim, aire(iip2), 1 )
         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
      ENDIF

!


      DO l = 1, llm
!
!   --------------------------------
!      CALCUL EN LATITUDE
!   --------------------------------

!   On commence par calculer la valeur du traceur moyenne sur le premier cercle
!   de latitude autour du pole (qpns pour le pole nord et qpsn pour
!    le pole nord) qui sera utilisee pour evaluer les pentes au pole.

      DO i = 1, iim
      airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
      ENDDO
      qpns   = SSUM( iim,  airescb ,1 ) / airej2
      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm

!   calcul des pentes aux points v

      DO ij=1,ip1jm
         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
         adyqv(ij)=abs(dyqv(ij))
      ENDDO

!   calcul des pentes aux points scalaires

      DO ij=iip2,ip1jm
         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
         dyqmax(ij)=pente_max*dyqmax(ij)
      ENDDO

!   calcul des pentes aux poles

      DO ij=1,iip1
         dyq(ij,l)=qpns-q(ij+iip1,l)
         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
      ENDDO

!   filtrage de la derivee
      dyn1=0.
      dys1=0.
      dyn2=0.
      dys2=0.
      DO ij=1,iim
         dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
         dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
         dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
         dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
      ENDDO
      DO ij=1,iip1
         dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
         dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
      ENDDO

!   calcul des pentes limites aux poles

      fn=1.
      fs=1.
      DO ij=1,iim
         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
         ENDIF
      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
         ENDIF
      ENDDO
      DO ij=1,iip1
         dyq(ij,l)=fn*dyq(ij,l)
         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
      ENDDO

!   calcul des pentes limitees

      DO ij=iip2,ip1jm
         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
         ELSE
            dyq(ij,l)=0.
         ENDIF
      ENDDO

      ENDDO

      DO l=1,llm
       DO ij=1,ip1jm
         IF( masse_adv_v(ij,l).GT.0. ) THEN
           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  + &
            dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
         ELSE
              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) * &
                         0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
         ENDIF
          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
       ENDDO
      ENDDO


      DO l=1,llm
         DO ij=iip2,ip1jm
            newmasse=masse(ij,l) &
            +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
               /newmasse
            masse(ij,l)=newmasse
         ENDDO
!.-. ancienne version
         convpn=SSUM(iim,qbyv(1,l),1)/apoln
         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
         DO ij = 1,iip1
            newmasse=masse(ij,l)+convmpn*aire(ij)
            q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/ &
                     newmasse
            masse(ij,l)=newmasse
         ENDDO
         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
         DO ij = ip1jm+1,ip1jmp1
            newmasse=masse(ij,l)+convmps*aire(ij)
            q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/ &
                     newmasse
            masse(ij,l)=newmasse
         ENDDO
!.-. fin ancienne version

!._. nouvelle version
!        convpn=SSUM(iim,qbyv(1,l),1)
!        convmpn=ssum(iim,masse_adv_v(1,l),1)
!        oldmasse=ssum(iim,masse(1,l),1)
!        newmasse=oldmasse+convmpn
!        newq=(q(1,l)*oldmasse+convpn)/newmasse
!        newmasse=newmasse/apoln
!        DO ij = 1,iip1
!           q(ij,l)=newq
!           masse(ij,l)=newmasse*aire(ij)
!        ENDDO
!        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
!        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
!        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
!        newmasse=oldmasse+convmps
!        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
!        newmasse=newmasse/apols
!        DO ij = ip1jm+1,ip1jmp1
!           q(ij,l)=newq
!           masse(ij,l)=newmasse*aire(ij)
!        ENDDO
!._. fin nouvelle version
      ENDDO

      RETURN
      END
1	guez	43	SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
2			!
3			! Auteurs: P.Le Van, F.Hourdin, F.Forget
4			!
5			! ********************************************************************
6			! Shema d'advection " pseudo amont " .
7			! ********************************************************************
8			! q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
9			! qsat est un argument de sortie pour le s-pg ....
10			!
11			!
12			! --------------------------------------------------------------------
13			!
14			use dimens_m
15			use paramet_m
16			use comconst
17	guez	66	use disvert_m
18	guez	57	use conf_gcm_m
19	guez	43	use comgeom
20			USE nr_util, ONLY : pi
21			IMPLICIT NONE
22			!
23			!
24			! Arguments:
25			! ----------
26			REAL masse(ip1jmp1,llm),pente_max
27			REAL masse_adv_v( ip1jm,llm)
28			REAL q(ip1jmp1,llm)
29			REAL qsat(ip1jmp1,llm)
30			!
31			! Local
32			! ---------
33			!
34			INTEGER i,ij,l
35			!
36			REAL airej2,airejjm,airescb(iim),airesch(iim)
37			REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
38			REAL adyqv(ip1jm),dyqmax(ip1jmp1)
39			REAL qbyv(ip1jm,llm)
40
41			REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
42			! REAL newq,oldmasse
43			Logical first,testcpu
44			REAL temps0,temps1,temps2,temps3,temps4,temps5
45			SAVE temps0,temps1,temps2,temps3,temps4,temps5
46			SAVE first,testcpu
47
48			REAL convpn,convps,convmpn,convmps
49			REAL sinlon(iip1),sinlondlon(iip1)
50			REAL coslon(iip1),coslondlon(iip1)
51			SAVE sinlon,coslon,sinlondlon,coslondlon
52			SAVE airej2,airejjm
53			!
54			!
55			REAL SSUM
56
57			DATA first,testcpu/.true.,.false./
58			DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
59
60			IF(first) THEN
61			PRINT*,'Shema Amont nouveau appele dans Vanleer '
62			first=.false.
63			do i=2,iip1
64			coslon(i)=cos(rlonv(i))
65			sinlon(i)=sin(rlonv(i))
66			coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
67			sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
68			ENDDO
69			coslon(1)=coslon(iip1)
70			coslondlon(1)=coslondlon(iip1)
71			sinlon(1)=sinlon(iip1)
72			sinlondlon(1)=sinlondlon(iip1)
73			airej2 = SSUM( iim, aire(iip2), 1 )
74			airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
75			ENDIF
76
77			!
78
79
80			DO l = 1, llm
81			!
82			! --------------------------------
83			! CALCUL EN LATITUDE
84			! --------------------------------
85
86			! On commence par calculer la valeur du traceur moyenne sur le premier cercle
87			! de latitude autour du pole (qpns pour le pole nord et qpsn pour
88			! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
89
90			DO i = 1, iim
91			airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
92			airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
93			ENDDO
94			qpns = SSUM( iim, airescb ,1 ) / airej2
95			qpsn = SSUM( iim, airesch ,1 ) / airejjm
96
97			! calcul des pentes aux points v
98
99			DO ij=1,ip1jm
100			dyqv(ij)=q(ij,l)-q(ij+iip1,l)
101			adyqv(ij)=abs(dyqv(ij))
102			ENDDO
103
104			! calcul des pentes aux points scalaires
105
106			DO ij=iip2,ip1jm
107			dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
108			dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
109			dyqmax(ij)=pente_max*dyqmax(ij)
110			ENDDO
111
112			! calcul des pentes aux poles
113
114			DO ij=1,iip1
115			dyq(ij,l)=qpns-q(ij+iip1,l)
116			dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
117			ENDDO
118
119			! filtrage de la derivee
120			dyn1=0.
121			dys1=0.
122			dyn2=0.
123			dys2=0.
124			DO ij=1,iim
125			dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
126			dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
127			dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
128			dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
129			ENDDO
130			DO ij=1,iip1
131			dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
132			dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
133			ENDDO
134
135			! calcul des pentes limites aux poles
136
137			fn=1.
138			fs=1.
139			DO ij=1,iim
140			IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
141			fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
142			ENDIF
143			IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
144			fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
145			ENDIF
146			ENDDO
147			DO ij=1,iip1
148			dyq(ij,l)=fn*dyq(ij,l)
149			dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
150			ENDDO
151
152			! calcul des pentes limitees
153
154			DO ij=iip2,ip1jm
155			IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
156			dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
157			ELSE
158			dyq(ij,l)=0.
159			ENDIF
160			ENDDO
161
162			ENDDO
163
164			DO l=1,llm
165			DO ij=1,ip1jm
166			IF( masse_adv_v(ij,l).GT.0. ) THEN
167			qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) + &
168			dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
169			ELSE
170			qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) * &
171			0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
172			ENDIF
173			qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
174			ENDDO
175			ENDDO
176
177
178			DO l=1,llm
179			DO ij=iip2,ip1jm
180			newmasse=masse(ij,l) &
181			+masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
182			q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
183			/newmasse
184			masse(ij,l)=newmasse
185			ENDDO
186			!.-. ancienne version
187			convpn=SSUM(iim,qbyv(1,l),1)/apoln
188			convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
189			DO ij = 1,iip1
190			newmasse=masse(ij,l)+convmpn*aire(ij)
191			q(ij,l)=(q(ij,l)masse(ij,l)+convpnaire(ij))/ &
192			newmasse
193			masse(ij,l)=newmasse
194			ENDDO
195			convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
196			convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
197			DO ij = ip1jm+1,ip1jmp1
198			newmasse=masse(ij,l)+convmps*aire(ij)
199			q(ij,l)=(q(ij,l)masse(ij,l)+convpsaire(ij))/ &
200			newmasse
201			masse(ij,l)=newmasse
202			ENDDO
203			!.-. fin ancienne version
204
205			!._. nouvelle version
206			! convpn=SSUM(iim,qbyv(1,l),1)
207			! convmpn=ssum(iim,masse_adv_v(1,l),1)
208			! oldmasse=ssum(iim,masse(1,l),1)
209			! newmasse=oldmasse+convmpn
210			! newq=(q(1,l)*oldmasse+convpn)/newmasse
211			! newmasse=newmasse/apoln
212			! DO ij = 1,iip1
213			! q(ij,l)=newq
214			! masse(ij,l)=newmasse*aire(ij)
215			! ENDDO
216			! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
217			! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
218			! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
219			! newmasse=oldmasse+convmps
220			! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
221			! newmasse=newmasse/apols
222			! DO ij = ip1jm+1,ip1jmp1
223			! q(ij,l)=newq
224			! masse(ij,l)=newmasse*aire(ij)
225			! ENDDO
226			!._. fin nouvelle version
227			ENDDO
228
229			RETURN
230			END