source: LMDZ6/branches/contrails/libf/phylmd/ecrad/ifs/yoe_spectral_planck.F90 @ 5450

! (C) Copyright 2019- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

MODULE YOE_SPECTRAL_PLANCK

! YOE_SPECTRAL_PLANCK
!
! PURPOSE
! -------
!   Calculate Planck function integrated across user-specified
!   spectral intervals, used in RADHEATN by approximate longwave
!   update scheme to modify longwave fluxes to account for the
!   spectral emissivity on the high-resolution model grid (rather than
!   the lower resolution grid seen by the radiation scheme).
!
! INTERFACE
! ---------
!   Call the INIT member routine to configure the look-up table of the
!   TSPECRALPLANCK type, followed by any number of CALC calls with the
!   temperatures at which the Planck function is required. FREE then
!   deallocates memory.
!
! AUTHOR
! ------
!   Robin Hogan, ECMWF
!   Original: 2019-02-04
!
! MODIFICATIONS
! -------------
!   A Dawson 2019-08-05 avoid single precision overflow in INIT

!-----------------------------------------------------------------------

USE PARKIND1,         ONLY :   JPRB,JPRD,JPIM
IMPLICIT NONE
SAVE

!-----------------------------------------------------------------------
! Type for storing Planck function look-up table
TYPE TSPECTRALPLANCK
  ! Number of intervals over which the integrated Planck function is
  ! required. Note that an interval need not be contiguous in
  ! wavelength.
  INTEGER(KIND=JPIM) :: NINTERVALS

  ! Number of temperatures in look-up table
  INTEGER(KIND=JPIM) :: NTEMPS

  ! Start temperature and temperature spacing of look-up table
  REAL(KIND=JPRB) :: TEMP1, DTEMP

  ! Integrated Planck functions in look-up table, dimensioned
  ! (NINTERVALS,NTEMPS)
  REAL(KIND=JPRB),    ALLOCATABLE :: PLANCK_LUT(:,:)

  ! Store interval data
  REAL(KIND=JPRB),    ALLOCATABLE :: WAVLEN_BOUND(:)
  INTEGER(KIND=JPIM), ALLOCATABLE :: INTERVAL_MAP(:)

CONTAINS
  PROCEDURE :: INIT
  PROCEDURE :: CALC
  PROCEDURE :: PRINT=> PRINT_SPECTRAL_PLANCK
  PROCEDURE :: FREE => FREE_SPECTRAL_PLANCK

END TYPE TSPECTRALPLANCK

CONTAINS

!-----------------------------------------------------------------------
! Generate a Planck function look-up table consisting of KINTERVALS
! spectral intervals (which need not be contiguous in wavelength),
! whose wavelength bounds are defined by PWAVLEN_BOUND and mapping
! on to KINTERVALS described by KINTERVAL_MAP.
SUBROUTINE INIT(SELF, KINTERVALS, PWAVLEN_BOUND, KINTERVAL_MAP)

  USE YOMCST,   ONLY : RPI, RKBOL, RHPLA, RCLUM
  USE YOMHOOK,  ONLY : LHOOK, DR_HOOK, JPHOOK
  USE YOMLUN,   ONLY : NULOUT

  CLASS(TSPECTRALPLANCK), INTENT(INOUT) :: SELF
  INTEGER(KIND=JPIM)    , INTENT(IN)    :: KINTERVALS
  REAL(KIND=JPRB)       , INTENT(IN)    :: PWAVLEN_BOUND(:)
  INTEGER(KIND=JPIM)    , INTENT(IN)    :: KINTERVAL_MAP(:)

  ! Current temperature (K)
  REAL(KIND=JPRB) :: ZTEMP

  ! Combinations of constants in the Planck function
  REAL(KIND=JPRB) :: ZCOEFF1, ZCOEFF2

  ! Wavelengths at start and end of range
  REAL(KIND=JPRB) :: ZWAVLEN1, ZWAVLEN2, DWAVLEN

  ! Wavelength, wavelength squared
  REAL(KIND=JPRB) :: ZWAVLEN, ZWAVLEN_SQR

  ! Sum of Planck values, integration weight
  REAL(KIND=JPRB) :: ZSUM, ZWEIGHT

  ! A double-precision temporary to hold the exponential term in Planck's law.
  ! A single precision float can overflow during this calculation.
  REAL(KIND=JPRD) :: ZPLANCKEXP

  ! Number of wavelength ranges represented by PWAVLEN_BOUND and
  ! KINTERVAL_MAP
  INTEGER(KIND=JPIM) :: NRANGES, NWAVLEN

  INTEGER(KIND=JPIM) :: JT, JI, JW, JR

  REAL(KIND=JPHOOK)    :: ZHOOK_HANDLE

