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

Annotation of /trunk/Sources/dyn3d/advtrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (hide annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
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 guez 40 module advtrac_m
2 guez 3
3 guez 40 IMPLICIT NONE
4 guez 3
5 guez 40 contains
6 guez 23
7 guez 40 SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, pk)
8 guez 3
9 guez 40 ! From dyn3d/advtrac.F, version 1.4 2005/04/13 08:58:34
10     ! Author: F. Hourdin
11 guez 3
12 guez 71 USE comconst, ONLY : dtvr
13     USE conf_gcm_m, ONLY : iapp_tracvl
14 guez 40 USE dimens_m, ONLY : iim, jjm, llm, nqmx
15 guez 71 USE iniadvtrac_m, ONLY : iadv
16 guez 108 use interpre_m, only: interpre
17 guez 91 use massbar_m, only: massbar
18 guez 40 USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
19     llmp1
20 guez 157 use vlsplt_m, only: vlsplt
21 guez 108 use vlspltqs_m, only: vlspltqs
22 guez 31
23 guez 71 REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
24     REAL, intent(in):: p(ip1jmp1, llmp1)
25     real, intent(in):: masse(ip1jmp1, llm)
26 guez 40 REAL, intent(inout):: q(ip1jmp1, llm, nqmx)
27 guez 71 INTEGER, intent(out):: iapptrac
28 guez 44 real, intent(in):: teta(ip1jmp1, llm)
29 guez 71 REAL, intent(in):: pk(ip1jmp1, llm)
30 guez 3
31 guez 40 ! Variables locales
32 guez 3
33 guez 40 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 guez 31
39 guez 40 INTEGER:: iadvtr = 0
40     INTEGER ij, l, iq
41     REAL zdpmin, zdpmax
42     EXTERNAL minmax
43 guez 31
44 guez 40 ! Rajouts pour PPM
45 guez 3
46 guez 40 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 guez 3
58 guez 40 !-----------------------------------------------------------
59 guez 3
60 guez 40 IF (iadvtr==0) THEN
61     CALL initial0(ijp1llm, pbaruc)
62     CALL initial0(ijmllm, pbarvc)
63     END IF
64 guez 3
65 guez 40 ! 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 guez 3
75 guez 40 ! selection de la masse instantannee des mailles avant le transport.
76 guez 71 IF (iadvtr==0) massem = masse
77 guez 3
78 guez 40 iadvtr = iadvtr + 1
79     iapptrac = iadvtr
80 guez 3
81 guez 40 ! Test pour savoir si on advecte a ce pas de temps
82 guez 71 IF (iadvtr == iapp_tracvl) THEN
83 guez 40 ! traitement des flux de masse avant advection.
84     ! 1. calcul de w
85     ! 2. groupement des mailles pres du pole.
86 guez 3
87 guez 150 CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)
88 guez 3
89 guez 40 ! 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 guez 3
100 guez 40 CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)
101 guez 3
102 guez 40 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 guez 3
107 guez 40 ! Advection proprement dite
108 guez 3
109 guez 150 ! Calcul des moyennes bas\'ees sur la masse
110 guez 40 CALL massbar(massem, massebx, masseby)
111 guez 31
112 guez 40 ! Appel des sous programmes d'advection
113 guez 31
114 guez 40 DO iq = 1, nqmx
115 guez 157 select case (iadv(iq))
116     case (10)
117     ! Schema de Van Leer I MUSCL
118 guez 40 CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
119 guez 157 case (14)
120 guez 40 ! 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 guez 157 case (12)
125 guez 40 ! Schema de Frederic Hourdin
126     ! Pas de temps adaptatif
127 guez 104 CALL adaptdt(dtbon, n, pbarug, massem)
128 guez 40 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 guez 157 case (13)
136 guez 40 ! Pas de temps adaptatif
137 guez 104 CALL adaptdt(dtbon, n, pbarug, massem)
138 guez 40 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 guez 157 case (20)
146 guez 40 ! Schema de pente SLOPES
147     CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
148 guez 157 case (30)
149 guez 40 ! Schema de Prather
150     ! Pas de temps adaptatif
151 guez 104 CALL adaptdt(dtbon, n, pbarug, massem)
152 guez 40 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 guez 157 case (11, 16:18)
158 guez 40 ! Schemas PPM Lin et Rood
159     ! Test sur le flux horizontal
160     ! Pas de temps adaptatif
161 guez 104 CALL adaptdt(dtbon, n, pbarug, massem)
162 guez 40 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 guez 3
180 guez 40 ! 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 guez 31
185 guez 40 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 guez 3
213 guez 40 ! Ss-prg interface PPM3d-LMDZ.4
214     CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
215 guez 157 END select
216 guez 40 END DO
217 guez 3
218 guez 40 ! 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