123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225 |
- ! First include the set of model-wide compiler flags
- #include "tm5.inc"
- module Mrweight
- !---------------------------------------------------------------------------
- ! Contains routine to compute the intensity at the top of the atmosphere
- ! for a Rayleigh scattering atmosphere with a Lambertian surface.
- !
- ! call rweight(wavelx,p_sfc,a_sfc,th0,th,dphix,refl)
- !
- ! Pepijn Veefkind, Henk Eskes, KNMI, 1999
- !---------------------------------------------------------------------------
- implicit none
- private
- public :: rweight
-
- contains
- subroutine csalbr(xtau,xalb)
- implicit none
- real,intent(in) :: xtau
- real,intent(out) :: xalb
- xalb=0
- xalb=(3*xtau-fintexp3(xtau)*(4+2*xtau)+2*exp(-xtau))
- xalb=xalb/(4.+3*xtau)
- end subroutine csalbr
- function fintexp3(xtau)
- implicit none
- real,intent(in) :: xtau
- real :: fintexp3
- real :: xx
- xx=(exp(-xtau)*(1.-xtau)+xtau*xtau*fint1exp(xtau))/2.
- fintexp3=xx
- end function fintexp3
- function fint1exp(xtau)
- ! accuracy 2e-07... for 0<xtau<1
- implicit none
- real,intent(in) :: xtau
- real :: fint1exp
- real :: xx,xftau
- real,dimension(0:5),parameter :: &
- a=(/-.57721566, 0.99999193,-0.24991055, &
- 0.05519968,-0.00976004, 0.00107857 /)
- integer :: i
- xx=a(0)
- xftau=1.
- do i=1,5
- xftau=xftau*xtau
- xx=xx+a(i)*xftau
- enddo
- fint1exp=xx-log(xtau)
- end function fint1exp
- subroutine chand (xphi,xmuv,xmus,xtau,xrray)
- !
- ! input parameters: xphi,xmus,xmuv,xtau
- ! xphi: azimuthal difference between sun and observation (xphi=0,
- ! in backscattering) and expressed in degree (0.:360.)
- ! xmus: cosine of the sun zenith angle
- ! xmuv: cosine of the observation zenith angle
- ! xtau: molecular optical depth
- ! output parameter:
- ! xrray : molecular reflectance (0.:1.)
- ! constant : xdep: depolarization factor (0.0279)
- !
- real,intent(in) :: xphi,xmuv,xmus,xtau
- real,intent(out) :: xrray
- !
- real :: xdep,pl(10)
- real :: fs0,fs1,fs2
- real :: as0(10),as1(2),as2(2)
- real :: fac,pi,phios,xcosf1,xcosf2
- real :: xcosf3,xbeta2,xfd,xph1,xph2,xph3,xitm, xp1, xp2, xp3
- real :: cfonc1,cfonc2,cfonc3,xlntau,xitot1,xitot2,xitot3
- integer :: i
- !
- as0 = (/.33243832,-6.777104e-02,.16285370 &
- ,1.577425e-03,-.30924818,-1.240906e-02,-.10324388 &
- ,3.241678e-02,.11493334,-3.503695e-02/)
- as1 = (/.19666292, -5.439061e-02/)
- as2 = (/.14545937, -2.910845e-02/)
- !
- pi=acos(-1.)
- fac=pi/180.
- phios=180.-xphi
- xcosf1=1.
- xcosf2=cos(phios*fac)
- xcosf3=cos(2*phios*fac)
- xbeta2=0.5
- xdep=0.0279
- xfd=xdep/(2-xdep)
- xfd=(1-xfd)/(1+2*xfd)
- xph1=1+(3*xmus*xmus-1)*(3*xmuv*xmuv-1)*xfd/8.
- xph2=-xmus*xmuv*sqrt(1-xmus*xmus)*sqrt(1-xmuv*xmuv)
- xph2=xph2*xfd*xbeta2*1.5
- xph3=(1-xmus*xmus)*(1-xmuv*xmuv)
- xph3=xph3*xfd*xbeta2*0.375
- xitm=(1-exp(-xtau*(1/xmus+1/xmuv)))*xmus/(4*(xmus+xmuv))
- xp1=xph1*xitm
- xp2=xph2*xitm
- xp3=xph3*xitm
- xitm=(1-exp(-xtau/xmus))*(1-exp(-xtau/xmuv))
- cfonc1=xph1*xitm
- cfonc2=xph2*xitm
- cfonc3=xph3*xitm
- xlntau=log(xtau)
- pl(1)=1.
- pl(2)=xlntau
- pl(3)=xmus+xmuv
- pl(4)=xlntau*pl(3)
- pl(5)=xmus*xmuv
- pl(6)=xlntau*pl(5)
- pl(7)=xmus*xmus+xmuv*xmuv
- pl(8)=xlntau*pl(7)
- pl(9)=xmus*xmus*xmuv*xmuv
- pl(10)=xlntau*pl(9)
- fs0=0.
- do i=1,10
- fs0=fs0+pl(i)*as0(i)
- enddo
- fs1=pl(1)*as1(1)+pl(2)*as1(2)
- fs2=pl(1)*as2(1)+pl(2)*as2(2)
- xitot1=xp1+cfonc1*fs0*xmus
- xitot2=xp2+cfonc2*fs1*xmus
- xitot3=xp3+cfonc3*fs2*xmus
- xrray=xitot1*xcosf1
- xrray=xrray+xitot2*xcosf2*2
- xrray=xrray+xitot3*xcosf3*2
- xrray=xrray/xmus
- !
- end subroutine chand
- subroutine rweight(wavelx,p_sfc,a_sfc,th0,th,dphix,refl)
- !-----------------------------------------------------------------------
- ! name :weight
- ! description:computes radiance weight as a function of cloud fraction
- ! date :18 August 1999
- ! version :0.0
- ! authors : J.P. Veefkind
- !-----------------------------------------------------------------------
- ! Variable declaration
- ! Input:
- ! Wavelength [nm] : wavel0
- ! Surface pressure [hPa] : p_sfc
- ! Surface albedo : a_sfc
- ! Solar zenith angle : th0
- ! Viewing zenith angle : th
- ! Relative azimuth angle : dphix
- ! Output:
- ! total reflectance : refl
- !-----------------------------------------------------------------------
- implicit none
- !
- real,intent(in) :: wavelx, p_sfc, a_sfc, th0, th, dphix
- real,intent(out) :: refl
- !
- real :: refl_a, t_mu, t_mu0, a_sph, a
- real :: mu, mu0, pi, dphi, wavel
- real :: tau_sfc, tau
- !
- pi=acos(-1.)
-
- !-----------------------------------------------------------------------
- ! Main
- !-----------------------------------------------------------------------
- ! Conversion of some parameters
- ! wavelength in micrometers:
- wavel=wavelx*1e-3
- ! compute mu and mu0
- mu0=cos(th0/180.*pi)
- mu =cos(th /180.*pi)
- !
- ! use symmetry of dphi
- dphi=abs(dphix)
- if( dphi > 180.0 ) dphi=360.0-dphi
- !
- ! redefine dphi
- dphi=180.-dphi
- !
- ! Compute the rayleigh optical depth for surface and cloud top pressure
- ! approximation formula by Hansen and Travis
- ! tau=0.008569 * wavel**(-4.) * (1. + 0.0113 * wavel**(-2.) + 0.00013 * wavel**(-4.))
-
- tau = (1. + 0.0113/(wavel*wavel) + 0.00013/(wavel*wavel*wavel*wavel))*0.008569/(wavel*wavel*wavel*wavel)
- tau_sfc= p_sfc / 1013.25 * tau
- tau=tau_sfc
- a=a_sfc
-
- ! compute the transmission, see Vermote and Tanre JQSRT 47, pp 305, 1992
- t_mu0=( (2./3.+mu0)+(2./3.-mu0)*exp(-tau/mu0) )/( (4./3.)+tau )
- t_mu =( (2./3.+mu )+(2./3.-mu )*exp(-tau/mu ) )/( (4./3.)+tau )
-
- ! compute sherical albedo
- call csalbr(tau, a_sph)
-
- ! compute atm reflectance
- call chand (dphi,mu,mu0,tau,refl_a)
- ! compute total reflectance
- refl=refl_a + t_mu0 * t_mu* a / (1. - a_sph * a)
-
- ! Write result of RT to screen
- ! print*,'dphix, dphi ',dphix, dphi
- ! print*,'mu0, mu ',mu0, mu
- ! write(*,'(a,f9.3)') 'Transmission mu0 ',t_mu0
- ! write(*,'(a,f9.3)') 'Transmission mu0 ',t_mu
- ! write(*,'(a,f9.3)') 'Spherical albedo ',a_sph
- ! write(*,'(a,f9.3)') 'Atmospheric reflectance ',refl_a
- ! write(*,'(a,f9.3)') 'Reflectance ',refl
- end subroutine rweight
- end module Mrweight
|