/[lmdze]/trunk/libf/dyn3d/advtrac.f90
ViewVC logotype

Contents of /trunk/libf/dyn3d/advtrac.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (show annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
File size: 11341 byte(s) 
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".
1 SUBROUTINE advtrac(pbaru,pbarv,p,masse,q,iapptrac,teta,pk)
2
3 ! From dyn3d/advtrac.F,v 1.4 2005/04/13 08:58:34
4 ! Auteur : F. Hourdin
5
6 ! Modif. P. Le Van (20/12/97)
7 ! F. Codron (10/99)
8 ! D. Le Croller (07/2001)
9 ! M.A Filiberti (04/2002)
10
11 USE dimens_m
12 USE paramet_m
13 USE comconst
14 USE comvert
15 USE conf_gcm_m
16 USE logic
17 USE comgeom
18 USE temps
19 USE ener
20 USE advtrac_m
21 USE comdissip
22 IMPLICIT NONE
23
24
25 !-------------------------------------------------------------------
26 ! Arguments
27 !-------------------------------------------------------------------
28 ! Ajout PPM
29 !--------------------------------------------------------
30 REAL massebx(ip1jmp1,llm), masseby(ip1jm,llm)
31 !--------------------------------------------------------
32 INTEGER iapptrac
33 REAL pbaru(ip1jmp1,llm), pbarv(ip1jm,llm)
34 REAL q(ip1jmp1,llm,nqmx), masse(ip1jmp1,llm)
35 REAL, intent(in):: p(ip1jmp1,llmp1)
36 real teta(ip1jmp1,llm)
37 REAL pk(ip1jmp1,llm)
38
39 !-------------------------------------------------------------
40 ! Variables locales
41 !-------------------------------------------------------------
42
43 REAL pbaruc(ip1jmp1,llm), pbarvc(ip1jm,llm)
44 REAL massem(ip1jmp1,llm), zdp(ip1jmp1)
45 REAL pbarug(ip1jmp1,llm), pbarvg(ip1jm,llm), wg(ip1jmp1,llm)
46 REAL cpuadv(nqmx)
47 COMMON /cpuadv/cpuadv
48
49 INTEGER iadvtr
50 INTEGER ij, l, iq
51 REAL zdpmin, zdpmax
52 EXTERNAL minmax
53 SAVE iadvtr, massem, pbaruc, pbarvc
54 DATA iadvtr/0/
55 !----------------------------------------------------------
56 ! Rajouts pour PPM
57 !----------------------------------------------------------
58 INTEGER indice, n
59 ! Pas de temps adaptatif pour que CFL<1
60 REAL dtbon
61 REAL cflmaxz, aaa, & ! CFL maximum
62 bbb
63 REAL psppm(iim,jjp1) ! pression au sol
64 REAL unatppm(iim,jjp1,llm), vnatppm(iim,jjp1,llm)
65 REAL qppm(iim*jjp1,llm,nqmx)
66 REAL fluxwppm(iim,jjp1,llm)
67 REAL apppm(llmp1), bpppm(llmp1)
68 LOGICAL dum, fill
69 DATA fill/ .TRUE./
70 DATA dum/ .TRUE./
71
72
73 IF (iadvtr==0) THEN
74 CALL initial0(ijp1llm,pbaruc)
75 CALL initial0(ijmllm,pbarvc)
76 END IF
77
78 ! accumulation des flux de masse horizontaux
79 DO l = 1, llm
80 DO ij = 1, ip1jmp1
81 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
82 END DO
83 DO ij = 1, ip1jm
84 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
85 END DO
86 END DO
87
88 ! selection de la masse instantannee des mailles avant le transport.
89 IF (iadvtr==0) THEN
90
91 CALL scopy(ip1jmp1*llm,masse,1,massem,1)
92 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
93
94 END IF
95
96 iadvtr = iadvtr + 1
97 iapptrac = iadvtr
98
99
100 ! Test pour savoir si on advecte a ce pas de temps
101 IF (iadvtr==iapp_tracvl) THEN
102
103 !c .. Modif P.Le Van ( 20/12/97 ) ....
104 !c
105
106 ! traitement des flux de masse avant advection.
107 ! 1. calcul de w
108 ! 2. groupement des mailles pres du pole.
109
110 CALL groupe(massem,pbaruc,pbarvc,pbarug,pbarvg,wg)
111
112
113 ! test sur l'eventuelle creation de valeurs negatives de la masse
114 DO l = 1, llm - 1
115 DO ij = iip2 + 1, ip1jm
116 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) - pbarvg(ij-iip1,l) + &
117 pbarvg(ij,l) + wg(ij,l+1) - wg(ij,l)
118 END DO
119 CALL scopy(jjm-1,zdp(iip1+iip1),iip1,zdp(iip2),iip1)
120 DO ij = iip2, ip1jm
121 zdp(ij) = zdp(ij)*dtvr/massem(ij,l)
122 END DO
123
124
125 CALL minmax(ip1jm-iip1,zdp(iip2),zdpmin,zdpmax)
126
127 IF (max(abs(zdpmin),abs(zdpmax))>0.5) THEN
128 PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
129 END IF
130
131 END DO
132
133 !-------------------------------------------------------------------
134 ! Advection proprement dite (Modification Le Croller (07/2001)
135 !-------------------------------------------------------------------
136
137 !----------------------------------------------------
138 ! Calcul des moyennes basées sur la masse
139 !----------------------------------------------------
140 CALL massbar(massem,massebx,masseby)
141
142 !-----------------------------------------------------------
143 ! Appel des sous programmes d'advection
144 !-----------------------------------------------------------
145 DO iq = 1, nqmx
146 ! call clock(t_initial)
147 IF (iadv(iq)==0) CYCLE
148 ! ----------------------------------------------------------------
149 ! Schema de Van Leer I MUSCL
150 ! ----------------------------------------------------------------
151 IF (iadv(iq)==10) THEN
152 CALL vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
153
154
155 ! ----------------------------------------------------------------
156 ! Schema "pseudo amont" + test sur humidite specifique
157 ! pour la vapeur d'eau. F. Codron
158 ! ----------------------------------------------------------------
159 ELSE IF (iadv(iq)==14) THEN
160
161 CALL vlspltqs(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta)
162 ! ----------------------------------------------------------------
163 ! Schema de Frederic Hourdin
164 ! ----------------------------------------------------------------
165 ELSE IF (iadv(iq)==12) THEN
166 ! Pas de temps adaptatif
167 CALL adaptdt(iadv(iq),dtbon,n,pbarug,massem)
168 IF (n>1) THEN
169 WRITE (*,*) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
170 'n=', n
171 END IF
172 DO indice = 1, n
173 CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
174 END DO
175 ELSE IF (iadv(iq)==13) THEN
176 ! Pas de temps adaptatif
177 CALL adaptdt(iadv(iq),dtbon,n,pbarug,massem)
178 IF (n>1) THEN
179 WRITE (*,*) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
180 'n=', n
181 END IF
182 DO indice = 1, n
183 CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
184 END DO
185 ! ----------------------------------------------------------------
186 ! Schema de pente SLOPES
187 ! ----------------------------------------------------------------
188 ELSE IF (iadv(iq)==20) THEN
189 CALL pentes_ini(q(1,1,iq),wg,massem,pbarug,pbarvg,0)
190
191 ! ----------------------------------------------------------------
192 ! Schema de Prather
193 ! ----------------------------------------------------------------
194 ELSE IF (iadv(iq)==30) THEN
195 ! Pas de temps adaptatif
196 CALL adaptdt(iadv(iq),dtbon,n,pbarug,massem)
197 IF (n>1) THEN
198 WRITE (*,*) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
199 'n=', n
200 END IF
201 CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
202 ! ----------------------------------------------------------------
203 ! Schemas PPM Lin et Rood
204 ! ----------------------------------------------------------------
205 ELSE IF (iadv(iq)==11 .OR. (iadv(iq)>=16 .AND. iadv(iq)<=18)) THEN
206
207 ! Test sur le flux horizontal
208 ! Pas de temps adaptatif
209 CALL adaptdt(iadv(iq),dtbon,n,pbarug,massem)
210 IF (n>1) THEN
211 WRITE (*,*) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
212 'n=', n
213 END IF
214 ! Test sur le flux vertical
215 cflmaxz = 0.
216 DO l = 2, llm
217 DO ij = iip2, ip1jm
218 aaa = wg(ij,l)*dtvr/massem(ij,l)
219 cflmaxz = max(cflmaxz,aaa)
220 bbb = -wg(ij,l)*dtvr/massem(ij,l-1)
221 cflmaxz = max(cflmaxz,bbb)
222 END DO
223 END DO
224 IF (cflmaxz>=1) THEN
225 WRITE (*,*) 'WARNING vertical', 'CFLmaxz=', cflmaxz
226 END IF
227
228 !-----------------------------------------------------------
229 ! Ss-prg interface LMDZ.4->PPM3d
230 !-----------------------------------------------------------
231
232 CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,apppm, &
233 bpppm,massebx,masseby,pbarug,pbarvg,unatppm,vnatppm,psppm)
234
235 DO indice = 1, n
236 !----------------------------------------------------------------
237 ! VL (version PPM) horiz. et PPM vert.
238 !---------------------------------------------------------------
239 IF (iadv(iq)==11) THEN
240 ! Ss-prg PPM3d de Lin
241 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm, &
242 fluxwppm,dtbon,2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01, &
243 6400000,fill,dum,220.)
244
245 !-----------------------------------------------------------
246 ! Monotonic PPM
247 !-------------------------------------------------------
248 ELSE IF (iadv(iq)==16) THEN
249 ! Ss-prg PPM3d de Lin
250 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm, &
251 fluxwppm,dtbon,3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01, &
252 6400000,fill,dum,220.)
253 ! Semi Monotonic PPM
254 ELSE IF (iadv(iq)==17) THEN
255 ! Ss-prg PPM3d de Lin
256 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm, &
257 fluxwppm,dtbon,4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01, &
258 6400000,fill,dum,220.)
259 ! Positive Definite PPM
260 ELSE IF (iadv(iq)==18) THEN
261 ! Ss-prg PPM3d de Lin
262 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm, &
263 fluxwppm,dtbon,5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01, &
264 6400000,fill,dum,220.)
265 END IF
266 END DO
267 !-----------------------------------------------------------------
268 ! Ss-prg interface PPM3d-LMDZ.4
269 !-----------------------------------------------------------------
270 CALL interpost(q(1,1,iq),qppm(1,1,iq))
271 END IF
272 !----------------------------------------------------------------------
273
274 !-----------------------------------------------------------------
275 ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
276 ! et Nord j=1
277 !-----------------------------------------------------------------
278
279 ! call traceurpole(q(1,1,iq),massem)
280
281 ! calcul du temps cpu pour un schema donne
282
283 ! call clock(t_final)
284 !ym tps_cpu=t_final-t_initial
285 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu
286
287 END DO
288
289
290 !------------------------------------------------------------------
291 ! on reinitialise a zero les flux de masse cumules
292 !---------------------------------------------------
293 iadvtr = 0
294
295
296 END IF
297 ! if iadvtr.EQ.iapp_tracvl
298 RETURN
299 END SUBROUTINE advtrac
  ViewVC Help
Powered by ViewVC 1.1.21