#include "abor1.intfb.h"

  IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',0,ZHOOK_HANDLE)

  IF (KINTERVALS == 1) THEN
    ! We can use Stefan-Boltzmann law without a look-up table
    WRITE(NULOUT,'(a)') 'YOE_SPECTRAL_PLANCK: Single-band look-up table requested: use Stefan-Boltzmann law'
    SELF%NINTERVALS = KINTERVALS
    CALL SELF%FREE
  ELSE
    ! Full look-up table required
    ZCOEFF1 = 2.0_JPRB * RPI * RHPLA * RCLUM * RCLUM
    ZCOEFF2 = RHPLA * RCLUM / RKBOL

    NRANGES = SIZE(KINTERVAL_MAP,1)
    IF (SIZE(PWAVLEN_BOUND,1) /= NRANGES-1) THEN
      CALL ABOR1('YOS_SPECTRAL_PLANCK:INIT: PWAVLEN_BOUND must have one fewer elements than KINTERVAL_MAP')
    ENDIF

    CALL SELF%FREE

    ALLOCATE(SELF%WAVLEN_BOUND(NRANGES-1))
    ALLOCATE(SELF%INTERVAL_MAP(NRANGES))
    SELF%WAVLEN_BOUND(1:NRANGES-1) = PWAVLEN_BOUND(1:NRANGES-1)
    SELF%INTERVAL_MAP(1:NRANGES)   = KINTERVAL_MAP(1:NRANGES)

    SELF%NINTERVALS = KINTERVALS
    ! Temperature in 1-K intervals from 150 K to 350 K
    SELF%TEMP1  = 150.0_JPRB
    SELF%DTEMP  = 1.0_JPRB
    SELF%NTEMPS = 1 + NINT((350.0_JPRB - SELF%TEMP1) / SELF%DTEMP)

    ALLOCATE(SELF%PLANCK_LUT(SELF%NINTERVALS,SELF%NTEMPS))
    SELF%PLANCK_LUT(:,:) = 0.0_JPRB

    ! Print the properties of the look-up table
    WRITE(NULOUT,'(a,i0,a,f5.1,a,f5.1,a)') &
         &  'YOE_SPECTRAL_PLANCK: Generating Planck look-up table with ', &
         &  SELF%NTEMPS, ' temperatures from ', &
         &  SELF%TEMP1, ' to ', SELF%TEMP1+SELF%DTEMP*(SELF%NTEMPS-1), ' K:'
    DO JI = 1,SELF%NINTERVALS
      WRITE(NULOUT,'(a,i0,a)',advance='no') '  Band ', JI, ':'
      DO JR = 1,NRANGES
        IF (KINTERVAL_MAP(JR) == JI) THEN
          IF (JR == 1) THEN
            WRITE(NULOUT,'(a,f0.2)',advance='no') ' 0.00-', &
                 &  PWAVLEN_BOUND(1)*1.0e6_JPRB
          ELSEIF (JR == NRANGES) THEN
            WRITE(NULOUT,'(a,f0.2,a)',advance='no') ' ', &
                 &  PWAVLEN_BOUND(JR-1)*1.0e6_JPRB, '-Inf'
          ELSE
            WRITE(NULOUT,'(a,f0.2,a,f0.2)',advance='no') ' ', &
                 &  PWAVLEN_BOUND(JR-1)*1.0e6_JPRB, '-', &
                 &  PWAVLEN_BOUND(JR)*1.0e6_JPRB
          ENDIF
        ENDIF
      ENDDO
      WRITE(NULOUT,'(a)') ' microns'
    ENDDO

    ! Create the look-up table
    DO JT = 1,SELF%NTEMPS

      ZTEMP = SELF%TEMP1 + (JT-1) * SELF%DTEMP

      DO JI = 1,NRANGES

        IF (JI == 1) THEN
          ZWAVLEN1 = MIN(1.0E-6_JPRB, 0.8_JPRB * PWAVLEN_BOUND(1))
          ZWAVLEN2 = PWAVLEN_BOUND(1)
        ELSEIF (JI == NRANGES) THEN
          ZWAVLEN1 = PWAVLEN_BOUND(NRANGES-1)
          ! Simulate up to at least 200 microns wavelength
          ZWAVLEN2 = MAX(200.0E-6_JPRB, PWAVLEN_BOUND(NRANGES-1)+20.0E-6_JPRB)
        ELSE
          ZWAVLEN1 = PWAVLEN_BOUND(JI-1)
          ZWAVLEN2 = PWAVLEN_BOUND(JI)
        ENDIF

        NWAVLEN = 100
        DWAVLEN = (ZWAVLEN2 - ZWAVLEN1) / NWAVLEN
        ZSUM = 0.0_JPRB
        DO JW = 0,NWAVLEN
          ZWAVLEN = ZWAVLEN1 + DWAVLEN*JW
          ! Weights for trapezoidal rule
          !IF (JW > 0 .AND. JW < NWAVLEN) THEN
          !  ZWEIGHT = 2.0_JPRB
          !ELSE
          !  ZWEIGHT = 1.0_JPRB
          !ENDIF
          ! Weights for Simpson's rule
          IF (JW > 0 .AND. JW < NWAVLEN) THEN
            ZWEIGHT = 2.0_JPRB + 2.0_JPRB * MOD(JW,2)
          ELSE
            ZWEIGHT = 1.0_JPRB
          ENDIF
          ! Planck's law
          !
          ! The exponential term is computed in double precision to avoid
          ! overflow. The final result should still be in the range of a single
          ! precision float.
          ZWAVLEN_SQR = ZWAVLEN*ZWAVLEN
          ZPLANCKEXP = EXP(REAL(ZCOEFF2, JPRD) &
                     &     / (REAL(ZWAVLEN, JPRD) * REAL(ZTEMP, JPRD)))
          ZSUM = ZSUM + ZWEIGHT / (ZWAVLEN_SQR*ZWAVLEN_SQR*ZWAVLEN &
               &  * (ZPLANCKEXP - 1.0_JPRB))
        ENDDO
        SELF%PLANCK_LUT(KINTERVAL_MAP(JI),JT) = SELF%PLANCK_LUT(KINTERVAL_MAP(JI),JT) &
             &  + ZCOEFF1 * ZSUM * DWAVLEN / 3.0_JPRB
      ENDDO

    ENDDO

  ENDIF

  IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',1,ZHOOK_HANDLE)

