! $Id: integrd_p.F 1299 2010-01-20 14:27:21Z fairhead $

SUBROUTINE integrd_loc &
        (  nq,vcovm1,ucovm1,tetam1,psm1,massem1, &
        dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
  USE parallel_lmdz
  USE control_mod
  USE lmdz_filtreg_p
  USE write_field_loc
  USE lmdz_write_field
  USE integrd_mod
  USE comconst_mod, ONLY: pi
  USE logic_mod, ONLY: leapf
  USE comvert_mod, ONLY: ap, bp
  USE temps_mod, ONLY: dt
  USE lmdz_strings, ONLY: int2str
  USE lmdz_iniprint, ONLY: lunout, prt_level
  USE lmdz_ssum_scopy, ONLY: ssum


  IMPLICIT NONE


  !=======================================================================
  !
  !   Auteur:  P. Le Van
  !   -------
  !
  !   objet:
  !   ------
  !
  !   Incrementation des tendances dynamiques
  !
  !=======================================================================
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------

  include "dimensions.h"
  include "paramet.h"
  include "comgeom.h"

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

  INTEGER,INTENT(IN) :: nq ! number of tracers to handle in this routine

  REAL,INTENT(INOUT) :: vcov(ijb_v:ije_v,llm) ! covariant meridional wind
  REAL,INTENT(INOUT) :: ucov(ijb_u:ije_u,llm) ! covariant zonal wind
  REAL,INTENT(INOUT) :: teta(ijb_u:ije_u,llm) ! potential temperature
  REAL,INTENT(INOUT) :: q(ijb_u:ije_u,llm,nq) ! advected tracers
  REAL,INTENT(INOUT) :: ps0(ijb_u:ije_u) ! surface pressure
  REAL,INTENT(INOUT) :: masse(ijb_u:ije_u,llm) ! atmospheric mass
  REAL,INTENT(INOUT) :: phis(ijb_u:ije_u) ! ground geopotential !!! unused
  ! values at previous time step
  REAL,INTENT(INOUT) :: vcovm1(ijb_v:ije_v,llm)
  REAL,INTENT(INOUT) :: ucovm1(ijb_u:ije_u,llm)
  REAL,INTENT(INOUT) :: tetam1(ijb_u:ije_u,llm)
  REAL,INTENT(INOUT) :: psm1(ijb_u:ije_u)
  REAL,INTENT(INOUT) :: massem1(ijb_u:ije_u,llm)
  ! the tendencies to add
  REAL,INTENT(INOUT) :: dv(ijb_v:ije_v,llm)
  REAL,INTENT(INOUT) :: du(ijb_u:ije_u,llm)
  REAL,INTENT(INOUT) :: dteta(ijb_u:ije_u,llm)
  REAL,INTENT(INOUT) :: dp(ijb_u:ije_u)
  REAL,INTENT(INOUT) :: dq(ijb_u:ije_u,llm,nq) !!! unused
   ! REAL,INTENT(INOUT) ::finvmaold(ijb_u:ije_u,llm) !!! unused

  !   Local:
  !   ------

  REAL :: vscr( ijb_v:ije_v ),uscr( ijb_u:ije_u )
  REAL :: hscr( ijb_u:ije_u ),pscr(ijb_u:ije_u)
  REAL :: massescr( ijb_u:ije_u,llm )
   ! REAL finvmasse(ijb_u:ije_u,llm)
  REAL :: tpn,tps,tppn(iim),tpps(iim)
  REAL :: qpn,qps,qppn(iim),qpps(iim)

  INTEGER :: l,ij,iq,i,j

  INTEGER :: ijb,ije,jjb,jje
  LOGICAL :: checksum
  LOGICAL,SAVE :: checksum_all=.TRUE.
  INTEGER :: stop_it
  INTEGER :: ierr

  !WRITE(*,*) 'integrd 88: entree, nq=',nq
  !-----------------------------------------------------------------------

!$OMP BARRIER
  IF (pole_nord) THEN
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO  l = 1,llm
      DO  ij = 1,iip1
       ucov(    ij    , l) = 0.
       uscr(     ij      ) = 0.
       ENDDO
    ENDDO
!$OMP END DO NOWAIT
  ENDIF

  IF (pole_sud) THEN
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO  l = 1,llm
      DO  ij = 1,iip1
       ucov( ij +ip1jm, l) = 0.
       uscr( ij +ip1jm   ) = 0.
      ENDDO
    ENDDO
!$OMP END DO NOWAIT
  ENDIF

  !    ............    integration  de       ps         ..............

   ! CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)

  ijb=ij_begin
  ije=ij_end
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO  l = 1,llm
    massescr(ijb:ije,l)=masse(ijb:ije,l)
  ENDDO
