! 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