!= Lacis and Hansen (1974) による短波放射モデル
!
!= Short wave radiation model described by Lacis and Hansen (1974)
!
! Authors::   Yoshiyuki O. Takahashi, Satoshi NODA
! Version::   $Id: radiation_LH74.f90,v 1.5 2010-04-12 02:29:06 noda Exp $
! Tag Name::  $Name: dcpam5-20101008-1 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module radiation_LH74
  !
  != Lacis and Hansen (1974) による短波放射モデル
  !
  != Short wave radiation model described by Lacis and Hansen (1974)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 短波放射モデル.
  !
  ! This is a model of short wave radiation. 
  !
  !== References
  !
  !  Lacis, A. A., and J. E. Hansen, A parameterization for the absorption of solar 
  !    radiation in the Earth's atmosphere, J. Atmos. Sci., 31, 118-133, 1974.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#radiation_DennouAGCM_nml
  !

  ! USE statements
  !

  ! 
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  ! Declaration statements
  !
  implicit none
  private

  ! 
  ! Public procedure
  !
  public:: RadiationLH74Flux, RadiationLH74AddO3Heating


  real(DP), save     :: H2OScaleIndex

  real(DP), save     :: DiffFactorH2O
  real(DP), save     :: DiffFactorO3

  real(DP), save     :: ShortAtmosAlbedo
                              ! 大気アルベド.
                              ! Albedo of air.
  logical , save     :: FlagSimpleMagFac
                              ! エアマス計算のためのフラグ
                              ! Flag for air-mass calculation
  logical , save     :: radiation_lh74_inited

  data radiation_lh74_inited /.false./


  character(*), parameter:: module_name = 'radiation_LH74'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20101008-1 $' // &
    & '$Id: radiation_LH74.f90,v 1.5 2010-04-12 02:29:06 noda Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  subroutine RadiationLH74Flux( &
    & xy_SurfAlbedo, xyz_Temp, xyz_QVap, xyr_Press, xyz_Press, & ! (in )
    & xyr_RadSFlux                                             & ! (out)
    & )


    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, & ! 
                               ! Number of grid points in longitude
      &                jmax, & ! 
                               ! Number of grid points in latitude
      &                kmax    ! 
                               ! Number of vertical level

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, & ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration
                         PI      ! $ \pi $ .
                                 ! Circular constant

    ! 短波入射 (太陽入射)
    ! Short wave (insolation) incoming
    !
    use radiation_short_income, only: ShortIncoming

    real(DP), intent(in ):: xy_SurfAlbedo   (0:imax-1, 1:jmax)
    real(DP), intent(in ):: xyz_Temp        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QVap        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press       (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyr_RadSFlux    (0:imax-1, 1:jmax, 0:kmax)


    !
    ! Work variables
    !
    real(DP):: RefPress
    real(DP):: RefTemp

    real(DP):: xyz_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyr_H2OColDen   (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_H2OAbsAmt   (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xy_IncomRadSFlux(0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス.
                              ! Short wave (insolation) flux
    real(DP):: xy_InAngle      (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP):: xy_MagFac       (0:imax-1, 1:jmax)

!!$    logical :: flag_dry_atmosphere

    integer :: i
    integer :: j
    integer :: k


    if ( .not. radiation_lh74_inited ) then
      call RadiationLH74Init
    end if


!!$    ! Check for dry atmosphere
!!$    !
!!$    if ( all( xyz_QVap <= 0.0d0 ) ) then
!!$      flag_dry_atmosphere = .true.
!!$      write( 6, * ) 'Dry atmosphere'
!!$    else
!!$      flag_dry_atmosphere = .false.
!!$    end if


    ! 短波入射の計算
    ! Calculate short wave (insolation) incoming radiation
    !
    call ShortIncoming( &
      & xy_IncomRadSFlux, xy_InAngle & ! (out)
      & )

    ! 大気アルベドの考慮
    ! Taking atmospheric albedo into consideration
    !
    xy_IncomRadSFlux = xy_IncomRadSFlux * ( 1.0d0 - ShortAtmosAlbedo )

!!$    if ( flag_dry_atmosphere ) then
!!$      do k = 0, kmax
!!$        xyr_RadSFlux(:,:,k) = - xy_IncomRadSFlux(:,:) + ...
!!$      end do
!!$      return
!!$    end if


    RefPress = 1.013d5
    RefTemp  = 273.0d0

    do k = 1, kmax
      xyz_H2ODelAbsAmt(:,:,k) = &
        &   ( xyz_Press(:,:,k) / RefPress        )**H2OScaleIndex &
        & * ( RefTemp          / xyz_Temp(:,:,k) )**0.5           &
        & * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
!!$      xyz_H2ODelAbsAmt(:,:,k) = &
!!$        &   ( xyz_Press(:,:,k) / RefPress ) &
!!$        & * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do


    xyr_H2OColDen(:,:,:) = 0.0d0

    do k = kmax-1, 0, -1
      xyr_H2OColDen(:,:,k) = xyr_H2OColDen(:,:,k+1) + xyz_H2ODelAbsAmt(:,:,k+1)
    end do


    if ( FlagSimpleMagFac ) then

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = xy_InAngle(i,j)
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    else

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = 35.0d0 &
              & / sqrt( 1224.0d0 * ( 1.0d0 / xy_InAngle(i,j) )**2 + 1.0d0 )
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    end if


    xyr_RadSFlux(:,:,:) = 0.0d0

    ! Downward flux
    do k = 0, kmax
      xyr_H2OAbsAmt(:,:,k) = &
        &   xyr_H2OColDen(:,:,k) * xy_MagFac(:,:)
    end do
    call RadiationLH74CalcAbsorptivity( &
      & xyr_H2OAbsAmt,    & ! (in )
      & xyr_Absorptivity  & ! (out)
      & )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
        & + xy_IncomRadSFlux(:,:)               &
        & * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do

    ! Upward flux
    do k = 0, kmax
      xyr_H2OAbsAmt(:,:,k) = &
        &   xyr_H2OColDen(:,:,0) * xy_MagFac(:,:) &
        & + ( xyr_H2OColDen(:,:,0) - xyr_H2OColDen(:,:,k) ) * DiffFactorH2O
    end do
    call RadiationLH74CalcAbsorptivity( &
      & xyr_H2OAbsAmt,    & ! (in )
      & xyr_Absorptivity  & ! (out)
      & )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
        & - xy_IncomRadSFlux(:,:)               &
        & * xy_SurfAlbedo(:,:)                  &
        & * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do


  end subroutine RadiationLH74Flux

  !**************************************************************************

  subroutine RadiationLH74AddO3Heating( &
    & xyz_QO3, xy_SurfAlbedo, & ! (in   )
    & xyr_Press, xyz_Press,   & ! (in   )
    & xyz_DTempDtRadS         & ! (inout)
    & )
    !
    ! Lacis and Hansen (1974) の計算方を用いて O3 の短波吸収加熱率を計算し, 
    ! 短波放射加熱率に加える. 
    !
    ! Calculate radiative heating rate due to absorption of short wave radiation of O3
    ! and add to the short wave radiative heating rate
    !


    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, & ! 
                               ! Number of grid points in longitude
      &                jmax, & ! 
                               ! Number of grid points in latitude
      &                kmax    ! 
                               ! Number of vertical level

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, & ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration
                         PI,   & ! $ \pi $ .
                                 ! Circular constant
                         CpDry   ! $ C_p $ [J kg-1 K-1].
                                 ! 乾燥大気の定圧比熱.
                                 ! Specific heat of air at constant pressure

    ! 短波入射 (太陽入射)
    ! Short wave (insolation) incoming
    !
    use radiation_short_income, only: ShortIncoming
!    use radiation_short_income_sr, only: ShortIncoming

!!$  ! 時刻管理
!!$  ! Time control
!!$  !
!!$  use timeset, only: &
!!$    & TimeN                   ! ステップ $ t $ の時刻.
!!$                              ! Time of step $ t $.

    use read_time_series, only: SetValuesFromTimeSeriesWrapper

!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

    real(DP), intent(in   ):: xyz_QO3        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ):: xy_SurfAlbedo  (0:imax-1, 1:jmax)
    real(DP), intent(in   ):: xyr_Press      (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in   ):: xyz_Press      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyz_DTempDtRadS(0:imax-1, 1:jmax, 1:kmax)


    !
    ! Work variables
    !
    real(DP):: xyz_O3DelAbsAmt (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyr_O3ColDen    (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_O3AbsAmt    (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_RadSFlux    (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xy_IncomRadSFlux(0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス.
                              ! Short wave (insolation) flux
    real(DP):: xy_InAngle      (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP):: xy_MagFac       (0:imax-1, 1:jmax)

!!$    logical :: flag_dry_atmosphere

    character(STRING):: O3FileName

    integer :: i
    integer :: j
    integer :: k


    if ( .not. radiation_lh74_inited ) then
      call RadiationLH74Init
    end if


!!$    write( O3FileName, '(a,i3.3,a)' ) &
!!$      & "../../../data_Earth/O3_CMIP5_climatology_zonalmean_T", (imax-1)/3, ".nc"
!!$    call SetValuesFromTimeSeriesWrapper( &
!!$      & O3FileName, "O3", &
!!$      & xyz_Press,         &               ! (in)
!!$      & xyz_QO3,           &               ! (inout)
!!$      & "O3"            &
!!$      & )
!!$
!!$    call HistoryAutoPut( TimeN, "O3", xyz_QO3 )


    ! 短波入射の計算
    ! Calculate short wave (insolation) incoming radiation
    !
    call ShortIncoming( &
      & xy_IncomRadSFlux, xy_InAngle & ! (out)
      & )


    do k = 1, kmax
      xyz_O3DelAbsAmt(:,:,k) = &
        &   xyz_QO3(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    xyr_O3ColDen(:,:,:) = 0.0d0
    do k = kmax-1, 0, -1
      xyr_O3ColDen(:,:,k) = xyr_O3ColDen(:,:,k+1) + xyz_O3DelAbsAmt(:,:,k+1)
    end do

    if ( FlagSimpleMagFac ) then

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = xy_InAngle(i,j)
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    else

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = 35.0d0 &
              & / sqrt( 1224.0d0 * ( 1.0d0 / xy_InAngle(i,j) )**2 + 1.0d0 )
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    end if


    xyr_RadSFlux(:,:,:) = 0.0d0

    ! Downward flux
    do k = 0, kmax
      xyr_O3AbsAmt(:,:,k) = &
        &   xyr_O3ColDen(:,:,k) * xy_MagFac(:,:)
    end do
    call RadiationLH74CalcO3Absorptivity( &
      & xyr_O3AbsAmt,    & ! (in )
      & xyr_Absorptivity  & ! (out)
      & )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
        & + xy_IncomRadSFlux(:,:)               &
        & * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do

    ! Upward flux
    do k = 0, kmax
      xyr_O3AbsAmt(:,:,k) = &
        &   xyr_O3ColDen(:,:,0) * xy_MagFac(:,:) &
        & + ( xyr_O3ColDen(:,:,0) - xyr_O3ColDen(:,:,k) ) * DiffFactorO3
    end do
    call RadiationLH74CalcO3Absorptivity( &
      & xyr_O3AbsAmt,    & ! (in )
      & xyr_Absorptivity  & ! (out)
      & )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
        & - xy_IncomRadSFlux(:,:)               &
        & * xy_SurfAlbedo(:,:)                  &
        & * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do

    ! 放射加熱率の計算
    ! Calculate radiation heating rate
    !
    do k = 1, kmax
      xyz_DTempDtRadS(:,:,k) =                                   &
        &    xyz_DTempDtRadS(:,:,k)                              &
        & + (     xyr_RadSFlux(:,:,k-1) - xyr_RadSFlux(:,:,k) )  &
        &     / ( xyr_Press(:,:,k-1)    - xyr_Press(:,:,k) )     &
        &     / CpDry * Grav
    end do

!!$    i = imax / 2
!!$    j = jmax / 2
!!$    j = jmax * 3 / 4
!!$    do k = 0, kmax
!!$      write( 93, * ) xyr_RadSFlux(i,j,k), xyr_Press(i,j,k)
!!$    end do
!!$    call flush( 93 )
!!$
!!$    do k = 1, kmax
!!$      write( 83, * ) &
!!$        & + (     xyr_RadSFlux(i,j,k-1) - xyr_RadSFlux(i,j,k) )  &
!!$        &     / ( xyr_Press(i,j,k-1)    - xyr_Press(i,j,k) )     &
!!$        &     / CpDry * Grav, &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 83 )
!!$
!!$
!!$    stop

  end subroutine RadiationLH74AddO3Heating

  !**************************************************************************

  subroutine RadiationLH74CalcO3Absorptivity( &
    & xyr_O3AbsAmt,    & ! (in )
    & xyr_Absorptivity  & ! (out)
    & )

    ! USE statements
    !

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: GasRUniv
                              ! $ R^{*} $ [J K-1 mol-1].
                              ! 普遍気体定数.  Universal gas constant

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, & ! 
                               ! Number of grid points in longitude
      &                jmax, & ! 
                               ! Number of grid points in latitude
      &                kmax    ! 
                               ! Number of vertical level

    real(DP), intent(in ):: xyr_O3AbsAmt    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)

    !
    ! Work variables
    !
    real(DP):: O3DensNTP
    real(DP):: xyr_O3cm(0:imax-1, 1:jmax, 0:kmax)


    O3DensNTP = 1013.25d2 / ( GasRUniv / 48.0d-3 * 273.15d0 )

    xyr_O3cm(:,:,:) = xyr_O3AbsAmt(:,:,:) / O3DensNTP * 1.0d2

    xyr_Absorptivity(:,:,:) =                                   &
      !
      &   0.02118d0 * xyr_O3cm(:,:,:)                           &
      & / (   1.0d0                                             &
      &     + 0.042d0 * xyr_O3cm(:,:,:)                         &
      &     + 0.000323d0 * xyr_O3cm(:,:,:) * xyr_O3cm(:,:,:) )  &
      !
      & + 1.082d0 * xyr_O3cm(:,:,:)                             &
      & / (   1.0d0                                             &
      &     + 138.6d0 * xyr_O3cm(:,:,:) )**0.805                &
      & + 0.0658d0 * xyr_O3cm(:,:,:)                            &
      & / (   1.0d0                                             &
      &     + ( 103.6d0 * xyr_O3cm(:,:,:) )**3 )


  end subroutine RadiationLH74CalcO3Absorptivity

  !**************************************************************************

  subroutine RadiationLH74CalcAbsorptivity( &
    & xyr_H2OAbsAmt,    & ! (in )
    & xyr_Absorptivity  & ! (out)
    & )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, & ! 
                               ! Number of grid points in longitude
      &                jmax, & ! 
                               ! Number of grid points in latitude
      &                kmax    ! 
                               ! Number of vertical level

    real(DP), intent(in ):: xyr_H2OAbsAmt   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)

    !
    ! Work variables
    !
    real(DP):: xyr_H2Oprcm(0:imax-1, 1:jmax, 0:kmax)


    xyr_H2Oprcm(:,:,:) = xyr_H2OAbsAmt(:,:,:) / ( 1.0d0 * 1.0d-3 * 1.0d6 ) * 1.0d2

    xyr_Absorptivity(:,:,:) =                                 &
      & 2.9d0 * xyr_H2Oprcm(:,:,:)                            &
      & / (                                                   &
      &       ( 1.0d0 + 141.5d0 * xyr_H2Oprcm(:,:,:) )**0.635 &
      &     + 5.925d0 * xyr_H2Oprcm(:,:,:)                    &
      &   )


  end subroutine RadiationLH74CalcAbsorptivity

  !**************************************************************************

  subroutine RadiationLH74Init


    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoAddVariable


    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /radiation_LH74_nml/ &
      & DiffFactorH2O, DiffFactorO3, ShortAtmosAlbedo, FlagSimpleMagFac
          !
          ! デフォルト値については初期化手続 "radiation_LH74#RadiationLH74Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_LH74#RadiationLH74Init" for the default values.
          !

    ! デフォルト値の設定
    ! Default values settings
    !
    DiffFactorH2O    = 5.0d0 / 3.0d0
    DiffFactorO3     = 1.9_DP

    ShortAtmosAlbedo = 0.2_DP

    FlagSimpleMagFac = .false.

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                     & ! (in)
        & nml = radiation_LH74_nml, & ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    H2OScaleIndex = 1.0d0


!!$    call HistoryAutoAddVariable( "O3",         & ! (in)
!!$      & (/ 'lon ', 'lat ', 'sig ', 'time' /),  & ! (in)
!!$      & "ozone", 'kg kg-1' )                     ! (in)


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'DiffFactorH2O    = %f', d = (/ DiffFactorH2O /) )
    call MessageNotify( 'M', module_name, 'DiffFactorO3     = %f', d = (/ DiffFactorO3 /) )
    call MessageNotify( 'M', module_name, 'ShortAtmosAlbedo = %f', d = (/ ShortAtmosAlbedo /) )
    call MessageNotify( 'M', module_name, 'FlagSimpleMagFac = %b', l = (/ FlagSimpleMagFac /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    radiation_lh74_inited = .true.
  end subroutine RadiationLH74Init


end module radiation_LH74
