MODULE flusv_mod
      
      IMPLICIT NONE
      
      CONTAINS
      
      SUBROUTINE flusv(KDLON,nsf,n,omega,g,tau,emis,bh,bsol,fah,fdh)
        use dimradmars_mod, only: ndlo2, ndlon, nflev
        IMPLICIT NONE
c.......................................................................
c
c  compute the upward and downward fluxes at the interface between n layers
c  * in the infrared
c  * B is a linear function of $\tau$ in each layer
c  * B at the surface can be different than what corresponds to the profile
c    in the n-th layer
c  * the work hypothes isthat we have two isotropic fluxes for each 
c    hemisphere ("hemispheric constant") + "technical source function"
c    (see Toon et al. 1988)
c  * the downward flux at the top of the atmosphere is zero
c  * layers are numbered from the top of the atmosphere to the ground
c
c  in :   * KDLON      ---> vectorisation dimension
c         * nsf        ---> nsf=0 ==> "hemispheric constant"
c                           nsf>0 ==> "hemispheric constant" + "source function"
c         * n          ---> number of layers
c         * omega(i)   ---> single scattering albedo for the i-th layer
c         * g(i)       ---> asymmetry parameter for the i-th layer
c         * tau(i)     ---> optical thickness of the i-th layer
c         * emis       ---> ground emissivity
c         * bh(i)      ---> black body luminance at the top of the i-th layer,
c                           bh(n+1) for the ground value which
c                           corresponds to the profile for the n-th layer
c         * bsol       ---> black body luminance of the ground
c
c  out :  * fah(i)     ---> upward flux at the top of the i-th layer,
c                           fah(n+1) for the ground
c         * fdh(i)     ---> downward flux at the top of the i-th layer,
c                           fdh(n+1) for the ground 
c
c.......................................................................
c  arguments 
c
      INTEGER,INTENT(IN) :: KDLON,nsf,n
      REAL,INTENT(IN) :: omega(NDLO2,n),g(NDLO2,n)
      REAL,INTENT(IN) :: tau(NDLO2,n),emis(NDLO2)
      REAL,INTENT(IN) :: bh(NDLO2,n+1),bsol(NDLO2)
      REAL,INTENT(OUT) :: fah(NDLO2,n+1),fdh(NDLO2,n+1)
c.......................................................................
c  local variables
c 
      REAL,PARAMETER :: pi=3.141592653589793E+0
      INTEGER iv,i,j
      REAL beta,gama1,gama2,amu1,grgama,b0,b1
      REAL a(NDLON,4*nflev),b(NDLON,4*nflev)
     &    ,d(NDLON,4*nflev),e(NDLON,4*nflev)
     &    ,y(NDLON,4*nflev)
     &    ,alambda(NDLON,2*nflev)
     &    ,e1(NDLON,2*nflev),e2(NDLON,2*nflev)
     &    ,e3(NDLON,2*nflev),e4(NDLON,2*nflev)
     &    ,cah(NDLON,2*nflev),cab(NDLON,2*nflev)
     &    ,cdh(NDLON,2*nflev),cdb(NDLON,2*nflev)
      REAL grg(NDLON,2*nflev),grh(NDLON,2*nflev)
     &    ,grj(NDLON,2*nflev),grk(NDLON,2*nflev)
     &    ,alpha1(NDLON,2*nflev),alpha2(NDLON,2*nflev)
     &    ,sigma1(NDLON,2*nflev),sigma2(NDLON,2*nflev)
      INTEGER,PARAMETER :: nq=8
      REAL,PARAMETER :: x(nq) = 
     &     (/1.9855071751231860E-2 , 0.1016667612931866E+0
     &     , 0.2372337950418355E+0 , 0.4082826787521751E+0
     &     , 0.5917173212478250E+0 , 0.7627662049581645E+0
     &     , 0.8983332387068134E+0 , 0.9801449282487682E+0/)
      REAL,PARAMETER :: w(nq) = 
     &     (/5.0614268145185310E-2 , 0.1111905172266872E+0
     &     , 0.1568533229389437E+0 , 0.1813418916891810E+0
     &     , 0.1813418916891810E+0 , 0.1568533229389437E+0
     &     , 0.1111905172266872E+0 , 5.0614268145185310E-2/)
      REAL :: gri(NDLON,nq)
c....................................................................... 
c
c....................................................................... 
      do i=1,n
        do iv=1,KDLON
      beta=(1.E+0-g(iv,i))/2.E+0
      gama1=2.E+0*(1.E+0-omega(iv,i)*(1.E+0-beta))
      gama2=2.E+0*omega(iv,i)*beta
      amu1=5.E-1
      alambda(iv,i)=sqrt(gama1**2-gama2**2)
      grgama=(gama1-alambda(iv,i))/gama2
