; $Id: AtmosTransGR.pro,v 6.0
;
;+
; NAME:
;       AtmosTransGR
;
; PURPOSE:
;       Computes direct & diffuse solar irradiance for clear skies  for the entire GRID
;
; CATEGORY:
;       Solar Radiation
;       Atmospheric Sciences
;
; CALLING SEQUENCE:
;       result = AtmosTransGR(zenith,day,DEM,visibility,RH,TempK,O3,alphag)
;
; INPUTS:
;       zenith: solar zenithal angle (degrees)
;               (relative to horizontal surface) (REF:1)
;       day: day of the year (1 : 365)
;       DEM: digital elevation model
;       visibility: in Km, for aerosol attenuation (5 < vis < 180 Km)
;               aerosol alternatives to be implemented
;       RH: 2D array of relative humidity (0:100)%
;       TempK: 2D array of  temperature, degrees Kelvin
;       O3: ozone layer thickness in cm (from REF:2, see notes)
;     alphaG: albedoground
;
; KEYWORDS:
;       MS: float variable.
;               set this value to ground albedo (0.0 to 1.0)
;               to compute multiple scattering between ground and Sky
;       TAU: named structure.  Will contain individual transmittances due to:
;               tauR: Rayleigh scattering
;               tauG: uniformly mixed gases other that water vapour
;               tauW: water vapour
;               tauO: ozone
;               tauA: aerosols
;
; OUTPUTS:
;       Result: 3-dimensional array  [2,m,n] where [0,m,n] is direct irradiance normal to surface
;       and [1,m,n] is diffuse irradiance
;
;
; SUBROUTINES:
;
; REQUIREMENTS:
;       FUNCTION wvapsat
;   FUNCTION  z2P
;
; NOTES:
;       1) 1 Dobson Unit (DU) is defined to be 0.01 mm thickness at stp.
;
; COMMENTS:
;
;
; REFERENCES:
;
;; REF:1
;; http://www.dict.org/bin/Dict?Form=Dict2&Database=devils&Query=ZENITH
;;
;; REF:2
;; TOMS/EP Total Ozone Mapping Spectrometer / Earth Probe
;; http://toms.gsfc.nasa.gov/ozone/ozone.html
;;
;; Spencer, J. W.: 1971, Fourier series representation of the position of
;; the sun, Search 2, 172.
;;
;; REFERENCES FOR COMPUTATION: link to URL (too many)Iqbal83 ParMod C
;
;
; EXAMPLE:
;       Insol = AtmosTransGR(45,180,DEM,50,RH,TempK,.321,.15)
;
;
; MODIFICATION HISTORY:
;       JGC, 24 NOvember, 2004
;
;
;  Copyright (c) 2002 Javier G Corripio.
;  corripio@members.ises.org corripio@ifu.baug.ethz.ch
;
;  This software may be used, copied, or redistributed as long as it is not
;  sold and this copyright notice is reproduced on each copy made.  This
;  routine is provided as is without any express or implied warranties
;  whatsoever.
;
;-

@wvapsat
@z2P

FUNCTION AtmosTransGR,zenith,day,DEM,visibility,RH,TempK,O3,alphag

;ON_ERROR, 2
IF (n_params() LT 8) THEN BEGIN
        message, 'USAGE: Result = AtmosTransGR(zenith,day,DEM,visibility,RH,'+$
                'TempK,O3,alphag)',/INFO
        RETURN,!values.F_NaN
ENDIF
        ;**night time
IF (zenith GT 90 ) THEN return, [0.0,0.0]

;----- PARAMETERS & Variables ------------------------------------------------
stlapse = -0.0065      ; standard lapse rate K/m
radeg = 180.0d/!dPI    ; radianes to degrees double precission
Theta = zenith/!RADEG  ; solar zenith angle radians
Isc = 1367.0   ; solar constant (Wm^(-2))
ssctalb = 0.9  ; single scattering albedo (aerosols)(Iqbal, 1983)
Fc = 0.84      ; ratio of forward to total energy scattered (Iqbal, 1983)
REarth = 6.3756766D6    ; Average earth's radius (m)
Md = 28.966d0     ; Molecular weight of dry air
atmos_P0 = 1.013250D5   ; Standard sea-level atmospheric pressure Pa
atmos_T0 = 0288.15d0    ; Standard sea-level Temperature
atmos_R = 8.31432D0     ; Universal gas constant (J deg-1 kmol-1)
earth_G = 9.80665D0     ; Acceleration due to gravity (m s-2)
stlapse = -0.0065d0 ; standard lapse rate K/m



; alpha_atmos   ; atmospheric albedo
; alphag        ; ground albedo
; In            ; direct normal irradiance (Wm^(-2))
; Id            ; total diffuse irradiance on a horizontal surface(Wm^(-2))
; Idr           ; diffuse irradiance on a h.s.-> molecular scattering (Wm^(-2))
; Ida           ; diffuse irradiance on a h.s.-> aerosols (Wm^(-2))
; Idm           ; diffuse irradiance on a h.s.-> multiplescattering (Wm^(-2))
; Mr            ; relative optical air mass
; Ma            ; relative optical air mass pressure corrected
; Pz            ; pressure at height z
; zenith        ; solar zenith angle degrees
; wvap_s        ; saturated vapour pressure (Leckner, 1978)
; Wprec         ; precipitable water (Leckner, 1978)
; TempK         ; temperature in Kelvin
;------------------------------------------------------------------------------



