! $Id: advtrac.F90 1907 2013-11-26 13:10:46Z fhourdin $ SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk) ! Auteur : F. Hourdin ! ! Modif. P. Le Van (20/12/97) ! F. Codron (10/99) ! D. Le Croller (07/2001) ! M.A Filiberti (04/2002) ! USE infotrac USE control_mod IMPLICIT NONE ! include "dimensions.h" include "paramet.h" include "comconst.h" include "comvert.h" include "comdissip.h" include "comgeom2.h" include "logic.h" include "temps.h" include "ener.h" include "description.h" include "iniprint.h" !------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------- ! Ajout PPM !-------------------------------------------------------- REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) !-------------------------------------------------------- INTEGER iapptrac REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) REAL pk(ip1jmp1,llm) REAL flxw(ip1jmp1,llm) !------------------------------------------------------------- ! Variables locales !------------------------------------------------------------- REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) REAL massem(ip1jmp1,llm),zdp(ip1jmp1) REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu INTEGER iadvtr INTEGER ij,l,iq,iiq REAL zdpmin, zdpmax EXTERNAL minmax SAVE iadvtr, massem, pbaruc, pbarvc DATA iadvtr/0/ !---------------------------------------------------------- ! Rajouts pour PPM !---------------------------------------------------------- INTEGER indice,n REAL dtbon ! Pas de temps adaptatif pour que CFL<1 REAL CFLmaxz,aaa,bbb ! CFL maximum REAL psppm(iim,jjp1) ! pression au sol REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) REAL qppm(iim*jjp1,llm,nqtot) REAL fluxwppm(iim,jjp1,llm) REAL apppm(llmp1), bpppm(llmp1) LOGICAL dum,fill DATA fill/.true./ DATA dum/.true./ integer,save :: countcfl=0 real cflx(ip1jmp1,llm) real cfly(ip1jm,llm) real cflz(ip1jmp1,llm) real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) IF(iadvtr.EQ.0) THEN CALL initial0(ijp1llm,pbaruc) CALL initial0(ijmllm,pbarvc) ENDIF ! accumulation des flux de masse horizontaux DO l=1,llm DO ij = 1,ip1jmp1 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) ENDDO DO ij = 1,ip1jm pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) ENDDO ENDDO ! selection de la masse instantannee des mailles avant le transport. IF(iadvtr.EQ.0) THEN CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) ! ENDIF iadvtr = iadvtr+1 iapptrac = iadvtr ! Test pour savoir si on advecte a ce pas de temps IF ( iadvtr.EQ.iapp_tracvl ) THEN !c .. Modif P.Le Van ( 20/12/97 ) .... !c ! traitement des flux de masse avant advection. ! 1. calcul de w ! 2. groupement des mailles pres du pole. CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) ! ... Flux de masse diaganostiques traceurs flxw = wg / REAL(iapp_tracvl) ! test sur l'eventuelle creation de valeurs negatives de la masse DO l=1,llm-1 DO ij = iip2+1,ip1jm zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & - pbarvg(ij-iip1,l) + pbarvg(ij,l) & + wg(ij,l+1) - wg(ij,l) ENDDO CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) DO ij = iip2,ip1jm zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) ENDDO CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, & ' MAX:', zdpmax ENDIF ENDDO !------------------------------------------------------------------- ! Calcul des criteres CFL en X, Y et Z !------------------------------------------------------------------- if (countcfl == 0. ) then cflxmax(:)=0. cflymax(:)=0. cflzmax(:)=0. endif countcfl=countcfl+iapp_tracvl cflx(:,:)=0. cfly(:,:)=0. cflz(:,:)=0. do l=1,llm do ij=iip2,ip1jm-1 if (pbarug(ij,l)>=0.) then cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) else cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) endif enddo enddo do l=1,llm do ij=iip2,ip1jm-1,iip1 cflx(ij+iip1,l)=cflx(ij,l) enddo enddo do l=1,llm do ij=1,ip1jm if (pbarvg(ij,l)>=0.) then cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) else cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) endif enddo enddo do l=2,llm do ij=1,ip1jm if (wg(ij,l)>=0.) then cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) else cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) endif enddo enddo do l=1,llm cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Par defaut, on sort le diagnostic des CFL tous les jours. ! Si on veut le sortir a chaque pas d'advection en cas de plantage ! if (countcfl==iapp_tracvl) then !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (countcfl==day_step) then do l=1,llm write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), & cflzmax(l) enddo countcfl=0 endif !------------------------------------------------------------------- ! Advection proprement dite (Modification Le Croller (07/2001) !------------------------------------------------------------------- !---------------------------------------------------- ! Calcul des moyennes basées sur la masse !---------------------------------------------------- call massbar(massem,massebx,masseby) !----------------------------------------------------------- ! Appel des sous programmes d'advection !----------------------------------------------------------- do iq=1,nqtot ! call clock(t_initial) if(iadv(iq) == 0) cycle ! ---------------------------------------------------------------- ! Schema de Van Leer I MUSCL ! ---------------------------------------------------------------- if(iadv(iq).eq.10) THEN call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) ! ---------------------------------------------------------------- ! Schema "pseudo amont" + test sur humidite specifique ! pour la vapeur d'eau. F. Codron ! ---------------------------------------------------------------- else if(iadv(iq).eq.14) then ! CALL vlspltqs( q(1,1,1), 2., massem, wg , & pbarug,pbarvg,dtvr,p,pk,teta ) ! ---------------------------------------------------------------- ! Schema de Frederic Hourdin ! ---------------------------------------------------------------- else if(iadv(iq).eq.12) then ! Pas de temps adaptatif call adaptdt(iadv(iq),dtbon,n,pbarug,massem) if (n.GT.1) then write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & dtvr,'n=',n endif do indice=1,n call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) end do else if(iadv(iq).eq.13) then ! Pas de temps adaptatif call adaptdt(iadv(iq),dtbon,n,pbarug,massem) if (n.GT.1) then write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & dtvr,'n=',n endif do indice=1,n call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) end do ! ---------------------------------------------------------------- ! Schema de pente SLOPES ! ---------------------------------------------------------------- else if (iadv(iq).eq.20) then call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) ! ---------------------------------------------------------------- ! Schema de Prather ! ---------------------------------------------------------------- else if (iadv(iq).eq.30) then ! Pas de temps adaptatif call adaptdt(iadv(iq),dtbon,n,pbarug,massem) if (n.GT.1) then write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & dtvr,'n=',n endif call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, & n,dtbon) ! ---------------------------------------------------------------- ! Schemas PPM Lin et Rood ! ---------------------------------------------------------------- else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. & iadv(iq).LE.18)) then ! Test sur le flux horizontal ! Pas de temps adaptatif call adaptdt(iadv(iq),dtbon,n,pbarug,massem) if (n.GT.1) then write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & dtvr,'n=',n endif ! Test sur le flux vertical CFLmaxz=0. do l=2,llm do ij=iip2,ip1jm aaa=wg(ij,l)*dtvr/massem(ij,l) CFLmaxz=max(CFLmaxz,aaa) bbb=-wg(ij,l)*dtvr/massem(ij,l-1) CFLmaxz=max(CFLmaxz,bbb) enddo enddo if (CFLmaxz.GE.1) then write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz endif !----------------------------------------------------------- ! Ss-prg interface LMDZ.4->PPM3d !----------------------------------------------------------- call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & apppm,bpppm,massebx,masseby,pbarug,pbarvg, & unatppm,vnatppm,psppm) do indice=1,n !---------------------------------------------------------------- ! VL (version PPM) horiz. et PPM vert. !---------------------------------------------------------------- if (iadv(iq).eq.11) then ! Ss-prg PPM3d de Lin call ppm3d(1,qppm(1,1,iq), & psppm,psppm, & unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, & iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & fill,dum,220.) !------------------------------------------------------------- ! Monotonic PPM !------------------------------------------------------------- else if (iadv(iq).eq.16) then ! Ss-prg PPM3d de Lin call ppm3d(1,qppm(1,1,iq), & psppm,psppm, & unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, & iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & fill,dum,220.) !------------------------------------------------------------- !------------------------------------------------------------- ! Semi Monotonic PPM !------------------------------------------------------------- else if (iadv(iq).eq.17) then ! Ss-prg PPM3d de Lin call ppm3d(1,qppm(1,1,iq), & psppm,psppm, & unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, & iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & fill,dum,220.) !------------------------------------------------------------- !------------------------------------------------------------- ! Positive Definite PPM !------------------------------------------------------------- else if (iadv(iq).eq.18) then ! Ss-prg PPM3d de Lin call ppm3d(1,qppm(1,1,iq), & psppm,psppm, & unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, & iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & fill,dum,220.) !------------------------------------------------------------- endif enddo !----------------------------------------------------------------- ! Ss-prg interface PPM3d-LMDZ.4 !----------------------------------------------------------------- call interpost(q(1,1,iq),qppm(1,1,iq)) endif !---------------------------------------------------------------------- !----------------------------------------------------------------- ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 ! et Nord j=1 !----------------------------------------------------------------- ! call traceurpole(q(1,1,iq),massem) ! calcul du temps cpu pour un schema donne ! call clock(t_final) !ym tps_cpu=t_final-t_initial !ym cpuadv(iq)=cpuadv(iq)+tps_cpu end DO !------------------------------------------------------------------ ! on reinitialise a zero les flux de masse cumules !--------------------------------------------------- iadvtr=0 ENDIF ! if iadvtr.EQ.iapp_tracvl END SUBROUTINE advtrac