MODULE advect_tracer_mod
  USE icosa
  PRIVATE
  INTEGER,PARAMETER::iapp_tracvl= 3
  REAL(rstd),SAVE :: dt

  TYPE(t_field),POINTER :: f_normal(:)
  TYPE(t_field),POINTER :: f_tangent(:)
  TYPE(t_field),POINTER :: f_gradq3d(:)

  PUBLIC init_advect_tracer, advect_tracer

CONTAINS

  SUBROUTINE init_advect_tracer(dt_in)
    USE advect_mod
    IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: dt_in
    REAL(rstd),POINTER :: tangent(:,:)
    REAL(rstd),POINTER :: normal(:,:)
    INTEGER :: ind

    dt=dt_in
    CALL allocate_field(f_normal,field_u,type_real,3)
    CALL allocate_field(f_tangent,field_u,type_real,3)
    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3)

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       normal=f_normal(ind)
       tangent=f_tangent(ind)
       CALL init_advect(normal,tangent)
    END DO

  END SUBROUTINE init_advect_tracer

  SUBROUTINE advect_tracer(f_ps,f_u,f_q)
    USE icosa
    USE advect_mod
    USE disvert_mod
    IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_q(:)
    REAL(rstd),POINTER :: q(:,:,:)
    REAL(rstd),POINTER :: u(:,:) 
    REAL(rstd),POINTER :: ps(:) 
    REAL(rstd),POINTER :: massflx(:,:)
    REAL(rstd),POINTER :: rhodz(:,:)
    TYPE(t_field),POINTER,SAVE :: f_massflxc(:)
    TYPE(t_field),POINTER,SAVE :: f_massflx(:)
    TYPE(t_field),POINTER,SAVE :: f_uc(:)
    TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:)
    TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
    REAL(rstd),POINTER,SAVE :: massflxc(:,:)
    REAL(rstd),POINTER,SAVE :: uc(:,:)
    REAL(rstd),POINTER,SAVE :: rhodzm1(:,:)
    REAL(rstd):: bigt
    INTEGER :: ind,it,i,j,l,ij
    INTEGER,SAVE :: iadvtr=0
    LOGICAL,SAVE:: first=.TRUE. 
