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

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

! #ifdef TROPNLL2DP

! #define IF_NOTOK_RETURN(notok_string) call fortranlog(notok_string,len(notok_string),2); status=1; return
! #define WRITE_WARNING(notok_string)   call fortranlog(notok_string,len(notok_string),2)

! #else

#define IF_NOTOK_RETURN(notok_string) write (*,'(a)') notok_string; status=1; return
#define WRITE_WARNING(notok_string)   write (*,'(a)') notok_string

! #endif


module MAmfGet
  !
  ! Module to return the height-dependent air-mass factor for a given geometry
  !
  !           Henk Eskes, KNMI, sept 2015
  !
  use physconstants,   only : pi

  ! --- Error message codes ---
  use pqf_module

  implicit none

  private

  public :: GetAmf, GetAmfGeo, GetAmf_vectorised
  public :: FitWindowCentre

  ! --------------------------------------------------------------
  ! FitWindowCentre: central wavelength of fit window
  ! --------------------------------------------------------------

  real, parameter  :: FitWindowCentre = 437.5    ! nm, used in DOMINO-1 & 2
  ! real, parameter  :: FitWindowCentre = 440.0    ! nm

  ! for error messages
  character(256) :: errmessage

contains


  subroutine GetAmfGeo ( cos_vza, cos_sza, amfGeo )
    !=======================================================================
    !
    ! GetAmfGeo:    Return simple geometric AMF
    !   
    ! input:
    !   cos_vza  : cos viewing angle (degrees) [-60..60]
    !   cos_sza  : cos solar zenith angle (degrees) [-90..90]
    !   amfGeo   : geometrical air mass factor
    !
    ! geometrical air mass factor
    !
    !                                             Henk Eskes, KNMI, aug 2016
    !=======================================================================

    implicit none

    ! in/out
    real, intent(in)      :: cos_vza, cos_sza
    real, intent(out)     :: amfGeo    

    ! local

    amfGeo = 1.0 / cos_vza + 1.0 / cos_sza

  end subroutine GetAmfGeo


  subroutine GetAmfGeo_Leue ( cos_vza, cos_sza, amfGeo )
    !=======================================================================
    !
    ! GetAmfGeo_Leue:    Return geometric AMF
    !   
    ! input:
    !   cos_vza  : cos viewing angle (degrees) [-60..60]
    !   cos_sza  : cos solar zenith angle (degrees) [-90..90]
    !   amfGeo   : geometrical air mass factor
    !
    ! geometrical air mass factor, based on a
    !    simple form with corrections for large SZAs,
    !    from Carsten Leue, thesis
    !
    !                                            Henk Eskes, KNMI, sept 2015
    !=======================================================================

    use physconstants,  	only : Earth_rad     !  pdiff2moleccm2

    implicit none

    ! in/out
    real, intent(in)      :: cos_vza, cos_sza
    real, intent(out)     :: amfGeo    

    ! local
    real, parameter   :: htopatm = 60.0   ! (km)  height of top of atmosphere
    real, parameter   :: eps = htopatm*1000.0/Earth_rad 
    real              :: xh

    xh = sqrt( cos_sza*cos_sza + eps*eps + 2.*eps )
    amfGeo = 1.0/cos_vza + ( xh - cos_sza ) / eps    

  end subroutine GetAmfGeo_Leue



  subroutine GetAmfFactors ( amfLutValue, amfLutArray, iLoc, fFactor )
    !==========================================================================
    !
    ! GetAmfIncreasing: Interpolate in increasing array
    !   
    ! input:
    !   amfLutValue   : value for which interpolation factors are determined
    !   amfLutArray   : interpolation array
    !
    ! input:
    !   iLoc          : index in amfLutArray where value just below amfLutValue
    !   fFactor       : interpolation factors
    !
    ! bug fix:          bordeline cases, e.g. azi == 0 (feb 2016)
    !
    !                                           Martien de Haan, KNMI, feb 2016
    !==========================================================================

    implicit none

    ! in/out
    real, intent(in)                :: amfLutValue
    real, dimension(:), intent(in)  :: amfLutArray
    integer, intent(out)            :: iLoc
    real, dimension(2), intent(out) :: fFactor

    integer           :: dum(1)
	
    ! --> Determine whether amfLutArray is increasing or decreasing
    if (amfLutArray(2) > amfLutArray(1)) then
       ! --> Increasing
       dum = maxloc ( amfLutArray, amfLutArray < amfLutValue ) 
    else
       dum = minloc ( amfLutArray, amfLutArray >= amfLutValue )
    end if
	
    iLoc = dum(1)
	
    ! --> Border case: amfLutValue is lowest value in array. Set to first index
    if (iLoc == 0) iLoc = 1
	
    ! --> Set interpolation factors
    fFactor(1) = ( amfLutArray(iLoc+1) - amfLutValue ) / ( amfLutArray(iLoc+1) - amfLutArray(iLoc) )
    fFactor(2) = 1.0 - fFactor(1)

  end subroutine GetAmfFactors


  subroutine GetAmf ( pres, azi1, view, sza, albedo, spres, amfLut, amf, status )
    !=======================================================================
    !
    ! GetAmf:    Return AMF (divided by the geometrical AMF) 
    !            interpolated from the lookup table
    ! input:
    !   pres   : model pressure level in hPa
    !   azi1   : azimuth angle (degrees) [-360..360]
    !   view   : viewing angle (degrees) [-60..60]
    !   sza    : solar zenith angle (degrees) [-90..90]
    !   albedo : ground albedo el. [0..1]
    !   spres  : surface pressure in hPa
    !   amfLut : a structure containing the AMF array and description of the axis
    !   status : 0 = OK, > 0 is error
    !
    !                                 Henk Eskes, KNMI, nov 1999 - sept 2015
    !=======================================================================

    use MTAmf,        only : TAmfLut

    implicit none

    ! in/out
    real,intent(in)            ::  pres, azi1, view, sza, albedo, spres
    type(TAmfLut), intent(in)  ::  amfLut    
    real,intent(out)           ::  amf
    integer, intent(out)       ::  status

    ! local
    character(*), parameter :: rname = "GetAmf"

    real              :: pres_tmp, spres_tmp
    integer           :: imu, imu0, iazi, ipres, ispres, ialbedo
    integer           :: i1, i2, i3, i4, i5, i6
    integer           :: dum(1)
    real,dimension(2) :: xmu, xmu0, xazi, xpres, xspres, xalbedo
    real              :: azi, mu, mu0
    !
    !    if ( pres > spres )   write(*,'(a,f,2x,f)')'WARNING GetAmf: pressure > surface pressure ',pres,spres
    status = 0
 
    ! the routine assumes the input varaibles have been range checked

    ! Azimuth fix (Ronald van der A), azi = abs(180.0-azi1)
    ! For the new AMF lookup table (August 2015, Alba Lorente) this fix is no longer needed (Henk Eskes) 
    azi = abs(azi1) 
    if ( azi > 180.0 )  azi = abs( 360.0 - azi ) 
    mu  = cos(pi*view*(1./180.0))
    mu0 = cos(pi*sza*(1./180.0))

    ! Additional range checks: is the point inside the LUT domain ?
    if ( (mu < amfLut%mu(1)) .or. (mu > amfLut%mu(amfLut%mmu)) ) then
       write ( errmessage, * ) 'WARNING GetAmf: cos(vza) out of range, value = ',mu
       IF_NOTOK_RETURN( trim(errmessage) ) 
    end if
    if ( (mu0 < amfLut%mu0(1)) .or. (mu0 > amfLut%mu0(amfLut%mmu0)) ) then
       write ( errmessage, * ) 'WARNING GetAmf: cos(sza) out of range, value = ',mu0
       IF_NOTOK_RETURN( trim(errmessage) ) 
    end if
    if ( (azi < amfLut%azi(1)) .or. (azi > amfLut%azi(amfLut%mazi)) ) then
       write ( errmessage, * ) 'WARNING GetAmf: azi out of range, value = ',azi
       IF_NOTOK_RETURN( trim(errmessage) ) 
    end if
    if ( (albedo < amfLut%albedo(1)) .or. (albedo > amfLut%albedo(amfLut%malbedo)) ) then
       write ( errmessage, * ) 'WARNING GetAmf: albedo out of range, value = ',albedo
       IF_NOTOK_RETURN( trim(errmessage) ) 
    end if

    ! force surface pressure to be within the LUT range 
    spres_tmp = spres
    if( (spres_tmp < amfLut%spres(amfLut%mspres)) ) then
       ! code die clipt, RuudDirksen 14 juni 2009        
       ! write ( errmessage,'(a,f8.2,a,f8.2)') 'WARNING: overwrite input surface pressure ',spres,'  with first AMF surface pressure ',amfLut%spres(amfLut%mspres)
       ! WRITE_WARNING( trim(errmessage) )
       spres_tmp = 1.00001 * amfLut%spres(amfLut%mspres)
    end if
    if ( spres_tmp > amfLut%spres(1) ) then
       spres_tmp = 0.99999 * amfLut%spres(1)       
       !write ( errmessage, * ) 'WARNING: Overwriting surface pressure ', spres, &
       !     ' with AMF first surf p level ', amfLut%spres(1) 
       !WRITE_WARNING( trim(errmessage) )
    end if

    ! force pressure to be within range
    pres_tmp = pres
    if ( pres < amfLut%pres(amfLut%mpres) ) then
       pres_tmp = 1.00001 * amfLut%pres(amfLut%mpres)
    end if
    if ( pres > amfLut%pres(1) ) then
       pres_tmp = 0.99999 * amfLut%pres(1)       
       write ( errmessage, * ) 'WARNING: Overwriting Model pressure ', pres, &
            ' with AMF first pressure level  ', amfLut%pres(1) 
       WRITE_WARNING( trim(errmessage) )
    end if

    !
    ! Determine reference index (imu,imu0,iazi,ipres,ispres,ialbedo)
    !
    ! --> Cos viewing angle mu: (increasing)
    call GetAmfFactors(mu, amfLut%mu, imu, xmu)

    ! --> Cos solar zenith angle: (increasing)
    call GetAmfFactors(mu0, amfLut%mu0, imu0, xmu0)

    ! --> Azimuth angle: (increasing)
    call GetAmfFactors(azi, amfLut%azi, iazi, xazi)

    ! --> Pressure: (decreasing)
    call GetAmfFactors(pres_tmp, amfLut%pres, ipres, xpres)
    
    ! --> Surface pressure: (decreasing)
    call GetAmfFactors(spres_tmp, amfLut%spres, ispres, xspres)

    ! --> Ground albedo: (increasing)
    call GetAmfFactors(albedo, amfLut%albedo, ialbedo, xalbedo)

    ! Linear interpolation
    amf = 0.0
    do i1 = 1, 2
       do i2 = 1, 2
          do i3 = 1, 2
             do i4 = 1, 2
                do i5 = 1, 2
                   do i6 = 1, 2
                      amf = amf + xmu(i1)*xmu0(i2)*xazi(i3)*xpres(i4)*xspres(i5)*xalbedo(i6)* &
                           amfLut%amf ( imu+i1-1, imu0+i2-1, iazi+i3-1, ipres+i4-1, ispres+i5-1, ialbedo+i6-1 )
                   end do
                end do
             end do
          end do
       end do
    end do
    !
  end subroutine GetAmf



  subroutine GetAmf_vectorised ( nObs, lm, pres1, azi1, view1, sza1, albedo, spres1, amfLut, amf, flag, debug ) 

    !=======================================================================
    !
    ! GetAmf_vectorised:    Return AMF (divided by the geometrical AMF) 
    !                       interpolated from the lookup table
    !                       for all observations at the same time
    !
    !   call GetAMF_vectorised(phybrid(:,:),aza(:),vza(:),sza(:),albclear(:),surfpres(:),amfclearrel(:,:))
    !
    ! input:
    !   nObs   : number of observations
    !   lm     : number of vertical layers
    !   pres1  : model pressure level in hPa, dimension: (nObs,lm)
    !   azi1   : azimuth angle (degrees) [-360..360]
    !   view1  : viewing angle (degrees) [-60..60]
    !   sza1   : solar zenith angle (degrees) [-90..90]
    !   albedo : ground albedo el. [0..1]
    !   spres1 : surface pressure in hPa
    !   amfLut : a structure containing the AMF array and description of the axis
    !   flag   : error flags for the individual observations
    !   debug  : provide debug info to screen 
    !
    ! output:
    !   amf    : air mass factor for all observations and on all layers
    !   flag   : error flags for the individual observations
    !
    !   note: all these quantities are arrays of size (nr of observations) or 
    !         (nr of observations, nr of levels)
    !
    !                                 Henk Eskes, KNMI, nov 1999 - sept 2015
    !                         Vectorized version by Ruud Dirksen 27 nov 2007
    !=======================================================================

    use MTAmf,       only : TAmfLut

    implicit none

    ! in/out
    integer, intent(in)                   :: nObs, lm
    real, dimension(:,:), intent(in)      :: pres1
    real, dimension(:), intent(in)        :: spres1, view1, azi1, sza1, albedo
    type(TAmfLut), intent(in)             :: amfLut    
    real, dimension(:,:), intent(out)     :: amf
    integer, dimension(:), intent(inout)  :: flag
    logical, intent(in)                   :: debug

    ! local
    character(*), parameter         :: rname = "GetAmf_vectorised"

    real,dimension(:),allocatable   :: azi, mu, mu0, spres
    real,dimension(:,:),allocatable :: pres

    integer                         :: i1, i2, i3, i4, i5, i6
    integer                         :: dum(1), iObs, l
    integer                         :: ipres, ispres, imu, iazi, imu0, ialbedo
    real,dimension(2)               :: xpres, xspres, xmu, xazi, xmu0, xalbedo

    ! start code

    if ( debug ) then
       print *, ' '
       print *, 'AMF LUT axes'
       print '(a,24f6.3)',  '  amf cvza   ',amfLut%mu
       print '(a,24f6.3)',  '  amf csza   ',amfLut%mu0
       print '(a,24f5.0)',  '  amf azi    ',amfLut%azi
       print '(a,240f6.0)', '  amf pres   ',amfLut%pres
       print '(a,24f6.0)',  '  amf spres  ',amfLut%spres
       print '(a,50f6.3)',  '  amf albedo ',amfLut%albedo
    end if

    ! initialise AMF
    amf(:,:) = 0.0

    ! helper arrays
    allocate ( azi(nObs) )
    allocate ( mu(nObs) )
    allocate ( mu0(nObs) )
    allocate ( spres(nObs) )
    allocate ( pres(nObs,lm) )

    ! make copies of the surface pressure and level pressures
    spres(:)   = spres1(1:nObs)
    pres(:,:)  = pres1(1:nObs,1:lm)
    !    if ( pres > spres )   write(*,*)'WARNING GetAmf: pressure > surface pressure'

    ! routine assumes the ranges of the input values have been checked

    ! Azimuth fix (Ronald van der A), azi = abs(180.0-azi1)
    ! For the new AMF lookup table (August 2015, Alba Lorente) this fix is no longer needed (Henk Eskes) 
    do iObs = 1, nObs
       azi(iObs) = abs(azi1(iObs))
       if ( azi(iObs) .gt. 180.0 ) azi(iObs) = abs(360.0 - azi(iObs))
    end do
    do iObs = 1, nObs
       mu(iObs) = cos(pi*view1(iObs)*(1./180.0))
       mu0(iObs) = cos(pi*sza1(iObs)*(1./180.0))
    end do

    !$OMP parallel do private(xazi,xview,xmu0,xalbedo,xspres,xpres) schedule(static)

    observationLoop: do iObs = 1, nObs

       ! skip to next iObs if (flag = error),  accept warnings
       if ( iand(flag(iObs), 255) /=  0 ) cycle

       ! Make sure the inputs are inside the LUT domain, 
       ! write warning and set error flag
       if ( (mu(iObs) < amfLut%mu(1)) .or. (mu(iObs) > amfLut%mu(amfLut%mmu)) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: cos(vza) out of range, value = ',mu(iObs),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
          flag(iObs) = PQF_E_LUT_RANGE_ERROR
       end if
       if ( (mu0(iObs) < amfLut%mu0(1)) .or. (mu0(iObs) > amfLut%mu0(amfLut%mmu0)) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: cos(sza) out of range, value = ',mu0(iObs),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
          flag(iObs) = PQF_E_LUT_RANGE_ERROR
       end if
       if ( (azi(iObs) < amfLut%azi(1)) .or. (azi(iObs) > amfLut%azi(amfLut%mazi)) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: azi out of range, value = ',azi(iObs),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
          flag(iObs) = PQF_E_LUT_RANGE_ERROR
       end if
       if ( (albedo(iObs) < amfLut%albedo(1)) .or. (albedo(iObs) > amfLut%albedo(amfLut%malbedo)) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: albedo out of range, value = ',albedo(iObs),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
          flag(iObs) = PQF_E_LUT_RANGE_ERROR
       end if
       if ( spres(iObs) < amfLut%spres(amfLut%mspres) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: surface pressure out of range, value = ',spres(iObs),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
          spres(iObs) = 1.00001 * amfLut%spres(amfLut%mspres)
       end if
       if ( spres(iObs) > amfLut%spres(1) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: surface pressure out of range, value = ',spres(iObs),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
          spres(iObs) = 0.99999 * amfLut%spres(1)
       end if
       if ( maxval(pres(iObs,:)) > amfLut%pres(1) ) then
          write ( errmessage, * ) 'WARNING GetAmf_vectorised: pressure out of range, value = ',maxval(pres1(iObs,:)),'  iObs = ',iObs
          WRITE_WARNING( trim(errmessage) ) 
       end if

       ! skip to next iObs if (flag = error),  accept warnings
       if ( iand(flag(iObs), 255) /=  0 ) cycle

       ! force surface pressure to be within the LUT range 
       ! spres = spres1
       ! where ( spres < amfLut%spres(amfLut%mspres) )  spres = 1.00001 * amfLut%spres(amfLut%mspres)
       ! where ( spres > amfLut%spres(1)             )  spres = 0.99999 * amfLut%spres(1)

       ! force pressure to be within the LUT range 
       where ( pres(iObs,:) < amfLut%pres(amfLut%mpres) )  pres(iObs,:) = 1.00001 * amfLut%pres(amfLut%mpres)
       where ( pres(iObs,:) > amfLut%pres(1)            )  pres(iObs,:) = 0.99999 * amfLut%pres(1)

       ! Determine reference index (ipres,iview,iazi,imu0,ialbedo)

       ! --> Cos viewing angle mu: (increasing)
       call GetAmfFactors(mu(iObs), amfLut%mu, imu, xmu)

       ! --> Cos solar zenith angle: (increasing)
       call GetAmfFactors(mu0(iObs), amfLut%mu0, imu0, xmu0)

       ! --> Azimuth angle: (increasing)
       call GetAmfFactors(azi(iObs), amfLut%azi, iazi, xazi)
	   
       ! --> Surface pressure: (decreasing)
       call GetAmfFactors(spres(iObs), amfLut%spres, ispres, xspres)

       ! --> Ground albedo: (increasing)
       call GetAmfFactors(albedo(iObs), amfLut%albedo, ialbedo, xalbedo)

       ! now loop over vertical levels
       levelLoop: do l = 1, lm

          ! Pressure: (decreasing)
          call GetAmfFactors(pres(iObs,l), amfLut%pres, ipres, xpres)

          ! Linear interpolation
          do i1 = 1, 2
             do i2 = 1, 2
                do i3 = 1, 2
                   do i4 = 1, 2
                      do i5 = 1, 2
                         do i6 = 1, 2
                            amf(iObs,l) = amf(iObs,l) + xmu(i1)*xmu0(i2)*xazi(i3)*xpres(i4)*xspres(i5)*xalbedo(i6)* &
                                 amfLut%amf ( imu+i1-1, imu0+i2-1, iazi+i3-1, ipres+i4-1, ispres+i5-1, ialbedo+i6-1 )
                         end do
                      end do
                   end do
                end do
             end do
          end do

       end do levelLoop  ! l


       if ( debug .and. (iObs == max(nObs/2,1)) ) then
          print '(a)', ' '
          print '(a,i6)', 'iobs = ',iObs
          print '(a,3f8.2)', '  vza,sza,aza1  = ',view1(iObs),sza1(iObs),azi1(iObs)
          print '(a,3f8.2)', '  cvza,csza,aza = ',mu(iObs),mu0(iObs),azi(iObs)
          print '(a,3f8.2)', '  spres1,spres,albclear = ',spres1(iObs),spres(iObs),albedo(iObs)
          print '(a,3i4,a,2i4)', 'imu, imu0, iazi, x, ispres, ialbedo = ', imu, imu0, iazi, ' X ', ispres, ialbedo
          print '(a,200f6.0)', '  amfpres = ',amfLut%pres
          print '(a,500f6.3)', '  amflut1 = ',amflut%amf(imu,  imu0,  iazi,  :,ispres,  ialbedo  )
          print '(a,500f6.3)', '  amflut2 = ',amflut%amf(imu+1,imu0+1,iazi+1,:,ispres+1,ialbedo+1)
          print '(a,50f6.0)', '  phybrid = ',pres(iObs,:)
          print '(a,50f6.3)', '  amf     = ',amf(iObs,:)
          print '(a)', ' '
       end if


    end do observationLoop   ! iObs

    !$OMP end parallel do

    if ( allocated(azi) )    deallocate ( azi )
    if ( allocated(mu) )     deallocate ( mu )
    if ( allocated(mu0) )    deallocate ( mu0 )
    if ( allocated(pres) )   deallocate ( pres )
    if ( allocated(spres) )  deallocate ( spres )

  end subroutine GetAmf_vectorised


end module MAmfGet