END SUBROUTINE INIT


!-----------------------------------------------------------------------
! Calculate Planck function in spectral intervals from temperature
SUBROUTINE CALC(SELF, KIDIA, KFDIA, KLON, PTEMPERATURE, PPLANCK)

  USE YOMCST,   ONLY : RSIGMA
  USE YOMHOOK,  ONLY : LHOOK, DR_HOOK, JPHOOK

  CLASS(TSPECTRALPLANCK), INTENT(IN)  :: SELF
  ! Process columns KIDIA-KFDIA from total of KLON columns
  INTEGER(KIND=JPIM)    , INTENT(IN)  :: KIDIA, KFDIA, KLON
  ! Temperature in Kelvin
  REAL(KIND=JPRB)       , INTENT(IN)  :: PTEMPERATURE(KLON)
  ! Integrated Planck function as an irradiance, in W m-2
  REAL(KIND=JPRB)       , INTENT(OUT) :: PPLANCK(KLON,SELF%NINTERVALS)

  ! Column loop counter, index to temperature interval
  INTEGER(KIND=JPRB) :: JL, ITEMP

  ! Interpolation weight, highest temperature in look-up table
  REAL(KIND=JPRB) :: ZWEIGHT, ZTEMP2

  REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

  IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',0,ZHOOK_HANDLE)

  IF (SELF%NINTERVALS == 1) THEN
    ! Stefan-Boltzmann law
    PPLANCK(KIDIA:KFDIA,1) = RSIGMA * PTEMPERATURE(KIDIA:KFDIA)**4

  ELSE
    ! Look-up table
    ZTEMP2 = SELF%TEMP1 + SELF%DTEMP * (SELF%NTEMPS - 1)

    DO JL = KIDIA,KFDIA

      IF (PTEMPERATURE(JL) <= SELF%TEMP1) THEN
        ! Cap the Planck function at the low end
        ITEMP   = 1
        ZWEIGHT = 0.0_JPRB
      ELSEIF (PTEMPERATURE(JL) < ZTEMP2) THEN
        ! Linear interpolation
        ZWEIGHT = 1.0_JPRB + (PTEMPERATURE(JL) - SELF%TEMP1) / SELF%DTEMP
        ITEMP   = NINT(ZWEIGHT)
        ZWEIGHT = ZWEIGHT - ITEMP
      ELSE
        ! Linear extrapolation at high temperatures off the scale
        ITEMP   = SELF%NTEMPS-1
        ZWEIGHT = 1.0_JPRB + (PTEMPERATURE(JL) - SELF%TEMP1) / SELF%DTEMP - ITEMP
      ENDIF

      PPLANCK(JL,:) = SELF%PLANCK_LUT(:,ITEMP) &
           &  + ZWEIGHT * (SELF%PLANCK_LUT(:,ITEMP+1) - SELF%PLANCK_LUT(:,ITEMP))

      ! Force sum to equal Stefan-Boltzmann law
      PPLANCK(JL,:) = PPLANCK(JL,:) * RSIGMA * PTEMPERATURE(JL)**4 / SUM(PPLANCK(JL,:),1)

    ENDDO

  ENDIF

  IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',1,ZHOOK_HANDLE)

