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 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), dq( ip1jmp1,llm)
!
!      Local
!   ---------
!
      INTEGER i,ij,l
!
      REAL airej2,airejjm,airescb(iim),airesch(iim)
      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
      REAL qbyv(ip1jm,llm)

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

      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,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

      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			!
11			!
12			! ----------------------------------------------------------------
13			use dimens_m
14			use paramet_m
15			use comconst
16	guez	66	use disvert_m
17	guez	57	use conf_gcm_m
18	guez	31	use comgeom
19	guez	39	USE nr_util, ONLY : pi
20	guez	31	IMPLICIT NONE
21			!
22			!
23			!
24			! Arguments:
25			! ----------
26			REAL masse(ip1jmp1,llm),pente_max
27			REAL masse_adv_v( ip1jm,llm)
28			REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
29			!
30			! Local
31			! ---------
32			!
33			INTEGER i,ij,l
34			!
35			REAL airej2,airejjm,airescb(iim),airesch(iim)
36			REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
37			REAL adyqv(ip1jm),dyqmax(ip1jmp1)
38			REAL qbyv(ip1jm,llm)
39
40			REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
41			! REAL newq,oldmasse
42			Logical extremum,first,testcpu
43			REAL temps0,temps1,temps2,temps3,temps4,temps5,second
44			SAVE temps0,temps1,temps2,temps3,temps4,temps5
45			SAVE first,testcpu
46
47			REAL convpn,convps,convmpn,convmps
48			real massepn,masseps,qpn,qps
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			DO l = 1, llm
79			!
80			! --------------------------------
81			! CALCUL EN LATITUDE
82			! --------------------------------
83
84			! On commence par calculer la valeur du traceur moyenne sur le
85			! premier cercle
86			! de latitude autour du pole (qpns pour le pole nord et qpsn pour
87			! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
88
89			DO i = 1, iim
90			airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
91			airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
92			ENDDO
93			qpns = SSUM( iim, airescb ,1 ) / airej2
94			qpsn = SSUM( iim, airesch ,1 ) / airejjm
95
96			! calcul des pentes aux points v
97
98			DO ij=1,ip1jm
99			dyqv(ij)=q(ij,l)-q(ij+iip1,l)
100			adyqv(ij)=abs(dyqv(ij))
101			ENDDO
102
103			! calcul des pentes aux points scalaires
104
105			DO ij=iip2,ip1jm
106			dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
107			dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
108			dyqmax(ij)=pente_max*dyqmax(ij)
109			ENDDO
110
111			! calcul des pentes aux poles
112
113			DO ij=1,iip1
114			dyq(ij,l)=qpns-q(ij+iip1,l)
115			dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
116			ENDDO
117
118			! filtrage de la derivee
119			dyn1=0.
120			dys1=0.
121			dyn2=0.
122			dys2=0.
123			DO ij=1,iim
124			dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
125			dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
126			dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
127			dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
128			ENDDO
129			DO ij=1,iip1
130			dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
131			dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
132			ENDDO
133
134			! calcul des pentes limites aux poles
135
136			goto 8888
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			8888 continue
152			DO ij=1,iip1
153			dyq(ij,l)=0.
154			dyq(ip1jm+ij,l)=0.
155			ENDDO
156
157			! calcul des pentes limitees
158
159			DO ij=iip2,ip1jm
160			IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
161			dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
162			ELSE
163			dyq(ij,l)=0.
164			ENDIF
165			ENDDO
166
167			ENDDO
168
169			DO l=1,llm
170			DO ij=1,ip1jm
171			IF(masse_adv_v(ij,l).gt.0) THEN
172			qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* &
173			& 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
174			ELSE
175			qbyv(ij,l)=q(ij,l)-dyq(ij,l)* &
176			& 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
177			ENDIF
178			qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
179			ENDDO
180			ENDDO
181
182
183			DO l=1,llm
184			DO ij=iip2,ip1jm
185			newmasse=masse(ij,l) &
186			& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
187			q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
188			& /newmasse
189			masse(ij,l)=newmasse
190			ENDDO
191			!.-. ancienne version
192			! convpn=SSUM(iim,qbyv(1,l),1)/apoln
193			! convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
194
195			convpn=SSUM(iim,qbyv(1,l),1)
196			convmpn=ssum(iim,masse_adv_v(1,l),1)
197			massepn=ssum(iim,masse(1,l),1)
198			qpn=0.
199			do ij=1,iim
200			qpn=qpn+masse(ij,l)*q(ij,l)
201			enddo
202			qpn=(qpn+convpn)/(massepn+convmpn)
203			do ij=1,iip1
204			q(ij,l)=qpn
205			enddo
206
207			! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
208			! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
209
210			convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
211			convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
212			masseps=ssum(iim, masse(ip1jm+1,l),1)
213			qps=0.
214			do ij = ip1jm+1,ip1jmp1-1
215			qps=qps+masse(ij,l)*q(ij,l)
216			enddo
217			qps=(qps+convps)/(masseps+convmps)
218			do ij=ip1jm+1,ip1jmp1
219			q(ij,l)=qps
220			enddo
221
222			!.-. fin ancienne version
223
224			!._. nouvelle version
225			! convpn=SSUM(iim,qbyv(1,l),1)
226			! convmpn=ssum(iim,masse_adv_v(1,l),1)
227			! oldmasse=ssum(iim,masse(1,l),1)
228			! newmasse=oldmasse+convmpn
229			! newq=(q(1,l)*oldmasse+convpn)/newmasse
230			! newmasse=newmasse/apoln
231			! DO ij = 1,iip1
232			! q(ij,l)=newq
233			! masse(ij,l)=newmasse*aire(ij)
234			! ENDDO
235			! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
236			! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
237			! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
238			! newmasse=oldmasse+convmps
239			! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
240			! newmasse=newmasse/apols
241			! DO ij = ip1jm+1,ip1jmp1
242			! q(ij,l)=newq
243			! masse(ij,l)=newmasse*aire(ij)
244			! ENDDO
245			!._. fin nouvelle version
246			ENDDO
247
248			RETURN
249			END