/[lmdze]/trunk/dyn3d/ADVN/advnx.f90
ViewVC logotype

Annotation of /trunk/dyn3d/ADVN/advnx.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 350 - (hide annotations)
Mon Dec 23 16:07:24 2019 UTC (4 years, 5 months ago) by guez
File size: 6879 byte(s) 
Add argument `itau_phy_redem` to phyredem0

Add argument `itau_phy_redem` to phyredem0. Motivation: being able to
call phyredem0 with `itau_phy_redem` = 0 in program ce0l, avoiding the
side effect on nday. We can now add attribute protected to variable
nday of module `conf_gcm`. And we do not need to set `itau_phy` to 0
in procedure `phyetat0_new`.

Remove variable `prt_level` of module `conf_gcm`, which was only used
in procedure advnx. Motivation: cluttering the output is not viable,
even for debug. Better use a debugger.

Forgot to add `dyn3d/ADVN/CMakeLists.txt` in previous revision.
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))*(zsigp*zqp+(zsig-zsigp)*zqm)
90     END IF
91     ELSE
92     IF (zsig<=zsigp) THEN
93     u_mq(ij, l) = u_m(ij, l)*(zqp-0.5*zsig/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.5*zsig/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.5*zsig/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
  ViewVC Help
Powered by ViewVC 1.1.21