source: LMDZ6/trunk/libf/dyn3dmem/advect_new_loc.f90 @ 5308

!
! $Header$
!
SUBROUTINE advect_new_loc(ucov,vcov,teta,w,massebx,masseby, &
        du,dv,dteta)
  USE comgeom_mod_h
  USE parallel_lmdz
  USE write_field_loc
  USE advect_new_mod
  USE comconst_mod, ONLY: daysec
  USE logic_mod, ONLY: conser
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
  USE paramet_mod_h
IMPLICIT NONE
  !=======================================================================
  !
  !   Auteurs:  P. Le Van , Fr. Hourdin  .
  !   -------
  !
  !   Objet:
  !   ------
  !
  !   *************************************************************
  !   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
  !   *************************************************************
  !    ces termes sont ajoutes a du,dv,dteta et dq .
  !  Modif F.Forget 03/94 : on retire q de advect
  !
  !=======================================================================
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------

  !   Arguments:
  !   ----------

  REAL :: vcov(ijb_v:ije_v,llm),ucov(ijb_u:ije_u,llm)
  REAL :: teta(ijb_u:ije_u,llm)
  REAL :: massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
  REAL :: w(ijb_u:ije_u,llm)
  REAL :: dv(ijb_v:ije_v,llm),du(ijb_u:ije_u,llm)
  REAL :: dteta(ijb_u:ije_u,llm)
  !   Local:
  !   ------

  REAL :: wsur2(ijb_u:ije_u)
  REAL :: unsaire2(ijb_u:ije_u), ge(ijb_u:ije_u)
  REAL :: deuxjour, ww, gt, uu, vv

  INTEGER :: ij,l,ijb,ije
  EXTERNAL  SSUM
  REAL :: SSUM



  !-----------------------------------------------------------------------
  !   2. Calculs preliminaires:
  !   -------------------------

  IF (conser.AND.1==0)  THEN
     deuxjour = 2. * daysec

     DO  ij   = 1, ip1jmp1
     unsaire2(ij) = unsaire(ij) * unsaire(ij)
     END DO
  END IF


  !------------------  -yy ----------------------------------------------
  !   .  Calcul de     u

!$OMP MASTER
  ijb=ij_begin
  ije=ij_end
  if (pole_nord) ijb=ijb+iip1
  if (pole_sud)  ije=ije-iip1

  DO ij=ijb,ije
    du2(ij,1)=0.
    du1(ij,llm)=0.
  ENDDO

  ijb=ij_begin
  ije=ij_end
  if (pole_sud)  ije=ij_end-iip1

  DO ij=ijb,ije
    dv2(ij,1)=0.
    dv1(ij,llm)=0.
  ENDDO

  ijb=ij_begin
  ije=ij_end

  DO ij=ijb,ije
    dteta2(ij,1)=0.
    dteta1(ij,llm)=0.
  ENDDO
!$OMP END MASTER
!$OMP BARRIER

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO  l=1,llm

     ijb=ij_begin
     ije=ij_end
     if (pole_nord) ijb=ijb+iip1
     if (pole_sud)  ije=ije-iip1

      ! DO    ij     = iip2, ip1jmp1
      !    uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
      ! ENDDO

      ! DO    ij     = iip2, ip1jm
      !    uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
      ! ENDDO

     DO    ij     = ijb, ije

       uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l)) &
             +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
     ENDDO

     if (pole_nord) then
       DO      ij         = 1, iip1
          uav(ij      ,l) = 0.
       ENDDO
     endif

     if (pole_sud) then
       DO      ij         = 1, iip1
          uav(ip1jm+ij,l) = 0.
       ENDDO
     endif

  ENDDO
!$OMP END DO
   ! call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))

  !------------------  -xx ----------------------------------------------
  !   .  Calcul de     v

  ijb=ij_begin
  ije=ij_end
  if (pole_sud)  ije=ij_end-iip1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO  l=1,llm

     DO    ij   = ijb+1, ije
       vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
     ENDDO

     DO    ij   = ijb,ije,iip1
      vav(ij,l) = vav(ij+iim,l)
     ENDDO


     DO    ij   = ijb, ije-1
      vav(ij,l) = vav(ij,l) + vav(ij+1,l)
     ENDDO

     DO    ij       = ijb, ije, iip1
      vav(ij+iim,l) = vav(ij,l)
     ENDDO

  ENDDO
!$OMP END DO
    ! call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))

  !-----------------------------------------------------------------------
!$OMP BARRIER

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llmm1


    ! ......   calcul de  - w/2.    au niveau  l+1   .......
  ijb=ij_begin
  ije=ij_end+iip1
  if (pole_sud)  ije=ij_end

  DO   ij   = ijb, ije
  wsur2( ij ) = - 0.5 * w( ij,l+1 )
  END DO


  ! .....................     calcul pour  du     ..................

  ijb=ij_begin
  ije=ij_end
  if (pole_nord) ijb=ijb+iip1
  if (pole_sud)  ije=ije-iip1

  DO ij = ijb ,ije-1
  ww        = wsur2 (  ij  )     + wsur2( ij+1 )
  uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
  du1(ij,l)  =  ww * ( uu - uav(ij, l ) )/massebx(ij, l )
  du2(ij,l+1)=  ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
  END DO

  ! .................    calcul pour   dv      .....................
  ijb=ij_begin
  ije=ij_end
  if (pole_sud)  ije=ij_end-iip1

  DO ij = ijb, ije
  ww        = wsur2( ij+iip1 )   + wsur2( ij )
  vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
  dv1(ij,l)  =  ww * (vv - vav(ij, l ) )/masseby(ij, l )
  dv2(ij,l+1)=  ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
  END DO

  !

  ! ............................................................
  ! ...............    calcul pour   dh      ...................
  ! ............................................................

  !                   ---z
  !   calcul de  - d( teta  * w )      qu'on ajoute a   dh
  !               ...............
    ijb=ij_begin
    ije=ij_end

    DO ij = ijb, ije
     ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
     dteta1(ij, l ) =   ww
     dteta2(ij,l+1) =   ww
    END DO

  ! ym ---> conser a voir plus tard

   ! IF( conser)  THEN
  !
  !    DO 17 ij = 1,ip1jmp1
  !    ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
  !  17    CONTINUE
  !    gt       = SSUM( ip1jmp1,ge,1 )
  !    gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
  !  END IF

  END DO
!$OMP END DO

  ijb=ij_begin
  ije=ij_end
  if (pole_nord) ijb=ijb+iip1
  if (pole_sud)  ije=ije-iip1
IF (CPPKEY_DEBUGIO) THEN
   CALL WriteField_u('du_bis',du)
END IF
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    DO ij=ijb,ije-1
      du(ij,l)=du(ij,l)+du2(ij,l)-du1(ij,l)
    ENDDO

    DO   ij   = ijb+iip1-1, ije, iip1
     du( ij, l  ) = du( ij -iim, l  )
    ENDDO
  ENDDO
!$OMP END DO NOWAIT
IF (CPPKEY_DEBUGIO) THEN
  CALL WriteField_u('du1',du1)
  CALL WriteField_u('du2',du2)
  CALL WriteField_u('du_bis',du)
END IF
  ijb=ij_begin
  ije=ij_end
  if (pole_sud)  ije=ij_end-iip1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    DO ij=ijb,ije
      dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l)
    ENDDO
  ENDDO
!$OMP END DO NOWAIT
  ijb=ij_begin
  ije=ij_end

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    DO ij=ijb,ije
      dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l)
    ENDDO
  ENDDO
!$OMP END DO NOWAIT

  RETURN
END SUBROUTINE advect_new_loc