source: trunk/LMDZ.PLUTO/libf/phypluto/sfluxi.F @ 3585

      module sfluxi_mod

      implicit none

      contains

      SUBROUTINE SFLUXI(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI,DWNI,
     *                  COSBI,WBARI,NFLUXTOPI,NFLUXTOPI_nu,
     *                  FMNETI,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,
     *                  FZEROI,TAUGSURF)

      use radinc_h, only: NTfac, NTstart, L_LEVELS, L_NSPECTI, L_NGAUSS
      use radinc_h, only: L_NLAYRAD, L_NLEVRAD
      use radcommon_h, only: planckir, tlimit,sigma, gweight
      use comcstfi_mod, only: pi
      use gfluxi_mod, only: gfluxi

      implicit none

      integer NLEVRAD, L, NW, NG, NTS, NTT

      real*8 TLEV(L_LEVELS), PLEV(L_LEVELS)
      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
      real*8 FMNETI(L_NLAYRAD)
      real*8 FMNETI_NU(L_NLAYRAD,L_NSPECTI)
      real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI)
      real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
      real*8 FMUPI(L_NLEVRAD), FMDI(L_NLEVRAD)
      real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
      real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
      real*8 NFLUXTOPI
      real*8 NFLUXTOPI_nu(L_NSPECTI)
      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)
      real*8 FTOPUP

      real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP
      real*8 PLANCK, PLTOP
      real*8 fluxupi(L_NLAYRAD), fluxdni(L_NLAYRAD)
      real*8 FZEROI(L_NSPECTI)
      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero

      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
      real*8 PLANCKSUM,PLANCKREF

! AB : variables for interpolation
      REAL*8 C1
      REAL*8 C2
      REAL*8 P1

!======================================================================C

      NLEVRAD = L_NLEVRAD

! ZERO THE NET FLUXES
      NFLUXTOPI = 0.0D0

      DO NW=1,L_NSPECTI
        NFLUXTOPI_nu(NW) = 0.0D0
        DO L=1,L_NLAYRAD
           FLUXUPI_nu(L,NW) = 0.0D0
           FMNETI_nu(L,NW)  = 0.0
           fup_tmp(nw)=0.0D0
           fdn_tmp(nw)=0.0D0
        END DO
      END DO

      DO L=1,L_NLAYRAD
        FMNETI(L)  = 0.0D0
        FLUXUPI(L) = 0.0D0
        FLUXDNI(L) = 0.0D0
      END DO

! WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
! TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL

      TTOP  = TLEV(2)  ! JL12 why not (1) ???
      TSURF = TLEV(L_LEVELS)

      NTS   = int(TSURF*NTfac)-NTstart+1
      NTT   = int(TTOP *NTfac)-NTstart+1

!JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
!JL12 this ensure that no flux is lost due to:
!JL12          -truncation of the planck function at high/low wavenumber
!JL12          -numerical error during first spectral integration
!JL12          -discrepancy between Tsurf and NTS/NTfac
      PLANCKSUM = 0.d0
      PLANCKREF = TSURF * TSURF
      PLANCKREF = sigma * PLANCKREF * PLANCKREF

      DO NW=1,L_NSPECTI
! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF
         C1 = TSURF * NTfac - int(TSURF * NTfac)
         P1 = (1.0D0 - C1) * PLANCKIR(NW,NTS) + C1 * PLANCKIR(NW,NTS+1)
         PLANCKSUM = PLANCKSUM + P1 * DWNI(NW)
      ENDDO

      PLANCKSUM = PLANCKREF / (PLANCKSUM * Pi)
!JL12

      DO 501 NW=1,L_NSPECTI
! SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF
! AB : idem for PLANCKIR(NW,NTT) and PLTOP
         C1 = TSURF * NTfac - int(TSURF * NTfac)
         C2 = TTOP  * NTfac - int(TTOP  * NTfac)
         P1 = (1.0D0 - C1) * PLANCKIR(NW,NTS) + C1 * PLANCKIR(NW,NTS+1)
         BSURF = (1. - RSFI) * P1 * PLANCKSUM
         PLTOP = (1.0D0 - C2) * PLANCKIR(NW,NTT) + C2*PLANCKIR(NW,NTT+1)

! If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
! special Gauss point at the end.
         FZERO = FZEROI(NW)

         IF(FZERO.ge.0.99) goto 40

         DO NG=1,L_NGAUSS-1

            if(TAUGSURF(NW,NG).lt. TLIMIT) then
               fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
               goto 30
            end if

! SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
! CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
! OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY

!            TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
            TAUTOP = TAUCUMI(2,NW,NG)
            BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP

! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER

            CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
     *                TAUCUMI(1,NW,NG),
     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
     *                BSURF,FTOPUP,FMUPI,FMDI)

! NOW CALCULATE THE CUMULATIVE IR NET FLUX
            NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)
     *                * (1.0D0-FZEROI(NW))

! and same thing by spectral band... (RDW)
            NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + FTOPUP * DWNI(NW)
     *                       * GWEIGHT(NG) * (1.0D0-FZEROI(NW))

            DO L=1,L_NLEVRAD-1
!           CORRECT FOR THE WAVENUMBER INTERVALS
               FMNETI(L)  = FMNETI(L) + (FMUPI(L)-FMDI(L)) * DWNI(NW)
     *                    * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
               FMNETI_NU(L,NW)  = FMNETI_NU(L,NW) + (FMUPI(L)-FMDI(L))
     *                    * DWNI(NW) * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
               FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)
     *                    * (1.0D0-FZEROI(NW))
               FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)
     *                    * (1.0D0-FZEROI(NW))
!         and same thing by spectral band... (RW)
               FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + FMUPI(L)*DWNI(NW)
     *                          * GWEIGHT(NG) * (1.0D0 - FZEROI(NW))
            END DO

   30       CONTINUE

         END DO       !End NGAUSS LOOP

   40    CONTINUE

! SPECIAL 17th Gauss point
         NG     = L_NGAUSS

!         TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
         TAUTOP = TAUCUMI(2,NW,NG)
         BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP

! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER

         CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
     *                TAUCUMI(1,NW,NG),
     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
     *                BSURF,FTOPUP,FMUPI,FMDI)

! NOW CALCULATE THE CUMULATIVE IR NET FLUX
         NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO

!         and same thing by spectral band... (RW)
         NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
     *      +FTOPUP*DWNI(NW)*FZERO

         DO L=1,L_NLEVRAD-1
! CORRECT FOR THE WAVENUMBER INTERVALS
            FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
            FMNETI_nu(L,NW)  = FMNETI_nu(L,NW)+(FMUPI(L)-FMDI(L))
     *       *DWNI(NW)*FZERO
            FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO
            FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO
! and same thing by spectral band... (RW)
            FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW)
     *                       + FMUPI(L) * DWNI(NW) * FZERO
         END DO

  501 CONTINUE      !End Spectral Interval LOOP
! *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****

      END SUBROUTINE SFLUXI

      end module sfluxi_mod