program ssssss c**********************************************************************c c c c c c c c c c c c ******************************************************** c c * Second Simulation of a Satellite Signal * c c * in the Solar Spectrum * c c * ... (6SV) ....... (6SV) ...... (6SV) ... * c c * Version 2.0 * c c * * c c * Vector Code * c c * * c c * This code predicts the satellite signal from 0.25 * c c * to 4.0 microns assuming cloudless atmosphere. The * c c * main atmospheric effects (gaseous absorption by * c c * by water vapor, carbon dioxyde, oxygen and ozone; * c c * scattering by molecules and aerosols) are taken * c c * into account. Non-uniform surfaces, as well as * c c * bidirectional reflectances may be considered as * c c * boundary conditions. * cf c * * c c * The following input parameters are needed: * c c * geometrical conditions * c c * atmospheric model for gaseous components * c c * aerosol model (type and concentration) * c c * spectral condition * c c * ground reflectance (type + spectral variation) * c c * At each step, you can either select some proposed * c c * standard conditions (for example, spectral bands of * c c * satellite for spectral conditions) or define your * c c * own conditions (in the example, you have to define * c c * the assumed spectral response). * c c * * c c * More details are given at each data input step. * c c * * c c ******************************************************** c c c c c c c c c c c c**********************************************************************c c**********************************************************************c c c c c c ******************************************************** c c * The authors of this code are * c c * * c c * (1) E.F. Vermote and S.Y. Kotchenova; * c c * (2) J.C. Roger; * c c * (3) D. Tanre, J.L. Deuze, M. Herman; * c c * (4) J.J. Morcrette; * c c * (5) T. Miura. * c c * * c c * affiliated with * c c * * c c * (1) Department of Geography, University of * c c * Maryland (4321 Hartwick Road, College Park, * c c * MD 20740, USA) and NASA Goddard Space * c c * Flight Center (code 614.5, Greenbelt, MD * c c * 20771, USA) * c c * * c c * (2) Observatoire de Physique du Globe de * c c * Glermont-Ferrand Universite Blaise Pascal * c c * (24 Avenue des Landais, 63177 Aubiere, * c c * France) * c c * * c c * (3) Laboratoire d'Optique Atmospherique, * c c * Universite des Sciences et Techniques de * c c * Lille (u.e.r. de Physique Fondamentale, * c c * 59655 Villeneuve d'Ascq Cedex, France) * c c * * c c * (4) European Center for Medium-Range Weather * c c * Forecasts (Shinfield Park, Reading, RG2 * c c * 9AX, United Kingdom) * c c * * c c * (5) University of Hawaii at Manoa * c c * (1910 East_West Road, Sherman Lab 101 * c c * Honolulu, HI 96822) * c c * * c c * * c c * * c c * * c c ******************************************************** c c c c c c**********************************************************************c c**********************************************************************c c ******************************************************** c c * limits of validity * c c * * c c * geometrical parameters no limitations * c c * * c c * atmospheric model no limitations * c c * * c c * aerosol model the visibility must be * c c * better than 5.0km * c c * for smaller values * c c * calculations might be * c c * no more valid. * c c * * c c * spectral conditions the gaseous transmittance* c c * and the scattering func * c c * tions are valid from 0.25* c c * to 4.0 micron. but the * c c * treatment of interaction * c c * between absorption and * c c * scattering is correct for* c c * not too large absorption * c c * if you want to compute * c c * signal within absorption * c c * bands,this interaction * c c * ought to be reconsidered * c c * * c c * ground reflectance (type) you can consider a patchy* c c * structure:that is a circu* c c * lar target of radius rad * c c * and of reflectance roc, * c c * within an environnement * c c * of reflectance roe. * c c * * c c * ground reflectance (type continued): for uniform * c c * surface conditions only, * c c * you may consider directio* c c * nal reflectance as bounda* c c * ry conditions. * c c * some analytical model are* c c * proposed, the user can * c c * specify his own values. * c c * the code assumes that the* c c * brdf is spectrally inde- * c c * pendent * c c * * c c * ground reflectance (spectral variation) four typi * c c * cal reflectances are pro * c c * posed, defined within * c c * given spectral range. * c c * this range differs accor * c c * ding to the selected case* c c * the reflectance is set to* c c * 0 outside this range,due * c c * to the deficiency of data* c c * user must verify these * c c * limits. that is obviously* c c * irrelevant for brdf * c c * * c c ******************************************************** c c**********************************************************************c c****************************************************************************c c for considering brdf< we have to compute the downward radiance in the c c whole hemisphere. to perform such computions, we selected the successive c c orders of scattering method. that method requires numerical integration c c over angles and optical depth. the integration method is the gauss method,c c mu is the number of angles nmu+1, nmu is settled to 24. the accuracy of c c the computations is obviously depending on the nmu value. this value c c can be easily changed as a parameter as well as the nt value which c c is the number of layers for performing the vertical integration. the c c downward radiance is computed for nmu values of the zenith angle and np c c values of the azimuth angle. the integration of the product of the c c radiance by the brdf is so performed over the nmu*np values. np is settledc c to 13, that value can be also changed. mu2 is equal to 2 times nmu. c c xlmus is the downward radiance, xf the downward irradiance, rm and gb c c the angles and the weights for the gauss integration over the zenith, rp c c and gp respectively for the azimuth integration. c c****************************************************************************c include "paramdef.inc" dimension anglem(mu2_p),weightm(mu2_p), s rm(-mu_p:mu_p),gb(-mu_p:mu_p),rp(np_p),gp(np_p) dimension xlmus(-mu_p:mu_p,np_p),xlmuv(-mu_p:mu_p,np_p) dimension angmu(10),angphi(13),brdfints(-mu_p:mu_p,np_p) s ,brdfdats(10,13),sbrdftmp(-1:1,1),sbrdf(1501), s srm(-1:1),srp(1), s brdfintv(-mu_p:mu_p,np_p),brdfdatv(10,13),robar(1501), s robarp(1501),robard(1501),xlm1(-mu_p:mu_p,np_p), s xlm2(-mu_p:mu_p,np_p),asolvec(10),phivvec(21),avisvec(15), s windvec(4) real romix_fi(nfi_p),rorayl_fi(nfi_p),ratm2_fi(nfi_p), s refet_fi(nfi_p),roatm_fi(3,20,nfi_p),xlphim(nfi_p) c*********************************************************************** c for including BRDF as ground boundary condition c in OSSUR (Vermote, 11/19/2010) c*********************************************************************** real rosur(0:mu_p,mu_p,83) real wfisur(83),fisur(83) real xlsurf(-mu_p:mu_p,np_p),rolutsurf(mu_p,61) real lddiftt,lddirtt,lsphalbt,ludiftt,ludirtt real lxtrans(-1:1) real rbar,rbarp,rbarc,rbarpc,rbard c*********************************************************************** real asolvec,phivvec,avisvec,windvec c************************************************************************* real rolut(mu_p,61),roluts(20,mu_p,61),roluti(mu_p,61) real rolutq(mu_p,61),rolutsq(20,mu_p,61),rolutiq(mu_p,61) real rolutu(mu_p,61),rolutsu(20,mu_p,61),rolutiu(mu_p,61) real filut(mu_p,61) integer aerod integer solind,visind,phivind,windind real its,lutmuv,luttv,iscama,iscami,scaa,cscaa,cfi integer nfilut(mu_p),nbisca real dtr real anglem,weightm,rm,gb,accu2,accu3 real rp,gp,xlmus,xlmuv,angmu,angphi,brdfints,brdfdats real brdfintv,brdfdatv,robar,robarp,robard,xlm1,xlm2 real c,wldisc,ani,anr,aini,ainr,rocl,roel,zpl,ppl,tpl,whpl real wopl,xacc,s,wlinf,wlsup,delta real nwlinf,nwlsup integer niinf,nisup real sigma,z,p,t,wh,wo,ext,ome,gasym,phase,qhase,roatm,dtdir real dtdif,utdir,utdif,sphal,wldis,trayl,traypl,pi,pi2,step real asol,phi0,avis,phiv,tu,xlon,xlat,xlonan,hna,dsol,campm real phi,phirad,xmus,xmuv,xmup,xmud,adif,uw,uo3,taer55 real taer,v,xps,uwus,uo3us,xpp,taer55p,puw,puo3,puwus real v_con,taer55_con,xps_con,xpp_con real pws_con,phi_wind_con,xsal_con,pcl_con real puo3us,wl,wlmoy,tamoy,tamoyp,pizmoy,pizmoyp,trmoy real trmoyp,fr,rad,spalt,sha,sham,uhase real albbrdf,par1,par2,par3,par4,robar1,xnorm1,rob,xnor,rodir real rdown,rdir,robar2,xnorm2,ro,roc,roe,rapp,rocave,roeave real seb,sbor,swl,sb,refet,refet1,refet2,refet3,alumet real refeti,pinst,ksiinst,ksirad real rpfet,rpfet1,rpfet2,rpfet3,plumet,plumeas real tgasm,rog,dgasm,ugasm,sdwava,sdozon,sddica,sdoxyg real sdniox,sdmoca,sdmeth,suwava,suozon,sudica,suoxyg real suniox,sumoca,sumeth,stwava,stozon,stdica,stoxyg,stniox real stmoca,stmeth,sodray,sodaer,sodtot,fophsr,fophsa,sroray real sroaer,srotot,ssdaer,sdtotr,sdtota,sdtott,sutotr,sutota real sutott,sasr,sasa,sast,dtozon,dtdica,dtoxyg real dtniox,dtmeth,dtmoca,utozon,utdica,utoxyg,utniox real utmeth,utmoca,attwava,ttozon,ttdica,ttoxyg,ttniox real ttmeth,ttmoca,dtwava,utwava,ttwava,coef,romix,rorayl real roaero,phaa,phar,tsca,tray,trayp,taerp,dtott,utott real rqmix,rqrayl,rqaero,qhaa,qhar,foqhsr,foqhsa,foqhst real rumix,rurayl,ruaero,uhaa,uhar,rpmix,rpaero,rprayl real srpray,srpaer,srptot,rpmeas1,rpmeas2,rpmeas3 real srqray,srqaer,srqtot,sruray,sruaer,srutot real astot,asray,asaer,utotr,utota,dtotr,dtota,dgtot,tgtot real tgp1,tgp2,rqatm,ruatm,fouhst,fouhsr,fouhsa,coefp real ugtot,edifr,edifa,tdird,tdiru,tdifd,tdifu,fra real fae,avr,romeas1,romeas2,romeas3,alumeas,sodrayp real sdppray,sdppaer,sdpptot,rop,sdpray,sdpaer,sdptot real spdpray,spdpaer,spdptot real ratm1,ratm2,ratm3,rsurf,rpatm1,rpatm2,rpatm3,rpsurf real sodaerp,sodtotp,tdir,tdif,etn,esn,es,ea0n,ea0,ee0n real ee0,tmdir,tmdif,xla0n,xla0,xltn,xlt,xlen,xle,pizera real fophst,pizerr,pizert,xrad,xa,xb,xc integer nt,mu,mu2,np,k,iwr,mum1,idatmp,ipol integer j,iread,l,igeom,month,jday,nc,nl,idatm,iaer,iaerp,n integer iaer_con,iaerp_con integer iwave,iinf,isup,ik,i,inhomo,idirec,ibrdf,igroun integer iwave_con,inhomo_con,idirec_con,ibrdf_con,irapp_con,igroun_con integer igrou1,igrou2,isort,irapp,ilut c variables used in the BRDF coupling correction process real robarstar,robarpstar,robarbarstar,tdd,tdu,tsd,tsu real coefa,coefb,coefc,discri,rogbrdf,roglamb,rbardest real romixatm,romixsur c variables related to surface polarization integer irop real ropq,ropu,pveg,wspd,azw,razw c*********************************************************************** c to vary the number of quadratures c*********************************************************************** integer nquad common /num_quad/ nquad c*********************************************************************** c the aerosol profile c*********************************************************************** integer iaer_prof,num_z real alt_z,taer_z,taer55_z,total_height,height_z(0:nt_p_max) common/aeroprof/num_z,alt_z(0:nt_p_max),taer_z(0:nt_p_max), &taer55_z(0:nt_p_max) character aer_model(15)*50 c*********************************************************************** c return to 6s c*********************************************************************** dimension c(4),wldisc(20),ani(2,3),anr(2,3),aini(2,3),ainr(2,3) dimension rocl(1501),roel(1501) real rfoaml(1501),rglitl(1501),rwatl(1501) real rn,ri,x1,x2,x3,cij,rsunph,nrsunph,rmax,rmin,cij_out(4) real rn_con(20,4),ri_con(20,4),x1_con(4),x2_con(4),cij_con(4),rmax_con,rmin_con, s cij_out_con(4) integer icp,irsunph,i1,i2 integer icp_con character etiq1(8)*60,nsat(206)*17,atmid(7)*51,reflec(8)*100 character FILE*80,FILE2*80 logical ier integer igmax common/sixs_ier/iwr,ier common /mie_in/ rmax,rmin,icp,rn(20,4),ri(20,4),x1(4),x2(4), s x3(4),cij(4),irsunph,rsunph(50),nrsunph(50) common /multorder/ igmax c*********************************************************************** c for considering pixel and sensor altitude c*********************************************************************** real pps,palt,ftray common /sixs_planesim/zpl(34),ppl(34),tpl(34),whpl(34),wopl(34) common /sixs_test/xacc c*********************************************************************** c for considering aerosol and ???? c*********************************************************************** integer options(5) integer pild,pihs real optics(3),struct(4) real pxLt,pc,pRl,pTl,pRs real pws,phi_wind,xsal,pcl,paw,rfoam,rwat,rglit real rfoamave,rwatave,rglitave real uli,eei,thmi,sli,cabi,cwi,vaii,rnci,rsl1i real p1,p2,p3,p1p,p2p,p3p c*********************************************************************** c return to 6s c*********************************************************************** common /sixs_ffu/s(1501),wlinf,wlsup common /sixs_del/ delta,sigma common /sixs_atm/z(34),p(34),t(34),wh(34),wo(34) common /sixs_aer/ext(20),ome(20),gasym(20),phase(20),qhase(20), suhase(20) common /sixs_disc/ roatm(3,20),dtdir(3,20),dtdif(3,20), s utdir(3,20),utdif(3,20),sphal(3,20),wldis(20),trayl(20), s traypl(20),rqatm(3,20),ruatm(3,20) c****************************************************************************c c angmu and angphi are the angles were the  is measured. these values c c can be changed as soon as they are well distributed over the whole space c c before the gauss integration, these values are interpolated to the gauss c c angles c c****************************************************************************c data angmu /85.0,80.0,70.0,60.0,50.0,40.0,30.0,20.0,10.0,0.00/ data angphi/0.00,30.0,60.0,90.0,120.0,150.0,180.0, s 210.0,240.0,270.0,300.0,330.0,360.0/ C****************************************************************************mino C 3/16/2020 C***** Looping change sol_zenith, Delta_phil , view zenith, wind speed data asolvec / 0.0, 12.0, 24.0, 36.0, 48.0, 54.0, 60.0, 66.0, s 72.0, 78.0/ C data avisvec / 0.0, 6.0, 12.0, 18.0, 24.0, 30.0, 36.0, 42.0, s 48.0, 54.0, 60.0, 66.0, 72.0, 78.0, 84.0/ data phivvec / 0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 72.0, 84.0, s 96.0, 108.0, 120.0, 132.0, 144.0, 156.0, 168.0, 170.0, s 172.0, 174.0, 176.0, 178.0, 180.0/ data windvec / 2.0, 6.0, 10.0, 14.0/ C*******INTERNAL FLAG********* integer iprtspr data iprtspr /0/ c*********************************************************************** c return to 6s c*********************************************************************** data wldisc /0.350,0.400,0.412,0.443,0.470,0.488,0.515,0.550, s 0.590,0.633,0.670,0.694,0.760,0.860,1.240,1.536, s 1.650,1.950,2.250,3.750/ data etiq1/ s '(1h*,22x,34h user defined conditions ,t79,1h*)', s '(1h*,22x,24h meteosat observation ,t79,1h*) ', s '(1h*,22x,25h goes east observation ,t79,1h*) ', s '(1h*,22x,25h goes west observation ,t79,1h*) ', s '(1h*,22x,30h avhrr (AM noaa) observation ,t79,1h*) ', s '(1h*,22x,30h avhrr (PM noaa) observation ,t79,1h*) ', s '(1h*,22x,24h h.r.v. observation ,t79,1h*) ', s '(1h*,22x,24h t.m. observation ,t79,1h*) '/ data nsat/ s ' constant ',' user s ', s ' meteosat ',' goes east ',' goes west ', s ' avhrr 1 (noaa6) ',' avhrr 2 (noaa6) ', s ' avhrr 1 (noaa7) ',' avhrr 2 (noaa7) ', s ' avhrr 1 (noaa8) ',' avhrr 2 (noaa8) ', s ' avhrr 1 (noaa9) ',' avhrr 2 (noaa9) ', s ' avhrr 1 (noaa10)',' avhrr 2 (noaa10)', s ' avhrr 1 (noaa11)',' avhrr 2 (noaa11)', s ' hrv1 1 ',' hrv1 2 ',' hrv1 3 ', s ' hrv1 pan ', s ' hrv2 1 ',' hrv2 2 ',' hrv2 3 ', s ' hrv2 pan ', s ' tm 1 ',' tm 2 ',' tm 3 ', s ' tm 4 ',' tm 5 ',' tm 7 ', s ' mss 4 ',' mss 5 ', s ' mss 6 ',' mss 7 ', s ' mas 1 ',' mas 2 ',' mas 3 ', s ' mas 4 ',' mas 5 ',' mas 6 ', s ' mas 7 ',' modis 1 ',' modis 2 ', s ' modis 3 ',' modis 4 ',' modis 5 ', s ' modis 6 ',' modis 7 ', s ' avhrr 1 (noaa12)',' avhrr 2 (noaa12)', s ' avhrr 1 (noaa14)',' avhrr 2 (noaa14)', s ' polder 1 ',' polder 2 ', s ' polder 3 ',' polder 4 ',' polder 5 ', s ' polder 6 ',' polder 7 ',' polder 8 ', s ' seawifs 1 ',' seawifs 2 ', s ' seawifs 3 ',' seawifs 4 ',' seawifs 5 ', s ' seawifs 6 ',' seawifs 7 ',' seawifs 8 ', s ' aatsr 1 ',' aatsr 2 ',' aatsr 3 ', s ' aatsr 4 ',' meris 1 ',' meris 2 ', s ' meris 3 ',' meris 4 ',' meris 5 ', s ' meris 6 ',' meris 7 ',' meris 8 ', s ' meris 9 ',' meris 10 ',' meris 11 ', s ' meris 12 ',' meris 13 ',' meris 14 ', s ' meris 15 ',' gli 1 ',' gli 2 ', s ' gli 3 ',' gli 4 ',' gli 5 ', s ' gli 6 ',' gli 7 ',' gli 8 ', s ' gli 9 ',' gli 10 ',' gli 11 ', s ' gli 12 ',' gli 13 ',' gli 14 ', s ' gli 15 ',' gli 16 ',' gli 17 ', s ' gli 18 ',' gli 19 ',' gli 20 ', s ' gli 21 ',' gli 22 ',' gli 23 ', s ' gli 24 ',' gli 25 ',' gli 26 ', s ' gli 27 ',' gli 28 ',' gli 29 ', s ' gli 30 ',' ali 1p ',' ali 1 ', s ' ali 2 ',' ali 3 ',' ali 4 ', s ' ali 4p ',' ali 5p ',' ali 5 ', s ' ali 7 ',' aster 1 ',' aster 2 ', s ' aster 3n ',' aster 3b ',' aster 4 ', s ' aster 5 ',' aster 6 ',' aster 7 ', s ' aster 8 ',' aster 9 ',' etm 1 ', s ' etm 2 ',' etm 3 ',' etm 4 ', s ' etm 5 ',' etm 7 ',' hypblue 1 ', s ' hypblue 2 ',' vgt 1 ',' vgt 2 ', s ' vgt 3 ',' vgt 4 ',' viirs m1 ', s ' viirs m2 ',' viirs m3 ',' viirs m4 ', s ' viirs m5 ',' viirs m6 ',' viirs m7 ', s ' viirs m8 ',' viirs m9 ',' viirs m10 ', s ' viirs m11 ',' viirs m12 ',' viirs i1 ', s ' viirs i2 ',' viirs i3 ',' viirs i4 ', s ' ldcm 1 ',' ldcm 2 ',' ldcm 3 ', s ' ldcm 4 ',' ldcm 5 ',' ldcm 6 ', s ' ldcm 7 ',' ldcm 8 ',' ldcm 9 ', s ' modis 8 ',' modis 9 ',' modis 10 ', s ' modis 11 ',' modis 12 ',' modis 13 ', s ' modis 14 ',' modis 15 ',' modis 16 ', s ' modis 17 ',' modis 18 ',' modis 19 ', s ' cavis 1 ',' cavis 2 ',' cavis 3 ', s ' cavis 4 ',' cavis 5 ',' cavis 6 ', s ' cavis 7 ',' cavis 8 ',' cavis 9 ', s ' cavis 10 ',' cavis 11 ', s ' dmc 1 ',' dmc 2 ',' dmc 3 ', s ' abi 1 ',' abi 2 ',' abi 3 ', s ' abi 4 ',' abi 5 ',' abi 6 '/ data atmid / s 'no absorption computed ', s 'tropical (uh2o=4.12g/cm2,uo3=.247cm-atm)', s 'midlatitude summer (uh2o=2.93g/cm2,uo3=.319cm-atm)', s 'midlatitude winter (uh2o=.853g/cm2,uo3=.395cm-atm)', s 'subarctic summer (uh2o=2.10g/cm2,uo3=.480cm-atm)', s 'subarctic winter (uh2o=.419g/cm2,uo3=.480cm-atm)', s 'us standard 1962 (uh2o=1.42g/cm2,uo3=.344cm-atm)'/ data reflec / & '(1h*,12x,39h user defined spectral reflectance ,f6.3,t79,1h*) ', & '(1h*,12x,27h monochromatic reflectance ,f6.3,t79,1h*)', & '(1h*,12x,39h constant reflectance over the spectra ,f6.3,t79,1h*) ', & '(1h*,12x,39h spectral vegetation ground reflectance,f6.3,t79,1h*) ', & '(1h*,12x,39h spectral clear water reflectance ,f6.3,t79,1h*) ', & '(1h*,12x,39h spectral dry sand ground reflectance ,f6.3,t79,1h*) ', & '(1h*,12x,39h spectral lake water reflectance ,f6.3,t79,1h*) ', & '(1h*,12x,39h spectral volcanic debris reflectance ,f6.3,t79,1h*) '/ FILE=' ' FILE2=' ' c*********************************************************************** c Parameters initialization c*********************************************************************** nt=nt_p mu=mu_p mu2=mu2_p np=np_p nfi=nfi_p iwr=6 ier=.FALSE. iinf=1 isup=1501 igmax=20 c*********************************************************************** c preliminary computations for gauss integration c*********************************************************************** pi=acos(-1.) pi2=2*pi accu2=1.E-03 accu3=1.E-07 do k=1,13 angphi(k)=angphi(k)*pi/180. enddo do k=1,10 angmu(k)=cos(angmu(k)*pi/180.) enddo call gauss(-1.,1.,anglem,weightm,mu2) call gauss(0.,pi2,rp,gp,np) mum1=mu-1 do 581 j=-mum1,-1 k=mu+j rm(-j-mu)=anglem(k) gb(-j-mu)=weightm(k) 581 continue do 582 j=1,mum1 k=mum1+j rm(mu-j)=anglem(k) gb(mu-j)=weightm(k) 582 continue gb(-mu)=0. gb(0)=0. gb(mu)=0. c*********************************************************************** c return to 6s c*********************************************************************** c constantes values sigma=0.056032 delta=0.0279 CCC pinst=0.02 CCC ksiinst=0. xacc=1.e-06 iread=5 step=0.0025 do 1111 l=1,20 wldis(l)=wldisc(l) 1111 continue c**********************************************************************c c c c * sun c c \ * / c c * * * * * c c z / * \ c c + /* c c satellite / + / c c o/ + / c c /.\ + /. c c / . \ _avis-_+_-asol_/ . c c . \- -+ / . north c c . \ + / . + c c . \ + / .+ c c . \ + / +. c c . \ + / + . c c . \ + / + . c c . \ +/ + . c c west + + + + + + + . + + + + +\+ + + + + . + + + + + + + + east c c . +.. . c c . + . . . c c . + . . . c c . + . .'. . c c . + .. . , ' .. c c .+ . \ . c c +. . \ . c c + . . \ . c c south . . (phiv-phi0) c c c c c c c c**********************************************************************c c**********************************************************************c c igeom geometrical conditions c c -------------------------------------- c c c c c c you choose your own conditions; igeom=0 c c 0 enter solar zenith angle (in degrees ) c c solar azimuth angle " c c satellite zenith angle " c c satellite azimuth angle " c c month c c day of the month c c c c or you select one of the following satellite conditions:igeom=1to7 c c 1 meteosat observation c c enter month,day,decimal hour (universal time-hh.ddd) c c n. of column,n. of line.(full scale 5000*2500) c c c c 2 goes east observation c c enter month,day,decimal hour (universal time-hh.ddd) c c n. of column,n. of line.(full scale 17000*12000)c c c c 3 goes west observation c c enter month,day,decimal hour (universal time-hh.ddd) c c n. of column,n. of line.(full scale 17000*12000)c c c c 4 avhrr ( PM noaa ) c c enter month,day,decimal hour (universal time-hh.ddd) c c n. of column(1-2048),xlonan,hna c c give long.(xlonan) and overpass hour (hna) at c c the ascendant node at equator c c c c 5 avhrr ( AM noaa ) c c enter month,day,decimal hour (universal time-hh.ddd) c c n. of column(1-2048),xlonan,hna c c give long.(xlonan) and overpass hour (hna) at c c the ascendant node at equator c c c c 6 hrv ( spot ) * enter month,day,hh.ddd,long.,lat.c c c c 7 tm ( landsat ) * enter month,day,hh.ddd,long.,lat.c c c c c c note: for hrv and tm experiments long. and lat. are the c c coordinates of the scene center. c c lat. must be > 0 for north lat., < 0 for south lat. c c long. must be > 0 for east long., <0 for west long. c c c c solar and viewing positions are computed c c c c**********************************************************************c C wind speed index C IAER=8 read(iread,*) iaer iaer_con=iaer if(iaer.eq.8)then read(iread,*) rmin,rmax,icp rmin_con=rmin rmax_con=rmax icp_con=icp do i=1,icp read(5,*)x1(i),x2(i),cij(i) x1_con=x1 x2_con=x2 cij_con=cij read(5,*)(rn(l,i),l=1,20) read(5,*)(ri(l,i),l=1,20) rn_con=rn ri_con=ri enddo do i=1,icp cij_out(i)=cij(i) enddo cij_out_con=cij_out read(5,*)iaerp iaerp_con=iaerp if (iaerp.eq.1)read(5,'(A80)')FILE i1=index(FILE,' ')-1 FILE2=FILE(1:I1)//'.mie' i2=index(FILE2,' ')-1 endif C MinOo this "end" if iaer.eq.12 if(iaer.eq.12)then read(5,'(A80)')FILE2 i2=index(FILE2,' ')-1 call aeroso(iaer,c,xmud,wldis,FILE2,ipol) endif read(iread,*) v v_con=v CC if(v) 71,10,11 CC if "v" is -1 goto 71, if v is 0 goto 10, if v is 1 goto 11 ( mino) CC 10 read(iread,*) taer55 read(iread,*) taer55 taer55_con=taer55 read(iread,*) xps xps_con=xps read(iread,*) xpp xpp_con=xpp read(iread,*) iwave iwave_con=iwave read(iread,*) inhomo CX if(inhomo) 30,30,31 ( is zero) inhomo_con=inhomo CX 30 read(iread,*) idirec ( is zero) read(iread,*) idirec idirec_con=idirec CX if(idirec)21,21,25 c**********************************************************************c c uniform conditions with brdf conditions c the following is for ocean , so comment them c c**********************************************************************c c cCX 25 read(iread,*) ibrdf cCO read(iread,*) ibrdf cCO ibrdf_con=ibrdf cCO if(ibrdf.eq.6) then cCO read(iread,*) pws,phi_wind,xsal,pcl cCO endif cCO pws_con=pws cCO phi_wind_con=phi_wind cCO xsal_con=xsal cCO pcl_con=pcl read(iread,*) igroun igroun_con=igroun read(iread,*) irapp irapp_con=irapp CX******************************************************************* CX LOOPING HERE C view phi (azimuth) index cCO do 1211 windind=1,4 do 1212 phivind=1,21 c 1,21 C view zenith index do 1213 visind=1,15 c 1,15 C solar zenith index do 1214 solind=1,10 c 1,10 CX read(iread,*) igeom C*****************************************8 igeom=0 iaer=iaer_con rmin=rmin_con rmax=rmax_con icp=icp_con x1=x1_con x2=x2_con cij=cij_con CX write(6,*)"cij",cij rn=rn_con ri=ri_con cij_out=cij_out_con CX write(6,*)"cij_out",cij_out iaerp=iaerp_con v=v_con taer55=taer55_con xps=xps_con xpp=xpp_con iwave=iwave_con inhomo=inhomo_con idirec=idirec_con igroun=igroun_con cCO ibrdf=ibrdf_con cCO if(ibrdf.eq.6) then cCO pws_con=windvec(windind) cCO pws=pws_con cCO phi_wind=phi_wind_con cCO xsal=xsal_con cCO pcl=pcl_con c endif irapp=irapp_con C*******************************************8 c goto 1311 CMin asol,phi0,avis,phiv,month,jday if (igeom.lt.0) then if (igeom.lt.-10) then igmax=int(abs(igeom/10)) igeom=igeom+igmax*10 endif ilut=0 igeom=0 endif ilut=0 goto(1001,1002,1003,1004,1005,1006,1007),igeom c igeom=0..... CX read(iread,*) asol,phi0,avis,phiv,month,jday asol=asolvec(solind) phi0=180.0 avis=avisvec(visind) phiv=phivvec(phivind) write(6,*)"solind,visind,phivind",solind ,visind ,phivind cx write(6,*)"igeom,iaer,rmin,rmax",igeom,iaer,rmin,rmax cx write(6,*)"rn",rn cx write(6,*)"ri",ri month=1 jday=23 goto 22 c 1001 read(iread,*) month,jday,tu,nc,nl call posmto(month,jday,tu,nc,nl, 1 asol,phi0,avis,phiv,xlon,xlat) goto 22 1002 read(iread,*) month,jday,tu,nc,nl call posge(month,jday,tu,nc,nl, 1 asol,phi0,avis,phiv,xlon,xlat) goto 22 1003 read(iread,*) month,jday,tu,nc,nl call posgw(month,jday,tu,nc,nl, 1 asol,phi0,avis,phiv,xlon,xlat) goto 22 1004 read(iread,*) month,jday,tu,nc,xlonan,hna campm=1.0 call posnoa(month,jday,tu,nc,xlonan,hna,campm, 1 asol,phi0,avis,phiv,xlon,xlat) goto 22 1005 read(iread,*) month,jday,tu,nc,xlonan,hna campm=-1.0 call posnoa(month,jday,tu,nc,xlonan,hna,campm, 1 asol,phi0,avis,phiv,xlon,xlat) goto 22 1006 read(iread,*) month,jday,tu,xlon,xlat call posspo(month,jday,tu,xlon,xlat, a asol,phi0,avis,phiv) goto 22 1007 read(iread,*) month,jday,tu,xlon,xlat call poslan(month,jday,tu,xlon,xlat, s asol,phi0,avis,phiv) 22 continue if(ier) stop dsol=1. call varsol(jday,month,dsol) c**********************************************************************c c c c / scattered direction c c / c c / c c / adif c c incident + + + + + + + + + + + + + + + c c direction c c c c**********************************************************************c phi=abs(phiv-phi0) phirad=(phi0-phiv)*pi/180. if (phirad.lt.0.) phirad=phirad+2.*pi if (phirad.gt.(2.*pi)) phirad=phirad-2.*pi xmus=cos(asol*pi/180.) xmuv=cos(avis*pi/180.) xmup=cos(phirad) xmud=-xmus*xmuv-sqrt(1.-xmus*xmus)*sqrt(1.-xmuv*xmuv)*xmup c test vermote bug if (xmud.gt.1.) xmud=1. if (xmud.lt.-1.) xmud=-1. adif=acos(xmud)*180./pi c**********************************************************************c c idatm atmospheric model c c -------------------- c c c c c c you select one of the following standard atmosphere: idatm=0 to 6 c c 0 no gaseous absorption c c 1 tropical ) c c 2 midlatitude summer ) c c 3 midlatitude winter ) c c 4 subarctic summer ) from lowtran c c 5 subarctic winter ) c c 6 us standard 62 ) c c c c or you define your own atmospheric model idatm=7 or 8 c c 7 user profile (radiosonde data on 34 levels) c c enter altitude ( in km ) c c pressure ( in mb ) c c temperature ( in k ) c c h2o density (in g/m3) c c o3 density (in g/m3) c c c c for example, altitudes are from 0 to 25km step of 1km c c from 25 to 50km step of 5km c c and two values at 70km and 100km c c so you have 34*5 values to input. c c 8 enter water vapor and ozone contents c c uw (in g/cm2 ) c c uo3 (in cm-atm) c c profil is taken from us62 c c c c**********************************************************************c uw=0. uo3=0. Cx read(iread,*) idatm idatm=0 if(idatm.eq.0) go to 5 if(idatm.eq.8) read(iread,*) uw,uo3 if(idatm.ne.7) go to 6 do 7 k=1,34 read(iread,*) z(k),p(k),t(k),wh(k),wo(k) 7 continue go to 5 6 if(idatm.eq.1) call tropic if(idatm.eq.2) call midsum if(idatm.eq.3) call midwin if(idatm.eq.4) call subsum if(idatm.eq.5) call subwin if(idatm.eq.6) call us62 c we have to define an atmosphere to compute rayleigh optical depth 5 if(idatm.eq.0.or.idatm.eq.8) call us62 c**********************************************************************c c THIS OPTION IS NOT AVAILABLE THE CODE RUNS WITH IPOL=1 c c ipol computation of the atmospheric polarization c c ------------------------------------------- c c c c**********************************************************************c c read(iread,*) ipol ipol=1 c write(6,*) "WARNING IPOL IS EQUAL 0" c**********************************************************************c c c c iaer aerosol model(type) and profile c c -------------- c c iaer = -1 The user-defined profile. You have to input the c c number of layers first, then the height (km), c c optical thickness (at 550 nm), and type of aerosol c c (see below) for each layer, starting from the c c ground. The present version of the program works c c only with the same type of aerosol for each layer. c c c c Example for iaer = -1: c c 4 c c 2.0 0.200 1 c c 10.0 0.025 1 c c 8.0 0.003 1 c c 80.0 0.000 1 c c c c The maximum total height of all layers cannot exceed 300 km. c c c c If you do not input iaer = -1, the program will use the default c c exponential profile. In this case, you need to select one of c c the following standard aerosol models: c c c c iaer = 0 no aerosols c c 1 continental ) c c 2 maritime ) according to d'Almeida's models c c 3 urban ) (see the manual) c c 5 background desert ) c c 6 biomass burning ) from AERONET measurements c c 7 stratospheric ) according to Russel's model c c c c or you define your own model using basic components: iaer=4 c c 4 enter the volumetric percentage of each component c c c(1) = volumetric % of dust-like c c c(2) = volumetric % of water-soluble c c c(3) = volumetric % of oceanic c c c(4) = volumetric % of soot c c between 0 to 1 c c c c or you define your own model using a size distribution function: c c 8 Multimodal Log-Normal distribution (up to 4 modes) c c 9 Modified Gamma distribution c c 10 Junge Power-Law distribution c c c c or you define a model using sun-photometer measurements: c c 11 Sun Photometer distribution (50 values max) c c you have to enter: r and dV/d(logr) c c where r is the radius (in micron), V is the volume, c c and dV/d(logr) is in (cm3/cm2/micron) c c then you have to enter: nr and ni for each wavelength c c where nr and ni are respectively the real and the c c imaginary parts of the refractive index c c c c or you can use the results computed and saved previously c c 12 Reading of data previously saved into FILE c c you have to enter the identification name FILE in the c c next line of inputs. c c c c c c iaerp and FILE aerosol model(type)-Printing of results c c --------------------------------------- c c c c For iaer=8,9,10,and 11: c c results from the MIE subroutine may be saved into the file c c FILE.mie (Extinction and scattering coefficients, single c c scattering albedo, asymmetry parameter, phase function at c c predefined wavelengths) and then can be re-used with the c c option iaer=12 where FILE is an identification name you c c have to enter. c c c c So, if you select iaer=8,9,10, or 11, you will have to enter c c iaerp after the inputs requested by options 8,9,10, or 11: c c c c iaerp=0 results will not be saved c c iaerp=1 results will be saved into the file FILE.mie c c next line enter FILE c c c c c c Example for iaer and iaerp c c 8 Multimodal Log-Normal distribution selected c c 0.001 20 3 Rmin, Rmax, 3 components c c 0.471 2.512 0.17 Rmean, Sigma, % density - 1st component c c 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.528 c c 1.52 1.462 1.4 1.368 1.276 1.22 1.2 nr for 20 wavelengths c c 0.008 0.008 0.008 0.008 0.008 0.008 0.008 0.008 0.008 0.008 0.008 c c 0.008 0.008 0.008 0.008 0.008 0.008 0.008 0.0085 0.011 ni c c 0.0285 2.239 0.61 Rmean, Sigma, % density - 2nd component c c 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.53 1.528 c c 1.52 1.51 1.42 1.42 1.42 1.42 1.452 nr for 20 wavelengths c c 0.005 0.005 0.005 0.005 0.005 0.005 0.0053 0.006 0.006 0.0067 0.007 c c 0.007 0.0088 0.0109 0.0189 0.0218 0.0195 0.0675 0.046 0.004 ni c c 0.0118 2.0 0.22 Rmean, Sigma, % density - 3rd component c c 1.75 1.75 1.75 1.75 1.75 1.75 1.75 1.75 1.75 1.75 1.75 1.75 1.75 c c 1.75 1.77 1.791 1.796 1.808 1.815 1.9 nr for 20 wavelengths c c 0.465 0.46 0.4588 0.4557 0.453 0.4512 0.447 0.44 0.436 0.435 0.433 c c 0.4306 0.43 0.433 0.4496 0.4629 0.472 0.488 0.5 0.57 ni c c 1 Results will be saved into FILE.mie c c Urban_Indust Identification of the output file (FILE) c c -> results will be saved into Urban_Indust.mie c c c c**********************************************************************c C comment the whole block with goto 499 goto 499 CX rmin=0. CX rmax=0. icp=1 do i=1,4 x1(i)=0.0 x2(i)=0.0 x3(i)=0.0 do l=1,20 rn(l,i)=0.0 ri(l,i)=0.0 enddo enddo do i=1,50 rsunph(i)=0. nrsunph(i)=0. enddo cij(1)=1.00 499 continue CX write(6,*)"cij",cij CX taer=0. CX taer55=0. CX iaer_prof=0 CMin read(iread,*) iaer c the user-defined aerosol profile if (iaer.lt.0) then total_height=0.0 iaer_prof=1 num_z=0 do i=0,50 alt_z(i)=0.0 taer55_z(i)=0.0 taer_z(i)=0.0 height_z(i)=0.0 enddo read(5,*) num_z do i=0,num_z-1 read(5,*) height_z(num_z-i),taer55_z(num_z-i),iaer alt_z(num_z-1-i)=total_height+height_z(num_z-i) total_height=total_height+height_z(num_z-i) taer55=taer55+taer55_z(num_z-i) enddo endif c the user-defined aerosol profile if (iaer.ge.0.and.iaer.le.7) nquad=nqdef_p if (iaer.ge.8.and.iaer.le.11) nquad=nquad_p if(iaer.eq.4) read(iread,*) (c(n),n=1,4) c(1)=0.25 c(2)=0.25 c(3)=0.25 c(4)=0.25 CX goto(49,40,41,42,49,49,49,49,43,44,45,46,47),iaer+1 goto(49,40,41,42,49,49,49,49,49,44,45,46,47),iaer+1 40 c(1)=0.70 c(2)=0.29 c(3)=0.00 c(4)=0.01 go to 49 41 c(1)=0.00 c(2)=0.05 c(3)=0.95 c(4)=0.00 go to 49 42 c(1)=0.17 c(2)=0.61 c(3)=0.00 c(4)=0.22 go to 49 CX 43 read(iread,*) rmin,rmax,icp CX do i=1,icp CX read(5,*)x1(i),x2(i),cij(i) CX read(5,*)(rn(l,i),l=1,20) CX read(5,*)(ri(l,i),l=1,20) CX enddo CX do i=1,icp CX cij_out(i)=cij(i) CX enddo go to 49 44 read(iread,*) rmin,rmax read(iread,*) x1(1),x2(1),x3(1) read(5,*)(rn(l,1),l=1,20) read(5,*)(ri(l,1),l=1,20) go to 49 45 read(iread,*) rmin,rmax read(iread,*) x1(1) read(5,*)(rn(l,1),l=1,20) read(5,*)(ri(l,1),l=1,20) go to 49 46 read(5,*)irsunph do i=1,irsunph read(5,*)rsunph(i),nrsunph(i) C nrsunph(i)=nrsunph(i)/(rsunph(i)**4.)/(4*3.1415/3) enddo rmin=rsunph(1) rmax=rsunph(irsunph)+1e-07 read(5,*)(rn(l,1),l=1,20) read(5,*)(ri(l,1),l=1,20) go to 49 C mino here is using in Dust_min.mie to print a line CX 47 write(6,*)"here" 47 continue CX 47 read(5,'(A80)')FILE2 CX i2=index(FILE2,' ')-1 go to 49 49 continue CX write(6,*)"herebefore" CX write(6,*)"rmax,rmin,rn",rmax,rmin,rn CX write(6,*)"ri",ri CX write(6,*)"x1,x2,x3,rsunph,nrsunph",x1,x2,x3,rsunph,nrsunph if (iaer.ge.8.and.iaer.le.11)then CX read(5,*)iaerp if (iaerp.eq.1)read(5,'(A80)')FILE i1=index(FILE,' ')-1 FILE2=FILE(1:I1)//'.mie' i2=index(FILE2,' ')-1 endif CX mino move aeroso to top for dustmie case CX comment the line below if IAER is NOT equal 12 call aeroso(iaer,c,xmud,wldis,FILE2,ipol) c**********************************************************************c c aerosol model (concentration) c c ---------------------------- c c (only for the default exponential profile) c c c c v if you have an estimate of the meteorological c c parameter: the visibility v, enter directly the c c value of v in km (the aerosol optical depth will c c be computed from a standard aerosol profile) c c c c v=0, taer55 if you have an estimate of aerosol optical depth , c c enter v=0 for the visibility and enter the aerosol c c optical depth at 550 c c c c v=-1 warning: if iaer=0, enter v=-1 c c c c**********************************************************************c if (iaer_prof.eq.0) then CX read(iread,*) v CX v=0 CX if(v) 71,10,11 CC if "v" is -1 goto 71, if v is 0 goto 10, if v is 1 goto 11 ( mino) CX 10 read(iread,*) taer55 v=exp(-log(taer55/2.7628)/0.79902) CX goto 71 CX 11 call oda550(iaer,v,taer55) CX 71 continue endif c**********************************************************************c c xps is the parameter to express the altitude of target c c c c c c xps >=0. the pressure is given in mb c c c c xps <0. means you know the altitude of the target c c expressed in km and you put that value as xps c c c c c c**********************************************************************c CX 771 read(iread,*) xps if (idatm.ne.8) then call pressure(uw,uo3,xps) else call pressure(uwus,uo3us,xps) endif c**********************************************************************c c c c xpp is the parameter to express the sensor altitude c c c c c c xpp= -1000 means that the sensor is a board a satellite c c xpp= 0 means that the sensor is at the ground level c c c c c c for aircraft simulations c c -100< xpp <0 means you know the altitude of the sensor expressed c c in kilometers units c c this altitude is relative to the target altitude c c c c for aircraft simulations only, you have to give c c puw,po3 (water vapor content,ozone content between the c c aircraft and the surface) c c taerp (the aerosol optical thickness at 550nm between the c c aircraft and the surface) c c if these data are not available, enter negative values for all c c of them, puw,po3 will then be interpolated from the us62 standard c C profile according to the values at ground level. Taerp will be c c computed according to a 2km exponential profile for aerosol. c c**********************************************************************c CX read(iread,*) xpp xpp=-xpp if (xpp.le.0.0) then c ground measurement option palt=0. pps=p(1) idatmp=0 taer55p=0. puw=0. puoz=0. else if (xpp.ge.100.) then c satellite case of equivalent palt=1000. pps=0. taer55p=taer55 ftray=1. idatmp=4 else c "real" plane case read(iread,*) puw,puo3 if (puw.lt.0.) then call presplane(puw,puo3,xpp,ftray) idatmp=2 if (idatm.eq.8) then puwus=puw puo3us=puo3 puw=puw*uw/uwus puo3=puo3*uo3/uo3us idatmp=8 endif else call presplane(puwus,puo3us,xpp,ftray) idatmp=8 endif if(ier) stop palt=zpl(34)-z(1) pps=ppl(34) read(iread,*) taer55p if ((taer55p.lt.0.).or.((taer55-taer55p).lt.accu2)) then c a scale heigh of 2km is assumed in case no value is given for taer55p taer55p=taer55*(1.-exp(-palt/2.)) else C compute effective scale heigh sham=exp(-palt/4.) sha=1.-(taer55p/taer55) if (sha.ge.sham) then taer55p=taer55*(1.-exp(-palt/4.)) else sha=-palt/log(sha) taer55p=taer55*(1.-exp(-palt/sha)) endif endif endif endif c**********************************************************************c c iwave input of the spectral conditions c c -------------------------------- c c c c you choose to define your own spectral conditions: iwave=-1,0 or 1 c c (three user s conditions ) c c -2 enter wlinf, wlsup, the filter function will be equal to 1c c over the whole band (as iwave=0) but step by step output c c will be printed c c -1 enter wl (monochr. cond, gaseous absorption is included) c c c c 0 enter wlinf, wlsup. the filter function will be equal to 1c c over the whole band. c c c c 1 enter wlinf, wlsup and user's filter function s(lambda) c c ( by step of 0.0025 micrometer). c c c c c c or you select one of the following satellite spectral bands: c c with indication in brackets of the band limits used in the code : c c iwave=2 to 60 c c 2 vis band of meteosat ( 0.350-1.110 ) c c 3 vis band of goes east ( 0.490-0.900 ) c c 4 vis band of goes west ( 0.490-0.900 ) c c 5 1st band of avhrr(noaa6) ( 0.550-0.750 ) c c 6 2nd " ( 0.690-1.120 ) c c 7 1st band of avhrr(noaa7) ( 0.500-0.800 ) c c 8 2nd " ( 0.640-1.170 ) c c 9 1st band of avhrr(noaa8) ( 0.540-1.010 ) c c 10 2nd " ( 0.680-1.120 ) c c 11 1st band of avhrr(noaa9) ( 0.530-0.810 ) c c 12 2nd " ( 0.680-1.170 ) c c 13 1st band of avhrr(noaa10 ( 0.530-0.780 ) c c 14 2nd " ( 0.600-1.190 ) c c 15 1st band of avhrr(noaa11 ( 0.540-0.820 ) c c 16 2nd " ( 0.600-1.120 ) c c 17 1st band of hrv1(spot1) ( 0.470-0.650 ) c c 18 2nd " ( 0.600-0.720 ) c c 19 3rd " ( 0.730-0.930 ) c c 20 pan " ( 0.470-0.790 ) c c 21 1st band of hrv2(spot1) ( 0.470-0.650 ) c c 22 2nd " ( 0.590-0.730 ) c c 23 3rd " ( 0.740-0.940 ) c c 24 pan " ( 0.470-0.790 ) c c 25 1st band of tm(landsat5) ( 0.430-0.560 ) c c 26 2nd " ( 0.500-0.650 ) c c 27 3rd " ( 0.580-0.740 ) c c 28 4th " ( 0.730-0.950 ) c c 29 5th " ( 1.5025-1.890 ) c c 30 7th " ( 1.950-2.410 ) c c 31 MSS band 1 (0.475-0.640) c c 32 MSS band 2 (0.580-0.750) c c 33 MSS band 3 (0.655-0.855) c c 34 MSS band 4 ( 0.785-1.100 ) c c 35 1st band of MAS (ER2) ( 0.5025-0.5875) c c 36 2nd " ( 0.6075-0.7000) c c 37 3rd " ( 0.8300-0.9125) c c 38 4th " ( 0.9000-0.9975) c c 39 5th " ( 1.8200-1.9575) c c 40 6th " ( 2.0950-2.1925) c c 41 7th " ( 3.5800-3.8700) c c 42 MODIS band 1 ( 0.6100-0.6850) c c 43 MODIS band 2 ( 0.8200-0.9025) c c 44 MODIS band 3 ( 0.4500-0.4825) c c 45 MODIS band 4 ( 0.5400-0.5700) c c 46 MODIS band 5 ( 1.2150-1.2700) c c 47 MODIS band 6 ( 1.6000-1.6650) c c 48 MODIS band 7 ( 2.0575-2.1825) c c 49 1st band of avhrr(noaa12 ( 0.500-1.000 ) c c 50 2nd " ( 0.650-1.120 ) c c 51 1st band of avhrr(noaa14 ( 0.500-1.110 ) c c 52 2nd " ( 0.680-1.100 ) c c 53 POLDER band 1 ( 0.4125-0.4775) c c 54 POLDER band 2 (non polar( 0.4100-0.5225) c c 55 POLDER band 3 (non polar( 0.5325-0.5950) c c 56 POLDER band 4 P1 ( 0.6300-0.7025) c c 57 POLDER band 5 (non polar( 0.7450-0.7800) c c 58 POLDER band 6 (non polar( 0.7000-0.8300) c c 59 POLDER band 7 P1 ( 0.8100-0.9200) c c 60 POLDER band 8 (non polar( 0.8650-0.9400) c c 61 SEAWIFS band 1 ( 0.3825-0.70) c c 62 SEAWIFS band 2 ( 0.3800-0.58) c c 63 SEAWIFS band 3 ( 0.3800-1.02) c c 64 SEAWIFS band 4 ( 0.3800-1.02) c c 65 SEAWIFS band 5 ( 0.3825-1.15) c c 66 SEAWIFS band 6 ( 0.3825-1.05) c c 67 SEAWIFS band 7 ( 0.3800-1.15) c c 68 SEAWIFS band 8 ( 0.3800-1.15) c c 69 AATSR band 1 ( 0.5250-0.5925) c c 70 AATSR band 2 ( 0.6275-0.6975) c c 71 AATSR band 3 ( 0.8325-0.9025) c c 72 AATSR band 4 ( 1.4475-1.7775) c c 73 MERIS band 01 ( 0.412) c c 74 MERIS band 02 ( 0.442) c c 75 MERIS band 03 ( 0.489) c c 76 MERIS band 04 ( 0.509) c c 77 MERIS band 05 ( 0.559) c c 78 MERIS band 06 ( 0.619) c c 79 MERIS band 07 ( 0.664) c c 80 MERIS band 08 ( 0.681) c c 81 MERIS band 09 ( 0.708) c c 82 MERIS band 10 ( 0.753) c c 83 MERIS band 11 ( 0.760) c c 84 MERIS band 12 ( 0.778) c c 85 MERIS band 13 ( 0.865) c c 86 MERIS band 14 ( 0.885) c c 87 MERIS band 15 ( 0.900) c c 88 GLI band 1 (0.380-1km) c c 89 GLI band 2 (0.400-1km) c c 90 GLI band 3 (0.412-1km) c c 91 GLI band 4 (0.443-1km) c c 92 GLI band 5 (0.460-1km) c c 93 GLI band 6 (0.490-1km) c c 94 GLI band 7 (0.520-1km) c c 95 GLI band 8 (0.545-1km) c c 96 GLI band 9 (0.565-1km) c c 97 GLI band 10 (0.625-1km) c c 98 GLI band 11 (0.666-1km) c c 99 GLI band 12 (0.680-1km) c c 100 GLI band 13 (0.678-1km) c c 101 GLI band 14 (0.710-1km) c c 102 GLI band 15 (0.710-1km) c c 103 GLI band 16 (0.749-1km) c c 104 GLI band 17 (0.763-1km) c c 105 GLI band 18 (0.865-1km) c c 106 GLI band 19 (0.865-1km) c c 107 GLI band 20 (0.460-0.25km) c c 108 GLI band 21 (0.545-0.25km) c c 109 GLI band 22 (0.660-0.25km) c c 110 GLI band 23 (0.825-0.25km) c c 111 GLI band 24 (1.050-1km) c c 112 GLI band 25 (1.135-1km) c c 113 GLI band 26 (1.240-1km) c c 114 GLI band 27 (1.338-1km) c c 115 GLI band 28 (1.640-1km) c c 116 GLI band 29 (2.210-1km) c c 117 GLI band 30 (3.715-1km) c c 118 ALI band 1p (0.4225-0.4625) c c 119 ALI band 1 (0.4325-0.550) c c 120 ALI band 2 (0.500-0.630) c c 121 ALI band 3 (0.5755-0.730) c c 122 ALI band 4 (0.7525-0.8375) c c 123 ALI band 4p (0.8025-0.935) c c 124 ALI band 5p (1.130-1.345) c c 125 ALI band 5 (1.470-1.820) c c 126 ALI band 7 (1.980-2.530) c c 127 ASTER band 1 (0.485-0.6425) c c 128 ASTER band 2 (0.590-0.730) c c 129 ASTER band 3n (0.720-0.9075) c c 130 ASTER band 3b (0.720-0.9225) c c 131 ASTER band 4 (1.570-1.7675) c c 132 ASTER band 5 (2.120-2.2825) c c 133 ASTER band 6 (2.150-2.295) c c 134 ASTER band 7 (2.210-2.390) c c 135 ASTER band 8 (2.250-2.440) c c 136 ASTER band 9 (2.2975-2.4875) c c 137 ETM band 1 (0.435-0.52) c c 138 ETM band 2 (0.5-0.6225) c c 139 ETM band 3 (0.615-0.7025) c c 140 ETM band 4 (0.74-0.9125) c c 141 ETM band 5 (1.51-1.7875) c c 142 ETM band 7 (2.015-2.3775) c c 143 HYPBLUE band 1 (0.4375-0.500) c c 144 HYPBLUE band 2 (0.435-0.52) c c 145 VGT band 1 (0.4175-0.500) c c 146 VGT band 2 (0.5975-0.7675) c c 147 VGT band 3 (0.7325-0.9575) c c 148 VGT band 4 (1.5225-1.800) c c 149 VIIRS band M1 (0.4025-0.4225) c c 150 VIIRS band M2 (0.4350-0.4550) c c 151 VIIRS band M3 (0.4775-0.4975) c c 152 VIIRS band M4 (0.5450-0.5650) c c 153 VIIRS band M5 (0.6625-0.6825) c c 154 VIIRS band M6 (0.7375-0.7525) c c 155 VIIRS band M7 (0.8450-0.8850) c c 156 VIIRS band M8 (1.2300-1.2500) c c 157 VIIRS band M9 (1.3700-1.3850) c c 158 VIIRS band M10 (1.5800-1.6400) c c 159 VIIRS band M11 (2.2250-2.2750) c c 160 VIIRS band M12 (3.6100-3.7900) c c 161 VIIRS band I1 (0.6000-0.6800) c c 162 VIIRS band I2 (0.8450-0.8850) c c 163 VIIRS band I3 (1.5800-1.6400) c c 164 VIIRS band I4 (3.5500-3.9300) c c 165 LDCM band 1 (0.4275-0.4575) c c 166 LDCM band 2 (0.4375-0.5275) c c 167 LDCM band 3 (0.5125-0.6000) c c 168 LDCM band 4 (0.6275-0.6825) c c 169 LDCM band 5 (0.8300-0.8950) c c 170 LDCM band 6 (1.5175-1.6950) c c 171 LDCM band 7 (2.0375-2.3500) c c 172 LDCM band 8 (0.4875-0.6925) c c 173 LDCM band 9 (1.3425-1.4025) c c 174 MODISkm band 8 (0.4025-0.4225) c c 175 MODISkm band 9 (0.4325-0.4500) c c 176 MODISkm band 10 (0.4775-0.4950) c c 177 MODISkm band 11 (0.5200-0.5400) c c 178 MODISkm band 12 (0.5375-0.5550) c c 179 MODISkm band 13 (0.6575-0.6750) c c 180 MODISkm band 14 (0.6675-0.6875) c c 181 MODISkm band 15 (0.7375-0.7575) c c 182 MODISkm band 16 (0.8525-0.8825) c c 183 MODISkm band 17 (0.8725-0.9375) c c 184 MODISkm band 18 (0.9225-0.9475) c c 185 MODISkm band 19 (0.8900-0.9875) c c 186 CAVIS band 1 (0.4275-0.4575) c c 187 CAVIS band 2 (0.4375-0.5275) c c 188 CAVIS band 3 (0.5125-0.6000) c c 189 CAVIS band 4 (0.6275-0.6825) c c 190 CAVIS band 5 (0.8300-0.8950) c c 191 CAVIS band 6 (1.3425-1.4025) c c 192 CAVIS band 7 (1.5175-1.6950) c c 193 CAVIS band 8 (2.0375-2.3500) c c 194 CAVIS band 9 (0.4875-0.6925) c c 195 CAVIS band 10 (0.4875-0.6925) c c 196 CAVIS band 11 (0.5100-0.6200) c c 197 DMC band 1 (0.4875-0.6925) c c 198 DMC band 2 (0.6100-0.7100) c c 199 DMC band 3 (0.7525-0.9275) c c 200 ABI band 1 (0.47) c c 201 ABI band 2 (0.67) c c 202 ABI band 3 (0.86) c c 203 ABI band 4 (1.38) c c 204 ABI band 5 (1.60) c c 205 ABI band 6 (2.30) c c note: wl has to be in micrometer c c**********************************************************************c do 38 l=iinf,isup s(l)=1. 38 continue CX read(iread,*) iwave if (iwave.eq.-2) goto 1600 if (iwave) 16,17,18 16 read(iread,*) wl wlinf=wl wlsup=wl go to 19 17 read(iread,*) wlinf,wlsup go to 19 1600 read(iread,*) wlinf,wlsup go to 19 c 110 c 111 band of meteosat (2) c 112 band of goes (3,4) c 114 band of avhr (5,16) c 118 band of hrv1 (17,24) c 121 band of tm (25,30) c 127 band of mss (31,34) c 128 band of MAS (35,41) c 129 MODIS band (42,49) c 130 band of avhrr (50,53) c 131 POLDER band (54,61) c 113 SEAWIFS band (62,69) c 150 AATSR band (70,73) c 151 MERIS band (74,88) c 152 GLI band (89,118) c 153 ALI band (119,127) c 154 ASTER band (128,137) c 155 ETM band (138,143) c 156 HYPBLUE band (144,145) c 157 VGT band (146,149) c 159 VIIRS band (149,164) c 161 LDCM band (165,173) c 162 MODIS1km band (174,185) c 163 CAVIS band (186,196) c 164 DMC band (197,199) c 165 ABI band (200,205) 18 goto (110, s 111, s 112,112, s 114,114,114,114,114,114,114,114,114,114,114,114, s 118,118,118,118,118,118,118,118, s 121,121,121,121,121,121, s 127,127,127,127, s 128,128,128,128,128,128,128, s 129,129,129,129,129,129,129, s 130,130,130,130, s 131,131,131,131,131,131,131,131, s 113,113,113,113,113,113,113,113, s 150,150,150,150, s 151,151,151,151,151,151,151,151, s 151,151,151,151,151,151,151, s 152,152,152,152,152,152,152,152,152,152, s 152,152,152,152,152,152,152,152,152,152, s 152,152,152,152,152,152,152,152,152,152, s 153,153,153,153,153,153,153,153,153, s 154,154,154,154,154,154,154,154,154,154, s 155,155,155,155,155,155, s 156,156, s 157,157,157,157, s 159,159,159,159,159,159,159,159,159,159, s 159,159,159,159,159,159, s 161,161,161,161,161,161,161,161,161, s 162,162,162,162,162,162,162,162,162,162,162,162, s 163,163,163,163,163,163,163,163,163,163,163, s 164,164,164, s 165,165,165,165,165,165),iwave 110 read(iread,*) wlinf,wlsup iinf=(wlinf-.25)/0.0025+1.5 isup=(wlsup-.25)/0.0025+1.5 do 1113 ik=iinf,isup s(ik)=0. 1113 continue read(iread,*) (s(i),i=iinf,isup) goto 20 111 call meteo go to 19 112 call goes(iwave-2) go to 19 114 call avhrr(iwave-4) go to 19 118 call hrv(iwave-16) go to 19 121 call tm(iwave-24) go to 19 127 call mss(iwave-30) goto 19 128 call mas(iwave-34) goto 19 129 call modis(iwave-41) goto 19 130 call avhrr(iwave-48) goto 19 131 call polder(iwave-52) goto 19 113 call seawifs(iwave-60) goto 19 150 call aatsr(iwave-68) goto 19 151 call meris(iwave-72) goto 19 152 call gli(iwave-87) goto 19 153 call ali(iwave-117) goto 19 154 call aster(iwave-126) goto 19 155 call etm(iwave-136) goto 19 156 call hypblue(iwave-142) goto 19 157 call vgt(iwave-144) goto 19 159 call viirs(iwave-148) goto 19 161 call ldcm(iwave-164) goto 19 162 call modis1km(iwave-173) goto 19 163 call cavis(iwave-185) goto 19 164 call dmc(iwave-196) goto 19 165 call abi(iwave-199) goto 19 19 iinf=(wlinf-.25)/0.0025+1.5 isup=(wlsup-.25)/0.0025+1.5 if (iprtspr.eq.1) then do i=iinf,isup write(6,*) "spres ",(i-1)*0.0025+0.25,s(i) enddo endif CX write(6,*) "wlinf,wlsup ****** ",wlinf,wlsup 20 continue C*********************************************************************** C LOOK UP TABLE INITIALIZATION C*********************************************************************** C initialization of look up table variable C Write(6,*) "TOTO THE HERO" do i=1,mu nfilut(i)=0 do j=1,41 rolut(i,j)=0. rolutq(i,j)=0. rolutu(i,j)=0. filut(i,j)=0. roluti(i,j)=0. rolutiq(i,j)=0. rolutiu(i,j)=0. enddo enddo xmus=cos(asol*pi/180.) its=acos(xmus)*180.0/pi C Case standart LUT if (ilut.eq.1) then do i=1,mu-1 lutmuv=rm(i) luttv=acos(lutmuv)*180./pi iscama=(180-abs(luttv-its)) iscami=(180-(luttv+its)) nbisca=int(0.01+(iscama-iscami)/4.0)+1 nfilut(i)=nbisca filut(i,1)=0.0 filut(i,nbisca)=180.0 scaa=iscama do j=2,nfilut(i)-1 scaa=scaa-4.0 cscaa=cos(scaa*pi/180.) cfi=-(cscaa+xmus*lutmuv)/(sqrt(1-xmus*xmus) S *sqrt(1.-lutmuv*lutmuv)) filut(i,j)=acos(cfi)*180.0/pi enddo enddo i=mu lutmuv=cos(avis*pi/180.) luttv=acos(lutmuv)*180./pi iscama=(180-abs(luttv-its)) iscami=(180-(luttv+its)) nbisca=int((iscama-iscami)/4)+1 nfilut(i)=nbisca filut(i,1)=0.0 filut(i,nbisca)=180.0 scaa=iscama do j=2,nfilut(i)-1 scaa=scaa-4.0 cscaa=cos(scaa*pi/180.) cfi=-(cscaa+xmus*lutmuv)/(sqrt(1-xmus*xmus) S *sqrt(1.-lutmuv*lutmuv)) filut(i,j)=acos(cfi)*180.0/pi enddo endif C END Case standart LUT C Case LUT for APS if (ilut.eq.3) then do i=1,mu-1 nbisca=2 nfilut(i)=nbisca filut(i,1)=(phi0-phiv) filut(i,nbisca)=(phi0-phiv)+180.0 enddo i=mu nbisca=1 nfilut(i)=nbisca filut(i,1)=(phi0-phiv) endif C END Case LUT for APS CCCC Check initialization (debug) do i=1,mu lutmuv=rm(i) luttv=acos(lutmuv)*180./pi do j=1,nfilut(i) cscaa=-xmus*lutmuv-cos(filut(i,j)*pi/180.)*sqrt(1.-xmus*xmus) S *sqrt(1.-lutmuv*lutmuv) scaa=acos(cscaa)*180./pi write(6,*) its,luttv,filut(i,j),scaa enddo enddo CCCC Check initialization (debug) C*********************************************************************** C END LOOK UP TABLE INITIALIZATION C*********************************************************************** c**********************************************************************c c here, we first compute an equivalent wavelenght which is the input c c value for monochromatic conditions or the integrated value for a c c filter functionr (call equivwl) then, the atmospheric properties are c c computed for that wavelength (call discom then call specinterp) c c molecular optical thickness is computed too (call odrayl). lastly c c the successive order of scattering code is called three times. c c first for a sun at thetas with the scattering properties of aerosols c c and molecules, second with a pure molecular atmosphere, then with thec c actual atmosphere for a sun at thetav. the iso code allows us to c c compute the scattering transmissions and the spherical albedo. all c c these computations are performed for checking the accuracy of the c c analytical expressions and in addition for computing the averaged c c directional reflectances c c**********************************************************************c if(iwave.ne.-1) then call equivwl(iinf,isup,step, s wlmoy) else wlmoy=wl endif call discom (idatmp,iaer,iaer_prof,xmus,xmuv,phi,taer55,taer55p, a palt,phirad,nt,mu,np,rm,gb,rp,ftray,ipol,xlm1,xlm2, a roatm_fi,nfi, a nfilut,filut,roluts,rolutsq,rolutsu) CX write(6,*) "iaer,wlmoy,aer55,taer55p,tamoy",iaer,wlmoy,aer55,taer55p,tamoy CX write(6,*) "iaer,wlmoy,aer55,taer55p,tamoy,tamoyp,pizmoy,pizmoyp,ipol",iaer,wlmoy,taer55,taer55p,tamoy,tamoyp,pizmoy,pizmoyp,ipol if(iaer.ne.0) then call specinterp(wlmoy,taer55,taer55p, s tamoy,tamoyp,pizmoy,pizmoyp,ipol) CX write(6,*) "iaer,wlmoy,aer55,taer55p,tamoy",iaer,wlmoy,aer55,taer55p,tamoy CX write(6,*) "iaer,wlmoy,aer55,taer55p,tamoy,tamoyp,pizmoy,pizmoyp,ipol",iaer,wlmoy,taer55,taer55p,tamoy,tamoyp,pizmoy,pizmoyp,ipol else tamoy=0. tamoyp=0. endif CX write(6,*) "wlmoy",wlmoy CX write(6,*) "tamoy",tamoy call odrayl(wlmoy, s trmoy) trmoyp=trmoy*ftray if (idatmp.eq.4) then trmoyp=trmoy tamoyp=tamoy endif if (idatmp.eq.0) then trmoyp=0. tamoyp=0. endif c*********************************************************************c c inhomo ground reflectance (type) c c ------------------ c c c c you consider an homogeneous surface: c c enter - inhomo=0 c c you may consider directional surface effects c c idirec=0 (no directional effect) c c you have to specify the surface reflectance:c c igroun (see note1) which is uniform and c c lambertian c c idirec=1 ( directional effect) c c you have to specify the brdf of the surface c c for the actual solar illumination you are c c considering as well as the  for a sun c c which would be at an angle thetav, in c c addition you have to give the surface c c albedo (spherical albedo). you can also c c select one of the selected model from the c c ibrdf value (see note2). 3 reflectances c c are computed, robar,robarp and robard c c c c you consider a non uniform surface, the surface is considered as a c c circular target with a reflectance roc and of radius r c c (expressed in km) within an environment of reflectance c c roe c c enter - inhomo=1, then c c igrou1,igrou2,rad c c - the target reflectance :igrou1 (see note1) c c - the envir. reflectance :igrou2 (see note1) c c - the target radius in km c c c c c c ****tree**** c c c c inhomo c c / \ c c / \ c c / \ c c / \ c c ------- 0 ------- -----1 ----- c c / / \ \ c c idirec / \ \ c c / \ / \ \ c c / \ / \ \ c c / \ igrou1 igrou2 rad c c 0 1 roc roe f(r) c c / \ c c / \ c c igroun ibrdf c c (roc = roe) (roc) c c (robar) c c (robarp) c c (robard) c c c c ground reflectance (spectral variation) c c --------------------------------------- c c note1: values of the reflectance selected by igroun,igrou1 or igrou2 c c may correspond to the following cases, c c 0 constant value of ro (or roc,or roe) whatever the wavelen c c gth. you enter this constant value of ro (or roc or roe). c c -1 you have to enter the value of ro (or roc,or roe) by step c c of 0.0025 micron from wlinf to wlsup (if you have used thec c satellite bands,see implicit values for these limits). c c 1 mean spectral value of green vegetation c c 2 mean spectral value of clear water c c 3 mean spectral value of sand c c 4 mean spectral value of lake water c c c c ground reflectance (brdf) c c ------------------------- c c note2: values of the directional reflectance is assumed spectrally c c independent, so you have to specify, the brdf at the c c wavelength for monochromatic condition of the mean value c c over the spectral band c c 0 you have to enter the value of ro for sun at thetas by c c step of 10 degrees for zenith view angles (from 0 to 80 c c and the value for 85) and by step of 30 degrees for c c azimuth view angles from 0 to 360 degrees, you have to do c c same for a sun which would be at thetav. in addition, the c c spherical albedo of the surface has to be specified ,as c C well as the observed reflectance in the selected geometry c c rodir(sun zenith,view zenith, relative azimuth). c c c c you also may select one of the following models c c 1 hapke model c c the parameters are: om,af,s0,h c c om= albedo c c af=assymetry parameter for the phase function c c s0=amplitude of hot spot c c h=width of the hot spot c c c c 2 verstraete et al. model c c the parameters are: c c there is three lines of parameters: c c line 1 (choice of options) c c line 2 (structural parameters) c c line 3 (optical parameters) c c line 1: opt3 opt4 opt5 c c opt1=1 parametrized model (see verstraete et al., c c JGR, 95, 11755-11765, 1990) c c opt2=1 reflectance factor (see pinty et al., JGR, c c 95, 11767-11775, 1990) c c opt3=0 for given values of kappa (see struc below)c c 1 for goudriaan's parameterization of kappa c c 2 for dickinson et al's correction to c c goudriaan's parameterization of kappa (see c c dickinson et al., agricultural and forest c c meteorology, 52, 109-131, 1990) c c ---see the manual for complete references---- c c opt4=0 for isotropic phase function c c 1 for heyney and greensteins' phase function c c 2 for legendre polynomial phase function c c opt5=0 for single scattering only c c 1 for dickinson et al. parameterization of c c multiple scattering c c line 2: str1 str2 str3 str4 c c str1='leaf area density', in m2 m-3 c c str2=radius of the sun flecks on the scatterer (m)c c str3=leaf orientation parameter: c c if opt3=0 then str3=kappa1 c c if opt3=1 or 2 then str3=chil c c str4=leaf orientation parameter (continued): c c if opt3=0 then str4=kappa2 c c if opt3=1 or 2 then str4 is not used c c line 3: optics1 optics2 optics3 c c optics1=single scattering albedo, n/d value c c between 0.0 and 1.0 c c optics2= phase function parameter: c c if opt4=0 then this input is not used c c if opt4=1 then asymmetry factor, n/d value c c between -1.0and 1.0 c c if opt4=2 then first coefficient of legendre c c polynomial c c optics3=second coefficient of legendre polynomial c c (if opt4=2) c c c c 3 Roujean et al. model c c the parameters are: k0,k1,k2 c c k0=albedo. c c k1=geometric parameter for hot spot effect c c k2=geometric parameter for hot spot effect c c c c 4 walthall et al. model c c the parameters are: a,ap,b,c c c a=term in square ts*tv c c ap=term in square ts*ts+tv*tv c c b=term in ts*tv*cos(phi) (limacon de pascal) c c c=albedo c c c c 5 minnaert model c c the parameters are: par1,par2 c c c c 6 Ocean c c the parameter are: pws,phi_wind,xsal,pcl c c pws=wind speed (in m/s) c c phi_wind=azim. of the wind (in degres) c c xsal=salinity (in ppt) xsal=34.3ppt if xsal<0 c c pcl=pigment concentration (in mg/m3) c c c c 7 Iaquinta and Pinty model c c the parameters are: c c there is 3 lines of parameters: c c line 1: choice of option (pild,pihs) c c line 2: structural parameters (pxLt,pc) c c line 3: optical parameters (pRl,pTl,pRs) c c Line 1: pild,pihs c c pild=1 planophile leaf distribution c c pild=2 erectophile leaf distribution c c pild=3 plagiophile leaf distribution c c pild=4 extremophile leaf distribution c c pild=5 uniform leaf distribution c c c c pihs=0 no hot spot c c pihs=1 hot spot c c Line 2: pxLt,pc c c pxLt=Leaf area index [1.,15.] c c pc=Hot spot parameter: 2*r*Lambda [0.,2.] c c Line 3: pRl,pTl,pRs c c pRl=Leaf reflectance [0.,0.99] c c pTl=Leaf transmitance [0.,0.99] c c pRs=Soil albedo [0.,0.99] c c NB: pRl+PTl <0.99 c c c c 8 Rahman et al. model c c the parameters are: rho0,af,xk c c rho0=Intensity of the reflectance of the surface c c cover, N/D value greater or equal to 0 c c af=Asymmetry factor, N/D value between -1.0 and 1.0 c c xk=Structural parameter of the medium c c 9 Kuusk's multispectral CR model c c Reference: c c Kuusk A. A multispectral canopy reflectance model. c c Remote Sens. Environ., 1994, 50:75-82 c c c c c c the parameters are: c c c c line 1: structural parameters (ul,eps,thm,sl) c c line 2: optical parameters (cAB,cW,N,cn,s1) c c c c ul=LAI [0.1...10] c c eps,thm - LAD parameters c c eps [0.0..0.9] thm [0.0..90.0] c c sl - relative leaf size [0.01..1.0] c c cAB - chlorophyll content, ug/cm^2 [30] c c cW - leaf water equivalent thickness [0.01..0.03] c c N - the effective number of elementary layers c c inside a leaf [1.225] c c cn - the ratio of refractive indices of the leaf c c surface wax and internal material [1.0] c c s1 - the weight of the 1st Price function for the c c soil reflectance [0.1..0.8] c c 10 MODIS operational BDRF c c the parameters are: p1,p2,p3 c c p1 weight for lambertian kernel c c p2 weight for Ross Thick kernel c c p3 weight for Li Sparse kernel c c 11 RossLiMaigan BDRF model c c the parameters are: p1,p2,p3 c c p1 weight for lambertian kernel c c p2 weight for Ross Thick with Hot Spot kernel c c p3 weight for Li Sparse kernel c c**********************************************************************c fr=0. rad=0. do 1116 ik=iinf,isup rocl(ik)=0. roel(ik)=0. 1116 continue c**********************************************************************c c uniform or non-uniform surface conditions c c**********************************************************************c CX read(iread,*) inhomo if(inhomo) 30,30,31 CX 30 read(iread,*) idirec 30 continue if(idirec)21,21,25 c**********************************************************************c c uniform conditions with brdf conditions c c**********************************************************************c c 25 read(iread,*) ibrdf c*********************************************************************c cCO land model does not need brdf if(ibrdf)23,23,24 c**********************************************************************c c brdf from in-situ measurements c c**********************************************************************c 23 do 900 k=1,13 read(iread,*) (brdfdats(10-j+1,k),j=1,10) 900 continue do 901 k=1,13 read(iread,*) (brdfdatv(10-j+1,k),j=1,10) 901 continue read(iread,*) albbrdf read(iread,*) rodir rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call brdfgrid(mu,np,rm,rp,brdfdats,angmu,angphi, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call brdfgrid(mu,np,rm,rp,brdfdatv,angmu,angphi, s brdfintv) brdfints(mu,1)=rodir do l=iinf,isup sbrdf(l)=rodir enddo go to 69 c**********************************************************************c c brdf from hapke's model c c**********************************************************************c 24 if(ibrdf.eq.1) then read(iread,*) par1,par2,par3,par4 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call hapkbrdf(par1,par2,par3,par4,1,1,srm,srp, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call hapkbrdf(par1,par2,par3,par4,mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call hapkbrdf(par1,par2,par3,par4,mu,np,rm,rp, s brdfintv) call hapkalbe(par1,par2,par3,par4, s albbrdf) go to 69 endif c**********************************************************************c c brdf from verstraete et al's model c c**********************************************************************c if(ibrdf.eq.2) then read(iread,*) (options(i),i=3,5) options(1)=1 options(2)=1 read(iread,*) (struct(i),i=1,4) read(iread,*) (optics(i),i=1,3) srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call versbrdf(options,optics,struct,1,1,srm,srp, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call versbrdf(options,optics,struct,mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call versbrdf(options,optics,struct,mu,np,rm,rp, s brdfintv) call versalbe(options,optics,struct, s albbrdf) go to 69 endif c**********************************************************************c c brdf from Roujean et al's model c c**********************************************************************c if(ibrdf.eq.3) then read(iread,*) par1,par2,par3 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call roujbrdf(par1,par2,par3,1,1,srm,srp, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call roujbrdf(par1,par2,par3,mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call roujbrdf(par1,par2,par3,mu,np,rm,rp, s brdfintv) call roujalbe(par1,par2,par3, s albbrdf) go to 69 endif c**********************************************************************c c brdf from walthall et al's model c**********************************************************************c if(ibrdf.eq.4) then read(iread,*) par1,par2,par3,par4 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call waltbrdf(par1,par2,par3,par4,1,1,srm,srp, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call waltbrdf(par1,par2,par3,par4,mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call waltbrdf(par1,par2,par3,par4,mu,np,rm,rp, s brdfintv) call waltalbe(par1,par2,par3,par4, s albbrdf) go to 69 endif c**********************************************************************c c brdf from minnaert's model c c**********************************************************************c if(ibrdf.eq.5) then read(iread,*) par1,par2 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call minnbrdf(par1,par2,1,1,srm, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call minnbrdf(par1,par2,mu,np,rm, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call minnbrdf(par1,par2,mu,np,rm, s brdfintv) call minnalbe(par1,par2, s albbrdf) go to 69 endif c**********************************************************************c c brdf from ocean condition c**********************************************************************c if(ibrdf.eq.6) then CX read(iread,*) pws,phi_wind,xsal,pcl CX write(6,*) iaer_prof,tamoy,trmoy if (xsal.lt.0.001)xsal=34.3 paw=phi0-phi_wind do l=iinf,isup srm(-1)=phirad srm(1)=xmuv srm(0)=xmus wl=.25+(l-1)*step call oceabrdf(pws,paw,xsal,pcl,wl,rfoam,rwat,rglit, s 1,1,srm,srp, s sbrdftmp) rfoaml(l)=rfoam rwatl(l)=rwat rglitl(l)=rglit sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call oceabrdf(pws,paw,xsal,pcl,wlmoy,rfoam,rwat,rglit, s mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call oceabrdf(pws,paw,xsal,pcl,wlmoy,rfoam,rwat,rglit, s mu,np,rm,rp, s brdfintv) call oceaalbe(pws,paw,xsal,pcl,wlmoy, s albbrdf) CX write(6,*) iaer_prof,tamoy,trmoy go to 69 endif c c**********************************************************************c c brdf from Iaquinta and Pinty model c**********************************************************************c if(ibrdf.eq.7) then read(iread,*) pild,pihs read(iread,*) pxLt,pc read(iread,*) pRl,pTl,pRs srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call iapibrdf(pild,pxlt,prl,ptl,prs,pihs,pc,1,1,srm,srp, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call iapibrdf(pild,pxlt,prl,ptl,prs,pihs,pc,mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call iapibrdf(pild,pxlt,prl,ptl,prs,pihs,pc,mu,np,rm,rp, s brdfintv) call iapialbe(pild,pxlt,prl,ptl,prs,pihs,pc, s albbrdf) go to 69 endif c c**********************************************************************c c brdf from Rahman model c**********************************************************************c if(ibrdf.eq.8) then read(iread,*) par1,par2,par3 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call rahmbrdf(par1,par2,par3,1,1,srm,srp, s sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call rahmbrdf(par1,par2,par3,mu,np,rm,rp, s brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call rahmbrdf(par1,par2,par3,mu,np,rm,rp, s brdfintv) call rahmalbe(par1,par2,par3, s albbrdf) c call for ground boundary condition in OSSURF rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call rahmbrdffos(par1,par2,par3,mu,rm,rosur, s wfisur,fisur) c write(6,*) "rosur ",rosur go to 69 endif c c**********************************************************************c c brdf from kuusk's msrm model c c**********************************************************************c if(ibrdf.eq.9) then read(iread,*) uli,eei,thmi,sli read(iread,*) cabi,cwi,vaii,rnci,rsl1i do l=iinf,isup srm(-1)=phirad srm(1)=xmuv srm(0)=xmus wl=.25+(l-1)*step call akbrdf(eei,thmi,uli,sli,rsl1i,wl,rnci,cabi,cwi,vaii s ,1,1,srm,srp,sbrdftmp) sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call akbrdf(eei,thmi,uli,sli,rsl1i,wlmoy,rnci,cabi,cwi,vaii & ,mu,np,rm,rp,brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call akbrdf(eei,thmi,uli,sli,rsl1i,wlmoy,rnci,cabi,cwi,vaii & ,mu,np,rm,rp,brdfintv) c call akalbe * & (eei,thmi,uli,sli,rsl1i,wlmoy,rnci,cabi,cwi,vaii,albbrdf) & (albbrdf) go to 69 endif c c**********************************************************************c c brdf from MODIS BRDF model c c**********************************************************************c if(ibrdf.eq.10) then read(iread,*)p1,p2,p3 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call modisbrdf(p1,p2,p3 s ,1,1,srm,srp,sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call modisbrdf(p1,p2,p3 & ,mu,np,rm,rp,brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call modisbrdf(p1,p2,p3 & ,mu,np,rm,rp,brdfintv) c call modisalbe(p1,p2,p3 & ,albbrdf) c call for ground boundary condition in OSSURF rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call modisbrdffos(p1,p2,p3,mu,rm, s rosur,wfisur,fisur) go to 69 endif c c**********************************************************************c c brdf from ROSSLIMAIGNAN BRDF model c c**********************************************************************c if(ibrdf.eq.11) then read(iread,*)p1,p2,p3 srm(-1)=phirad srm(1)=xmuv srm(0)=xmus call rlmaignanbrdf(p1,p2,p3 s ,1,1,srm,srp,sbrdftmp) do l=iinf,isup sbrdf(l)=sbrdftmp(1,1) enddo c stop c rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call rlmaignanbrdf(p1,p2,p3 & ,mu,np,rm,rp,brdfints) rm(-mu)=2.*pi-phirad rm(mu)=xmus rm(0)=xmuv call rlmaignanbrdf(p1,p2,p3 & ,mu,np,rm,rp,brdfintv) c call rlmaignanalbe(p1,p2,p3 & ,albbrdf) c write(6,*) "GOT TILL HERE " c call for ground boundary condition in OSSURF rm(-mu)=phirad rm(mu)=xmuv rm(0)=xmus call rlmaignanbrdffos(p1,p2,p3,mu,rm, s rosur,wfisur,fisur) c do i=0,mu c do j=1,mu c do k=1,83 c write(6,*) i,j,k,rosur(i,j,k),acos(rm(i))*180./pi,acos(rm(j))*180./pi,fisur(k)*180./pi+180. c enddo c enddo c enddo go to 69 endif c c 69 continue c**********************************************************************c c compute the downward irradiance for a sun at thetas and then at c c tetav c c**********************************************************************c c call os to compute downward radiation field for robar CX IRON-MAN rm(-mu)=-xmuv rm(mu)=xmuv rm(0)=-xmus spalt=1000. CX write(6,*) "tamoy,trmoy",tamoy,trmoy c write(6,*) iaer_prof,tamoy,trmoy,pizmoy,tamoyp,trmoyp,spalt, c s phirad,nt,mu,np,rm,gb,rp, c s xlmus,xlphim,nfi,rolut call os(iaer_prof,tamoy,trmoy,pizmoy,tamoyp,trmoyp,spalt, s phirad,nt,mu,np,rm,gb,rp, s xlmus,xlphim,nfi,rolut) c write(6,*) "xlmus",xlmus romixatm=(xlmus(-mu,1)/xmus) c write(6,*) "romix atm", romix,tamoy,trmoy,phirad c call os to compute downward radiation field for robarp if (idatmp.ne.0) then rm(-mu)=-xmus rm(mu)=xmus rm(0)=-xmuv call os(iaer_prof,tamoyp,trmoyp,pizmoy,tamoyp,trmoyp,spalt, s phirad,nt,mu,np,rm,gb,rp, s xlmuv,xlphim,nfi,rolut) endif c call ossurf to compute the actual brdf coupling rm(-mu)=-xmuv rm(mu)=xmuv rm(0)=-xmus spalt=1000. call ossurf(iaer_prof,tamoyp,trmoyp,pizmoy,tamoyp,trmoyp,spalt, s phirad,nt,mu,np,rm,gb,rp,rosur,wfisur,fisur, s xlsurf,xlphim,nfi,rolutsurf) romixsur=(xlsurf(-mu,1)/xmus)-romixatm c write(6,*) "romix surf", romix c call ISO (twice) to compute the spherical albedo for the equivalent wavelength c and diffuse and direct transmission at equivalent vavelength rm(-mu)=-xmuv rm(mu)=xmuv rm(0)=xmus call iso(iaer_prof,tamoyp,trmoyp,pizmoy,tamoyp,trmoyp,spalt, a nt,mu,rm,gb,lxtrans) ludiftt=lxtrans(1)-exp(-(tamoyp+trmoyp)/xmuv) ludirtt=exp(-(tamoyp+trmoyp)/xmuv) rm(-mu)=-xmus rm(mu)=xmus rm(0)=xmus call iso(iaer_prof,tamoyp,trmoyp,pizmoy,tamoyp,trmoyp,spalt, a nt,mu,rm,gb,lxtrans) lddiftt=lxtrans(1)-exp(-(tamoyp+trmoyp)/xmus) lddirtt=exp(-(tamoyp+trmoyp)/xmus) lsphalbt=lxtrans(0)*2. c write(6,*) "sphalbt ddiftt ddirtt udiftt udirtt", c a lsphalbt,lddiftt,lddirtt,ludiftt,ludirtt,xmus,xmuv c stop c**********************************************************************c c the downward irradiance was computed for a sun at thetas and c c several viewing directions (mu zenith times np azimuth). then, the c c code computes the product of ldown*brdf integrated over the total c c hemisphere and gives the averaged directional reflectance after the c c normalization. the resulting reflectance is named robar c c**********************************************************************c robar1=0. xnorm1=0. CX write(6,*) "xlmus",xlmus do 83 j=1,np rob=0. xnor=0. do 84 k=1,mu-1 rdown=xlmus(-k,j) rdir=brdfintv(k,j) rob=rob+rdown*rdir*rm(k)*gb(k) xnor=xnor+rdown*rm(k)*gb(k) 84 continue robar1=robar1+rob*gp(j) xnorm1=xnorm1+xnor*gp(j) 83 continue c**********************************************************************c c the downward irradiance was computed for a sun at thetav and c c several viewing directions (mu zenith times np azimuth). then, the c c code computes the product of ldown*brdf integrated over the total c c hemisphere and gives the averaged directional reflectance after the c c normalization. the resulting reflectance is named robarp c c**********************************************************************c robar2=0. xnorm2=0. do 85 j=1,np rob=0. xnor=0. do 86 k=1,mu-1 rdown=xlmuv(-k,j) rdir=brdfints(k,j) rob=rob+rdown*rdir*rm(k)*gb(k) xnor=xnor+rdown*rm(k)*gb(k) 86 continue robar2=robar2+rob*gp(j) xnorm2=xnorm2+xnor*gp(j) 85 continue Write(6,*) "ROBAR",robar1,robar2,xnorm1,xnorm2,romix c robard is assumed equal to albbrdf c print 301,brdfints(mu,1),robar1,xnorm1, c s robar2,xnorm2,albbrdf c print 301,robar1/xnorm1,robar2/xnorm2 c print 301,betal(0)/3,pizmoy c301 format(6(f10.4,2x)) c501 format(5(i10,2x)) rbar=robar1/xnorm1 rbarp=robar2/xnorm2 rbarc=rbar*lddiftt*ludirtt rbarpc=rbarp*ludiftt*lddirtt rdirc=sbrdftmp(1,1)*ludirtt*lddirtt CX write(6,*) "romixsur,rbarc,rbarpc,rdirc",romixsur,rbarc,rbarpc,rdirc coefc=-(romixsur-rbarc-rbarpc-rdirc) c write(6,*) " lddiftt,ludiftt ", lddiftt,ludiftt coefb=lddiftt*ludiftt coefa=(lddiftt+lddirtt)*(ludiftt+ludirtt)*lsphalbt a /(1.-lsphalbt*albbrdf) CX write(6,*) "a,b,c",coefa,coefb,coefc CX write(6,*) "discri2 ",(coefb*coefb-4*coefa*coefc) c ADD by JC to avoid negative discriminant to determine the rhobarbar xtestdiscri=(coefb*coefb-4*coefa*coefc) if (xtestdiscri.lt.0.0)then rbard=albbrdf goto 487 Endif discri=sqrt(coefb*coefb-4*coefa*coefc) rbard=(-coefb+discri)/(2*coefa) CX Write(6,*) "rbard albbrdf 1rst iteration", rbard,albbrdf coefa=(lddiftt+lddirtt)*(ludiftt+ludirtt)*lsphalbt a /(1.-lsphalbt*rbard) CX write(6,*) "a,b,c",coefa,coefb,coefc CX write(6,*) "discri2 ",(coefb*coefb-4*coefa*coefc) discri=sqrt(coefb*coefb-4*coefa*coefc) rbard=(-coefb+discri)/(2*coefa) c Write(6,*) "rbard albbrdf 2nd iteration", rbard,albbrdf 487 continue do 335 l=iinf,isup rocl(l)=sbrdf(l) roel(l)=sbrdf(l) robar(l)=robar1/xnorm1 if (idatmp.ne.0) then robarp(l)=robar2/xnorm2 else robarp(l)=0. xnorm2=1. robar2=0. endif robard(l)=albbrdf C test version with Istvan (NOAA) CX robard(l)=rbard 335 continue go to 34 c**********************************************************************c c uniform surface with lambertian conditions c c**********************************************************************c CX 21 read(iread,*) igroun 21 if(igroun) 29,32,33 29 read(iread,*) nwlinf,nwlsup niinf=(nwlinf-.25)/0.0025+1.5 nisup=(nwlsup-.25)/0.0025+1.5 read(iread,*) (rocl(i),i=niinf,nisup) goto 36 32 read(iread,*) ro do 35 l=iinf,isup rocl(l)=ro 35 continue goto 36 igroun=igroun_con 33 if(igroun.eq.1) call vegeta(rocl) if(igroun.eq.2) call clearw(rocl) if(igroun.eq.3) call sand (rocl) if(igroun.eq.4) call lakew (rocl) 36 do 39 l=iinf,isup roel(l)=rocl(l) 39 continue go to 34 c**********************************************************************c c non-uniform conditions with lambertian conditions c c**********************************************************************c 31 read(iread,*) igrou1,igrou2,rad igrou1=2 igrou2=1 rad=0.5 if(igrou1) 59,60,63 59 read(iread,*) (rocl(i),i=iinf,isup) goto 61 60 read(iread,*) roc do 64 l=iinf,isup rocl(l)=roc 64 continue go to 61 63 if(igrou1.eq.1) call vegeta(rocl) if(igrou1.eq.2) call clearw(rocl) if(igrou1.eq.3) call sand (rocl) if(igrou1.eq.4) call lakew (rocl) 61 if(igrou2) 66,62,65 66 read(iread,*) (roel(i),i=iinf,isup) goto 34 62 read(iread,*) roe do 67 l=iinf,isup roel(l)=roe 67 continue go to 34 65 if(igrou2.eq.1) call vegeta(roel) if(igrou2.eq.2) call clearw(roel) if(igrou2.eq.3) call sand (roel) if(igrou2.eq.4) call lakew (roel) 34 continue c**********************************************************************c c c c irapp that input parameter allows to activate atmospheric c c correction mode c c c c -1: No atmospheric Correction is performed c c 0,1: Atmospheric Correction with Lambertian assumption c c and with the assumption that c c target BRDF is proportional to the input BRDF (see c c case idirec=1) c c c c rapp parameter that contains the reflectance/radiance c c to be corrected. c c c c if rapp >0. : the code retrieve the value of the c c surface reflectance (rog) that will produce a radiance c c equal to rapp [w/m2/str/mic] in the atmospheric c c conditions described by user before c c c c if -1.