!------------------------------------------------------sarvesh
       CALL transfert_request(f_ps,req_i1)
       CALL transfert_request(f_u,req_e1)
       CALL transfert_request(f_u,req_e1)
       CALL transfert_request(f_q,req_i1)
	
    IF ( first ) THEN 
       CALL allocate_field(f_rhodz,field_t,type_real,llm) 
       CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 
       CALL allocate_field(f_massflxc,field_u,type_real,llm) 
       CALL allocate_field(f_massflx,field_u,type_real,llm) 
       CALL allocate_field(f_uc,field_u,type_real,llm) 
       first = .FALSE.
    END IF

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       rhodz=f_rhodz(ind)
       massflx=f_massflx(ind)
       ps=f_ps(ind)
       u=f_u(ind)

       DO l = 1, llm
          DO j=jj_begin-1,jj_end+1
             DO i=ii_begin-1,ii_end+1
                ij=(j-1)*iim+i
                rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
             ENDDO
          ENDDO
       ENDDO

       DO l = 1, llm
          DO j=jj_begin-1,jj_end+1
             DO i=ii_begin-1,ii_end+1
                ij=(j-1)*iim+i
                massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
                massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
                massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    IF ( iadvtr == 0 ) THEN 
       DO ind=1,ndomain         
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
          rhodz=f_rhodz(ind)  
          rhodzm1 = f_rhodzm1(ind) 
          massflxc = f_massflxc(ind)  ! accumulated mass fluxes
          uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura)
          rhodzm1 = rhodz
          massflxc = 0.0 
          uc = 0.0 
       END DO
       CALL transfert_request(f_rhodzm1,req_i1)   !---->
       CALL transfert_request(f_massflxc,req_e1) !---->
       CALL transfert_request(f_massflxc,req_e1) !------>
       CALL transfert_request(f_uc,req_e1) !---->
       CALL transfert_request(f_uc,req_e1) 
    END IF

    iadvtr = iadvtr + 1 

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       massflx  = f_massflx(ind)
       massflxc = f_massflxc(ind) 
       uc = f_uc(ind) 
       u  = f_u(ind) 
       massflxc = massflxc + massflx 
       uc = uc + u 
    END DO

    IF ( iadvtr == iapp_tracvl ) THEN 
       PRINT *, 'Advection scheme'
       bigt = dt*iapp_tracvl
       DO ind=1,ndomain
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
          uc = f_uc(ind) 
          uc = uc/real(iapp_tracvl) 
          massflxc = f_massflxc(ind)
          massflxc = massflxc*dt
          ! now massflx is time-integrated
       END DO

       CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt) 
       iadvtr = 0
    END IF
  END SUBROUTINE advect_tracer

  SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)
    USE field_mod
    USE domain_mod
    USE dimensions
    USE grid_param
    USE geometry
    USE metric
    USE advect_mod
    IMPLICIT NONE

    TYPE(t_field),POINTER :: f_q(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_rhodz(:) 
    TYPE(t_field),POINTER :: f_massflx(:)

    TYPE(t_field),POINTER,SAVE :: f_wg(:)
    TYPE(t_field),POINTER,SAVE :: f_zm(:)
    TYPE(t_field),POINTER,SAVE :: f_zq(:)

    REAL(rstd)::bigt,dt
    REAL(rstd),POINTER :: q(:,:,:)
    REAL(rstd),POINTER :: u(:,:)
    REAL(rstd),POINTER :: rhodz(:,:) 
    REAL(rstd),POINTER :: massflx(:,:) 
    REAL(rstd),POINTER :: wg(:,:)
    REAL(rstd),POINTER :: zq(:,:,:) 
    REAL(rstd),POINTER :: zm(:,:) 
    REAL(rstd),POINTER :: tangent(:,:)
    REAL(rstd),POINTER :: normal(:,:)
    REAL(rstd),POINTER :: gradq3d(:,:,:)

    REAL(rstd)::pente_max
    LOGICAL,SAVE::first = .true. 
    INTEGER :: i,ij,l,j,ind,k
    REAL(rstd) :: zzpbar, zzw 
    REAL::qvmax,qvmin 

    IF ( first ) THEN 
       CALL allocate_field(f_wg,field_t,type_real,llm)
       CALL allocate_field(f_zm,field_t,type_real,llm)
       CALL allocate_field(f_zq,field_t,type_real,llm,nqtot) 
       first = .FALSE.
    END IF

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       q=f_q(ind) 
       rhodz=f_rhodz(ind) 
       zq=f_zq(ind)
       zm=f_zm(ind) 
       zm = rhodz ; zq = q  
       wg = f_wg(ind)
       wg = 0.0
       massflx=f_massflx(ind) 
       CALL advtrac(massflx,wg) ! compute vertical mass fluxes
    END DO

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       zq=f_zq(ind)
       zm=f_zm(ind) 
       wg=f_wg(ind)
       wg=wg*0.5
       DO k = 1, nqtot
          CALL vlz(zq(:,:,k),2.,zm,wg)
       END DO
    END DO

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       zq=f_zq(ind)
       zq = f_zq(ind)
       zm = f_zm(ind)
       massflx =f_massflx(ind)
       u = f_u(ind)
       tangent = f_tangent(ind)
       normal = f_normal(ind)
       gradq3d = f_gradq3d(ind)

       DO k = 1,nqtot
          CALL compute_gradq3d(zq(:,:,k),gradq3d)
          CALL compute_advect_horiz(tangent,normal,zq(:,:,k),gradq3d,zm,u,massflx,bigt) 
       END DO
    END DO

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       q = f_q(ind)
       zq = f_zq(ind) 
       zm = f_zm(ind) 
       wg = f_wg(ind)
       DO k = 1,nqtot
          CALL vlz(zq(:,:,k),2.,zm,wg)
       END DO
       q = zq 
    END DO

  END SUBROUTINE vlsplt

  !==============================================================================  
  SUBROUTINE advtrac(massflx,wgg)
    USE domain_mod
    USE dimensions
    USE grid_param
    USE geometry
    USE metric
    USE disvert_mod
    IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm)
    REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm)

    INTEGER :: i,j,ij,l
    REAL(rstd) :: convm(iim*jjm,llm) 

    DO l = 1, llm
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             ! divergence of horizontal flux
             convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  &
                  ne(ij,rup)*massflx(ij+u_rup,l)       +  &
                  ne(ij,lup)*massflx(ij+u_lup,l)       +  &
                  ne(ij,left)*massflx(ij+u_left,l)     +  &
                  ne(ij,ldown)*massflx(ij+u_ldown,l)   +  &
                  ne(ij,rdown)*massflx(ij+u_rdown,l))
          ENDDO
       ENDDO
    ENDDO

    ! accumulate divergence from top of model
    DO  l = llm-1, 1, -1
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             convm(ij,l) = convm(ij,l) + convm(ij,l+1)
          ENDDO
       ENDDO
    ENDDO

