/[lmdze]/trunk/dyn3d/advtrac.f
ViewVC logotype

Contents of /trunk/dyn3d/advtrac.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 91 - (show annotations)
Wed Mar 26 17:18:58 2014 UTC (10 years, 2 months ago) by guez
File size: 8131 byte(s) 
Removed unused variables lock_startdate and time_stamp of module
calendar.

Noticed that physiq does not change the surface pressure. So removed
arguments ps and dpfi of subroutine addfi. dpfi was always 0. The
computation of ps in addfi included some averaging at the poles. In
principle, this does not change ps but in practice it does because of
finite numerical precision. So the results of the simulation are
changed. Removed arguments ps and dpfi of calfis. Removed argument
d_ps of physiq.

du at the poles is not computed by dudv1, so declare only the
corresponding latitudes in dudv1. caldyn passes only a section of the
array dudyn as argument.

Removed variable niadv of module iniadvtrac_m.

Declared arguments of exner_hyb as assumed-shape arrays and made all
other horizontal sizes in exner_hyb dynamic. This allows the external
program test_disvert to use exner_hyb at a single horizontal position.
1 module advtrac_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, pk)
8
9 ! From dyn3d/advtrac.F, version 1.4 2005/04/13 08:58:34
10 ! Author: F. Hourdin
11
12 USE comconst, ONLY : dtvr
13 USE conf_gcm_m, ONLY : iapp_tracvl
14 USE dimens_m, ONLY : iim, jjm, llm, nqmx
15 USE iniadvtrac_m, ONLY : iadv
16 use massbar_m, only: massbar
17 USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
18 llmp1
19
20 REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
21 REAL, intent(in):: p(ip1jmp1, llmp1)
22 real, intent(in):: masse(ip1jmp1, llm)
23 REAL, intent(inout):: q(ip1jmp1, llm, nqmx)
24 INTEGER, intent(out):: iapptrac
25 real, intent(in):: teta(ip1jmp1, llm)
26 REAL, intent(in):: pk(ip1jmp1, llm)
27
28 ! Variables locales
29
30 REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)
31 REAL, save:: pbaruc(ip1jmp1, llm), pbarvc(ip1jm, llm)
32 REAL, save:: massem(ip1jmp1, llm)
33 real zdp(ip1jmp1)
34 REAL pbarug(ip1jmp1, llm), pbarvg(ip1jm, llm), wg(ip1jmp1, llm)
35 REAL cpuadv(nqmx)
36
37 INTEGER:: iadvtr = 0
38 INTEGER ij, l, iq
39 REAL zdpmin, zdpmax
40 EXTERNAL minmax
41
42 ! Rajouts pour PPM
43
44 INTEGER indice, n
45 ! Pas de temps adaptatif pour que CFL < 1
46 REAL dtbon
47 REAL cflmaxz ! CFL maximum
48 real aaa, bbb
49 REAL psppm(iim, jjp1) ! pression au sol
50 REAL unatppm(iim, jjp1, llm), vnatppm(iim, jjp1, llm)
51 REAL qppm(iim*jjp1, llm, nqmx)
52 REAL fluxwppm(iim, jjp1, llm)
53 REAL apppm(llmp1), bpppm(llmp1)
54 LOGICAL:: dum = .TRUE., fill = .TRUE.
55
56 !-----------------------------------------------------------
57
58 IF (iadvtr==0) THEN
59 CALL initial0(ijp1llm, pbaruc)
60 CALL initial0(ijmllm, pbarvc)
61 END IF
62
63 ! accumulation des flux de masse horizontaux
64 DO l = 1, llm
65 DO ij = 1, ip1jmp1
66 pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
67 END DO
68 DO ij = 1, ip1jm
69 pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
70 END DO
71 END DO
72
73 ! selection de la masse instantannee des mailles avant le transport.
74 IF (iadvtr==0) massem = masse
75
76 iadvtr = iadvtr + 1
77 iapptrac = iadvtr
78
79 ! Test pour savoir si on advecte a ce pas de temps
80 IF (iadvtr == iapp_tracvl) THEN
81 ! traitement des flux de masse avant advection.
82 ! 1. calcul de w
83 ! 2. groupement des mailles pres du pole.
84
85 CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg)
86
87 ! test sur l'eventuelle creation de valeurs negatives de la masse
88 DO l = 1, llm - 1
89 DO ij = iip2 + 1, ip1jm
90 zdp(ij) = pbarug(ij-1, l) - pbarug(ij, l) - pbarvg(ij-iip1, l) + &
91 pbarvg(ij, l) + wg(ij, l+1) - wg(ij, l)
92 END DO
93 CALL scopy(jjm-1, zdp(iip1+iip1), iip1, zdp(iip2), iip1)
94 DO ij = iip2, ip1jm
95 zdp(ij) = zdp(ij)*dtvr/massem(ij, l)
96 END DO
97
98 CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)
99
100 IF (max(abs(zdpmin), abs(zdpmax))>0.5) THEN
101 PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
102 END IF
103 END DO
104
105 ! Advection proprement dite
106
107 ! Calcul des moyennes basées sur la masse
108 CALL massbar(massem, massebx, masseby)
109
110 ! Appel des sous programmes d'advection
111
112 DO iq = 1, nqmx
113 IF (iadv(iq)==0) CYCLE
114
115 ! Schema de Van Leer I MUSCL
116
117 IF (iadv(iq)==10) THEN
118 CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
119 ! Schema "pseudo amont" + test sur humidite specifique
120 ! pour la vapeur d'eau. F. Codron
121 ELSE IF (iadv(iq)==14) THEN
122 CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, &
123 p, pk, teta)
124 ! Schema de Frederic Hourdin
125 ELSE IF (iadv(iq)==12) THEN
126 ! Pas de temps adaptatif
127 CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
128 IF (n>1) THEN
129 WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
130 'n=', n
131 END IF
132 DO indice = 1, n
133 CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
134 END DO
135 ELSE IF (iadv(iq)==13) THEN
136 ! Pas de temps adaptatif
137 CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
138 IF (n>1) THEN
139 WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
140 'n=', n
141 END IF
142 DO indice = 1, n
143 CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
144 END DO
145 ! Schema de pente SLOPES
146 ELSE IF (iadv(iq)==20) THEN
147 CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
148 ! Schema de Prather
149 ELSE IF (iadv(iq)==30) THEN
150 ! Pas de temps adaptatif
151 CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
152 IF (n>1) THEN
153 WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
154 'n=', n
155 END IF
156 CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
157 ! Schemas PPM Lin et Rood
158 ELSE IF (iadv(iq)==11 .OR. (iadv(iq)>=16 .AND. iadv(iq)<=18)) THEN
159 ! Test sur le flux horizontal
160 ! Pas de temps adaptatif
161 CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
162 IF (n>1) THEN
163 WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
164 'n=', n
165 END IF
166 ! Test sur le flux vertical
167 cflmaxz = 0.
168 DO l = 2, llm
169 DO ij = iip2, ip1jm
170 aaa = wg(ij, l)*dtvr/massem(ij, l)
171 cflmaxz = max(cflmaxz, aaa)
172 bbb = -wg(ij, l)*dtvr/massem(ij, l-1)
173 cflmaxz = max(cflmaxz, bbb)
174 END DO
175 END DO
176 IF (cflmaxz>=1) THEN
177 WRITE (*, *) 'WARNING vertical', 'CFLmaxz=', cflmaxz
178 END IF
179
180 ! Ss-prg interface LMDZ.4->PPM3d
181 CALL interpre(q(1, 1, iq), qppm(1, 1, iq), wg, fluxwppm, massem, &
182 apppm, bpppm, massebx, masseby, pbarug, pbarvg, unatppm, &
183 vnatppm, psppm)
184
185 DO indice = 1, n
186 ! VL (version PPM) horiz. et PPM vert.
187 IF (iadv(iq)==11) THEN
188 ! Ss-prg PPM3d de Lin
189 CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
190 vnatppm, fluxwppm, dtbon, 2, 2, 2, 1, iim, jjp1, 2, &
191 llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
192 ! Monotonic PPM
193 ELSE IF (iadv(iq)==16) THEN
194 ! Ss-prg PPM3d de Lin
195 CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
196 vnatppm, fluxwppm, dtbon, 3, 3, 3, 1, iim, jjp1, 2, &
197 llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
198 ! Semi Monotonic PPM
199 ELSE IF (iadv(iq)==17) THEN
200 ! Ss-prg PPM3d de Lin
201 CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
202 vnatppm, fluxwppm, dtbon, 4, 4, 4, 1, iim, jjp1, 2, &
203 llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
204 ! Positive Definite PPM
205 ELSE IF (iadv(iq)==18) THEN
206 ! Ss-prg PPM3d de Lin
207 CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
208 vnatppm, fluxwppm, dtbon, 5, 5, 5, 1, iim, jjp1, 2, &
209 llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
210 END IF
211 END DO
212
213 ! Ss-prg interface PPM3d-LMDZ.4
214 CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
215 END IF
216 END DO
217
218 ! on reinitialise a zero les flux de masse cumules
219 iadvtr = 0
220 END IF
221
222 END SUBROUTINE advtrac
223
224 end module advtrac_m
  ViewVC Help
Powered by ViewVC 1.1.21