advtrac_loc.F

Go to the documentation of this file.

!
! $Id$
!
c
c
#define DEBUG_IO
#undef DEBUG_IO
      SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
     *                   p,  massem,q,teta,
     *                   pk   )

c     Auteur :  F. Hourdin
c
c     Modif. P. Le Van     (20/12/97)
c            F. Codron     (10/99)
c            D. Le Croller (07/2001)
c            M.A Filiberti (04/2002)
c
      USE parallel_lmdz
      USE write_field_loc
      USE write_field
      USE bands
      USE mod_hallo
      USE vampir
      USE times
      USE infotrac, ONLY: nqtot, iadv, ok_iso_verif
      USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
      USE advtrac_mod, ONLY: finmasse
      IMPLICIT NONE
c
#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"

c-------------------------------------------------------------------
c     Arguments
c-------------------------------------------------------------------
c     Ajout PPM
c--------------------------------------------------------
      REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
c--------------------------------------------------------
      INTEGER iapptrac
      REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
      REAL wg(ijb_u:ije_u,llm)
      REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
      REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
      REAL pk(ijb_u:ije_u,llm)

c-------------------------------------------------------------
c     Variables locales
c-------------------------------------------------------------

      REAL zdp(ijb_u:ije_u)
      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
      INTEGER,SAVE :: iadvtr=0
c$OMP THREADPRIVATE(iadvtr)
      INTEGER ij,l,iq,iiq
      REAL zdpmin, zdpmax
c----------------------------------------------------------
c     Rajouts pour PPM
c----------------------------------------------------------
      INTEGER indice,n
      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
      REAL CFLmaxz,aaa,bbb ! CFL maximum
      REAL psppm(iim,jjb_u:jje_u) ! pression  au sol
      REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
      REAL qppm(iim*jjnb_u,llm,nqtot)
      REAL fluxwppm(iim,jjb_u:jje_u,llm)
      REAL apppm(llmp1), bpppm(llmp1)
      LOGICAL dum,fill
      DATA fill/.true./
      DATA dum/.true./
      integer ijb,ije,ijbu,ijbv,ijeu,ijev,j
      type(request),SAVE :: testRequest
!$OMP THREADPRIVATE(testRequest)

c  test sur l''eventuelle creation de valeurs negatives de la masse
         ijb=ij_begin
         ije=ij_end
         if (pole_nord) ijb=ij_begin+iip1
         if (pole_sud) ije=ij_end-iip1
         
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
         DO l=1,llm-1
            DO ij = ijb+1,ije
              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
     s                  +       wg(ij,l+1)  - wg(ij,l)
            ENDDO
            
c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
            
            do ij=ijb,ije-iip1+1,iip1
              zdp(ij)=zdp(ij+iip1-1)
            enddo
            
            DO ij = ijb,ije
               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
            ENDDO 


c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
c  ym ---> eventuellement a revoir
            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
            
            IF(max(abs(zdpmin),abs(zdpmax)).GT.0.5) THEN
            print*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
     s        '   MAX:', zdpmax
            ENDIF

         ENDDO
c$OMP END DO NOWAIT

c-------------------------------------------------------------------
c   Advection proprement dite (Modification Le Croller (07/2001)
c-------------------------------------------------------------------

c----------------------------------------------------
c        Calcul des moyennes basées sur la masse
c----------------------------------------------------

cym      ----> Normalement, inutile pour les schémas classiques
cym      ----> Revérifier lors de la parallélisation des autres schemas
   
cym          call massbar_p(massem,massebx,masseby)          

#ifdef DEBUG_IO    
          CALL writefield_u('massem',massem)
          CALL writefield_u('wg',wg)
          CALL writefield_u('pbarug',pbarug)
          CALL writefield_v('pbarvg',pbarvg)
          CALL writefield_u('p_tmp',p)
          CALL writefield_u('pk_tmp',pk)
          CALL writefield_u('teta_tmp',teta)
          do j=1,nqtot
            call writefield_u('q_adv'//trim(int2str(j)),
     .                q(:,:,j))
          enddo
#endif

!          
!       call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
!
!       call SendRequest(TestRequest)
!c$OMP BARRIER
!       call WaitRequest(TestRequest)
c$OMP BARRIER
                 
          !write(*,*) 'advtrac 157: appel de vlspltgen_loc'
          call vlspltgen_loc( q,iadv, 2., massem, wg ,
     *                        pbarug,pbarvg,dtvr,p,
     *                        pk,teta )

          !write(*,*) 'advtrac 162: apres appel vlspltgen_loc'
      if (ok_iso_verif) then
           call check_isotopes(q,ijb_u,ije_u,'advtrac 162')
      endif !if (ok_iso_verif) then

#ifdef DEBUG_IO     
          do j=1,nqtot
            call writefield_u('q_adv'//trim(int2str(j)),
     .                q(:,:,j))
          enddo
#endif
         
                  GOTO 1234     
c-----------------------------------------------------------
c     Appel des sous programmes d'advection
c-----------------------------------------------------------
      do iq=1,nqtot
c        call clock(t_initial)
        if(iadv(iq) == 0) cycle 
c   ----------------------------------------------------------------
c   Schema de Van Leer I MUSCL
c   ----------------------------------------------------------------
        if(iadv(iq).eq.10) THEN
     
!LF                 call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)

