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

Contents of /trunk/dyn3d/advtrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
Original Path: trunk/Sources/dyn3d/advtrac.f
File size: 8016 byte(s)
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.

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 interpre_m, only: interpre
17 use massbar_m, only: massbar
18 USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
19 llmp1
20 use vlsplt_m, only: vlsplt
21 use vlspltqs_m, only: vlspltqs
22
23 REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
24 REAL, intent(in):: p(ip1jmp1, llmp1)
25 real, intent(in):: masse(ip1jmp1, llm)
26 REAL, intent(inout):: q(ip1jmp1, llm, nqmx)
27 INTEGER, intent(out):: iapptrac
28 real, intent(in):: teta(ip1jmp1, llm)
29 REAL, intent(in):: pk(ip1jmp1, llm)
30
31 ! Variables locales
32
33 REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)
34 REAL, save:: pbaruc(ip1jmp1, llm), pbarvc(ip1jm, llm)
35 REAL, save:: massem(ip1jmp1, llm)
36 real zdp(ip1jmp1)
37 REAL pbarug(ip1jmp1, llm), pbarvg(ip1jm, llm), wg(ip1jmp1, llm)
38
39 INTEGER:: iadvtr = 0
40 INTEGER ij, l, iq
41 REAL zdpmin, zdpmax
42 EXTERNAL minmax
43
44 ! Rajouts pour PPM
45
46 INTEGER indice, n
47 ! Pas de temps adaptatif pour que CFL < 1
48 REAL dtbon
49 REAL cflmaxz ! CFL maximum
50 real aaa, bbb
51 REAL psppm(iim, jjp1) ! pression au sol
52 REAL unatppm(iim, jjp1, llm), vnatppm(iim, jjp1, llm)
53 REAL qppm(iim*jjp1, llm, nqmx)
54 REAL fluxwppm(iim, jjp1, llm)
55 REAL apppm(llmp1), bpppm(llmp1)
56 LOGICAL:: dum = .TRUE., fill = .TRUE.
57
58 !-----------------------------------------------------------
59
60 IF (iadvtr==0) THEN
61 CALL initial0(ijp1llm, pbaruc)
62 CALL initial0(ijmllm, pbarvc)
63 END IF
64
65 ! accumulation des flux de masse horizontaux
66 DO l = 1, llm
67 DO ij = 1, ip1jmp1
68 pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
69 END DO
70 DO ij = 1, ip1jm
71 pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
72 END DO
73 END DO
74
75 ! selection de la masse instantannee des mailles avant le transport.
76 IF (iadvtr==0) massem = masse
77
78 iadvtr = iadvtr + 1
79 iapptrac = iadvtr
80
81 ! Test pour savoir si on advecte a ce pas de temps
82 IF (iadvtr == iapp_tracvl) THEN
83 ! traitement des flux de masse avant advection.
84 ! 1. calcul de w
85 ! 2. groupement des mailles pres du pole.
86
87 CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)
88
89 ! test sur l'eventuelle creation de valeurs negatives de la masse
90 DO l = 1, llm - 1
91 DO ij = iip2 + 1, ip1jm
92 zdp(ij) = pbarug(ij-1, l) - pbarug(ij, l) - pbarvg(ij-iip1, l) + &
93 pbarvg(ij, l) + wg(ij, l+1) - wg(ij, l)
94 END DO
95 CALL scopy(jjm-1, zdp(iip1+iip1), iip1, zdp(iip2), iip1)
96 DO ij = iip2, ip1jm
97 zdp(ij) = zdp(ij)*dtvr/massem(ij, l)
98 END DO
99
100 CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)
101
102 IF (max(abs(zdpmin), abs(zdpmax))>0.5) THEN
103 PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
104 END IF
105 END DO
106
107 ! Advection proprement dite
108
109 ! Calcul des moyennes bas\'ees sur la masse
110 CALL massbar(massem, massebx, masseby)
111
112 ! Appel des sous programmes d'advection
113
114 DO iq = 1, nqmx
115 select case (iadv(iq))
116 case (10)
117 ! Schema de Van Leer I MUSCL
118 CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
119 case (14)
120 ! Schema "pseudo amont" + test sur humidite specifique
121 ! pour la vapeur d'eau. F. Codron
122 CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, &
123 p, pk, teta)
124 case (12)
125 ! Schema de Frederic Hourdin
126 ! Pas de temps adaptatif
127 CALL adaptdt(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 case (13)
136 ! Pas de temps adaptatif
137 CALL adaptdt(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 case (20)
146 ! Schema de pente SLOPES
147 CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
148 case (30)
149 ! Schema de Prather
150 ! Pas de temps adaptatif
151 CALL adaptdt(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 case (11, 16:18)
158 ! Schemas PPM Lin et Rood
159 ! Test sur le flux horizontal
160 ! Pas de temps adaptatif
161 CALL adaptdt(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 select
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