!
! $Id: add_phys_tend.F90 2611 2016-08-03 15:41:26Z jyg $
!
SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql, &
        zdqi,paprs,text,abortphy,flag_inhib_tend &
#ifdef ISO
        ,zdxt,zdxtl,zdxti &
#endif
        )
!======================================================================
! Ajoute les tendances des variables physiques aux variables 
! d'etat de la dynamique t_seri, q_seri ...
! On en profite pour faire des tests sur les tendances en question.
!======================================================================


!======================================================================
! Declarations
!======================================================================

USE dimphy, ONLY: klon, klev
USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, &
#ifdef ISO
                              xtl_seri, xts_seri, xt_seri, &
#endif
                              t_seri
USE phys_state_var_mod, ONLY: ftsol
USE geometry_mod, ONLY: longitude_deg, latitude_deg
USE print_control_mod, ONLY: prt_level
USE cmp_seri_mod
#ifdef ISO
  USE infotrac_phy, ONLY: ntraciso ! ajout C Risi pour isos  
#ifdef ISOVERIF
  USE isotopes_mod, ONLY: iso_eau
  USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
        iso_verif_egalite_vect2D
#endif    
#endif

IMPLICIT none
  include "YOMCST.h"
  include "clesphys.h"

! Arguments :
!------------
REAL zdu(klon,klev),zdv(klon,klev)
REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev),zdqi(klon,klev)
#ifdef ISO
REAL zdxt(ntraciso,klon,klev),zdxtl(ntraciso,klon,klev) &
        ,zdxti(ntraciso,klon,klev)
integer ixt
#endif
REAL paprs(klon,klev+1)
CHARACTER*(*) text
INTEGER abortphy
INTEGER flag_inhib_tend ! if flag_inhib_tend != 0, tendencies are not added

! Local :
!--------
REAL zt,zq
#ifdef ISO
REAL zxt(ntraciso)
real qprec(klon,klev)
REAL zxt_int(ntraciso), zxtp_int(ntraciso), zxtp(ntraciso,klev),zxt_new
#endif
REAL zq_int, zqp_int, zq_new

REAL zqp(klev)

INTEGER i, k,j
INTEGER jadrs(klon*klev), jbad
INTEGER jqadrs(klon*klev), jqbad
INTEGER kadrs(klon*klev)
INTEGER kqadrs(klon*klev)

LOGICAL done(klon)

integer debug_level
logical, save :: first=.true.
!$OMP THREADPRIVATE(first)
INTEGER, SAVE :: itap
!$OMP THREADPRIVATE(itap)
!======================================================================
! Initialisations

     IF (prt_level >= 5) then
        write (*,*) "In add_phys_tend, after ",text
        call flush(0) ! C Risi: attention, on a rajouté (0)
     end if

     ! if flag_inhib_tend != 0, tendencies are not added
     IF (flag_inhib_tend /= 0) then
        ! If requiered, diagnostics are shown
        IF (flag_inhib_tend > 0) then
           ! print some diagnostics if xxx_seri have changed
           call cmp_seri(flag_inhib_tend,text)
        END IF
        RETURN ! on n ajoute pas les tendance
     END IF

     IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
                              ! a deja plante.

     debug_level=10
     if (first) then
        itap=0
        first=.false.
     endif
! Incrementer le compteur de la physique
     itap   = itap + 1
!======================================================================
! Ajout des tendances sur le vent et l'eau liquide
!======================================================================

     u_seri(:,:)=u_seri(:,:)+zdu(:,:)
     v_seri(:,:)=v_seri(:,:)+zdv(:,:)
     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
#ifdef ISO 
     xtl_seri(:,:,:)=xtl_seri(:,:,:)+zdxtl(:,:,:) 
     xts_seri(:,:,:)=xts_seri(:,:,:)+zdxti(:,:,:)  
#endif

!======================================================================
! On ajoute les tendances de la temperature et de la vapeur d'eau
! en verifiant que ca ne part pas dans les choux
!======================================================================

      jbad=0
      jqbad=0
      DO k = 1, klev
         DO i = 1, klon
            zt=t_seri(i,k)+zdt(i,k)
            zq=q_seri(i,k)+zdq(i,k)
