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