c
c small hack here : if the optical depth of a layer is too small,
c $dB \over d\tau$ becomes very large and the scheme fails.
c In those cases we assume an isothermal layer.
c
      if (tau(iv,i).gt.1.E-3) then
        b0=bh(iv,i)
        b1=(bh(iv,i+1)-b0)/tau(iv,i)
      else
        b0=(bh(iv,i)+bh(iv,i+1))/2.E+0
        b1=0.E+0
      endif
c
      e1(iv,i)=1.E+0+grgama*exp(-alambda(iv,i)*tau(iv,i))
      e2(iv,i)=1.E+0-grgama*exp(-alambda(iv,i)*tau(iv,i))
      e3(iv,i)=grgama+exp(-alambda(iv,i)*tau(iv,i))
      e4(iv,i)=grgama-exp(-alambda(iv,i)*tau(iv,i))
      cah(iv,i)=2.E+0*pi*amu1*(b0+b1/(gama1+gama2))
      cab(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)+1.E+0/(gama1+gama2)))
      cdh(iv,i)=2.E+0*pi*amu1*(b0-b1/(gama1+gama2))
      cdb(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)-1.E+0/(gama1+gama2)))  
c
      grg(iv,i)=(1.E+0/amu1-alambda(iv,i))
      grh(iv,i)=grgama*(alambda(iv,i)+1.E+0/amu1)
      grj(iv,i)=grh(iv,i)
      grk(iv,i)=grg(iv,i)
      alpha1(iv,i)=2.E+0*pi*(b0+b1*(1.E+0/(gama1+gama2)-amu1))
      alpha2(iv,i)=2.E+0*pi*b1
      sigma1(iv,i)=2.E+0*pi*(b0-b1*(1.E+0/(gama1+gama2)-amu1))
      sigma2(iv,i)=alpha2(iv,i)
c
        enddo ! of do iv=1,KDLON
      enddo ! of do i=1,n
c.......................................................................
      do iv=1,KDLON
        a(iv,1)=0.E+0
        b(iv,1)=e1(iv,1)
        d(iv,1)=-e2(iv,1)
        e(iv,1)=-cdh(iv,1)
      enddo
c
      do i=1,n-1 
        j=2*i+1
        do iv=1,KDLON
          a(iv,j)=e2(iv,i)*e3(iv,i)-e4(iv,i)*e1(iv,i)
          b(iv,j)=e1(iv,i)*e1(iv,i+1)-e3(iv,i)*e3(iv,i+1)
          d(iv,j)=e3(iv,i)*e4(iv,i+1)-e1(iv,i)*e2(iv,i+1)
          e(iv,j)=e3(iv,i)*(cah(iv,i+1)-cab(iv,i))
     &           +e1(iv,i)*(cdb(iv,i)-cdh(iv,i+1))
        enddo
      enddo ! of do i=1,n-1
c
      do i=1,n-1
        j=2*i
        do iv=1,KDLON
          a(iv,j)=e2(iv,i+1)*e1(iv,i)-e3(iv,i)*e4(iv,i+1)
          b(iv,j)=e2(iv,i)*e2(iv,i+1)-e4(iv,i)*e4(iv,i+1)
          d(iv,j)=e1(iv,i+1)*e4(iv,i+1)-e2(iv,i+1)*e3(iv,i+1)
          e(iv,j)=e2(iv,i+1)*(cah(iv,i+1)-cab(iv,i))
     &           +e4(iv,i+1)*(cdb(iv,i)-cdh(iv,i+1))
        enddo
      enddo ! of do i=1,n-1
c     
      j=2*n
      do iv=1,KDLON
        a(iv,j)=e1(iv,n)-(1.E+0-emis(iv))*e3(iv,n)
        b(iv,j)=e2(iv,n)-(1.E+0-emis(iv))*e4(iv,n)
        d(iv,j)=0.E+0
        e(iv,j)=emis(iv)*pi*bsol(iv)-cab(iv,n)
     &          +(1.E+0-emis(iv))*cdb(iv,n)
      enddo
c.......................................................................
      call sys3v(KDLON,2*n,a,b,d,e,y)
