source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_gas.F90 @ 4013

! radiation_gas.F90 - Derived type to store the gas mixing ratios
!
! (C) Copyright 2014- 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.
!
! Author:  Robin Hogan
! Email:   r.j.hogan@ecmwf.int
!
! Modifications
!   2019-01-14  R. Hogan  Added out_of_physical_bounds routine

module radiation_gas

  use parkind1, only : jprb

  implicit none
  public

  ! Gas codes; these indices match those of RRTM-LW up to 7
  integer, parameter :: IGasNotPresent = 0
  integer, parameter :: IH2O   = 1
  integer, parameter :: ICO2   = 2
  integer, parameter :: IO3    = 3
  integer, parameter :: IN2O   = 4
  integer, parameter :: ICO    = 5
  integer, parameter :: ICH4   = 6
  integer, parameter :: IO2    = 7
  integer, parameter :: ICFC11 = 8
  integer, parameter :: ICFC12 = 9
  integer, parameter :: IHCFC22= 10
  integer, parameter :: ICCl4  = 11 
  integer, parameter :: INO2   = 12
  integer, parameter :: NMaxGases = 12

  ! Molar masses (g mol-1) of dry air and the various gases above
  real(jprb), parameter :: IAirMolarMass = 28.970
  real(jprb), parameter, dimension(0:NMaxGases) :: IGasMolarMass = (/ &
       & 0.0_jprb,        & ! Gas not present
       & 18.0152833_jprb, & ! H2O
       & 44.011_jprb,     & ! CO2
       & 47.9982_jprb,    & ! O3
       & 44.013_jprb,     & ! N2O
       & 28.0101_jprb,    & ! CO
       & 16.043_jprb,     & ! CH4
       & 31.9988_jprb,    & ! O2
       & 137.3686_jprb,   & ! CFC11
       & 120.914_jprb,    & ! CFC12
       & 86.469_jprb,     & ! HCFC22
       & 153.823_jprb,    & ! CCl4    
       & 46.0055_jprb /)    ! NO2

  ! The corresponding names of the gases in upper and lower case, used
  ! for reading variables from the input file
  character*6, dimension(NMaxGases), parameter :: GasName &
       &  = (/'H2O   ','CO2   ','O3    ','N2O   ','CO    ','CH4   ', &
       &      'O2    ','CFC11 ','CFC12 ','HCFC22','CCl4  ','NO2   '/)
  character*6, dimension(NMaxGases), parameter :: GasLowerCaseName &
       &  = (/'h2o   ','co2   ','o3    ','n2o   ','co    ','ch4   ', &
       &      'o2    ','cfc11 ','cfc12 ','hcfc22','ccl4  ','no2   '/)

  ! Available units
  enum, bind(c)
    enumerator IMassMixingRatio, IVolumeMixingRatio
  end enum

  !---------------------------------------------------------------------
  ! This derived type describes the gaseous composition of the
  ! atmosphere; gases may be stored as part of a 3D array (if their
  ! variation with height/column is to be represented) or one 1D array
  ! (if they are to be assumed globally well-mixed).
  type gas_type
    ! Units of each stored gas (or 0 if not present)
    integer :: iunits(NMaxGases) = 0

    ! Scaling factor that should be applied to each stored gas to get
    ! a dimensionless result, e.g. if iunits=IVolumeMixingRatio then
    ! 1.0e-6 is used to indicate the units are actually PPMV: need to
    ! multiply by 1e-6 to get mol/mol.
    real(jprb) :: scale_factor(NMaxGases) = 1.0_jprb
    
    ! Mixing ratios of variable gases, dimensioned (ncol, nlev,
    ! NMaxGases)
    real(jprb), allocatable, dimension(:,:,:) :: mixing_ratio

    ! Flag to indicate whether a gas is present
    logical :: is_present(NMaxGases) = .false.

    ! Flag to indicate whether a gas is well mixed
    logical :: is_well_mixed(NMaxGases) = .false.

    integer :: ntype          = 0 ! Number of gas types described

    integer :: ncol           = 0 ! Number of columns in mixing_ratio
    integer :: nlev           = 0 ! Number of levels  in mixing_ratio

    ! A list of length ntype of gases whose volume mixing ratios have
    ! been provided
    integer :: icode(NMaxGases) = 0
    
   contains
     procedure :: allocate   => allocate_gas
     procedure :: deallocate => deallocate_gas
     procedure :: put        => put_gas
     procedure :: put_well_mixed => put_well_mixed_gas
     procedure :: scale      => scale_gas
     procedure :: set_units  => set_units_gas
     procedure :: assert_units => assert_units_gas
     procedure :: get        => get_gas
     procedure :: reverse    => reverse_gas
     procedure :: out_of_physical_bounds
  end type gas_type