!$OMP END DO NOWAIT

!$OMP DO SCHEDULE(STATIC)
  DO ij = ijb,ije
   pscr (ij)    = ps0(ij)
   ps (ij)      = psm1(ij) + dt * dp(ij)

  END DO

!$OMP END DO
!$OMP BARRIER
  ! --> ici synchro OPENMP pour ps

  checksum=.TRUE.
  stop_it=0

!$OMP MASTER
  !c$OMP DO SCHEDULE(STATIC)
  DO ij = ijb,ije
     IF( ps(ij)<0. ) THEN
       IF (checksum) stop_it=ij
       checksum=.FALSE.
     ENDIF
   ENDDO
  !c$OMP END DO NOWAIT

   ! CALL MPI_ALLREDUCE(checksum,checksum_all,1,
  ! &                   MPI_LOGICAL,MPI_LOR,COMM_LMDZ,ierr)
  IF( .NOT. checksum ) THEN
     WRITE(lunout,*) "integrd: ps = ", ps(stop_it)
     WRITE(lunout,*) " at node ij =", stop_it
     ! since ij=j+(i-1)*jjp1 , we have
      j=modulo(stop_it,jjp1)
      i=1+(stop_it-j)/jjp1
      WRITE(lunout,*) " lon = ",rlonv(i)*180./pi, " deg", &
            " lat = ",rlatu(j)*180./pi, " deg"
     CALL abort_gcm("integrd_loc", "negative surface pressure", 1)
  ENDIF

!$OMP END MASTER
!$OMP BARRIER
    !WRITE(*,*) 'integrd 170'
  IF (.NOT. Checksum_all) THEN
    CALL WriteField_v('int_vcov',vcov)
    CALL WriteField_u('int_ucov',ucov)
    CALL WriteField_u('int_teta',teta)
    CALL WriteField_u('int_ps0',ps0)
    CALL WriteField_u('int_masse',masse)
    CALL WriteField_u('int_phis',phis)
    CALL WriteField_v('int_vcovm1',vcovm1)
    CALL WriteField_u('int_ucovm1',ucovm1)
    CALL WriteField_u('int_tetam1',tetam1)
    CALL WriteField_u('int_psm1',psm1)
    CALL WriteField_u('int_massem1',massem1)

    CALL WriteField_v('int_dv',dv)
    CALL WriteField_u('int_du',du)
    CALL WriteField_u('int_dteta',dteta)
    CALL WriteField_u('int_dp',dp)
     ! CALL WriteField_u('int_finvmaold',finvmaold)
    do j=1,nq
      CALL WriteField_u('int_q'//trim(int2str(j)), &
            q(:,:,j))
      CALL WriteField_u('int_dq'//trim(int2str(j)), &
            dq(:,:,j))
    enddo
    CALL abort_gcm("integrd_loc", "", 1)
  ENDIF


  !
  !   !WRITE(*,*) 'integrd 200'
!$OMP MASTER
  IF (pole_nord) THEN

    DO  ij    = 1, iim
     tppn(ij) = aire(   ij   ) * ps(  ij    )
    ENDDO
     tpn      = SSUM(iim,tppn,1)/apoln
    DO ij   = 1, iip1
     ps(   ij   )  = tpn
    ENDDO

  ENDIF

  IF (pole_sud) THEN

    DO  ij    = 1, iim
     tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
    ENDDO
     tps      = SSUM(iim,tpps,1)/apols
    DO ij   = 1, iip1
     ps(ij+ip1jm)  = tps
    ENDDO

  ENDIF
!$OMP END MASTER
!$OMP BARRIER
  !WRITE(*,*) 'integrd 217'
  !
  !  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
  !

  CALL pression_loc ( ip1jmp1, ap, bp, ps, p )

