ecearth3/sources/tm5mp/proj/domino/MAmfRead.F90

! First include the set of model-wide compiler flags 
#include "tm5.inc"

! Macro
#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN if (status/=0) then; TRACEBACK; return; end if
#define IF_ERROR_RETURN if (status> 0) then; TRACEBACK; return; end if

module MAmfRead

  !
  ! Module to load the height-dependent air-mass factor lookup table
  !
  !  subroutine ReadAmfLut
  !
  !                             Henk Eskes, KNMI, Sept 2015

  use GO, only : gol, goErr

  implicit none

  private

  public :: ReadAmfLut

  ! The LUT netcdf file contains:
  !
  ! float mu(mu)             ; long_name = "Cosine of viewing zenith angle" ; units = "1"
  ! float mu0(mu0)           ; long_name = "Cosine of solar zenith angle"   ; units = "1"
  ! int dphi(dphi)           ; long_name = "Relative azimuth angle"         ; units = "degree"
  ! float vza(mu)            ; long_name = "Viewing zenith angle"           ; units = "degree" (not used)
  ! float sza(mu0)           ; long_name = "Solar zenith angle"             ; units = "degree" (not used)
  ! float p(p)               ; longname = "Pressure Levels"                 ; units = "hPa" 
  ! int p_surface(p_surface) ; long_name = "Surface pressure levels"        ; units = "hPa" 
  ! float albedo(albedo)     ; long_name = "Surface Albedo"                 ; units = "1"
  !
  ! float amf(albedo, p_surface, p, dphi, mu0, mu) ;
  !                            (note: indices in reverse order as in F90)
  !                            long_name = "Box air mass factor"            ; units = "1" 
  !
  ! Note: 
  !     albedo, dphi, mu, mu0 : increasing
  !     p, p_surface          : decreasing

  character(len=*), parameter  ::  mname = 'MAmfRead'

contains


  subroutine ReadAmfLut ( amf_filename, divideByAMFGeo, amfLut, status )
    !=======================================================================
    !
    ! ReadAmfLut:    read lookup table for AMF calculations
    !                updated netcdf format lookup table, 7 Sept 2015 
    !                     from Alba Lorente, Folkert Boersma
    !
    ! input:
    !   amf_filename = path+name of the file containing the AMF lookup table
    !   divideByAMFGeo = provide inital division by AMF_geo ? (t/f)  
    !      (The latest AMF files contain AMF / AMF_geo, and the flag should be put to .false.)
    ! output:
    !   amfLut = structure with AMF lookup table
    !   status = status ( 0 means OK )
    !
    ! The routine will allocate the arrays in "anfLut",
    !     using the dimensions in the file  
    !                                            Henk Eskes, KNMI, sept 2015
    !=======================================================================

    use MDF,      only : MDF_Open,      MDF_NETCDF,  MDF_READ 
    use MDF,      only : MDF_Inq_VarID, MDF_Get_Var, MDF_Close 
    use MDF,      only : MDF_Inquire_Dimension, MDF_Inq_DimID

    use MTAmf,    only : TAmfLut, AllocateAmfLut

    use MAmfGet,  only : GetAmfGeo

    implicit none

    ! in/out
    character(*), intent(in)     ::  amf_filename
    logical, intent(in)          ::  divideByAMFGeo
    type(TAmfLut), intent(inout) ::  amfLut    
    integer, intent(out)         ::  status

    ! local
    character(len=*), parameter ::  rname = mname//'/ReadAmfLut'

    integer      :: fid
    integer      :: dimid, id_mu, id_mu0, id_dphi, id_p, id_psurf, id_alb, id_amf 
    integer      :: nmu_file, nmu0_file, ndphi_file, np_file, npsurf_file, nalbedo_file 
    integer      :: imu, imu0

    integer, dimension(:), allocatable :: dphi_int, psurf_int

    real :: amfGeo

    ! start code

    call MDF_Open( trim(amf_filename), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN

    ! get the 6 dimensions
    call MDF_Inq_DimID( fid, 'mu', dimid, status )                                                       
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=nmu_file )                                  
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'mu0', dimid, status )                                                       
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=nmu0_file )                                  
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'dphi', dimid, status )                                                       
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=ndphi_file )                                  
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'p', dimid, status )                                                       
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=np_file )                                  
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'p_surface', dimid, status )                                                       
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=npsurf_file )                                  
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'albedo', dimid, status )                                                       
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=nalbedo_file )                                  
    IF_NOTOK_RETURN

    ! Allocate storage space
    call AllocateAmfLut ( nmu_file, nmu0_file, ndphi_file, np_file, npsurf_file, nalbedo_file, amfLut )

    ! Define variables
    CALL MDF_Inq_VarID( fid, 'mu', id_mu, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'mu0', id_mu0, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'dphi', id_dphi, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'p', id_p, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'p_surface', id_psurf, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'albedo' , id_alb, status )
    IF_ERROR_RETURN

    CALL MDF_Inq_VarID( fid, 'amf', id_amf, status )
    IF_ERROR_RETURN

    ! Read variables
    CALL MDF_Get_Var( fid, id_mu, amfLut%mu, status )      ! cos(vza)
    IF_NOTOK_RETURN

    CALL MDF_Get_Var( fid, id_mu0, amfLut%mu0, status )    ! cos(sza)
    IF_NOTOK_RETURN

    allocate ( dphi_int(ndphi_file) )
    CALL MDF_Get_Var( fid, id_dphi, dphi_int, status )
    IF_NOTOK_RETURN
    amfLut%azi(:) = real ( dphi_int(:) )                   ! azi/dphi in degree
    deallocate ( dphi_int )

    CALL MDF_Get_Var( fid, id_p, amfLut%pres, status )     ! layer pressures
    IF_NOTOK_RETURN

    allocate ( psurf_int(npsurf_file) )
    CALL MDF_Get_Var( fid, id_psurf, psurf_int, status )
    IF_NOTOK_RETURN
    amfLut%spres(:) = real ( psurf_int(:) )                ! surface pressures
    deallocate ( psurf_int )

    CALL MDF_Get_Var( fid, id_alb, amfLut%albedo, status ) ! albedo values
    IF_NOTOK_RETURN

    CALL MDF_Get_Var( fid, id_amf, amfLut%amf, status )    ! the AMF lookup table
    IF_NOTOK_RETURN

    ! Close file
    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN

    ! The LUT file contains either AMF or AMF/AMFgeo
    if ( divideByAMFGeo ) then
       ! The AMF is pre-divided by the geometrical AMF, which makes it more smooth 
       !    in order to reduce interpolation errors
       do imu = 1, amfLut%mmu 
          do imu0 = 1, amfLut%mmu0 
             call GetAmfGeo ( amfLut%mu(imu), amfLut%mu0(imu0), amfGeo )
             amfLut%amf(imu,imu0,:,:,:,:) = amfLut%amf(imu,imu0,:,:,:,:) / amfGeo
             ! print*,' mu, mu0, amfgeo ',amfLut%mu(imu), amfLut%mu0(imu0), amfGeo, amfLut%amf(imu,imu0,4,150,1,7)
          end do
       end do
    end if

  end subroutine ReadAmfLut


end module MAmfRead