dyn3d/Vlsplt/vly.f

      SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
!
!     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 ....
!     dq              sont des arguments de sortie pour le s-pg ....

      use comconst
      use comgeom, only: aire
      use conf_gcm_m
      use dimensions
      use disvert_m
      USE dynetat0_m, only: rlonv, rlonu
      USE nr_util, ONLY : pi
      use paramet_m

      IMPLICIT NONE
!
!
!
!   Arguments:
!   ----------
      REAL masse(ip1jmp1,llm)
      real, intent(in):: pente_max
      REAL masse_adv_v( ip1jm,llm)
      REAL q(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
      SAVE first

      REAL convpn,convps,convmpn,convmps
      real massepn,masseps,qpn,qps
      REAL sinlon(iip1),sinlondlon(iip1)
      REAL coslon(iip1),coslondlon(iip1)
      SAVE sinlon,coslon,sinlondlon,coslondlon
      SAVE airej2,airejjm
!
!
      REAL      SSUM

      DATA first/.true./

      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

      goto 8888
      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
 8888   continue
      DO ij=1,iip1
         dyq(ij,l)=0.
         dyq(ip1jm+ij,l)=0.
      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)=q(ij+iip1,l)+dyq(ij+iip1,l)*                   &
     &                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
          ELSE
              qbyv(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

         convpn=SSUM(iim,qbyv(1,l),1)
         convmpn=ssum(iim,masse_adv_v(1,l),1)
         massepn=ssum(iim,masse(1,l),1)
         qpn=0.
         do ij=1,iim
            qpn=qpn+masse(ij,l)*q(ij,l)
         enddo
         qpn=(qpn+convpn)/(massepn+convmpn)
         do ij=1,iip1
            q(ij,l)=qpn
         enddo

!        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
!        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols

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