dyn3d/ADVN/advnx.f90

SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)

  ! Auteur : F. Hourdin

  ! ********************************************************************
  ! Shema  d'advection " pseudo amont " .
  ! ********************************************************************
  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....


  ! --------------------------------------------------------------------
  USE dimensions
  USE paramet_m
  USE comconst
  USE disvert_m
  USE conf_gcm_m
  IMPLICIT NONE


  ! Arguments:
  ! ----------
  INTEGER mode
  REAL masse(ip1jmp1, llm)
  REAL u_m(ip1jmp1, llm)
  REAL q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)

  ! Local
  ! ---------

  INTEGER i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
  INTEGER n0, nl(llm)

  REAL new_m, zu_m, zdq, zz
  REAL zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
  REAL u_mq(ip1jmp1, llm)

  REAL zm, zq, zsigm, zsigp, zqm, zqp, zu

  LOGICAL ladvplus(ip1jmp1, llm)

  REAL prec
  SAVE prec

  DATA prec/1.E-15/

  DO l = 1, llm
    DO ij = iip2, ip1jm
      zdq = qd(ij, l) - qg(ij, l)
      IF (abs(zdq)>prec) THEN
        zsigd(ij, l) = (q(ij,l)-qg(ij,l))/zdq
        zsigg(ij, l) = 1. - zsigd(ij, l)
      ELSE
        zsigd(ij, l) = 0.5
        zsigg(ij, l) = 0.5
        qd(ij, l) = q(ij, l)
        qg(ij, l) = q(ij, l)
      END IF
    END DO
  END DO

  ! calcul de la pente maximum dans la maille en valeur absolue

  DO l = 1, llm
    DO ij = iip2, ip1jm - 1
      IF (u_m(ij,l)>=0.) THEN
        zsigp = zsigd(ij, l)
        zsigm = zsigg(ij, l)
        zqp = qd(ij, l)
        zqm = qg(ij, l)
        zm = masse(ij, l)
        zq = q(ij, l)
      ELSE
        zsigm = zsigd(ij+1, l)
        zsigp = zsigg(ij+1, l)
        zqm = qd(ij+1, l)
        zqp = qg(ij+1, l)
        zm = masse(ij+1, l)
        zq = q(ij+1, l)
      END IF
      zu = abs(u_m(ij,l))
      ladvplus(ij, l) = zu > zm
      zsig = zu/zm
      IF (zsig==0.) zsigp = 0.1
      IF (mode==1) THEN
        IF (zsig<=zsigp) THEN
          u_mq(ij, l) = u_m(ij, l)*zqp
        ELSE IF (mode==1) THEN
          u_mq(ij, l) = sign(zm, u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
        END IF
      ELSE
        IF (zsig<=zsigp) THEN
          u_mq(ij, l) = u_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
        ELSE
          zz = 0.5*(zsig-zsigp)/zsigm
          u_mq(ij, l) = sign(zm, u_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)* &
            (zq+zz*(zqm-zq)))
        END IF
      END IF
    END DO
  END DO

  DO l = 1, llm
    DO ij = iip1 + iip1, ip1jm, iip1
      u_mq(ij, l) = u_mq(ij-iim, l)
      ladvplus(ij, l) = ladvplus(ij-iim, l)
    END DO
  END DO

  ! =================================================================
  ! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
  ! =================================================================
  ! tris des regions a traiter
  n0 = 0
  DO l = 1, llm
    nl(l) = 0
    DO ij = iip2, ip1jm
      IF (ladvplus(ij,l)) THEN
        nl(l) = nl(l) + 1
        u_mq(ij, l) = 0.
      END IF
    END DO
    n0 = n0 + nl(l)
  END DO

  IF (n0>1) THEN
    IF (prt_level>9) PRINT *, &
      'Nombre de points pour lesquels on advect plus que le', &
      'contenu de la maille : ', n0

    DO l = 1, llm
      IF (nl(l)>0) THEN
        iju = 0
        ! indicage des mailles concernees par le traitement special
        DO ij = iip2, ip1jm
          IF (ladvplus(ij,l) .AND. mod(ij,iip1)/=0) THEN
            iju = iju + 1
            indu(iju) = ij
          END IF
        END DO
        niju = iju

        ! traitement des mailles
        DO iju = 1, niju
          ij = indu(iju)
          j = (ij-1)/iip1 + 1
          zu_m = u_m(ij, l)
          u_mq(ij, l) = 0.
          IF (zu_m>0.) THEN
            ijq = ij
            i = ijq - (j-1)*iip1
            ! accumulation pour les mailles completements advectees
            DO WHILE (zu_m>masse(ijq,l))
              u_mq(ij, l) = u_mq(ij, l) + q(ijq, l)*masse(ijq, l)
              zu_m = zu_m - masse(ijq, l)
              i = mod(i-2+iim, iim) + 1
              ijq = (j-1)*iip1 + i
            END DO
            ! MODIFS SPECIFIQUES DU SCHEMA
            ! ajout de la maille non completement advectee
            zsig = zu_m/masse(ijq, l)
            IF (zsig<=zsigd(ijq,l)) THEN
              u_mq(ij, l) = u_mq(ij, l) + zu_m*(qd(ijq,l)-0.5*zsig/zsigd(ijq, &
                l)*(qd(ijq,l)-q(ijq,l)))
            ELSE
              ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
              ! goto 8888
              zz = 0.5*(zsig-zsigd(ijq,l))/zsigg(ijq, l)
              IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
                PRINT *, 'probleme2 au point ij=', ij, '  l=', l
                PRINT *, 'zz=', zz
                STOP
              END IF
              u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l)*(0.5*(q(ijq, &
                l)+qd(ijq,l))*zsigd(ijq,l)+(zsig-zsigd(ijq,l))*(q(ijq, &
                l)+zz*(qg(ijq,l)-q(ijq,l))))
            END IF
          ELSE
            ijq = ij + 1
            i = ijq - (j-1)*iip1
            ! accumulation pour les mailles completements advectees
            DO WHILE (-zu_m>masse(ijq,l))
              u_mq(ij, l) = u_mq(ij, l) - q(ijq, l)*masse(ijq, l)
              zu_m = zu_m + masse(ijq, l)
              i = mod(i, iim) + 1
              ijq = (j-1)*iip1 + i
            END DO
            ! ajout de la maille non completement advectee
            ! 2eme MODIF SPECIFIQUE
            zsig = -zu_m/masse(ij+1, l)
            IF (zsig<=zsigg(ijq,l)) THEN
              u_mq(ij, l) = u_mq(ij, l) + zu_m*(qg(ijq,l)-0.5*zsig/zsigg(ijq, &
                l)*(qg(ijq,l)-q(ijq,l)))
            ELSE
              ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
              ! goto 9999
              zz = 0.5*(zsig-zsigg(ijq,l))/zsigd(ijq, l)
              IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
                PRINT *, 'probleme22 au point ij=', ij, '  l=', l
                PRINT *, 'zz=', zz
                STOP
              END IF
              u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l)*(0.5*(q(ijq, &
                l)+qg(ijq,l))*zsigg(ijq,l)+(zsig-zsigg(ijq,l))*(q(ijq, &
                l)+zz*(qd(ijq,l)-q(ijq,l))))
            END IF
            ! fin de la modif
          END IF
        END DO
      END IF
    END DO
  END IF ! n0.gt.0

  ! bouclage en latitude
  DO l = 1, llm
    DO ij = iip1 + iip1, ip1jm, iip1
      u_mq(ij, l) = u_mq(ij-iim, l)
    END DO
  END DO

  ! =================================================================
  ! CALCUL DE LA CONVERGENCE DES FLUX
  ! =================================================================

  DO l = 1, llm
    DO ij = iip2 + 1, ip1jm
      new_m = masse(ij, l) + u_m(ij-1, l) - u_m(ij, l)
      q(ij, l) = (q(ij,l)*masse(ij,l)+u_mq(ij-1,l)-u_mq(ij,l))/new_m
      masse(ij, l) = new_m
    END DO
    ! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    DO ij = iip1 + iip1, ip1jm, iip1
      q(ij-iim, l) = q(ij, l)
      masse(ij-iim, l) = masse(ij, l)
    END DO
  END DO

  RETURN