;Check for NaN values:
;nan = WHERE (finite([zenith,day,height,visibility,RH,TempK,O3,alphag],/NaN),count)
;IF count GT 0 THEN BEGIN
;       message,'data contain NaN values',/INFO
 ;       return, !values.F_NaN
;ENDIF



;------ ATMOSPHERIC PARAMETERS---------------------------------------------

;Pz = 1013.25*EXP(-0.0001184*height)
;*** pressure as a function of altitude US standard atmosphere p. 12 & 3
;; stlapse*1000 as in table 4, p. 3  temperature gradient is given in K/km

;   H1 = (REarth * height) /(REarth + height)  ;Geopotential height
;   HB = 0.0d0

;Pz = atmos_P0*( atmos_T0/(atmos_T0+stlapse*(H1-HB)) )^( (earth_G*Md)/(atmos_R*stlapse*1000) )
;Pz = Pz/100.0d0    ;Pa to mb

sz = size(DEM,/dim)
cols = sz[0]
rows = sz[1]

InId = fltarr(2,cols,rows)

Pz = z2P(DEM)

Mr = 1.0/(cos(theta)+0.15*((93.885-zenith)^(-1.253)))

Ma = Mr*Pz/1013.25

;** Use Lowe(1977) Lowe's polynomials for vapor pressure
wvap_s =  wVapSat(tempK)

;wvap_s = wvap_s *100.0     ;; mb to Pa    ;***************************
;Wprec = 0.493*(RH/100.0)*wvap_s/TempK   ;precipitable water in cm Leckner (1978)

Wprec = 46.5*(RH/100.0)*wvap_s/TempK  ;Prata 1996


;------- SUN EARTH DISTANCE: Spencer(1971) ---------------------------------
        ;eccentricity correction
 dayang_S = 2.0*!PI*(day-1)/365.0
 rho2 = 1.00011 + $
                0.034221*COS(dayang_S)+0.00128*SIN(dayang_S)+ $
                0.000719*COS(2*dayang_S)+0.000077*SIN(2*dayang_S)

;--------------- DIRECT IRRADIANCE ------------------------------------------

TauR = EXP((-.09030*(Ma^0.84) )*(1.0+Ma-(Ma^1.01)) )

TauO = 1.0-( ( 0.1611*(O3*Mr)*(1.0+139.48*(O3*Mr))^(-0.3035) ) - $
        0.002715*(O3*Mr)*( 1.0+0.044*(O3*Mr)+0.0003*(O3*Mr)^2 )^(-1))

TauG = EXP(-0.0127*(Ma^0.26))

TauW = 1.0-2.4959*(Wprec*Mr)*( (1.0+79.034*(Wprec*Mr))^0.6828 + $
        6.385*(Wprec*Mr) )^(-1)
        ;** Wprec is pressure & temperature corrected,
        ;** Mr need not to be correcred (Iqbal83, p. 176)

TauA = ( 0.97-1.265*(visibility^(-0.66)) )^(Ma^0.9)   ;Mächler, 1983
        ;** uses parameterization model A (visibility) (Iqbal83)

;;;**************************************************************************

;;; parameters are corrected for sea leval.  Apply correction for altitude:
;;; Bintanja (1996)
B_z = 2.2*10.0^(-5) ; m^-1
    ;; Bz increases linearly linear up to 3000m, then constant up to 5-6000m

Bintcorr = DEM * B_z
Bintcorr = Bintcorr < B_z*3000.0

;;;**************************************************************************

TauTotal = TauR*TauO*TauG*TauW*TauA + B_z   ;*******************

InId[0,*,*] = 0.9751*rho2*Isc*TauTotal



;--------------- DIFUSSE IRRADIANCE ------------------------------------------

tauaa = 1.0-(1.0-ssctalb)*(1.0-Ma+Ma^1.06)*(1.0-TauA)
Idr = 0.79*rho2*Isc*COS(theta)*TauO*TauG*TauW*tauaa*0.5*(1.0-TauR)/(1.0-Ma+Ma^(1.02))

tauas = (TauA)/tauaa
Ida = 0.79*rho2*Isc*COS(theta)*TauO*TauG*TauW*tauaa*Fc*(1.0-tauas)/(1.0-Ma+Ma^1.02)

;this one do without MS and incorporate fraction snow ground...
alpha_atmos = 0.0685+(1.0-Fc)*(1.0-tauas)
Idm = (InId[0,*,*]*COS(theta)+Idr+Ida)*alphag*alpha_atmos/(1.0-alphag*alpha_atmos)

InId[1,*,*] = Idr+Ida+Idm

Return, [InId]


END