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
    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	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	DO l = 1, llm
128	IF (nl(l)>0) THEN
129	iju = 0
130	! indicage des mailles concernees par le traitement special
131	DO ij = iip2, ip1jm
132	IF (ladvplus(ij,l) .AND. mod(ij,iip1)/=0) THEN
133	iju = iju + 1
134	indu(iju) = ij
135	END IF
136	END DO
137	niju = iju
138
139	! traitement des mailles
140	DO iju = 1, niju
141	ij = indu(iju)
142	j = (ij-1)/iip1 + 1
143	zu_m = u_m(ij, l)
144	u_mq(ij, l) = 0.
145	IF (zu_m>0.) THEN
146	ijq = ij
147	i = ijq - (j-1)*iip1
148	! accumulation pour les mailles completements advectees
149	DO WHILE (zu_m>masse(ijq,l))
150	u_mq(ij, l) = u_mq(ij, l) + q(ijq, l)*masse(ijq, l)
151	zu_m = zu_m - masse(ijq, l)
152	i = mod(i-2+iim, iim) + 1
153	ijq = (j-1)*iip1 + i
154	END DO
155	! MODIFS SPECIFIQUES DU SCHEMA
156	! ajout de la maille non completement advectee
157	zsig = zu_m/masse(ijq, l)
158	IF (zsig<=zsigd(ijq,l)) THEN
159	u_mq(ij, l) = u_mq(ij, l) + zu_m(qd(ijq,l)-0.5zsig/zsigd(ijq, &
160	l)*(qd(ijq,l)-q(ijq,l)))
161	ELSE
162	! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
163	! goto 8888
164	zz = 0.5*(zsig-zsigd(ijq,l))/zsigg(ijq, l)
165	IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
166	PRINT *, 'probleme2 au point ij=', ij, ' l=', l
167	PRINT *, 'zz=', zz
168	STOP
169	END IF
170	u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l)(0.5(q(ijq, &
171	l)+qd(ijq,l))zsigd(ijq,l)+(zsig-zsigd(ijq,l))(q(ijq, &
172	l)+zz*(qg(ijq,l)-q(ijq,l))))
173	END IF
174	ELSE
175	ijq = ij + 1
176	i = ijq - (j-1)*iip1
177	! accumulation pour les mailles completements advectees
178	DO WHILE (-zu_m>masse(ijq,l))
179	u_mq(ij, l) = u_mq(ij, l) - q(ijq, l)*masse(ijq, l)
180	zu_m = zu_m + masse(ijq, l)
181	i = mod(i, iim) + 1
182	ijq = (j-1)*iip1 + i
183	END DO
184	! ajout de la maille non completement advectee
185	! 2eme MODIF SPECIFIQUE
186	zsig = -zu_m/masse(ij+1, l)
187	IF (zsig<=zsigg(ijq,l)) THEN
188	u_mq(ij, l) = u_mq(ij, l) + zu_m(qg(ijq,l)-0.5zsig/zsigg(ijq, &
189	l)*(qg(ijq,l)-q(ijq,l)))
190	ELSE
191	! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
192	! goto 9999
193	zz = 0.5*(zsig-zsigg(ijq,l))/zsigd(ijq, l)
194	IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
195	PRINT *, 'probleme22 au point ij=', ij, ' l=', l
196	PRINT *, 'zz=', zz
197	STOP
198	END IF
199	u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l)(0.5(q(ijq, &
200	l)+qg(ijq,l))zsigg(ijq,l)+(zsig-zsigg(ijq,l))(q(ijq, &
201	l)+zz*(qd(ijq,l)-q(ijq,l))))
202	END IF
203	! fin de la modif
204	END IF
205	END DO
206	END IF
207	END DO
208	END IF ! n0.gt.0
209
210	! bouclage en latitude
211	DO l = 1, llm
212	DO ij = iip1 + iip1, ip1jm, iip1
213	u_mq(ij, l) = u_mq(ij-iim, l)
214	END DO
215	END DO
216
217	! =================================================================
218	! CALCUL DE LA CONVERGENCE DES FLUX
219	! =================================================================
220
221	DO l = 1, llm
222	DO ij = iip2 + 1, ip1jm
223	new_m = masse(ij, l) + u_m(ij-1, l) - u_m(ij, l)
224	q(ij, l) = (q(ij,l)*masse(ij,l)+u_mq(ij-1,l)-u_mq(ij,l))/new_m
225	masse(ij, l) = new_m
226	END DO
227	! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
228	DO ij = iip1 + iip1, ip1jm, iip1
229	q(ij-iim, l) = q(ij, l)
230	masse(ij-iim, l) = masse(ij, l)
231	END DO
232	END DO
233
234	RETURN
235	END SUBROUTINE advnx