c   ----------------------------------------------------------------
c   Schema "pseudo amont" + test sur humidite specifique
C    pour la vapeur d'eau. F. Codron
c   ----------------------------------------------------------------
        else if(iadv(iq).eq.14) then
c
cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
!LF           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
!LF     *                 pbarug,pbarvg,dtvr,p,pk,teta )
c   ----------------------------------------------------------------
c   Schema de Frederic Hourdin
c   ----------------------------------------------------------------
        else if(iadv(iq).eq.12) then
          stop 'advtrac : schema non parallelise'
c            Pas de temps adaptatif
           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
           if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             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
          stop 'advtrac : schema non parallelise'
c            Pas de temps adaptatif
           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
           if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             dtvr,'n=',n
           endif
          do indice=1,n
            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
          end do
c   ----------------------------------------------------------------
c   Schema de pente SLOPES
c   ----------------------------------------------------------------
        else if (iadv(iq).eq.20) then
          stop 'advtrac : schema non parallelise'

            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)

c   ----------------------------------------------------------------
c   Schema de Prather
c   ----------------------------------------------------------------
        else if (iadv(iq).eq.30) then
          stop 'advtrac : schema non parallelise'
c            Pas de temps adaptatif
           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
           if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             dtvr,'n=',n
           endif
           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
     s                     n,dtbon)
c   ----------------------------------------------------------------
c   Schemas PPM Lin et Rood
c   ----------------------------------------------------------------
       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
     s                     iadv(iq).LE.18)) then

           stop 'advtrac : schema non parallelise'

c        Test sur le flux horizontal
c        Pas de temps adaptatif
         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
         if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             dtvr,'n=',n
         endif
c        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

c-----------------------------------------------------------
c        Ss-prg interface LMDZ.4->PPM3d
c-----------------------------------------------------------

          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
     s                 unatppm,vnatppm,psppm)

          do indice=1,n
c---------------------------------------------------------------------
c                         VL (version PPM) horiz. et PPM vert.
c---------------------------------------------------------------------
                if (iadv(iq).eq.11) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)

c----------------------------------------------------------------------
c                           Monotonic PPM
c----------------------------------------------------------------------
               else if (iadv(iq).eq.16) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c                           Semi Monotonic PPM
c---------------------------------------------------------------------
               else if (iadv(iq).eq.17) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c                         Positive Definite PPM
c---------------------------------------------------------------------
                else if (iadv(iq).eq.18) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)
c---------------------------------------------------------------------
                endif
            enddo
c-----------------------------------------------------------------
c               Ss-prg interface PPM3d-LMDZ.4
c-----------------------------------------------------------------
                  call interpost(q(1,1,iq),qppm(1,1,iq))
            endif
c----------------------------------------------------------------------

c-----------------------------------------------------------------
c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
c et Nord j=1
c-----------------------------------------------------------------

c                  call traceurpole(q(1,1,iq),massem)

c calcul du temps cpu pour un schema donne


       end DO

1234  CONTINUE
c$OMP BARRIER

      if (planet_type=="earth") then

      ijb=ij_begin
      ije=ij_end

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
       DO l = 1, llm
         DO ij = ijb, ije
           finmasse(ij,l) =  p(ij,l) - p(ij,l+1) 
         ENDDO
       ENDDO
c$OMP END DO

        ! CRisi: on passe nqtot et non nq
       CALL qminimum_loc( q, nqtot, finmasse )

      endif ! of if (planet_type=="earth")

       RETURN
       END