program main ! ***** ! Run gisscoxmunk_vfunction_plus.f90, breon.f90, and guerin.f90 side by side ! Results are placed in the output f.out ascii file ! Specify the following below by hand ! winddpeedbreon,wmic,nreal,nim,phisun,phiview,windaz ! Have NSTOKES=4 below ! For ione = 1 will write out from the .f90 subroutines for one latitude ! For ione = 0 will loop over num=16 latitudes and write out summary data ! ***** integer, parameter :: numlatmax=100 integer, parameter :: ndatguerin=25 integer :: idat,nopr,iout,num,i,j integer :: ione,nnval ! Write flags integer :: noprzen,noprgiss,noprnew,noprbreon ! character(len=5) :: fout ! real :: wave(nwavemax) ! real :: rn(nwavemax),ri(nwavemax) ! see rdzen.f90 for Lite file averages integer :: numqf0(numlatmax),numqf1(numlatmax) real :: lons(numlatmax),lats(numlatmax) real :: zensun(numlatmax),zensensor(numlatmax) real :: azsun(numlatmax),azsensor(numlatmax) real :: phidiffs(numlatmax) ! See rdguerin.f90 and guerin.f90 real :: u10(ndatguerin) real :: mu(ndatguerin),mc(ndatguerin),mus(ndatguerin),muc(ndatguerin) real :: c40(ndatguerin),c04(ndatguerin),c22(ndatguerin),c12(ndatguerin) real :: c30(ndatguerin),thetapmax(ndatguerin) ! Added here DOUBLE PRECISION suntheta,sensorview DOUBLE PRECISION windspeed,SIGMA2,CN2,lambertian DOUBLE PRECISION SHADOW double precision wmic,wavenumber double precision mre,mim double precision m11,m12,m21,m22 double precision m11div,m12div,m21div,m22div double precision ratiomtm21m11 double precision exvalsign,term4,term5 double precision factor62,dmodmt,factor62mult double precision r62mtm11,r62mtm12 double precision theta86,shadow87 ! ****************************** ! Subroutine arguments for gisscpxmunk_vfunction_plus.f90 INTEGER, PARAMETER :: NSTOKES=4, NPARS=4 double precision PARS(NPARS) double precision XI, SXI, XJ, SXJ double precision CKPHI_REF, SKPHI_REF double precision R1(NSTOKES) double precision Ls_R1(NSTOKES,NPARS) ! useful for conversions double precision radius,pi,circum,convr,conv ! For breon.f90 and guerin.f90 real :: nreal,nim real :: m11h,m12h,m33h,m44h real :: reflm11,reflm12,reflm33,reflm44 ! Fro recalcuation real :: a1,a2,a3,a4,a5,a6 real :: b1,b2,b3,b4,b5,b6 real :: c1,c2,c3,c4,c5,c6 ! For gisscoxmunk_vfunction_plus.f90 real :: rspgiss(numlatmax) real :: sunglint11giss(numlatmax),sunglint12giss(numlatmax) real :: ratiom12m11giss(numlatmax) real :: shadowgiss(numlatmax),shadowbreon(numlatmax) ! For breon.f90 real :: rspbreon(numlatmax),rspdivbreon(numlatmax) real :: gcbreon(numlatmax),exvalbreon(numlatmax),cs4breon(numlatmax) real :: betadegbreon(numlatmax),fresnelbreon(numlatmax) real :: ratiom12m11breon(numlatmax) real :: m11breon(numlatmax),m12breon(numlatmax) real :: sunglint11breon(numlatmax),sunglint12breon(numlatmax) real :: sunglint33breon(numlatmax),sunglint44breon(numlatmax) ! For coxmunk_new.f90 real :: rspnew(numlatmax),m11new(numlatmax),m12new(numlatmax) real :: m11divnew(numlatmax),m12divnew(numlatmax) real :: shadownew(numlatmax),exvalnew(numlatmax) real :: term4new(numlatmax),term5new(numlatmax) real :: factor62new(numlatmax),factor62multnew(numlatmax) real :: comparenew(numlatmax),comparebreon(numlatmax) real :: theta86new(numlatmax),ratiom12m11mt(numlatmax) real :: sunglint11new(numlatmax),sunglint12new(numlatmax) ! For guerin.f90 real :: rspguerin(numlatmax) real :: m11guerin(numlatmax),m12guerin(numlatmax) real :: sunglint11guerin(numlatmax),sunglint12guerin(numlatmax) real :: sunglint33guerin(numlatmax),sunglint44guerin(numlatmax) real :: shadowguerin(numlatmax) real :: m11g,m12g,m33g,m44g ! ****************************** ! to run one of the latitudes put ione = 1 and specify nnval ! nnval=1 for agreement (sza=17, sensor zenith=13) ! nnval=12 for disagreement (sza=57, sensor zenith=44) ! nnval=15 for disagreement (sza=71, sensor zenith=54) ione=0 nnval=1 ! Set print flags if (ione .eq. 1) then noprzen=1 noprgiss=1 noprnew=1 noprbreon=1 noprguerin=1 endif if (ione .ne. 1) then noprzen=0 noprgiss=0 noprnew=0 noprbreon=0 noprguerin=0 endif ! ***** ! Output ascii file f.out iout=6 open(iout,form='formatted',file='f.out',status='unknown') write(iout,fmt=10) 10 format(/,2x,"main: run brdf routines side by side") ! ****************************** ! Helpful for conversions radius=6370.0d0 ! pi=3.14159265d0 pi=atan(1.)*4. twopi=2.0d0*pi circum=2.00*pi*radius convr=pi/180.0d0 ! km per degree (110.0) conv=circum/360.0d0 ! ****************************** ! Read in the Lite file values ! There are num latitudes ! with longitude,latitude,solar zenith, sensor ! zenith, solar azimuth, sensor azimuth ! equal to lons(i),lats(i),zensun(i),zensensor(i() ! azsun(i),azsensor(i) for i=1,num call rdzen(noprzen,iout,& num,numqf0,numqf1,lons,lats,& zensun,zensensor,azsun,azsensor,& phidiffs) ! ****************************** ! Read in the guerin coefficients (Table 1 values) ! See Guerin, Capelle, Jean-Michael Hartmann ! Revisiting the Cox and Munk wave-slope statistics using IASI observations ! of the sea surface ! October 11, 2022 preprint https://arxiv.org/abs/2210.05456 ! These coefficients are used in guerin.f90 call rdguerin(noprguerin,iout,& u10,mu,mc,mus,muc,c40,c04,c22,c12,c30,thetapmax) ! ****************************** ! User: specify the windspeed in m/sec ! windspeedbreon=15.0 windspeedbreon=8.0 ! windspeedbreon=4.0 windspeedguerin=windspeedbreon if (windspeedguerin .lt. 3.0) then windspeedguerin=3.0 endif if (windspeedguerin .gt. 15.0) then windspeedguerin=15.0 endif ! For gisscoxmunk_vfunction_plus windspeed=1.0d0*windspeedbreon ! ****************************** ! User: specify the wavelength, and the nreal and ! nim refractive indices ! The wavelength in microns wmic=0.76d0 ! have in units near 1 for convenicnce wavenumber=twopi/wmic ! The complex index of refraction nreal=1.33 nim=-1.3e-7 mre=1.0d0*nreal mim=1.0d0*nim write(iout,fmt=15) wmic,wavenumber,nreal,nim,windspeed 15 format(/,2x,"main: wmic,wavenumber (2pi/wmic)",/,& 2x,1(1x,f10.4),1p,1(1x,e10.3),/,& 2x,"main: nreal,nim",/,& 2x,1p,2(1x,e10.3),/,& 2x,"main: windspeed (m/sec)",/,& 2x,0p,1(1x,f10.4)) ! ****************************** ! Loop over the latitudes nn1=1 nn2=num if (ione .eq. 1) then nn1=nnval nn2=nn1 endif do nn=nn1,nn2 ! **************** ! Will specify the inputs to the gisscoxmunk_vfunction.f90, ! breon.f90, and gruin.f90 routines ! *************** ! For gisscoxmunk_vfunction_plus.f90 ! Just use the real part ! CN2=CMPLX(mre,mim) CN2=mre ! Equation 18 0f Mishchenko and Travis (1997) ! See Mishchenko and Travis, "Satellite retrieval ! of aerosol properties over the ocean using polarization as ! well as intensity sunlight" ! JGR, vol 102, No. D14, pages 16989-17013, July 27, 1997 ! Refer to this paper as MT below SIGMA2=0.5d0*(0.003d0+(0.00512d0*windspeed)) ! Don't add lambertian to the surface reflectance lambertian=0.0d0 ! Used by gisscoxmunk_vfunction_plus.f90 PARS(1)=SIGMA2 PARS(2)=CN2 PARS(3)=lambertian ! set to 1.d0 to include shadowing PARS(4)=1.d0 suntheta=1.0d0*zensun(nn) sensorview=1.0d0*zensensor(nn) ! Note use of abs XI=ABS(cos(convr*suntheta)) XJ=ABS(cos(convr*sensorview)) SXI=ABS(sin(convr*suntheta)) SXJ=ABS(sin(convr*sensorview)) ! Careful ! routine has incident azimuth as zero ! For glint then (check) for reflected azimuth be 180 (and ! cos(180)=-1.0d0) ! this gave 0.02 (for sza=17, view=13 case) ! CKPHI_REF=-1.0d0 ! SKPHI_REF=0.0d0 ! this gave 0.127 (for sza=17, view=13 case) CKPHI_REF=1.0d0 SKPHI_REF=0.0d0 ! Note that incident solar azimuth in Mishchenko and Travis ! is 0.0 (with cos=1 and sin=0) ! The reflected ray is in the upper half of the x-y-z grid ! while the incident solar beam vector is in the lower half of ! the x-y-z grid ! See gisscoxmunk_vfunction_plus.f90 ! CKPHI_INC = 1.d0 ! SKPHI_INC = 0.d0 ! ********** ! For breon.f90 ! Input ! thetasun solar zenith angle (degrees) ! thetaview sensor zenith angle (degrees) ! phisun solar azimuth angle (degrees) ! phiview ssensor azimuth angle (degrees) ! windaz wind azimuth angle ! windspeedbreon m/sec ! nreal,nim real and imaginary indices ! phidiff relative azimuth angle between solar and sensor azimuths, ! for specular reflection phidiff is near 180 degrees ! windazdiff difference in wind and solar azimuths, windaz+180.0-phisun thetasun=zensun(nn) thetaview=zensensor(nn) ! obs point is at x=y=0 ! User: specify phisun in degrees ! sun is along the -x axis (to the left of the origin) with azimuth angle 180.0 phisun=180.0 ! User: specify phiview in degrees ! The sensor is along the +x axis (to the right of origin) at azimuth angle 0.0 ! phisun=180 and phiview=0 corresponds to angles near specular reflection ! (Note that phiview=0.0 ! hmmm, shouldn't phiview be offset by 30 degrees? i30=0 if (i30 .eq. 1) then phiview=30.0 endif ! User: specify windaz in degrees ! The windaz=0 for wind going from west to east along the y=0, x axis ! The windaz=180 for wind going from east to west along the y=0, x axis windaz=0.0 ! windaz=45.0 ! windaz=90.0 ! windaz=135.0 ! windaz=180.0 ! For guerin.f90 ! You have a separate variable for guerin.f90 ! in case there is a different corrdinate convention for guerin.f90 phiwind=windaz ! Real and complex indices nreal=mre nim=mim ! Carefull ! For breon.f90 phidiff should be near 180.0 ! for angles near specular reflection phidiff=phisun-phiview ! Carefull ! From Breon email phi(wind)=AZW(downwind)+180.0-AZSUN windazdiff=windaz+180.0-phisun if (windazdiff .gt. 360.0) then windazdiff=windazdiff-360.0 endif if (windazdiff .lt. 0.0) then windazdiff=windazdiff+360.0 endif ! ****************************** if (nn .eq. nn1) then write(iout,fmt=20) mre,mim,suntheta,sensorview,& phisun,phiview,windaz,windazdiff 20 format(/,2x,"main: mre,mim ",2x,1p,2(1x,d10.3),/,& 2x,"main:(for nn=nn1) suntheta,sensorview ",2x,1p,2(1x,d10.3),/,& 2x,"main: for Breon geometry, azimuth of sun and sensor",/,& 2x,"main: phisun ",0p,f10.4,/,& 2x,"main: phiview ",f10.4,/,& 2x,"main: windaz ",f10.4,/,& 2x,"main: windazdiff (windaz+180-phisun) ",f10.4) endif ! User spedcifications are done ! ****************************** ! Note that a print flag and SAHDOW is added to the original routine ! argument list call gisscoxmunk_vfunction_plus & (NSTOKES, NPARS, PARS,& !I XJ, SXJ, XI, SXI,& !I CKPHI_REF, SKPHI_REF,& !I R1, Ls_R1, noprgiss,iout,SHADOW) ! Store results for summary write rspgiss(nn)=r1(1) ! The final sun glint reflectance matrix germs sunglint11giss(nn)=r1(1) sunglint12giss(nn)=r1(2) ratiom12m11giss(nn)=r1(2)/r1(1) shadowgiss(nn)=1.0*SHADOW ! ****************************** ! Purpose of coxmunk_new.f90 is to calculate ! the reflection matrix ratios from scratch ! using Appendix A of MT (see getmtmatrix.f90 called ! in coxmunk_new.f90) ! The ratio of m12/m11 is ratiomtm21m11 ! Other values from gisscoxmunk_vfunction_plus can be compared ! to terms calculated from breon.f90 ! This helps to see why breon.f90 and guerin.f90 sun glint values ! differs from gisscoxmunk_vfundion_plus.f90 call coxmunk_new & (NSTOKES, NPARS, PARS,& !I XJ, SXJ, XI, SXI,& !I CKPHI_REF, SKPHI_REF,& !I R1, Ls_R1,& !I noprnew,iout,& !I mre,mim,& !I wmic,wavenumber,shadow87,comparen,exvalsign,term4,term5,& factor62,dmodmt,factor62mult,& m11,m12,m21,m22,& m11div,m12div,m21div,m22div,& theta86,r62mtm11,r62mtm12,ratiomtm21m11) ! Store results for summary write rspnew(nn)=R1(1) ! equations 64 and 65 of MT (the M11 and M12 terms) m11new(nn)=real(m11) m12new(nn)=real(m12) ! m11divnew is m11/dmodmt ! so you can caompre to m11breon m11divnew(nn)=real(m11div) m12divnew(nn)=real(m12div) ! from mtshadow.f90 (equation 87 of MT) shadownew(nn)=real(shadow87) comparenew(nn)=real(comparen) ! the exp(-xx)/sigma2 term exvalnew(nn)=real(exvalsign) term4new(nn)=real(term4) term5new(nn)=real(term5) factor62new(nn)=real(factor62) factor62multnew(nn)=real(factor62mult) ! equation 86 theta angle theta86new(nn)=real(theta86) ! since m12=m21 we can use Mishchenko and Traivs ratiomtm21m11 for ratiom1211mt ratiom12m11mt(nn)=real(ratiomtm21m11) ! The final sun glint reflectance matrix terms sunglint11new(nn)=real(r62mtm11) sunglint12new(nn)=real(r62mtm12) ! ****************************** ! Breon brdf, using Breon and Henriot (2006) paper equations ! Breon does not specify the reflection matrix ! Will use getmtmatrix.f90 (that calculates Appendix A values ! from the MT paper) to specify the m11 and m12 reflection matix values) call breon(noprbreon,iout,& thetasun,thetaview,phisun,phiview,windaz,windspeedbreon,nreal,nim,fresnel,& phidiff,windazdiff,rsp,rspdiv,shadowb,& m11h,m12h,m33h,m44h,& reflm11,reflm12,reflm33,reflm44,& reflm11b,reflm12b,reflm33b,reflm44b,& compareb,gc,exvalsig,cs4,betadeg,ratiom12m11) ! Store results for summary write rspbreon(nn)=rsp ! rsp divided by fresnel rspdivbreon(nn)=rspdiv ! from hafermanmatrix.f90 m11breon(nn)=m11h m12breon(nn)=m12h ! from mtshadow.f90 shadowbreon(nn)=shadowb ! from breon.f90 gcbreon(nn)=gc comparebreon(nn)=compareb exvalbreon(nn)=exvalsig cs4breon(nn)=cs4 betadegbreon(nn)=betadeg fresnelbreon(nn)=fresnel ratiom12m11breon(nn)=ratiom12m11 ! from breon.f90 the final sun glint reflection matrix terms sunglint11breon(nn)=reflm11 sunglint12breon(nn)=reflm12 sunglint33breon(nn)=reflm33 sunglint44breon(nn)=reflm44 ! ****************************** ! Guerin brdf ! THis uses the guerin etal 2022 preprint Table 1 in conjunction with ! Breon equations (since Guerin only discusses the Cox-Munk probability term) ! The guerin paper just discusses the Cox Munk probability term sig2=real(SIGMA2) call guerin(noprguerin,iout,& ndatguerin,u10,mu,mc,mus,muc,c40,c04,c22,c12,c30,thetapmax,& thetasun,thetaview,phisun,phiview,phiwind,windspeedguerin,nreal,nim,fresnel,& phidiff,& rsp,m11g,m12g,m33g,m44g,reflm11,reflm12,reflm33,reflm44,sig2,shadowg) ! Store results for summary write rspguerin(nn)=rsp ! from mtshadow.f90 shadowguerin(nn)=shadowg ! Based on hfermanmatrix.f90 toutine m11guerin(nn)=m11g m12guerin(nn)=m12g ! The final sunglint reflectance matrix terms sunglint11guerin(nn)=reflm11 sunglint12guerin(nn)=reflm12 enddo ! Loop over latitudes is done ! ****************************** ! Write out the results of the calculations write(iout,fmt=100) 100 format(/,2x,"main: the final (1,1) reflectance term",/,& 2x,"main: rspgiss from coxmunk_vfunction_plus.f90 ",/,& 2x,"main: rspnew from coxmunk_new.f90 ",/,& 2x,"main: rspbreon from breon.f90 ",/,& 2x,"main: rspguerin from guerin.f90 ",/,& 2x,"main: nn,lats,zensun,rspgiss,rspnew,rspbreon,rspguerin") do nn=nn1,nn2 write(iout,fmt=105) nn,lats(nn),zensun(nn),rspgiss(nn),rspnew(nn),rspbreon(nn),& rspguerin(nn) enddo 105 format(2x,i3,6(1x,f12.5)) write(iout,fmt=110) 110 format(/,2x,"main: the M matrix (1,1) term ",/,& 2x,"main: m11divnew is m11/dmodmt from coxmunk_new.f90 ",/,& 2x,"main: m11breon is from breon.f90 (uses hafermanmatrix.f90) ",/,& 2x,"main: m11guerin is from guerin.f90 ",/,& 2x,"main: nn,lats,zensun,m11divnew,m11breon,m11guerin") do nn=nn1,nn2 write(iout,fmt=115) nn,lats(nn),zensun(nn),m11divnew(nn),m11breon(nn),& m11guerin(nn) enddo 115 format(2x,i3,5(1x,f12.5)) write(iout,fmt=120) 120 format(/,2x,"main: the M matrix (1,2) term ",/,& 2x,"main: m12divnew is m11/dmodmt from coxmunk_new.f90 ",/,& 2x,"main: m12breon is from breon.f90 (uses hafermanmatrix.f90) ",/,& 2x,"main: m12guerin is from guerin.f90 ",/,& 2x,"main: nn,lats,zensun,m12divnew,m12breon,m12guerin") do nn=nn1,nn2 write(iout,fmt=125) nn,lats(nn),zensun(nn),m12divnew(nn),m12breon(nn),& m12guerin(nn) enddo 125 format(2x,i3,5(1x,f12.5)) write(iout,fmt=130) 130 format(/,2x,"main: the final sun glint reflectance matrix (1,1) term ",/,& 2x,"main: sunglint11giss is from gisscoxmunk_vfunction_plus.f90 ",/,& 2x,"main: sunglint11new is from coxmunk_new.f90 ",/,& 2x,"main: sunglint11breon is from breon.f90 (uses hafermanmatrix.f90) ",/,& 2x,"main: suunglint11guerin is from guerin.f90 ",/,& 2x,"main: nn,lats,zensun, singlint11new,sunglint11breon,sunglint11guerin") do nn=nn1,nn2 write(iout,fmt=135) nn,lats(nn),zensun(nn),sunglint11giss(nn),sunglint11new(nn),& sunglint11breon(nn),sunglint11guerin(nn) enddo 135 format(2x,i3,6(1x,f12.5)) write(iout,fmt=140) 140 format(/,2x,"main: the final sun glint reflectance matrix (1,2) term ",/,& 2x,"main: sunglint12giss is from gisscoxmunk_vfunction_plus.f90 ",/,& 2x,"main: sunglint12new is from coxmunk_new.f90 ",/,& 2x,"main: sunglint12breon is from breon.f90 (uses hafermanmatrix.f90) ",/,& 2x,"main: suunglint12guerin is from guerin.f90 ",/,& 2x,"main: nn,lats,zensun, singlint12new,sunglint12breon,sunglint12guerin") do nn=nn1,nn2 write(iout,fmt=145) nn,lats(nn),zensun(nn),sunglint12giss(nn),sunglint12new(nn),& sunglint12breon(nn),sunglint12guerin(nn) enddo 145 format(2x,i3,6(1x,f12.5)) write(iout,fmt=150) 150 format(/,2x,"main: nn,lats,zensun,shadowgiss,shadowbreon,shadowguerin") do nn=nn1,nn2 write(iout,fmt=155) nn,lats(nn),zensun(nn),shadowgiss(nn),shadowbreon(nn),& shadowguerin(nn) enddo 155 format(2x,i3,5(1x,f10.4)) write(iout,fmt=160) 160 format(/,2x,"main: fresnelbreon is the fresnel value",/,& 2x,"main: gcbreon is the breon Gram-Charlier series sum",/,& 2x,"main: comparebreon and comparenew should be similar",/,& 2x,"main: nn,lats,zensun,gcbreon,comparebreon,comparenew") do nn=nn1,nn2 write(iout,fmt=165) nn,lats(nn),zensun(nn),fresnelbreon(nn),gcbreon(nn),& comparebreon(nn),comparenew(nn) enddo 165 format(2x,i3,6(1x,f10.4)) write(iout,fmt=170) 170 format(/,2x,"main: exvalbreon and exvalnew are exp(-xx)/sigma2 values",/,& 2x,"main: theta86new is the eqn 86 theta from getmtmatrix",/,& 2x,"main: betadegbreon is the surface slope in degrees",/,& 2x,"main: cs4breon is the cos(suface slope)**4 term",/,& 2x,"main: term4new is RDZ4/(FACT1*FACT1) value",/,& 2x,"main: term5new is dmodmt value",/,& 2x,"main: nn,lats,zensun,exvalbreon,exvalnew,theta86new,betadegbreon,cs4breon",/,& 2x,"main: term4new,term5new") do nn=nn1,nn2 write(iout,fmt=175) nn,lats(nn),zensun(nn),exvalbreon(nn),exvalnew(nn),& theta86new(nn),betadegbreon(nn),cs4breon(nn),term4new(nn),term5new(nn) enddo 175 format(2x,i3,9(1x,f10.4)) write(iout,fmt=180) 180 format(/,2x,"main: compare refl(1,2) / refl(1,1) ratios",/,& 2x,"main: ratiom12m11giss is from gisscoxmunk_vfuncrtion_plus.f90",/,& 2x,"main: ratio of R1(2)/R1(1) ",/,& 2x,"main: ratiom12m11breon is from breon.f90 ((which uses",/,& 2x,"main: the hafermanmatrix.f90 routine)",/,& 2x,"main: ratiom12m11mt is from Mishchenko and Travis MT paper",/,& 2x,"main: nn,lats,zensun,ratiom12m11giss,ratiom12m11breon,ratiom12m11mt") do nn=nn1,nn2 write(iout,fmt=185) nn,lats(nn),zensun(nn),ratiom12m11giss(nn),ratiom12m11breon(nn),ratiom12m11mt(nn) enddo 185 format(2x,i3,5(1x,f10.4)) write(iout,fmt=190) 190 format(/,2x,"main: factor62mult is MT equation 62 scalar part * dmodmt",/,& 2x,"main: rspdivbreon is rsp / fresnel from breon.f90 ",/,& 2x,"main: ratio62=factor62multnew(nn)/rspdivbreon(nn) ",/,& 2x,"main: for perfect agreement ratio62 should be 1.00",/,& 2x,"main: nn,lats,zensun,factor62multnew,rspdivbreon") do nn=nn1,nn2 ratio62=factor62multnew(nn)/rspdivbreon(nn) write(iout,fmt=195) nn,lats(nn),zensun(nn),factor62multnew(nn),& rspdivbreon(nn),ratio62 enddo 195 format(2x,i3,5(1x,f10.4)) write(iout,fmt=200) windspeedbreon,windaz 200 format(/,2x,"main: recalculate from scalar and matrix terms",/,& 2x,"main: the fianl sun glint reflection matrix (1,1) terms",/,& 2x,"main: from coxmunk_new, breon.f90, guerin.f90",/,& 2x,"main: windspeedbreon,windaz ",/,2(1x,f10.4),/,& 2x,"main: a1=factor62new,b1=m11new,c1=a1*b1 ",/,& 2x,"main: a2=factor62multnew,b2=m11divnew,c2=a2*b2 ",/,& 2x,"main: a3=rspdivbreon,b3=m11breon,c3=a3*b3 ",/,& 2x,"main: c4=rspguerin",/,& 2x,"main: nn,lats,zensun, c1,c2,c3,c4 ") do nn=nn1,nn2 a1=factor62new(nn) b1=m11new(nn) c1=a1*b1 a2=factor62multnew(nn) b2=m11divnew(nn) c2=a2*b2 a3=rspdivbreon(nn) b3=m11breon(nn) c3=a3*b3 c4=rspguerin(nn) write(iout,fmt=205) nn,lats(nn),zensun(nn),c1,c2,c3,c4 enddo 205 format(2x,i3,11(1x,f10.4)) ! ****************************** ! Close the f.out ascii file close (iout) ! ***** stop end