END SUBROUTINE advnx
1	guez	349	SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)
2
3			! Auteur : F. Hourdin
4
5			! ********************************************************************
6			! Shema d'advection " pseudo amont " .
7			! ********************************************************************
8			! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
9
10
11			! --------------------------------------------------------------------
12			USE dimensions
13			USE paramet_m
14			USE comconst
15			USE disvert_m
16			USE conf_gcm_m
17			IMPLICIT NONE
18
19
20
21			! Arguments:
22			! ----------
23			INTEGER mode
24			REAL masse(ip1jmp1, llm)
25			REAL u_m(ip1jmp1, llm)
26			REAL q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)
27
28			! Local
29			! ---------
30
31			INTEGER i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
32			INTEGER n0, nl(llm)
33
34			REAL new_m, zu_m, zdq, zz
35			REAL zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
36			REAL u_mq(ip1jmp1, llm)
37
38			REAL zm, zq, zsigm, zsigp, zqm, zqp, zu
39
40			LOGICAL ladvplus(ip1jmp1, llm)
41
42			REAL prec
43			SAVE prec
44
45			DATA prec/1.E-15/
46
47			DO l = 1, llm
48			DO ij = iip2, ip1jm
49			zdq = qd(ij, l) - qg(ij, l)
50			IF (abs(zdq)>prec) THEN
51			zsigd(ij, l) = (q(ij,l)-qg(ij,l))/zdq
52			zsigg(ij, l) = 1. - zsigd(ij, l)
53			ELSE
54			zsigd(ij, l) = 0.5
55			zsigg(ij, l) = 0.5
56			qd(ij, l) = q(ij, l)
57			qg(ij, l) = q(ij, l)
58			END IF
59			END DO
60			END DO
61
62			! calcul de la pente maximum dans la maille en valeur absolue
63
64			DO l = 1, llm
65			DO ij = iip2, ip1jm - 1
66			IF (u_m(ij,l)>=0.) THEN
67			zsigp = zsigd(ij, l)
68			zsigm = zsigg(ij, l)
69			zqp = qd(ij, l)
70			zqm = qg(ij, l)
71			zm = masse(ij, l)
72			zq = q(ij, l)
73			ELSE
74			zsigm = zsigd(ij+1, l)
75			zsigp = zsigg(ij+1, l)
76			zqm = qd(ij+1, l)
77			zqp = qg(ij+1, l)
78			zm = masse(ij+1, l)
79			zq = q(ij+1, l)
80			END IF
81			zu = abs(u_m(ij,l))
82			ladvplus(ij, l) = zu > zm
83			zsig = zu/zm
84			IF (zsig==0.) zsigp = 0.1
85			IF (mode==1) THEN
86			IF (zsig<=zsigp) THEN
87			u_mq(ij, l) = u_m(ij, l)*zqp
88			ELSE IF (mode==1) THEN
89			u_mq(ij, l) = sign(zm, u_m(ij,l))(zsigpzqp+(zsig-zsigp)*zqm)
90			END IF
91			ELSE
92			IF (zsig<=zsigp) THEN
93			u_mq(ij, l) = u_m(ij, l)(zqp-0.5zsig/zsigp*(zqp-zq))
94			ELSE
95			zz = 0.5*(zsig-zsigp)/zsigm
96			u_mq(ij, l) = sign(zm, u_m(ij,l))(0.5(zq+zqp)zsigp+(zsig-zsigp) &
97			(zq+zz*(zqm-zq)))
98			END IF
99			END IF
100			END DO
101			END DO
102
103			DO l = 1, llm
104			DO ij = iip1 + iip1, ip1jm, iip1
105			u_mq(ij, l) = u_mq(ij-iim, l)
106			ladvplus(ij, l) = ladvplus(ij-iim, l)
107			END DO
108			END DO
109
110			! =================================================================
111			! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
112			! =================================================================
113			! tris des regions a traiter
114			n0 = 0
115			DO l = 1, llm
116			nl(l) = 0
117			DO ij = iip2, ip1jm
118			IF (ladvplus(ij,l)) THEN
119			nl(l) = nl(l) + 1
120			u_mq(ij, l) = 0.
121			END IF
122			END DO
123			n0 = n0 + nl(l)
124			END DO
125
126			IF (n0>1) THEN
127			IF (prt_level>9) PRINT *, &
128			'Nombre de points pour lesquels on advect plus que le', &
129			'contenu de la maille : ', n0
130
131			DO l = 1, llm
132			IF (nl(l)>0) THEN
133			iju = 0
134			! indicage des mailles concernees par le traitement special
135			DO ij = iip2, ip1jm
136			IF (ladvplus(ij,l) .AND. mod(ij,iip1)/=0) THEN
137			iju = iju + 1
138			indu(iju) = ij
139			END IF
140			END DO
141			niju = iju
142
143			! traitement des mailles
144			DO iju = 1, niju
145			ij = indu(iju)
146			j = (ij-1)/iip1 + 1
147			zu_m = u_m(ij, l)
148			u_mq(ij, l) = 0.
149			IF (zu_m>0.) THEN
150			ijq = ij
151			i = ijq - (j-1)*iip1
152			! accumulation pour les mailles completements advectees
153			DO WHILE (zu_m>masse(ijq,l))
154			u_mq(ij, l) = u_mq(ij, l) + q(ijq, l)*masse(ijq, l)
155			zu_m = zu_m - masse(ijq, l)
156			i = mod(i-2+iim, iim) + 1
157			ijq = (j-1)*iip1 + i
158			END DO
159			! MODIFS SPECIFIQUES DU SCHEMA
160			! ajout de la maille non completement advectee
161			zsig = zu_m/masse(ijq, l)
162			IF (zsig<=zsigd(ijq,l)) THEN
163			u_mq(ij, l) = u_mq(ij, l) + zu_m(qd(ijq,l)-0.5zsig/zsigd(ijq, &
164			l)*(qd(ijq,l)-q(ijq,l)))
165			ELSE
166			! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
167			! goto 8888
168			zz = 0.5*(zsig-zsigd(ijq,l))/zsigg(ijq, l)
169			IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
170			PRINT *, 'probleme2 au point ij=', ij, ' l=', l
171			PRINT *, 'zz=', zz
172			STOP
173			END IF
174			u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l)(0.5(q(ijq, &
175			l)+qd(ijq,l))zsigd(ijq,l)+(zsig-zsigd(ijq,l))(q(ijq, &
176			l)+zz*(qg(ijq,l)-q(ijq,l))))
177			END IF
178			ELSE
179			ijq = ij + 1
180			i = ijq - (j-1)*iip1
181			! accumulation pour les mailles completements advectees
182			DO WHILE (-zu_m>masse(ijq,l))
183			u_mq(ij, l) = u_mq(ij, l) - q(ijq, l)*masse(ijq, l)
184			zu_m = zu_m + masse(ijq, l)
185			i = mod(i, iim) + 1
186			ijq = (j-1)*iip1 + i
187			END DO
188			! ajout de la maille non completement advectee
189			! 2eme MODIF SPECIFIQUE
190			zsig = -zu_m/masse(ij+1, l)
191			IF (zsig<=zsigg(ijq,l)) THEN
192			u_mq(ij, l) = u_mq(ij, l) + zu_m(qg(ijq,l)-0.5zsig/zsigg(ijq, &
193			l)*(qg(ijq,l)-q(ijq,l)))
194			ELSE
195			! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
196			! goto 9999
197			zz = 0.5*(zsig-zsigg(ijq,l))/zsigd(ijq, l)
198			IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
199			PRINT *, 'probleme22 au point ij=', ij, ' l=', l
200			PRINT *, 'zz=', zz
201			STOP
202			END IF
203			u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l)(0.5(q(ijq, &
204			l)+qg(ijq,l))zsigg(ijq,l)+(zsig-zsigg(ijq,l))(q(ijq, &
205			l)+zz*(qd(ijq,l)-q(ijq,l))))
206			END IF
207			! fin de la modif
208			END IF
209			END DO
210			END IF
211			END DO
212			END IF ! n0.gt.0
213
214			! bouclage en latitude
215			DO l = 1, llm
216			DO ij = iip1 + iip1, ip1jm, iip1
217			u_mq(ij, l) = u_mq(ij-iim, l)
218			END DO
219			END DO
220
221			! =================================================================
222			! CALCUL DE LA CONVERGENCE DES FLUX
223			! =================================================================
224
225			DO l = 1, llm
226			DO ij = iip2 + 1, ip1jm
227			new_m = masse(ij, l) + u_m(ij-1, l) - u_m(ij, l)
228			q(ij, l) = (q(ij,l)*masse(ij,l)+u_mq(ij-1,l)-u_mq(ij,l))/new_m
229			masse(ij, l) = new_m
230			END DO
231			! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
232			DO ij = iip1 + iip1, ip1jm, iip1
233			q(ij-iim, l) = q(ij, l)
234			masse(ij-iim, l) = masse(ij, l)
235			END DO
236			END DO
237
238			RETURN
239			END SUBROUTINE advnx