Changeset 4171 for LMDZ6/branches/LMDZ-ECRAD/libf/dyn3d/advtrac.F90
- Timestamp:
- Jun 17, 2022, 4:24:49 PM (2 years ago)
- Location:
- LMDZ6/branches/LMDZ-ECRAD
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/dyn3d/advtrac.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/LMDZ-ECRAD
- Property svn:mergeinfo changed
-
LMDZ6/branches/LMDZ-ECRAD/libf/dyn3d/advtrac.F90
r2622 r4171 1 1 ! $Id$ 2 2 3 SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk) 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 infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif 12 USE control_mod, ONLY: iapp_tracvl, day_step 13 USE comconst_mod, ONLY: dtvr 14 15 IMPLICIT NONE 16 ! 17 include "dimensions.h" 18 include "paramet.h" 19 include "comdissip.h" 20 include "comgeom2.h" 21 include "description.h" 22 include "iniprint.h" 23 24 !------------------------------------------------------------------- 25 ! Arguments 26 !------------------------------------------------------------------- 27 INTEGER,INTENT(OUT) :: iapptrac 28 REAL,INTENT(IN) :: pbaru(ip1jmp1,llm) 29 REAL,INTENT(IN) :: pbarv(ip1jm,llm) 30 REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) 31 REAL,INTENT(IN) :: masse(ip1jmp1,llm) 32 REAL,INTENT(IN) :: p( ip1jmp1,llmp1 ) 33 REAL,INTENT(IN) :: teta(ip1jmp1,llm) 34 REAL,INTENT(IN) :: pk(ip1jmp1,llm) 35 REAL,INTENT(OUT) :: flxw(ip1jmp1,llm) 36 !------------------------------------------------------------------- 37 ! Ajout PPM 38 !-------------------------------------------------------- 39 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 40 !------------------------------------------------------------- 41 ! Variables locales 42 !------------------------------------------------------------- 43 44 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 45 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 46 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 47 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 48 INTEGER iadvtr 49 INTEGER ij,l,iq,iiq 50 REAL zdpmin, zdpmax 51 EXTERNAL minmax 52 SAVE iadvtr, massem, pbaruc, pbarvc 53 DATA iadvtr/0/ 54 !---------------------------------------------------------- 55 ! Rajouts pour PPM 56 !---------------------------------------------------------- 57 INTEGER indice,n 58 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 59 REAL CFLmaxz,aaa,bbb ! CFL maximum 60 REAL psppm(iim,jjp1) ! pression au sol 61 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 62 REAL qppm(iim*jjp1,llm,nqtot) 63 REAL fluxwppm(iim,jjp1,llm) 64 REAL apppm(llmp1), bpppm(llmp1) 65 LOGICAL dum,fill 66 DATA fill/.true./ 67 DATA dum/.true./ 68 69 integer,save :: countcfl=0 70 real cflx(ip1jmp1,llm) 71 real cfly(ip1jm,llm) 72 real cflz(ip1jmp1,llm) 73 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 74 75 IF(iadvtr.EQ.0) THEN 76 pbaruc(:,:)=0 77 pbarvc(:,:)=0 78 ENDIF 79 80 ! accumulation des flux de masse horizontaux 81 DO l=1,llm 82 DO ij = 1,ip1jmp1 83 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 84 ENDDO 85 DO ij = 1,ip1jm 86 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 87 ENDDO 88 ENDDO 89 90 ! selection de la masse instantannee des mailles avant le transport. 91 IF(iadvtr.EQ.0) THEN 92 3 #define DEBUG_IO 4 #undef DEBUG_IO 5 SUBROUTINE advtrac(pbaru, pbarv, p, masse,q,iapptrac,teta, flxw, pk) 6 ! Auteur : F. Hourdin 7 ! 8 ! Modif. P. Le Van (20/12/97) 9 ! F. Codron (10/99) 10 ! D. Le Croller (07/2001) 11 ! M.A Filiberti (04/2002) 12 ! 13 USE infotrac, ONLY: nqtot, tracers, isoCheck 14 USE control_mod, ONLY: iapp_tracvl, day_step 15 USE comconst_mod, ONLY: dtvr 16 17 IMPLICIT NONE 18 ! 19 include "dimensions.h" 20 include "paramet.h" 21 include "comdissip.h" 22 include "comgeom2.h" 23 include "description.h" 24 include "iniprint.h" 25 26 !--------------------------------------------------------------------------- 27 ! Arguments 28 !--------------------------------------------------------------------------- 29 INTEGER, INTENT(OUT) :: iapptrac 30 REAL, INTENT(IN) :: pbaru(ip1jmp1,llm) 31 REAL, INTENT(IN) :: pbarv(ip1jm, llm) 32 REAL, INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) 33 REAL, INTENT(IN) :: masse(ip1jmp1,llm) 34 REAL, INTENT(IN) :: p(ip1jmp1,llmp1 ) 35 REAL, INTENT(IN) :: teta(ip1jmp1,llm) 36 REAL, INTENT(IN) :: pk(ip1jmp1,llm) 37 REAL, INTENT(OUT) :: flxw(ip1jmp1,llm) 38 !--------------------------------------------------------------------------- 39 ! Ajout PPM 40 !--------------------------------------------------------------------------- 41 REAL :: massebx(ip1jmp1,llm), masseby(ip1jm,llm) 42 !--------------------------------------------------------------------------- 43 ! Variables locales 44 !--------------------------------------------------------------------------- 45 INTEGER :: ij, l, iq, iadv 46 ! REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu 47 REAL :: zdp(ip1jmp1), zdpmin, zdpmax 48 INTEGER, SAVE :: iadvtr=0 49 REAL, DIMENSION(ip1jmp1,llm) :: pbaruc, pbarug, massem, wg 50 REAL, DIMENSION(ip1jm, llm) :: pbarvc, pbarvg 51 EXTERNAL minmax 52 SAVE massem, pbaruc, pbarvc 53 !--------------------------------------------------------------------------- 54 ! Rajouts pour PPM 55 !--------------------------------------------------------------------------- 56 INTEGER indice, n 57 REAL :: dtbon ! Pas de temps adaptatif pour que CFL<1 58 REAL :: CFLmaxz, aaa, bbb ! CFL maximum 59 REAL, DIMENSION(iim,jjp1,llm) :: unatppm, vnatppm, fluxwppm 60 REAL :: qppm(iim*jjp1,llm,nqtot) 61 REAL :: psppm(iim,jjp1) ! pression au sol 62 REAL, DIMENSION(llmp1) :: apppm, bpppm 63 LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE. 64 65 INTEGER, SAVE :: countcfl=0 66 REAL, DIMENSION(ip1jmp1,llm) :: cflx, cflz 67 REAL, DIMENSION(ip1jm ,llm) :: cfly 68 REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax 69 70 IF(iadvtr == 0) THEN 71 pbaruc(:,:)=0 72 pbarvc(:,:)=0 73 END IF 74 75 !--- Accumulation des flux de masse horizontaux 76 DO l=1,llm 77 DO ij = 1,ip1jmp1 78 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 79 END DO 80 DO ij = 1,ip1jm 81 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 82 END DO 83 END DO 84 85 !--- Selection de la masse instantannee des mailles avant le transport. 86 IF(iadvtr == 0) THEN 93 87 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 94 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 95 ! 96 ENDIF 97 98 iadvtr = iadvtr+1 99 iapptrac = iadvtr 100 101 102 ! Test pour savoir si on advecte a ce pas de temps 103 IF ( iadvtr.EQ.iapp_tracvl ) THEN 104 105 !c .. Modif P.Le Van ( 20/12/97 ) .... 106 !c 107 108 ! traitement des flux de masse avant advection. 109 ! 1. calcul de w 110 ! 2. groupement des mailles pres du pole. 111 112 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 113 114 ! ... Flux de masse diaganostiques traceurs 115 flxw = wg / REAL(iapp_tracvl) 116 117 ! test sur l'eventuelle creation de valeurs negatives de la masse 118 DO l=1,llm-1 119 DO ij = iip2+1,ip1jm 120 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 121 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 122 + wg(ij,l+1) - wg(ij,l) 123 ENDDO 124 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 125 DO ij = iip2,ip1jm 126 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 127 ENDDO 128 129 130 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 131 132 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 133 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, & 134 ' MAX:', zdpmax 135 ENDIF 136 137 ENDDO 138 139 140 !------------------------------------------------------------------- 141 ! Calcul des criteres CFL en X, Y et Z 142 !------------------------------------------------------------------- 143 144 if (countcfl == 0. ) then 145 cflxmax(:)=0. 146 cflymax(:)=0. 147 cflzmax(:)=0. 148 endif 149 150 countcfl=countcfl+iapp_tracvl 151 cflx(:,:)=0. 152 cfly(:,:)=0. 153 cflz(:,:)=0. 154 do l=1,llm 155 do ij=iip2,ip1jm-1 156 if (pbarug(ij,l)>=0.) then 157 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 158 else 159 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 160 endif 161 enddo 162 enddo 163 do l=1,llm 164 do ij=iip2,ip1jm-1,iip1 165 cflx(ij+iip1,l)=cflx(ij,l) 166 enddo 167 enddo 168 169 do l=1,llm 170 do ij=1,ip1jm 171 if (pbarvg(ij,l)>=0.) then 172 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 173 else 174 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 175 endif 176 enddo 177 enddo 178 179 do l=2,llm 180 do ij=1,ip1jm 181 if (wg(ij,l)>=0.) then 182 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 183 else 184 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 185 endif 186 enddo 187 enddo 188 189 do l=1,llm 190 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 191 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 192 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 193 enddo 194 195 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 196 ! Par defaut, on sort le diagnostic des CFL tous les jours. 197 ! Si on veut le sortir a chaque pas d'advection en cas de plantage 198 ! if (countcfl==iapp_tracvl) then 199 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 200 if (countcfl==day_step) then 201 do l=1,llm 202 write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), & 203 cflzmax(l) 204 enddo 205 countcfl=0 206 endif 207 208 !------------------------------------------------------------------- 209 ! Advection proprement dite (Modification Le Croller (07/2001) 210 !------------------------------------------------------------------- 211 212 !---------------------------------------------------- 213 ! Calcul des moyennes basées sur la masse 214 !---------------------------------------------------- 215 call massbar(massem,massebx,masseby) 216 217 !----------------------------------------------------------- 218 ! Appel des sous programmes d'advection 219 !----------------------------------------------------------- 220 221 if (ok_iso_verif) then 222 write(*,*) 'advtrac 227' 223 call check_isotopes_seq(q,ip1jmp1,'advtrac 162') 224 endif !if (ok_iso_verif) then 225 226 do iq=1,nqperes 227 ! call clock(t_initial) 228 if(iadv(iq) == 0) cycle 229 ! ---------------------------------------------------------------- 230 ! Schema de Van Leer I MUSCL 231 ! ---------------------------------------------------------------- 232 if(iadv(iq).eq.10) THEN 233 ! CRisi: on fait passer tout q pour avoir acces aux fils 234 235 !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) 236 call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 237 238 ! ---------------------------------------------------------------- 239 ! Schema "pseudo amont" + test sur humidite specifique 240 ! pour la vapeur d'eau. F. Codron 241 ! ---------------------------------------------------------------- 242 else if(iadv(iq).eq.14) then 243 ! 244 !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) 245 CALL vlspltqs( q, 2., massem, wg , & 246 pbarug,pbarvg,dtvr,p,pk,teta,iq) 247 248 ! ---------------------------------------------------------------- 249 ! Schema de Frederic Hourdin 250 ! ---------------------------------------------------------------- 251 else if(iadv(iq).eq.12) then 252 ! Pas de temps adaptatif 253 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 254 if (n.GT.1) then 255 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 256 dtvr,'n=',n 257 endif 258 do indice=1,n 259 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 260 end do 261 else if(iadv(iq).eq.13) then 262 ! Pas de temps adaptatif 263 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 264 if (n.GT.1) then 265 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 266 dtvr,'n=',n 267 endif 268 do indice=1,n 269 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 270 end do 271 ! ---------------------------------------------------------------- 272 ! Schema de pente SLOPES 273 ! ---------------------------------------------------------------- 274 else if (iadv(iq).eq.20) then 275 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 276 277 ! ---------------------------------------------------------------- 278 ! Schema de Prather 279 ! ---------------------------------------------------------------- 280 else if (iadv(iq).eq.30) then 281 ! Pas de temps adaptatif 282 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 283 if (n.GT.1) then 284 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 285 dtvr,'n=',n 286 endif 287 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, & 288 n,dtbon) 289 290 ! ---------------------------------------------------------------- 291 ! Schemas PPM Lin et Rood 292 ! ---------------------------------------------------------------- 293 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. & 294 iadv(iq).LE.18)) then 295 296 ! Test sur le flux horizontal 297 ! Pas de temps adaptatif 298 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 299 if (n.GT.1) then 300 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 301 dtvr,'n=',n 302 endif 303 ! Test sur le flux vertical 304 CFLmaxz=0. 305 do l=2,llm 306 do ij=iip2,ip1jm 307 aaa=wg(ij,l)*dtvr/massem(ij,l) 308 CFLmaxz=max(CFLmaxz,aaa) 309 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 310 CFLmaxz=max(CFLmaxz,bbb) 311 enddo 312 enddo 313 if (CFLmaxz.GE.1) then 314 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 315 endif 316 317 !----------------------------------------------------------- 318 ! Ss-prg interface LMDZ.4->PPM3d 319 !----------------------------------------------------------- 320 321 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 322 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 323 unatppm,vnatppm,psppm) 324 325 do indice=1,n 326 !---------------------------------------------------------------- 327 ! VL (version PPM) horiz. et PPM vert. 328 !---------------------------------------------------------------- 329 if (iadv(iq).eq.11) then 330 ! Ss-prg PPM3d de Lin 331 call ppm3d(1,qppm(1,1,iq), & 332 psppm,psppm, & 333 unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, & 334 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 335 fill,dum,220.) 336 337 !------------------------------------------------------------- 338 ! Monotonic PPM 339 !------------------------------------------------------------- 340 else if (iadv(iq).eq.16) then 341 ! Ss-prg PPM3d de Lin 342 call ppm3d(1,qppm(1,1,iq), & 343 psppm,psppm, & 344 unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, & 345 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 346 fill,dum,220.) 347 !------------------------------------------------------------- 348 349 !------------------------------------------------------------- 350 ! Semi Monotonic PPM 351 !------------------------------------------------------------- 352 else if (iadv(iq).eq.17) then 353 ! Ss-prg PPM3d de Lin 354 call ppm3d(1,qppm(1,1,iq), & 355 psppm,psppm, & 356 unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, & 357 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 358 fill,dum,220.) 359 !------------------------------------------------------------- 360 361 !------------------------------------------------------------- 362 ! Positive Definite PPM 363 !------------------------------------------------------------- 364 else if (iadv(iq).eq.18) then 365 ! Ss-prg PPM3d de Lin 366 call ppm3d(1,qppm(1,1,iq), & 367 psppm,psppm, & 368 unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, & 369 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 370 fill,dum,220.) 371 !------------------------------------------------------------- 372 endif 373 enddo 374 !----------------------------------------------------------------- 375 ! Ss-prg interface PPM3d-LMDZ.4 376 !----------------------------------------------------------------- 377 call interpost(q(1,1,iq),qppm(1,1,iq)) 378 endif 379 !---------------------------------------------------------------------- 380 381 !----------------------------------------------------------------- 382 ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 383 ! et Nord j=1 384 !----------------------------------------------------------------- 385 386 ! call traceurpole(q(1,1,iq),massem) 387 388 ! calcul du temps cpu pour un schema donne 389 390 ! call clock(t_final) 391 !ym tps_cpu=t_final-t_initial 392 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 393 394 end DO 395 396 if (ok_iso_verif) then 397 write(*,*) 'advtrac 402' 398 call check_isotopes_seq(q,ip1jmp1,'advtrac 397') 399 endif !if (ok_iso_verif) then 400 401 !------------------------------------------------------------------ 402 ! on reinitialise a zero les flux de masse cumules 403 !--------------------------------------------------- 404 iadvtr=0 405 406 ENDIF ! if iadvtr.EQ.iapp_tracvl 88 ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 89 END IF 90 91 iadvtr = iadvtr+1 92 iapptrac = iadvtr 93 94 !--- Test pour savoir si on advecte a ce pas de temps 95 IF(iadvtr /= iapp_tracvl) RETURN 96 97 ! .. Modif P.Le Van ( 20/12/97 ) .... 98 ! 99 ! traitement des flux de masse avant advection. 100 ! 1. calcul de w 101 ! 2. groupement des mailles pres du pole. 102 103 CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg) 104 105 !--- Flux de masse diaganostiques traceurs 106 flxw = wg / REAL(iapp_tracvl) 107 108 !--- Test sur l'eventuelle creation de valeurs negatives de la masse 109 DO l=1,llm-1 110 DO ij = iip2+1,ip1jm 111 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 112 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 113 + wg(ij,l+1) - wg(ij,l) 114 END DO 115 ! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ? 116 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 117 DO ij = iip2,ip1jm 118 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 119 END DO 120 121 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 122 123 IF(MAX(ABS(zdpmin),ABS(zdpmax)) > 0.5) & 124 WRITE(*,*)'WARNING DP/P l=',l,' MIN:',zdpmin,' MAX:', zdpmax 125 126 END DO 127 128 !------------------------------------------------------------------------- 129 ! Calcul des criteres CFL en X, Y et Z 130 !------------------------------------------------------------------------- 131 IF(countcfl == 0. ) then 132 cflxmax(:)=0. 133 cflymax(:)=0. 134 cflzmax(:)=0. 135 END IF 136 137 countcfl=countcfl+iapp_tracvl 138 cflx(:,:)=0. 139 cfly(:,:)=0. 140 cflz(:,:)=0. 141 DO l=1,llm 142 DO ij=iip2,ip1jm-1 143 IF(pbarug(ij,l)>=0.) then 144 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 145 ELSE 146 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 147 END IF 148 END DO 149 END DO 150 151 DO l=1,llm 152 DO ij=iip2,ip1jm-1,iip1 153 cflx(ij+iip1,l)=cflx(ij,l) 154 END DO 155 END DO 156 157 DO l=1,llm 158 DO ij=1,ip1jm 159 IF(pbarvg(ij,l)>=0.) then 160 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 161 ELSE 162 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 163 END IF 164 END DO 165 END DO 166 167 DO l=2,llm 168 DO ij=1,ip1jm 169 IF(wg(ij,l) >= 0.) THEN 170 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 171 ELSE 172 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 173 END IF 174 END DO 175 END DO 176 177 DO l=1,llm 178 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 179 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 180 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 181 END DO 182 183 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 184 ! Par defaut, on sort le diagnostic des CFL tous les jours. 185 ! Si on veut le sortir a chaque pas d'advection en cas de plantage 186 ! IF(countcfl==iapp_tracvl) then 187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 188 IF(countcfl==day_step) then 189 DO l=1,llm 190 WRITE(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l) 191 END DO 192 countcfl=0 193 END IF 194 195 !--------------------------------------------------------------------------- 196 ! Advection proprement dite (Modification Le Croller (07/2001) 197 !--------------------------------------------------------------------------- 198 199 !--------------------------------------------------------------------------- 200 ! Calcul des moyennes basees sur la masse 201 !--------------------------------------------------------------------------- 202 CALL massbar(massem,massebx,masseby) 203 204 #ifdef DEBUG_IO 205 CALL WriteField_u('massem',massem) 206 CALL WriteField_u('wg',wg) 207 CALL WriteField_u('pbarug',pbarug) 208 CALL WriteField_v('pbarvg',pbarvg) 209 CALL WriteField_u('p_tmp',p) 210 CALL WriteField_u('pk_tmp',pk) 211 CALL WriteField_u('teta_tmp',teta) 212 DO iq=1,nqtot 213 CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq)) 214 END DO 215 #endif 216 217 IF(isoCheck) WRITE(*,*) 'advtrac 227' 218 CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162') 219 220 !------------------------------------------------------------------------- 221 ! Appel des sous programmes d'advection 222 !------------------------------------------------------------------------- 223 DO iq = 1, nqtot 224 ! CALL clock(t_initial) 225 IF(tracers(iq)%parent /= 'air') CYCLE 226 iadv = tracers(iq)%iadv 227 !----------------------------------------------------------------------- 228 SELECT CASE(iadv) 229 !----------------------------------------------------------------------- 230 CASE(0); CYCLE 231 !-------------------------------------------------------------------- 232 CASE(10) !--- Schema de Van Leer I MUSCL 233 !-------------------------------------------------------------------- 234 ! WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) 235 CALL vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 236 237 !-------------------------------------------------------------------- 238 CASE(14) !--- Schema "pseuDO amont" + test sur humidite specifique 239 !--- pour la vapeur d'eau. F. Codron 240 !-------------------------------------------------------------------- 241 ! WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) 242 CALL vlspltqs(q,2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta,iq) 243 244 !-------------------------------------------------------------------- 245 CASE(12) !--- Schema de Frederic Hourdin 246 !-------------------------------------------------------------------- 247 CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif 248 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 249 DO indice=1,n 250 CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 251 END DO 252 253 !-------------------------------------------------------------------- 254 CASE(13) !--- Pas de temps adaptatif 255 !-------------------------------------------------------------------- 256 CALL adaptdt(iadv,dtbon,n,pbarug,massem) 257 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 258 DO indice=1,n 259 CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 260 END DO 261 262 !-------------------------------------------------------------------- 263 CASE(20) !--- Schema de pente SLOPES 264 !-------------------------------------------------------------------- 265 CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 266 267 !-------------------------------------------------------------------- 268 CASE(30) !--- Schema de Prather 269 !-------------------------------------------------------------------- 270 ! Pas de temps adaptatif 271 CALL adaptdt(iadv,dtbon,n,pbarug,massem) 272 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 273 CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon) 274 275 !-------------------------------------------------------------------- 276 CASE(11,16,17,18) !--- Schemas PPM Lin et Rood 277 !-------------------------------------------------------------------- 278 ! Test sur le flux horizontal 279 CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif 280 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 281 ! Test sur le flux vertical 282 CFLmaxz=0. 283 DO l=2,llm 284 DO ij=iip2,ip1jm 285 aaa=wg(ij,l)*dtvr/massem(ij,l) 286 CFLmaxz=max(CFLmaxz,aaa) 287 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 288 CFLmaxz=max(CFLmaxz,bbb) 289 END DO 290 END DO 291 IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 292 !---------------------------------------------------------------- 293 ! Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin) 294 !---------------------------------------------------------------- 295 CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 296 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 297 unatppm,vnatppm,psppm) 298 299 !---------------------------------------------------------------- 300 DO indice=1,n !--- VL (version PPM) horiz. et PPM vert. 301 !---------------------------------------------------------------- 302 SELECT CASE(iadv) 303 !---------------------------------------------------------- 304 CASE(11) 305 !---------------------------------------------------------- 306 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 307 2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) 308 !---------------------------------------------------------- 309 CASE(16) !--- Monotonic PPM 310 !---------------------------------------------------------- 311 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 312 3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) 313 !---------------------------------------------------------- 314 CASE(17) !--- Semi monotonic PPM 315 !---------------------------------------------------------- 316 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 317 4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.) 318 !---------------------------------------------------------- 319 CASE(18) !--- Positive Definite PPM 320 !---------------------------------------------------------- 321 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 322 5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) 323 END SELECT 324 !---------------------------------------------------------------- 325 END DO 326 !---------------------------------------------------------------- 327 ! Ss-prg interface PPM3d-LMDZ.4 328 !---------------------------------------------------------------- 329 CALL interpost(q(1,1,iq),qppm(1,1,iq)) 330 !---------------------------------------------------------------------- 331 END SELECT 332 !---------------------------------------------------------------------- 333 334 !---------------------------------------------------------------------- 335 ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1 336 !---------------------------------------------------------------------- 337 ! CALL traceurpole(q(1,1,iq),massem) 338 339 !--- Calcul du temps cpu pour un schema donne 340 ! CALL clock(t_final) 341 !ym tps_cpu=t_final-t_initial 342 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 343 344 END DO 345 346 IF(isoCheck) WRITE(*,*) 'advtrac 402' 347 CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397') 348 349 !------------------------------------------------------------------------- 350 ! on reinitialise a zero les flux de masse cumules 351 !------------------------------------------------------------------------- 352 iadvtr=0 407 353 408 354 END SUBROUTINE advtrac
Note: See TracChangeset
for help on using the changeset viewer.