c....................................................................... 
      do i=1,n
        do iv=1,KDLON
          grg(iv,i)=grg(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
          grh(iv,i)=grh(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
          grj(iv,i)=grj(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
          grk(iv,i)=grk(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
        enddo
      enddo
c.......................................................................
c values of "hemispheric constant" fluxes
c 
      IF (nsf.eq.0) THEN
        do i=1,n
          do iv=1,KDLON
            fah(iv,i)=e3(iv,i)*y(iv,2*i-1)-e4(iv,i)*y(iv,2*i)+cah(iv,i)
            fdh(iv,i)=e1(iv,i)*y(iv,2*i-1)-e2(iv,i)*y(iv,2*i)+cdh(iv,i)
          enddo
        enddo
        do iv=1,KDLON
          fah(iv,n+1)=e1(iv,n)*y(iv,2*n-1)+e2(iv,n)*y(iv,2*n)+cab(iv,n)
          fdh(iv,n+1)=e3(iv,n)*y(iv,2*n-1)+e4(iv,n)*y(iv,2*n)+cdb(iv,n)
        enddo
      ELSE
c....................................................................... 
c going to the "source function" 
c 
c apply a quadrature over nq (fixed parameter) points
c x is the vector of the \mu of the quadrature
c w is the vector of corresponding weights
c x() et w() are fixed parameters
c 
c....................................................................... 
c start from the top and go down along the nq angles to compute all
c downward fluxes
c
      do j=1,nq 
        do iv=1,KDLON
          gri(iv,j)=0.E+0
        enddo
      enddo
      do iv=1,KDLON
        fdh(iv,1)=0.E+0
      enddo
      do i=1,n
        do j=1,nq
          do iv=1,KDLON
            gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
     &      +grj(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
     &      *(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
     &      +grk(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
     &      *(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
     &      +sigma1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
     &      +sigma2(iv,i)*(x(j)*exp(-tau(iv,i)/x(j))+tau(iv,i)-x(j))
          enddo ! of do iv=1,KDLON
        enddo ! of do j=1,nq
        do iv=1,KDLON
          fdh(iv,i+1)=0.E+0
        enddo
        do j=1,nq 
          do iv=1,KDLON
            fdh(iv,i+1)=fdh(iv,i+1)+w(j)*x(j)*gri(iv,j)
          enddo ! of do iv=1,KDLON
        enddo ! of do j=1,nq
      enddo ! of do i=1,n
c.......................................................................
c apply the reflexion condition on the ground
c
      do iv=1,KDLON
        fah(iv,n+1)=(1.E+0-emis(iv))*fdh(iv,n+1)+pi*emis(iv)*bsol(iv)
      enddo
      do j=1,nq
        do iv=1,KDLON
          gri(iv,j)=2.E+0*fah(iv,n+1)
        enddo
      enddo
c.......................................................................
c going back up to compute all the upward fluxes
c  
      do i=n,1,-1 
        do j=1,nq 
          do iv=1,KDLON
            gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
     &      +grg(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
     &      *(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
     &      +grh(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
     &      *(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
     &      +alpha1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
     &      +alpha2(iv,i)*(x(j)-(tau(iv,i)+x(j))*exp(-tau(iv,i)/x(j)))
          enddo ! of do iv=1,KDLON
        enddo ! of do j=1,nq
        do iv=1,KDLON
          fah(iv,i)=0.E+0
        enddo
        do j=1,nq 
          do iv=1,KDLON
            fah(iv,i)=fah(iv,i)+w(j)*x(j)*gri(iv,j)
          enddo
        enddo ! of do j=1,nq
      enddo ! of do i=n,1,-1
c.......................................................................
      ENDIF ! of IF (nsf.eq.0)
c.......................................................................
c
c.......................................................................

      END SUBROUTINE flusv

c ***************************************************************


      SUBROUTINE sys3v(KDLON,n,a,b,d,e,y)
      use dimradmars_mod, only: ndlon, ndlo2, nflev
      IMPLICIT NONE
c.......................................................................
c
c  solve a tridiagonal linear system such that:
c
c  |  b1  d1               |   | y1 |   | e1 |
c  |  a2  b2  d2           |   | y2 |   | e2 |
c  |      a3  b3  d3       | * | y3 | = | e3 |
c  |             ....      |   |    |   |    |
c  |              an  bn   |   | yn |   | en |
c
c  in :   * KDLON        --> vectorisation dimension
c         * n            --> system size
c         * a,b,d,e      --> coefficients as shown above
c
c  out :  * y            --> see above
c
c.......................................................................
c  arguments 
c
      INTEGER,INTENT(IN) :: KDLON,n
      REAL,INTENT(IN) :: a(NDLO2,n),b(NDLO2,n),d(NDLO2,n),e(NDLO2,n)
      REAL,INTENT(OUT) :: y(NDLO2,n)
c.......................................................................
c  local variables
c 
      INTEGER :: iv,i
      REAL :: as(NDLON,4*nflev),ds(NDLON,4*nflev)
     &    ,x(NDLON,4*nflev)
c.......................................................................
c
c.......................................................................
      do iv=1,KDLON
        as(iv,n)=a(iv,n)/b(iv,n)
        ds(iv,n)=e(iv,n)/b(iv,n)
      enddo
      do i=n-1,1,-1
        do iv=1,KDLON
          x(iv,i)=1.E+0/(b(iv,i)-d(iv,i)*as(iv,i+1))
          as(iv,i)=a(iv,i)*x(iv,i)
          ds(iv,i)=(e(iv,i)-d(iv,i)*ds(iv,i+1))*x(iv,i)
        enddo
      enddo 
      do iv=1,KDLON
        y(iv,1)=ds(iv,1)
      enddo
      do i=2,n 
        do iv=1,KDLON
          y(iv,i)=ds(iv,i)-as(iv,i)*y(iv,i-1)
        enddo
      enddo
c.......................................................................
c
c.......................................................................

      END SUBROUTINE sys3v

      END MODULE flusv_mod