!!! Compute vertical velocity
    DO l = 1,llm-1
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )) 
          ENDDO
       ENDDO
    ENDDO

    DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          wgg(ij,1)  = 0.
       ENDDO
    ENDDO
  END SUBROUTINE advtrac

  SUBROUTINE vlz(q,pente_max,masse,wgg)
    !c
    !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
    !c
    !c    ********************************************************************
    !c     Shema  d'advection " pseudo amont " .
    !c    ********************************************************************
    USE icosa
    IMPLICIT NONE
    !c
    !c   Arguments:
    !c   ----------
    REAL masse(iim*jjm,llm),pente_max
    REAL q(iim*jjm,llm)
    REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1) 
    REAL dq(iim*jjm,llm)
    INTEGER i,ij,l,j,ii
    !c
    REAL wq(iim*jjm,llm+1),newmasse

    REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax
    REAL sigw

    REAL      SSUM


    w(:,1:llm) = -wgg(:,:)  ! w>0 for downward transport
    w(:,llm+1) = 0.0 

    !c    On oriente tout dans le sens de la pression c'est a dire dans le
    !c    sens de W

    DO l=2,llm
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
             adzqw(ij,l)=abs(dzqw(ij,l))
          ENDDO
       ENDDO
    ENDDO

    DO l=2,llm-1
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
             ELSE
                dzq(ij,l)=0.
             ENDIF
             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
          ENDDO
       ENDDO
    ENDDO

    DO l=2,llm-1
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             dzq(ij,1)=0.
             dzq(ij,llm)=0.
          ENDDO
       ENDDO
    ENDDO
    !c ---------------------------------------------------------------
    !c   .... calcul des termes d'advection verticale  .......
    !c ---------------------------------------------------------------

    !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq

    DO l = 1,llm-1
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             IF(w(ij,l+1).gt.0.) THEN
                ! upwind only if downward transport
                sigw=w(ij,l+1)/masse(ij,l+1)
                wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
             ELSE
                ! upwind only if upward transport
                sigw=w(ij,l+1)/masse(ij,l)
                wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
             ENDIF
          ENDDO
       ENDDO
    END DO

    DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          wq(ij,llm+1)=0.
          wq(ij,1)=0.
       ENDDO
    END DO

    DO l=1,llm
       DO j=jj_begin,jj_end
          DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
             ! masse -= dw/dz but w>0 <=> downward
             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l)) 
!             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>>
             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse    
!             masse(ij,l)=newmasse
          ENDDO
       ENDDO
    END DO
    RETURN
  END SUBROUTINE vlz

END MODULE advect_tracer_mod