!$OMP BARRIER
  CALL massdair_loc (     p  , masse         )

  ! Ehouarn : we don't use/need finvmaold and finvmasse,
        ! so might as well not compute them
  !c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
   ! ijb=ij_begin
   ! ije=ij_end

  !c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
   ! DO  l = 1,llm
   !   finvmasse(ijb:ije,l)=masse(ijb:ije,l)
   ! ENDDO
  !c$OMP END DO NOWAIT

   ! jjb=jj_begin
   ! jje=jj_end
   ! CALL filtreg_p( finvmasse,jjb_u,jje_u,jjb,jje, jjp1, llm,
  ! &                -2, 2, .TRUE., 1  )
  !

  !    ............   integration  de  ucov, vcov,  h     ..............

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

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

  DO ij = ijb,ije
  uscr( ij )   =  ucov( ij,l )
  ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
  END DO

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

  DO ij = ijb,ije
  vscr( ij )   =  vcov( ij,l )
  vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
  END DO

  ijb=ij_begin
  ije=ij_end

  DO ij = ijb,ije
  hscr( ij )    =  teta(ij,l)
  teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l) &
        + dt * dteta(ij,l) / masse(ij,l)
  END DO

  !   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
  !
  !
  !   !WRITE(*,*) 'integrd 291'
  IF (pole_nord) THEN

    DO  ij   = 1, iim
      tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
    ENDDO
      tpn      = SSUM(iim,tppn,1)/apoln

    DO ij   = 1, iip1
      teta(   ij   ,l)  = tpn
    ENDDO

  ENDIF

  IF (pole_sud) THEN

    DO  ij   = 1, iim
      tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
    ENDDO
      tps      = SSUM(iim,tpps,1)/apols

    DO ij   = 1, iip1
      teta(ij+ip1jm,l)  = tps
    ENDDO

  ENDIF
  !

  IF(leapf)  THEN
      ! CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
      ! CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
      ! CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
    ijb=ij_begin
    ije=ij_end
    ucovm1(ijb:ije,l)=uscr(ijb:ije)
    tetam1(ijb:ije,l)=hscr(ijb:ije)
    IF (pole_sud) ije=ij_end-iip1
    vcovm1(ijb:ije,l)=vscr(ijb:ije)

  END IF

  END DO
!$OMP END DO NOWAIT

  !
  !   .......  integration de   q   ......
  !
  ijb=ij_begin
  ije=ij_end

     IF (planet_type=="earth") THEN
  ! Earth-specific treatment of first 2 tracers (water)
!$OMP BARRIER
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
       DO ij = ijb, ije
        deltap(ij,l) =  p(ij,l) - p(ij,l+1)
       ENDDO
      ENDDO

!$OMP END DO NOWAIT
!$OMP BARRIER

    CALL check_isotopes(q,ijb,ije,'integrd 342')

    !WRITE(*,*) 'integrd 341'
    CALL qminimum_loc( q, nq, deltap )
    !WRITE(*,*) 'integrd 343'

    CALL check_isotopes(q,ijb,ije,'integrd 346')
  !
  !    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
  !
!$OMP BARRIER
  IF (pole_nord) THEN

    DO iq = 1, nq

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

         DO ij = 1, iim
           qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
         ENDDO
           qpn  =  SSUM(iim,qppn,1)/apoln

         DO ij = 1, iip1
           q(   ij   ,l,iq)  = qpn
         ENDDO

      ENDDO
!$OMP END DO NOWAIT

    ENDDO

  ENDIF

  IF (pole_sud) THEN

    DO iq = 1, nq

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

         DO ij = 1, iim
           qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
         ENDDO
           qps  =  SSUM(iim,qpps,1)/apols

         DO ij = 1, iip1
           q(ij+ip1jm,l,iq)  = qps
         ENDDO

      ENDDO
!$OMP END DO NOWAIT

    ENDDO

  ENDIF

  CALL check_isotopes(q,ijb,ije,'integrd 409')

  ! Ehouarn: forget about finvmaold
  !c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )

  !c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
   ! DO l = 1, llm
   !   finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)
   ! ENDDO
  !c$OMP END DO NOWAIT

  ENDIF ! of if (planet_type.EQ."earth")

  !
  !
  ! .....   FIN  de l'integration  de   q    .......

      !WRITE(*,*) 'integrd 410'

!$OMP DO SCHEDULE(STATIC)
  DO ij=ijb,ije
    ps0(ij)=ps(ij)
  ENDDO
!$OMP END DO NOWAIT

  !    .................................................................


  IF( leapf )  THEN
    ! CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
    ! CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
!$OMP DO SCHEDULE(STATIC)
  DO ij=ijb,ije
    psm1(ij)=pscr(ij)
  ENDDO
!$OMP END DO NOWAIT

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        massem1(ijb:ije,l)=massescr(ijb:ije,l)
      ENDDO
!$OMP END DO NOWAIT
  END IF
!$OMP BARRIER

END SUBROUTINE integrd_loc