! 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