END SUBROUTINE CALC


!-----------------------------------------------------------------------
! Print look-up table to a unit
SUBROUTINE PRINT_SPECTRAL_PLANCK(SELF, IUNIT)

  CLASS(TSPECTRALPLANCK), INTENT(IN) :: SELF
  INTEGER(KIND=JPIM),     INTENT(IN) :: IUNIT

  INTEGER(KIND=JPIM) :: JT

  CHARACTER(len=24)  :: MY_FORMAT

  IF (SELF%NINTERVALS == 1) THEN

    WRITE(IUNIT,'(A)') 'Spectral Planck in only one interval: using Stefan-Boltzmann law'

  ELSE

    WRITE(IUNIT,'(A,I0,A)') 'Spectral Planck look-up table defined in ', &
         &  SELF%NINTERVALS, ' intervals:'
    WRITE(MY_FORMAT,'(A,I0,A)') '(f7.2,', SELF%NINTERVALS, 'e15.5)'
    DO JT = 1,SELF%NTEMPS
      WRITE(IUNIT,TRIM(MY_FORMAT)) SELF%TEMP1 + (JT-1) * SELF%DTEMP, &
           &  SELF%PLANCK_LUT(:,JT)
    ENDDO

  ENDIF

END SUBROUTINE PRINT_SPECTRAL_PLANCK


!-----------------------------------------------------------------------
! Free allocated memory
SUBROUTINE FREE_SPECTRAL_PLANCK(SELF)

  CLASS(TSPECTRALPLANCK), INTENT(INOUT) :: SELF

  IF (ALLOCATED(SELF%PLANCK_LUT)) THEN
    DEALLOCATE(SELF%PLANCK_LUT)
  ENDIF
  IF (ALLOCATED(SELF%WAVLEN_BOUND)) THEN
    DEALLOCATE(SELF%WAVLEN_BOUND)
  ENDIF
  IF (ALLOCATED(SELF%INTERVAL_MAP)) THEN
    DEALLOCATE(SELF%INTERVAL_MAP)
  ENDIF

END SUBROUTINE FREE_SPECTRAL_PLANCK


END MODULE YOE_SPECTRAL_PLANCK