source: trunk/LMDZ.PLUTO/libf/phypluto/sfluxi_pluto_mod.F @ 3586

      module sfluxi_PLUTO_mod

      implicit none

      contains

      SUBROUTINE SFLUXI_PLUTO(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 BSURFtest ! by RW for test

      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
      real*8 XX,YY   ! calcul intermediate


C======================================================================C

      NLEVRAD = L_NLEVRAD


C     ZERO THE NET FLUXES

      NFLUXTOPI = 0.0
      DO NW=1,L_NSPECTI
        NFLUXTOPI_nu(NW) = 0.0
        DO L=1,L_NLAYRAD
           FLUXUPI_nu(L,NW) = 0.0
           FMNETI_nu(L,NW)  = 0.0

              fup_tmp(nw)=0.0
              fdn_tmp(nw)=0.0

        END DO
      END DO

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

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

      TTOP  = TLEV(2)
      TSURF = TLEV(L_LEVELS)
      !PRINT*, 'TB18 NTstart: ',NTstart,TSURF
      NTS   = int(TSURF*10.0D0)-NTstart + 1
      !PRINT*, 'TB18 NTS: ',NTS
      !PRINT*, 'TB18 TLEV sfluxi: ',TLEV
      NTT   = int(TTOP *10.0D0)-NTstart +1
      !PRINT*, 'TB18 NTT: ',NTT

! if temperatures superior to 50K
!      NTS   = TSURF*10.0D0-499
!      NTT   = TTOP *10.0D0-499

! analyse
c      write(*,*)'TTOP',TTOP
c      write(*,*)'TSURF',TSURF
c      write(*,*)'NTS',NTS
c      write(*,*)'NTT',NTT
c     DO L=1,L_LEVELS
c         write(*,*)'TLEV',TLEV
c      ENDDO
      BSURFtest=0.0

      DO 501  NW=1,L_NSPECTI

C       SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
        BSURF = (1.-RSFI)*PLANCKIR(NW,NTS) ! interpolate for accuracy??
        PLTOP = PLANCKIR(NW,NTT)

        !BSURFtest=BSURFtest+BSURF*DWNI(NW)
        !if(NW.eq.L_NSPECTI)then
        !     print*,'eps*sigma*T^4',5.67e-8*(1.-RSFI)*TSURF**4
        !     print*,'BSURF in sfluxi=',pi*BSURFtest
        !endif


C  If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
C  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.0-FZEROI(NW))*GWEIGHT(NG)
            goto 30
          end if

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

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

C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
C         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)


C         NOW CALCULATE THE CUMULATIVE IR NET FLUX
          XX= FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
          NFLUXTOPI =        NFLUXTOPI+ XX
          NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + XX !same by spectral band.


          DO L=1,L_NLEVRAD-1

C           CORRECT FOR THE WAVENUMBER INTERVALS
            YY = DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))

            FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))* YY
            FMNETI_nu(L,NW)  = FMNETI_nu(L,NW)+(FMUPI(L)-FMDI(L))* YY !same by spectral band
            FLUXDNI(L) = FLUXDNI(L) + FMDI(L)* YY

            FLUXUPI(L) = FLUXUPI(L) +  FMUPI(L)* YY
            FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) +  FMUPI(L)* YY  !same by spectral band.
          END DO

         !fup_tmp(nw)=fup_tmp(nw)+FMUPI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG)
         !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG)
         !fup_tmp(nw)=fup_tmp(nw)+FMUPI(1)*DWNI(NW)*GWEIGHT(NG)
         !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(1)*DWNI(NW)*GWEIGHT(NG)

   30     CONTINUE

       END DO       !End NGAUSS LOOP

           ! print*,'-----------------------------------'
            !print*,'FMDI(',nw,')=',FMDI(L_NLEVRAD-1)
            !print*,'FMUPI(',nw,')=',FMUPI(L_NLEVRAD-1)
            !print*,'DWNI(',nw,')=',DWNI(nw)
           ! print*,'-----------------------------------'

   40  CONTINUE

C      SPECIAL 17th Gauss point

    !   print*,'fzero for ng=17',fzero


       NG     = L_NGAUSS

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

C      WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
C      CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
C      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)

C      NOW CALCULATE THE CUMULATIVE IR NET FLUX

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

       DO L=1,L_NLEVRAD-1

C        CORRECT FOR THE WAVENUMBER INTERVALS

         FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO

         XX = (FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
         FMNETI(L)  = FMNETI(L)+ XX
         FMNETI_nu(L,NW)  = FMNETI_nu(L,NW)+XX !same by spectral band.

         XX=  FMUPI(L)*DWNI(NW)*FZERO
         FLUXUPI(L) = FLUXUPI(L) + XX
         FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + XX !same by spectral band.

       END DO

       !     print*,'-----------------------------------'
       !     print*,'nw=',nw
       !     print*,'ng=',ng
       !     print*,'FMDI=',FMDI(L_NLEVRAD-1)
       !     print*,'FMUPI=',FMUPI(L_NLEVRAD-1)
       !     print*,'-----------------------------------'

  501 CONTINUE      !End Spectral Interval LOOP

C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****

            !print*,'-----------------------------------'
            !print*,'gweight=',gweight
            !print*,'FLUXDNI=',FLUXDNI(L_NLEVRAD-1)
            !print*,'FLUXUPI=',FLUXUPI(L_NLEVRAD-1)
            !print*,'-----------------------------------'

            !do nw=1,L_NSPECTI
            !   print*,'fup_tmp(',nw,')=',fup_tmp(nw)
            !   print*,'fdn_tmp(',nw,')=',fdn_tmp(nw)
            !enddo


      END SUBROUTINE SFLUXI_PLUTO

      end module sfluxi_pluto_mod