#ifdef ISO
            do ixt=1,ntraciso    
              zxt(ixt)=xt_seri(ixt,i,k)+zdxt(ixt,i,k)
            enddo !do ixt=1,ntraciso 
            qprec(i,k)=q_seri(i,k)
#endif
            IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
            jbad = jbad + 1
            jadrs(jbad) = i
            kadrs(jbad) = k
            ENDIF
            IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
            jqbad = jqbad + 1
            jqadrs(jqbad) = i
            kqadrs(jqbad) = k
            ENDIF
            t_seri(i,k)=zt
            q_seri(i,k)=zq
#ifdef ISO            
            do ixt=1,ntraciso    
              xt_seri(ixt,i,k)=zxt(ixt)
            enddo !do ixt=1,ntraciso 
#endif
         ENDDO
      ENDDO

#ifdef ISO
#ifdef ISOVERIF
        if (iso_eau.gt.0) then
          call iso_verif_egalite_vect2D( &
                xt_seri,q_seri, &
                'add_phys_tend 122'//text,ntraciso,klon,klev)
          call iso_verif_egalite_vect2D( &
                xtl_seri,ql_seri, &
                'add_phys_tend 138'//text,ntraciso,klon,klev)
         endif !if (iso_eau.gt.0) then
#endif
#endif

!=====================================================================================
! Impression et stop en cas de probleme important
!=====================================================================================

IF (jbad .GT. 0) THEN
      DO j = 1, jbad
         i=jadrs(j)
         if(prt_level.ge.debug_level) THEN
          print*,'PLANTAGE POUR LE POINT i lon lat =',&
                 i,longitude_deg(i),latitude_deg(i),text
          print*,'l    T     dT       Q     dQ    '
          DO k = 1, klev
             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
          ENDDO
          call print_debug_phys(i,debug_level,text)
         endif
      ENDDO
ENDIF
!
!=====================================================================================
! Impression, warning et correction en cas de probleme moins important
!=====================================================================================
IF (jqbad .GT. 0) THEN
      done(:) = .false.                         !jyg
      DO j = 1, jqbad
        i=jqadrs(j)
          if(prt_level.ge.debug_level) THEN
           print*,'WARNING  : EAU POUR LE POINT i lon lat =',&
                  i,longitude_deg(i),latitude_deg(i),text
           print*,'l    T     dT       Q     dQ    '
           DO k = 1, klev
              write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
           ENDDO
          endif
          IF (ok_conserv_q) THEN
!jyg<20140228 Corrections pour conservation de l'eau
            IF (.NOT.done(i)) THEN                  !jyg
              DO k = 1, klev
                zqp(k) = max(q_seri(i,k),1.e-15)
#ifdef ISO
                do ixt=1,ntraciso
                  zxtp(ixt,k)=xt_seri(ixt,i,k)*max(1.0,1.e-15/qprec(i,k))
                enddo
#endif
              ENDDO
              zq_int  = 0.
              zqp_int = 0.
              DO k = 1, klev
                zq_int  = zq_int  + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
                zqp_int = zqp_int + zqp(k)     *(paprs(i,k)-paprs(i,k+1))/Rg
              ENDDO
#ifdef ISO
              do ixt=1,ntraciso
              zxt_int(ixt)  = 0.
              zxtp_int(ixt) = 0.
              DO k = 1, klev
                zxt_int(ixt)  = zxt_int(ixt)  + xt_seri(ixt,i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
                zxtp_int(ixt) = zxtp_int(ixt) + zxtp(ixt,k)     *(paprs(i,k)-paprs(i,k+1))/Rg
              ENDDO
              enddo
#endif

              if(prt_level.ge.debug_level) THEN
               print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
                                    i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
              endif
              DO k = 1, klev
                zq_new = zqp(k)*zq_int/zqp_int
                zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
                q_seri(i,k) = zq_new
#ifdef ISO
                do ixt=1,ntraciso
                  zxt_new = zxtp(ixt,k)*zxt_int(ixt)/zxtp_int(ixt)
                  zdxt(ixt,i,k) = zdxt(ixt,i,k) + zxt_new - xt_seri(ixt,i,k)
                  xt_seri(ixt,i,k) = zxt_new
                enddo !ixt=1,ntraciso
#endif
              ENDDO
              done(i) = .true.
            ENDIF !(.NOT.done(i))
          ELSE 
!jyg>
            DO k = 1, klev
              zq=q_seri(i,k)+zdq(i,k)
              if (zq.lt.1.e-15) then
                 if (q_seri(i,k).lt.1.e-15) then
                  if(prt_level.ge.debug_level) THEN
                   print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
                  endif
                  q_seri(i,k)=1.e-15
                  zdq(i,k)=(1.e-15-q_seri(i,k))
#ifdef ISO
               do ixt=1,ntraciso     
                xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
                  ! zdxt(ixt,i,k)=0.0
                zdxt(ixt,i,k)=xt_seri(ixt,i,k)*(1.e-15/qprec(i,k)-1)
               enddo !do ixt=1,ntraciso   
#endif
                 endif
              endif
!              zq=q_seri(i,k)+zdq(i,k)
!              if (zq.lt.1.e-15) then
!                 zdq(i,k)=(1.e-15-q_seri(i,k))
!              endif
            ENDDO
!jyg<20140228
          ENDIF   !  (ok_conserv_q)
!jyg>
      ENDDO ! j = 1, jqbad
ENDIF

#ifdef ISO
#ifdef ISOVERIF
     if (iso_eau.gt.0) then
        call iso_verif_egalite_vect2D( &
                xt_seri,q_seri, &
                        'add_phys_tend 173'//text,ntraciso,klon,klev)
        call iso_verif_egalite_vect2D( &
                zdxt,zdq, &
                        'add_phys_tend 175'//text,ntraciso,klon,klev)
      endif !if (iso_eau.gt.0) then
#endif
#endif
!

!IM ajout memes tests pour reverifier les jbad, jqbad beg
      jbad=0
      jqbad=0
      DO k = 1, klev
         DO i = 1, klon
            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
            jbad = jbad + 1
            jadrs(jbad) = i
!            if(prt_level.ge.debug_level) THEN
!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
!            endif
            ENDIF
            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
            jqbad = jqbad + 1
            jqadrs(jqbad) = i
            kqadrs(jqbad) = k
!            if(prt_level.ge.debug_level) THEN
!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
!            endif
            ENDIF
         ENDDO
      ENDDO
IF (jbad .GT. 0) THEN
      DO j = 1, jbad
         i=jadrs(j)
         k=kadrs(j)
         if(prt_level.ge.debug_level) THEN
          print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
                 i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN 
          print*,'l    T     dT       Q     dQ    '
          DO k = 1, klev
             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
          ENDDO
          call print_debug_phys(i,debug_level,text)
         endif
      ENDDO
ENDIF 
!
IF (jqbad .GT. 0) THEN
      DO j = 1, jqbad
         i=jqadrs(j)
         k=kqadrs(j)
         if(prt_level.ge.debug_level) THEN
          print*,'WARNING  : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
                 i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN 
          print*,'l    T     dT       Q     dQ    '
          DO k = 1, klev
            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
          ENDDO
          call print_debug_phys(i,debug_level,text)
         endif
      ENDDO
ENDIF


!======================================================================
! Contrôle des min/max pour arrêt du modèle
! Si le modele est en mode abortphy, on retire les tendances qu'on
! vient d'ajouter. Pas exactement parce qu'on ne tient pas compte des
! seuils.
!======================================================================

      CALL hgardfou(t_seri,ftsol,text,abortphy) 
      IF (abortphy==1) THEN
        Print*,'ERROR ABORT hgardfou dans ',text
        u_seri(:,:)=u_seri(:,:)-zdu(:,:)
        v_seri(:,:)=v_seri(:,:)-zdv(:,:)
        ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
        qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
      ENDIF


  RETURN
END SUBROUTINE add_phys_tend