MODULE radiative_lw

#include "use_logging.h"

  IMPLICIT NONE
  SAVE

  PRIVATE

  PUBLIC :: lw

  LOGICAL, PARAMETER :: lstrong=.TRUE.
  REAL, PARAMETER :: stephan=5.67e-08

CONTAINS

  SUBROUTINE lw(ngrid,nlayer,coefir,emissiv, &
       pp,ps_rad,ptsurf,pt,              &
       pfluxir,pdtlw,                    &
       lverbose, lwrite)
    USE phys_const, ONLY : cpp, g
    USE writefield_mod, ONLY : writefield
    !=======================================================================
    !
    !   calcul de l evolution de la temperature sous l effet du rayonnement
    !   infra-rouge.
    !   Pour simplifier, les transmissions sont precalculees et ne
    !   dependent que de l altitude.
    !
    !   arguments:
    !   ----------
    !
    !   entree:
    !   -------
    !      ngrid             nombres de points de la grille horizontale
    !      nlayer              nombre de couches
    !      ptsurf(ngrid)     temperature de la surface
    !      pt(ngrid,nlayer)    temperature des couches
    !      pp(ngrid,nlayer+1)  pression entre les couches
    !      lwrite            variable logique pour sorties
    !
    !   sortie:
    !   -------
    !      pdtlw(ngrid,nlayer) taux de refroidissement
    !      pfluxir(ngrid)    flux infrarouge sur le sol
    !
    !=======================================================================

    !   declarations:
    !   -------------

    !   arguments:
    !   ----------

    INTEGER, INTENT(IN)  ::  ngrid,nlayer
    REAL,    INTENT(IN)  :: coefir,emissiv(ngrid),ps_rad
    REAL,    INTENT(IN)  ::  ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
    REAL,    INTENT(OUT) ::  pdtlw(ngrid,nlayer),pfluxir(ngrid)
    LOGICAL, INTENT(IN)  :: lwrite, lverbose

    !   variables locales:
    !   ------------------

    INTEGER nlevel,ilev,ig,i,il
    REAL zplanck(ngrid,nlayer+1),zcoef
    REAL zfluxup(ngrid,nlayer+1),zfluxdn(ngrid,nlayer+1)
    REAL zflux(ngrid,nlayer+1)
    REAL zlwtr1(ngrid),zlwtr2(ngrid)
    REAL zup(ngrid,nlayer+1),zdup(ngrid)

    CHARACTER(6), PARAMETER :: tag='rad/lw'
    !-----------------------------------------------------------------------
    !   initialisations:
    !   ----------------

    nlevel=nlayer+1

    !-----------------------------------------------------------------------
    !   2. calcul des quantites d absorbants:
    !   -------------------------------------

    !   absorption forte
    IF(lstrong) THEN
       DO ilev=1,nlevel
          DO ig=1,ngrid
             zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
          ENDDO
       ENDDO
       IF(lverbose) THEN
          DO ilev=1,nlayer
             WRITELOG(*,*) ' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
          ENDDO
          LOG_DBG(tag)
       ENDIF
       zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))

       !   absorption faible
    ELSE
       DO ilev=1,nlevel
          DO ig=1,ngrid
             zup(ig,ilev)=pp(ig,ilev)
          ENDDO
       ENDDO
       zcoef=-log(coefir)/ps_rad
    ENDIF


    !-----------------------------------------------------------------------
    !   2. calcul de la fonction de corps noir:
    !   ---------------------------------------

    DO ilev=1,nlayer
       DO ig=1,ngrid
          zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
          zplanck(ig,ilev)=stephan* &
               zplanck(ig,ilev)*zplanck(ig,ilev)
       ENDDO
    ENDDO

    !-----------------------------------------------------------------------
    !   4. flux descendants:
    !   --------------------

    DO ilev=1,nlayer
       DO ig=1,ngrid
          zfluxdn(ig,ilev)=0.
       ENDDO
       DO ig=1,ngrid
          zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
       ENDDO
       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)

       DO il=nlayer,ilev,-1
          zlwtr2(:)=zlwtr1(:)
          DO ig=1,ngrid
             zdup(ig)=zup(ig,ilev)-zup(ig,il)
          ENDDO
          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
          DO ig=1,ngrid
             zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+ &
                  zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
          ENDDO
       ENDDO
    ENDDO

    DO ig=1,ngrid
       zfluxdn(ig,nlevel)=0.
       pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
    ENDDO

    DO ig=1,ngrid
       zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
       zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1) &
            +(1.-emissiv(ig))*zfluxdn(ig,1)
    ENDDO

    !-----------------------------------------------------------------------
    !   3. flux montants:
    !   ------------------

    DO ilev=1,nlayer
       DO ig=1,ngrid
          zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
       ENDDO
       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
       DO ig=1,ngrid
          zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
       ENDDO
       DO il=1,ilev
          zlwtr2(:)=zlwtr1(:)
          DO ig=1,ngrid
             zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
          ENDDO
          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
          DO ig=1,ngrid
             zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+ &
                  zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
          ENDDO
       ENDDO

    ENDDO

    !-----------------------------------------------------------------------
    !   5. calcul des flux nets:
    !   ------------------------

    DO ilev=1,nlevel
       DO ig=1,ngrid
          zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
       ENDDO
    ENDDO

    !-----------------------------------------------------------------------
    !   6. Calcul des taux de refroidissement:
    !   --------------------------------------

    DO ilev=1,nlayer
       DO ig=1,ngrid
          pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))* &
               g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
       ENDDO
    ENDDO

    !-----------------------------------------------------------------------
    !   10. sorties eventuelles:
    !   ------------------------

    IF (lverbose) THEN
       WRITELOG(*,*) 'Diagnostique rayonnement thermique'
       WRITELOG(*,*) 'temperature     ', &
            'flux montant    flux desc.     taux de refroid.'
       i=ngrid/2+1
       WRITELOG(6,'(4e18.4)') ptsurf(i)
       DO ilev=1,nlayer
          WRITELOG(6,'(i4,4e18.4)') ilev,pt(i,ilev), &
               zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
       ENDDO
       WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)
       LOG_DBG(tag)
    ENDIF

    IF(lwrite) THEN
       CALL writefield('lwflux_up', 'Upward LW flux', 'W/m2', zfluxup)
       CALL writefield('lwflux_down', 'Downward LW flux', 'W/m2', zfluxdn)
    END IF
    !-----------------------------------------------------------------------

  END SUBROUTINE lw

  PURE SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
    INTEGER, INTENT(IN) :: ngrid
    REAL,    INTENT(IN) :: coef
    LOGICAL, INTENT(IN) :: lstrong
    REAL,    INTENT(IN) :: dup(ngrid)
    REAL,    INTENT(OUT) :: transm(ngrid)
    INTEGER ig
    IF(lstrong) THEN
       DO ig=1,ngrid
          transm(ig)=exp(-coef*sqrt(dup(ig)))
       ENDDO
    ELSE
       DO ig=1,ngrid
          transm(ig)=exp(-coef*dup(ig))
       ENDDO
    ENDIF

  END SUBROUTINE lwtr

END MODULE radiative_lw