contains


  !---------------------------------------------------------------------
  subroutine allocate_gas(this, ncol, nlev)

    use yomhook, only : lhook, dr_hook

    class(gas_type), intent(inout) :: this
    integer,         intent(in)    :: ncol, nlev

    real(jprb)          :: hook_handle

    if (lhook) call dr_hook('radiation_gas:allocate',0,hook_handle)

    call this%deallocate()

    allocate(this%mixing_ratio(ncol, nlev, NMaxGases))
    this%mixing_ratio = 0.0_jprb

    this%ncol = ncol
    this%nlev = nlev

    if (lhook) call dr_hook('radiation_gas:allocate',1,hook_handle)

  end subroutine allocate_gas


  !---------------------------------------------------------------------
  ! Deallocate memory and reset arrays
  subroutine deallocate_gas(this)

    use yomhook, only : lhook, dr_hook

    class(gas_type), intent(inout) :: this

    real(jprb)          :: hook_handle

    if (lhook) call dr_hook('radiation_gas:deallocate',0,hook_handle)

    if (allocated(this%mixing_ratio)) then
       deallocate(this%mixing_ratio)
    end if
    
    this%iunits = 0
    this%scale_factor = 0.0_jprb
    this%is_present = .false.
    this%is_well_mixed = .false.
    this%ntype = 0
    this%ncol = 0
    this%nlev = 0
    this%icode = 0

    if (lhook) call dr_hook('radiation_gas:deallocate',1,hook_handle)

  end subroutine deallocate_gas


  !---------------------------------------------------------------------
  ! Put gas mixing ratio corresponding to gas ID "igas" with units
  ! "iunits"
  subroutine put_gas(this, igas, iunits, mixing_ratio, scale_factor, &
       istartcol)

    use yomhook,        only : lhook, dr_hook
    use radiation_io,   only : nulerr, radiation_abort

    class(gas_type),      intent(inout) :: this
    integer,              intent(in)    :: igas
    integer,              intent(in)    :: iunits
    real(jprb),           intent(in)    :: mixing_ratio(:,:)
    real(jprb), optional, intent(in)    :: scale_factor
    integer,    optional, intent(in)    :: istartcol

    integer :: i1, i2

    real(jprb)                          :: hook_handle

    if (lhook) call dr_hook('radiation_gas:put',0,hook_handle)

    ! Check inputs
    if (igas <= IGasNotPresent .or. iunits > NMaxGases) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: provided gas ID (', &
           &   igas, ') must be in the range ', IGasNotPresent+1, ' to ', &
           &   NMaxGases
      call radiation_abort()
    end if
    if (iunits < IMassMixingRatio .or. iunits > IVolumeMixingRatio) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: provided gas units (', &
           &   iunits, ') must be in the range ', IMassMixingRatio, ' to ', &
           &   IVolumeMixingRatio
      call radiation_abort()
    end if

    if (.not. allocated(this%mixing_ratio)) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: attempt to put data to unallocated radiation_gas object'
      call radiation_abort()
    end if

    if (present(istartcol)) then
      i1 = istartcol
    else
      i1 = 1
    end if

    i2 = i1 + size(mixing_ratio,1) - 1

    if (i1 < 1 .or. i2 < 1 .or. i1 > this%ncol .or. i2 > this%ncol) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: attempt to put columns indexed ', &
           &   i1, ' to ', i2, ' to array indexed 1 to ', this%ncol
      call radiation_abort()
    end if

    if (size(mixing_ratio,2) /= this%nlev) then
      write(nulerr,'(a,i0,a)') &
           &  '*** Error: gas mixing ratio expected to have ', this%nlev, &
           &  ' levels'
      call radiation_abort()
    end if

    if (.not. this%is_present(igas)) then
      ! Gas not present until now
      this%ntype = this%ntype + 1
      this%icode(this%ntype) = igas
    end if
    this%is_present(igas) = .true.
    this%iunits(igas) = iunits
    this%is_well_mixed(igas) = .false.
    this%mixing_ratio(i1:i2,:,igas) = mixing_ratio

    if (present(scale_factor)) then
      this%scale_factor(igas) = scale_factor
    else
      this%scale_factor(igas) = 1.0_jprb
    end if

    if (lhook) call dr_hook('radiation_gas:put',1,hook_handle)

  end subroutine put_gas


  !---------------------------------------------------------------------
  ! Put well-mixed gas mixing ratio corresponding to gas ID "igas"
  ! with units "iunits"
  subroutine put_well_mixed_gas(this, igas, iunits, mixing_ratio, &
       scale_factor, istartcol, iendcol)

    use yomhook,        only : lhook, dr_hook
    use radiation_io,   only : nulerr, radiation_abort

    class(gas_type),      intent(inout) :: this
    integer,              intent(in)    :: igas
    integer,              intent(in)    :: iunits
    real(jprb),           intent(in)    :: mixing_ratio
    real(jprb), optional, intent(in)    :: scale_factor
    integer,    optional, intent(in)    :: istartcol, iendcol

    real(jprb)                          :: hook_handle

    integer :: i1, i2

    if (lhook) call dr_hook('radiation_gas:put_well_mixed',0,hook_handle)

    ! Check inputs
    if (igas <= IGasNotPresent .or. igas > NMaxGases) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: provided gas ID (', &
           &   igas, ') must be in the range ', IGasNotPresent+1, ' to ', &
           &   NMaxGases
      call radiation_abort()
    end if
    if (iunits < IMassMixingRatio .or. iunits > IVolumeMixingRatio) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: provided gas units (', &
           &   iunits, ') must be in the range ', IMassMixingRatio, ' to ', &
           &   IVolumeMixingRatio
      call radiation_abort()
    end if

    if (.not. allocated(this%mixing_ratio)) then
      write(nulerr,'(a)') '*** Error: attempt to put well-mixed gas data to unallocated radiation_gas object'
      call radiation_abort()
    end if

    if (present(istartcol)) then
      i1 = istartcol
    else
      i1 = 1
    end if

    if (present(iendcol)) then
      i2 = iendcol
    else
      i2 = this%ncol
    end if

    if (i1 < 1 .or. i2 < 1 .or. i1 > this%ncol .or. i2 > this%ncol) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: attempt to put columns indexed ', &
           &   i1, ' to ', i2, ' to array indexed 1 to ', this%ncol
      call radiation_abort()
    end if

    if (.not. this%is_present(igas)) then
      ! Gas not present until now
      this%ntype = this%ntype + 1
      this%icode(this%ntype) = igas
    end if
    ! Map uses a negative value to indicate a well-mixed value
    this%is_present(igas)              = .true.
    this%iunits(igas)                  = iunits
    this%is_well_mixed(igas)           = .true.
    this%mixing_ratio(i1:i2,:,igas)    = mixing_ratio

    if (present(scale_factor)) then
      this%scale_factor(igas) = scale_factor
    else
      this%scale_factor(igas) = 1.0_jprb
    end if

    if (lhook) call dr_hook('radiation_gas:put_well_mixed',1,hook_handle)

  end subroutine put_well_mixed_gas


  !---------------------------------------------------------------------
  ! Scale gas concentrations, e.g. igas=ICO2 and set scale_factor=2 to
  ! double CO2.  Note that this does not perform the scaling
  ! immediately, but changes the scale factor for the specified gas,
  ! ready to be used in set_units_gas.

  subroutine scale_gas(this, igas, scale_factor, lverbose)

    use radiation_io, only : nulout

    class(gas_type),      intent(inout) :: this
    integer,              intent(in)    :: igas
    real(jprb),           intent(in)    :: scale_factor
    logical,    optional, intent(in)    :: lverbose

    if (scale_factor /= 1.0_jprb) then
      this%scale_factor(igas) = this%scale_factor(igas) * scale_factor
      if (present(lverbose)) then
        if (lverbose) then
          write(nulout,'(a,a,a,f0.3)') '  Scaling ', trim(GasName(igas)), &
               &  ' concentration by ', scale_factor
        end if
      end if
    end if

  end subroutine scale_gas


  !---------------------------------------------------------------------
  ! Scale the gas concentrations so that they have the units "iunits"
  ! and are therefore ready to be used by the gas optics model within
  ! ecRad with no further scaling.  The existing scale_factor for each
  ! gas is applied.  If "igas" is present then apply only to gas with
  ! ID "igas", otherwise to all gases. Optional argument scale_factor
  ! specifies scaling that any subsequent access would need to apply
  ! to get a dimensionless result (consistent with definition of
  ! gas_type). So say that your gas optics model requires gas
  ! concentrations in PPMV, specify iunits=IVolumeMixingRatio and
  ! scale_factor=1.0e-6. If the gas concentrations were currently
  ! dimensionless volume mixing ratios, then the values would be
  ! internally divided by 1.0e-6.
  recursive subroutine set_units_gas(this, iunits, igas, scale_factor)
    class(gas_type),      intent(inout) :: this
    integer,              intent(in)    :: iunits
    integer,    optional, intent(in)    :: igas
    real(jprb), optional, intent(in)    :: scale_factor    

    integer :: ig

    ! Scaling factor to convert from old to new
    real(jprb) :: sf

    ! New scaling factor to store inside the gas object
    real(jprb) :: new_sf

    if (present(scale_factor)) then
      ! "sf" is the scaling to be applied now to the numbers (and may
      ! be modified below), "new_sf" is the value to be stored along
      ! with the numbers, informing subsequent routines how much you
      ! would need to multiply the numbers by to get a dimensionless
      ! result.
      sf     = 1.0_jprb / scale_factor
      new_sf = scale_factor
    else
      sf     = 1.0_jprb
      new_sf = 1.0_jprb
    end if

    if (present(igas)) then
      if (this%is_present(igas)) then
        if (iunits == IMassMixingRatio &
             &   .and. this%iunits(igas) == IVolumeMixingRatio) then
          sf = sf * IGasMolarMass(igas) / IAirMolarMass
        else if (iunits == IVolumeMixingRatio &
             &   .and. this%iunits(igas) == IMassMixingRatio) then
          sf = sf * IAirMolarMass / IGasMolarMass(igas)
        end if
        sf = sf * this%scale_factor(igas)
        
        if (sf /= 1.0_jprb) then
          this%mixing_ratio(:,:,igas) = this%mixing_ratio(:,:,igas) * sf
        end if
        ! Store the new units and scale factor for this gas inside the
        ! gas object
        this%iunits(igas) = iunits
        this%scale_factor(igas) = new_sf
      end if
    else
      do ig = 1,this%ntype
        call this%set_units(iunits, igas=this%icode(ig), scale_factor=new_sf)
      end do
    end if

  end subroutine set_units_gas


  !---------------------------------------------------------------------
  ! Assert that gas mixing ratio units are "iunits", applying to gas
  ! with ID "igas" if present, otherwise to all gases. Otherwise the
  ! program will exit. Otional argument scale factor specifies any
  ! subsequent multiplication to apply; for PPMV one would use
  ! iunits=IVolumeMixingRatio and scale_factor=1.0e6.
  recursive subroutine assert_units_gas(this, iunits, igas, scale_factor)
    
    use radiation_io,   only : nulerr, radiation_abort    

    class(gas_type),      intent(in) :: this
    integer,              intent(in) :: iunits
    integer,    optional, intent(in) :: igas
    real(jprb), optional, intent(in) :: scale_factor    

    integer :: ig

    real(jprb) :: sf

    if (present(scale_factor)) then
      sf = scale_factor
    else
      sf = 1.0_jprb
    end if

    if (present(igas)) then
      if (this%is_present(igas)) then
        if (iunits /= this%iunits(igas)) then
          write(nulerr,'(a,a,a)') '*** Error: ', trim(GasName(igas)), &
               &  ' is not in the required units'
          call radiation_abort()
        else if (sf /= this%scale_factor(igas)) then
          write(nulerr,'(a,a,a,e12.4,a,e12.4)') '*** Error: ', GasName(igas), &
               &  ' scaling of ', this%scale_factor(igas), &
               &  ' does not match required ', sf
          call radiation_abort()
        end if
      end if
    else
      do ig = 1,this%ntype
        call this%assert_units(iunits, igas=this%icode(ig), scale_factor=sf)
      end do
    end if

  end subroutine assert_units_gas


  !---------------------------------------------------------------------
  ! Get gas mixing ratio corresponding to gas ID "igas" with units
  ! "iunits" and return as a 2D array of dimensions (ncol,nlev).  The
  ! array will contain zeros if the gas is not stored.
  subroutine get_gas(this, igas, iunits, mixing_ratio, scale_factor, &
       &   istartcol)

    use yomhook,        only : lhook, dr_hook
    use radiation_io,   only : nulerr, radiation_abort

    class(gas_type),      intent(in)  :: this
    integer,              intent(in)  :: igas
    integer,              intent(in)  :: iunits
    real(jprb),           intent(out) :: mixing_ratio(:,:)
    real(jprb), optional, intent(in)  :: scale_factor
    integer,    optional, intent(in)  :: istartcol

    real(jprb)                        :: sf
    integer                           :: i1, i2

    real(jprb)                        :: hook_handle

    if (lhook) call dr_hook('radiation_gas:get',0,hook_handle)

    if (present(scale_factor)) then
      sf = scale_factor
    else
      sf = 1.0_jprb
    end if

    if (present(istartcol)) then
      i1 = istartcol
    else
      i1 = 1
    end if

    i2 = i1 + size(mixing_ratio,1) - 1

    if (i1 < 1 .or. i2 < 1 .or. i1 > this%ncol .or. i2 > this%ncol) then
      write(nulerr,'(a,i0,a,i0,a,i0)') '*** Error: attempt to get columns indexed ', &
           &   i1, ' to ', i2, ' from array indexed 1 to ', this%ncol
      call radiation_abort()
    end if

    if (size(mixing_ratio,2) /= this%nlev) then
      write(nulerr,'(a,i0,a)') &
           &  '*** Error: gas destination array expected to have ', this%nlev, &
           &  ' levels'
      call radiation_abort()
    end if

    if (.not. this%is_present(igas)) then
      mixing_ratio = 0.0_jprb
    else 
      if (iunits == IMassMixingRatio &
           &   .and. this%iunits(igas) == IVolumeMixingRatio) then
        sf = sf * IGasMolarMass(igas) / IAirMolarMass
      else if (iunits == IVolumeMixingRatio &
           &   .and. this%iunits(igas) == IMassMixingRatio) then
        sf = sf * IAirMolarMass / IGasMolarMass(igas)
      end if
      sf = sf * this%scale_factor(igas)
        
      if (sf /= 1.0_jprb) then
        mixing_ratio = this%mixing_ratio(i1:i2,:,igas) * sf
      else
        mixing_ratio = this%mixing_ratio(i1:i2,:,igas)
      end if
    end if

    if (lhook) call dr_hook('radiation_gas:get',1,hook_handle)

  end subroutine get_gas


  !---------------------------------------------------------------------
  ! Copy data to "gas_rev", reversing the height ordering of the gas
  ! data
  subroutine reverse_gas(this, istartcol, iendcol, gas_rev)

    class(gas_type)             :: this
    integer,        intent(in)  :: istartcol, iendcol
    type(gas_type), intent(out) :: gas_rev

    gas_rev%iunits = this%iunits
    gas_rev%scale_factor = this%scale_factor
    gas_rev%is_present = this%is_present
    gas_rev%is_well_mixed = this%is_well_mixed
    gas_rev%ntype = this%ntype
    gas_rev%ncol = this%ncol
    gas_rev%nlev = this%nlev
    gas_rev%icode = this%icode

    if (allocated(gas_rev%mixing_ratio)) deallocate(gas_rev%mixing_ratio)

    if (allocated(this%mixing_ratio)) then
      allocate(gas_rev%mixing_ratio(istartcol:iendcol,this%nlev,NMaxGases))
      gas_rev%mixing_ratio(istartcol:iendcol,:,:) &
           &  = this%mixing_ratio(istartcol:iendcol,this%nlev:1:-1,:)
    end if

  end subroutine reverse_gas

  !---------------------------------------------------------------------
  ! Return .true. if variables are out of a physically sensible range,
  ! optionally only considering columns between istartcol and iendcol
  function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad)

    use yomhook,          only : lhook, dr_hook
    use radiation_config, only : out_of_bounds_3d

    class(gas_type),   intent(inout) :: this
    integer,  optional,intent(in) :: istartcol, iendcol
    logical,  optional,intent(in) :: do_fix
    logical                       :: is_bad

    logical    :: do_fix_local

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_gas:out_of_physical_bounds',0,hook_handle)

    if (present(do_fix)) then
      do_fix_local = do_fix
    else
      do_fix_local = .false.
    end if

    is_bad = out_of_bounds_3d(this%mixing_ratio, 'gas%mixing_ratio', &
         &                    0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol)

    if (lhook) call dr_hook('radiation_gas:out_of_physical_bounds',1,hook_handle)

  end function out_of_physical_bounds
  
end module radiation_gas