.. _program_listing_file_src_modules_Iqbal_iswr.cpp:

Program Listing for File Iqbal_iswr.cpp
=======================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_modules_Iqbal_iswr.cpp>` (``src/modules/Iqbal_iswr.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   //
   // Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel
   // modular unstructured mesh based approach for hydrological modelling
   // Copyright (C) 2018 Christopher Marsh
   //
   // This file is part of Canadian Hydrological Model.
   //
   // Canadian Hydrological Model is free software: you can redistribute it and/or
   // modify
   // it under the terms of the GNU General Public License as published by
   // the Free Software Foundation, either version 3 of the License, or
   // (at your option) any later version.
   //
   // Canadian Hydrological Model is distributed in the hope that it will be useful,
   // but WITHOUT ANY WARRANTY; without even the implied warranty of
   // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   // GNU General Public License for more details.
   //
   // You should have received a copy of the GNU General Public License
   // along with Canadian Hydrological Model.  If not, see
   // <http://www.gnu.org/licenses/>.
   //
   
   #include "Iqbal_iswr.hpp"
   REGISTER_MODULE_CPP(Iqbal_iswr);
   
   Iqbal_iswr::Iqbal_iswr(config_file cfg)
           :module_base("Iqbal_iswr", parallel::data, cfg)
   {
   
       depends("t");
       depends("rh");
       depends("cloud_frac");
       depends("solar_el");
   
       provides("iswr_direct_no_slope");
       provides("iswr_diffuse_no_slope");
   
       provides("atm_trans");
   }
   
   Iqbal_iswr::~Iqbal_iswr()
   {
   
   }
   
   void Iqbal_iswr::run(mesh_elem &face)
   {
       double pressure = mio::Atmosphere::stdAirPressure(face->get_z());//101325.0;
       double altitude = face->get_z();
       double sun_elevation = (*face)["solar_el"_s];
       sun_elevation = sun_elevation < 0? 0. : sun_elevation;
   
       if (sun_elevation < 3)
       {
           (*face)["iswr_direct_no_slope"_s]=0;
           (*face)["iswr_diffuse_no_slope"_s]=0;
           return;
       }
   
       double ta = (*face)["t"_s]+273.15;
       double rh = (*face)["rh"_s]/100.0;
       double R_toa = 1375;
       double R_direct=0;
       double R_diffuse=0;
       double ground_albedo = 0.1;
   
   
       //these pow cost us a lot here, but replacing them by fastPow() has a large impact on accuracy (because of the exp())
       const double olt = 0.32;   //ozone layer thickness (cm) U.S.standard = 0.34 cm
       const double w0 = 0.9;     //fraction of energy scattered to total attenuation by aerosols (Bird and Hulstrom(1981))
       const double fc = 0.84;    //fraction of forward scattering to total scattering (Bird and Hulstrom(1981))
       const double alpha = 1.3;  //wavelength exponent (Iqbal(1983) p.118). Good average value: 1.3+/-0.5. Related to the size distribution of the particules
       const double beta = 0.03;  //amount of particules index (Iqbal(1983) p.118). Between 0 & .5 and above.
       const double zenith = 90. - sun_elevation; //this is the TRUE zenith because the elevation is the TRUE elevation
       const double cos_zenith = cos(zenith*mio::Cst::to_rad); //this uses true zenith angle
   
       // relative optical air mass Young (1994), see http://en.wikipedia.org/wiki/Airmass
       //const double mr = 1. / (cos_zenith + 0.50572 * pow( 96.07995-zenith , -1.6364 )); //pbl: this should use apparent zenith angle, and we only get true zenith angle here...
       // relative optical air mass, Young, A. T. 1994. Air mass and refraction. Applied Optics. 33:1108–1110.
       const double mr = ( 1.002432*cos_zenith*cos_zenith + 0.148386*cos_zenith + 0.0096467) /
                         ( cos_zenith*cos_zenith*cos_zenith + 0.149864*cos_zenith*cos_zenith
                           + 0.0102963*cos_zenith +0.000303978);
   
       // actual air mass: because mr is applicable for standard pressure
       // it is modified for other pressures (in Iqbal (1983), p.100)
       // pressure in Pa
       const double ma = mr * (pressure/101325.);
   
       // the equations for all the transmittances of the individual atmospheric constituents
       // are from Bird and Hulstrom (1980, 1981) and can be found summarized in Iqbal (1983)
       // on the quoted pages
   
       // broadband transmittance by Rayleigh scattering (Iqbal (1983), p.189)
       const double taur = exp( -0.0903 * pow(ma,0.84) * (1. + ma - pow(ma,1.01)) );
   
       // broadband transmittance by ozone (Iqbal (1983), p.189)
       const double u3 = olt * mr; // ozone relative optical path length
       const double alpha_oz = 0.1611 * u3 * pow(1. + 139.48 * u3, -0.3035) -
                               0.002715 * u3 / ( 1. + 0.044  * u3 + 0.0003 * u3 * u3); //ozone absorbance
       const double tauoz = 1. - alpha_oz;
   
       // broadband transmittance by uniformly mixed gases (Iqbal (1983), p.189)
       const double taug = exp( -0.0127 * pow(ma, 0.26) );
   
       // saturation vapor pressure in Pa
       //const double Ps = exp( 26.23 - 5416./ta ); //as used for the parametrization
       const double Ps = mio::Atmosphere::vaporSaturationPressure(ta);
   
       // Leckner (1978) (in Iqbal (1983), p.94), reduced precipitable water
       const double w = 0.493 * rh * Ps / ta;
   
       // pressure corrected relative optical path length of precipitable water (Iqbal (1983), p.176)
       // pressure and temperature correction not necessary since it is included in its numerical constant
       const double u1 = w * mr;
   
       // broadband transmittance by water vapor (in Iqbal (1983), p.189)
       const double tauw = 1. - 2.4959 * u1  / (pow(1.0 + 79.034 * u1, 0.6828) + 6.385 * u1);
   
       // broadband total transmittance by aerosols (in Iqbal (1983), pp.189-190)
       // using Angstroem's turbidity formula Angstroem (1929, 1930) for the aerosol thickness
       // in Iqbal (1983), pp.117-119
       // aerosol optical depth at wavelengths 0.38 and 0.5 micrometer
       const double ka1 = beta * pow(0.38, -alpha);
       const double ka2 = beta * pow(0.5, -alpha);
   
       // broadband aerosol optical depth:
       const double ka  = 0.2758 * ka1 + 0.35 * ka2;
   
       // total aerosol transmittance function for the two wavelengths 0.38 and 0.5 micrometer:
       const double taua = exp( -pow(ka, 0.873) * (1. + ka - pow(ka, 0.7088)) * pow(ma, 0.9108) );
   
       // broadband transmittance by aerosols due to absorption only (Iqbal (1983) p. 190)
       const double tauaa = 1. - (1. - w0) * (1. - ma + pow(ma, 1.06)) * (1. - taua);
   
       // broadband transmittance function due to aerosols scattering only
       // Iqbal (1983) p. 146 (Bird and Hulstrom (1981))
       const double tauas = taua / tauaa;
   
       // direct normal solar irradiance in range 0.3 to 3.0 micrometer (Iqbal (1983) ,p.189)
       // 0.9751 is for this wavelength range.
       // Bintanja (1996) (see Corripio (2002)) introduced a correction beta_z for increased
       // transmittance with altitude that is linear up to 3000 m and than fairly constant up to 5000 - 6000 m
       const double beta_z = (altitude<3000.)? 2.2*1.e-5*altitude : 2.2*1.e-5*3000.;
   
       //Now calculating the radiation
       //Top of atmosphere radiation (it will always be positive, because we check for sun elevation before)
       const double tau_commons = tauoz * taug * tauw * taua;
   
       // Diffuse radiation from the sky
       const double factor = 0.79 * R_toa * tau_commons / (1. - ma + pow( ma,1.02 ));  //avoid recomputing pow() twice
       // Rayleigh-scattered diffuse radiation after the first pass through atmosphere (Iqbal (1983), p.190)
       const double Idr = factor * 0.5 * (1. - taur );
   
       // aerosol scattered diffuse radiation after the first pass through atmosphere (Iqbal (1983), p.190)
       const double Ida = factor * fc  * (1. - tauas);
   
       // cloudless sky albedo Bird and Hulstrom (1980, 1981) (in Iqbal (1983) p. 190)
       //in Iqbal, it is recomputed with ma=1.66*pressure/101325.; and alb_sky=0.0685+0.17*(1.-taua_p)*w0;
       const double alb_sky = 0.0685 + (1. - fc) * (1. - tauas);
   
   
       //Now, we compute the direct and diffuse radiation components
       //Direct radiation. All transmitances, including Rayleigh scattering (Iqbal (1983), p.189)
       R_direct = 0.9751*( taur * tau_commons + beta_z ) * R_toa ;
   
       // multiple reflected diffuse radiation between surface and sky (Iqbal (1983), p.154)
       const double Idm = (Idr + Ida + R_direct) * ground_albedo * alb_sky / (1. - ground_albedo * alb_sky);
       R_diffuse = (Idr + Ida + Idm)*cos_zenith; //Iqbal always "project" diffuse radiation on the horizontal
   
   
       double elevation_threshold = 2.0 * mio::Cst::to_rad;
   
       if( sun_elevation < elevation_threshold ) {
           //if the Sun is too low on the horizon, we put all the radiation as diffuse
           //the splitting calculation that might take place later on will reflect this
           //instead point radiation, it becomes the radiation of a horizontal sky above the domain
           R_diffuse += R_direct*cos_zenith; //HACK
          // R_diffuse =0.0;
           R_direct = 0.;
       }
   
   
       double cf = (*face)["cloud_frac"_s];
       double dir = R_direct  * (0.6 + 0.2*cos_zenith) * (1.0-cf);
   
       dir = std::max(0.0,dir);
       R_diffuse = std::max(0.0,R_diffuse);
   
       (*face)["iswr_direct_no_slope"_s]=dir;
       (*face)["iswr_diffuse_no_slope"_s]=R_diffuse;
   
       (*face)["atm_trans"_s]=(dir+R_diffuse/1375.);
   
   }