subroutine guerin(nopr,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,shadow) ! 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 ! The guerin paper just discusses the Cox Munk probability term ! ****************************** ! See rdguerin.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) real :: fresnel,fresnelsayerval real :: omega,omegaprime real :: nair,nwater real :: shadow ! Used here integer :: i1 real :: nreal,nim real :: muval,mcval,musval,mucval real :: rsp,reflm11,reflm12,reflm33,reflm44 real :: m11g,m12g,m33g,m44g,ratiom12m11,ratiom33m11 ! ****************************** ! Input ! thetasun solar zenith angle (degrees) ! thetaview sensor zenith angle (degrees) ! phisun solar azimuth angle (degrees) ! phiview sensor azimuth angle (degrees) ! phiwind wind azimuth angle (degrees) ! windspeedguerin m/sec ! nreal,nim real and imaginary indices ! Output ! fresnel fresnel reflection coefficient ! rsp sun glint reflection value ! ****************************** ! Helpful for conversions radius=6370.0 ! pi=3.14159265 pi=atan(1.)*4. twopi=2.0*pi circum=2.00*pi*radius ! radians per degree convr=pi/180.0 ! km per degree (110.0) conv=circum/360.0 ! ****************************** ! Define cos and sin terms s1=convr*thetasun snthetasun=sin(s1) csthetasun=cos(s1) s1=convr*thetaview snthetaview=sin(s1) csthetaview=cos(s1) s1=convr*phisun snphisun=sin(s1) csphisun=cos(s1) s1=convr*phiview snphiview=sin(s1) csphiview=cos(s1) s1=convr*(phiview-phiwind) snphiwind=sin(s1) csphiwind=cos(s1) s1=convr*(phisun-phiwind) snsunwind=sin(s1) cssunwind=cos(s1) ! For the specular direction phidiff is near 180.0 degrees s1=convr*phidiff snphidiff=sin(s1) csphidiff=cos(s1) ! ****************************** ! Calculate the Fresnel coefficient (should be near 0.02) ! oops, this gives fresnel values only near 0.02 ! phir=thetasun-thetaview ! s1=convr*phir ! cos2chi=(csthetasun*csthetaview)+(snthetaview*snthetasun*cos(s1)) ! if (cos2chi .gt. 1.0) then ! cos2chi=0.99999999999 ! endif ! if (cos2chi .lt. -1.0) then ! cos2chi=-0.99999999999 ! endif ! coschi=sqrt(0.5*(1+cos2chi)) ! sinchi=sqrt(0.5*(1-cos2chi)) noprf=0 if (nopr .eq. 1) then noprf=1 endif ! hmmm, fresnel does not vary as a function of solar zemith angle (booh) ! call fresnelval(noprf,iout,nreal,nim,coschi,sinchi,fresnel) ! By Sayer Oxford AOPP doc nair=1.00029 nwater=nreal call fresnelsayer(nopr,iout,nair,nwater,csthetasun,snthetasun,& csthetaview,snthetaview,csphidiff,snphidiff,& thetasun,thetaview,phidiff,& omega,omegaprime,fresnelsayerval) fresnel=fresnelsayerval ! ****************************** ! Eequation (3) ! su = dz/dx a1=-((snthetasun*cssunwind)+(snthetaview*csphiwind)) a2=csthetasun+csthetaview su=a1/a2 ! Equation (4) ! sc = dz/dy b1=-((snthetasun*snsunwind)+(snthetaview*snphiwind)) sc=b1/a2 ! ****************************** ! Determine mu and mc du10=u10(2)-u10(1) a1=(windspeedguerin-u10(1))/du10 i1=1+int(a1) if (i1 .lt. 1) then i1=1 endif if (i1 .gt. ndatguerin) then i1=ndatguerin endif muval=mu(i1) mcval=mc(i1) musval=mus(i1) mucval=muc(i1) uval=su/sqrt(muval) cval=sc/sqrt(mcval) u1=uval u2=u1*u1 u3=u2*u1 u4=u3*u1 c1=cval c2=c1*c1 c3=c2*c1 c4=c3*c1 aterm=1.00/(twopi*sqrt((muval*mcval))) a1=0.5*uval*uval e1=exp(-a1) a2=0.5*cval*cval e2=exp(-a2) term1=-(c12(i1)*u1*(c2-1.00))/2.00 term2=-(c30(i1)*(u3-(3.0*u1)))/6.0 term3=(c40(i1)*(u4-(6.0*u2)+3.0))/24.0 term4=(c22(i1)*(u2-1.00)*(c2-1.00))/4.0 term5=(c04(i1)*(c4-(6.0*c2)+3.0))/24.0 ! Gram-Charlier series gc=1.00+term1+term2+term3+term4+term5 ! part of Equation (9) p0=aterm*e1*e2 ! Equation (9) pval=p0*gc ! **************************** if (nopr .eq. 1) then write(iout,100) thetasun,thetaview,phisun,phiview 100 format(/,2x,'guerin: thetasun,thetaview,phisun,phiview',/,& 2x,4(1x,f12.5)) write(iout,110) phiwind,windspeedguerin 110 format(2x,'guerin: phiwind,windspeedguerin',/,& 2x,2(1x,f12.5)) write(iout,120) nreal,nim,fresnel 120 format(2x,'guerin: nreal,nim,fresnel',/,& 2x,3(1x,f12.5)) write(iout,125) muval,mcval,musval,mcval 125 format(2x,'guerin: muval,mcval,musval,mcval',/,& 2x,4(1x,f12.5)) write(iout,137) su,sc 137 format(2x,'guerin: su,sc',/,& 2x,2(1x,f12.5)) write(iout,130) uval,cval 130 format(2x,'guerin: uval,cval',/,& 2x,4(1x,f12.5)) write(iout,135) u1,u2,u3,u4 135 format(2x,'guerin: u1,u2,u3,u4',/,& 2x,4(1x,f12.5)) write(iout,140) c1,c2,c3,c4 140 format(2x,'guerin: c1,c2,c3,c4',/,& 2x,4(1x,f12.5)) write(iout,165) term1,term2,term3,term4,term5 165 format(2x,'guerin: term1,term2,term3,term4,term5',/,& 2x,5(1x,f12.5)) write(iout,170) gc 170 format(2x,'guerin: gc ',/,& 2x,1(1x,f12.5)) write(iout,175) aterm,e1,e2 175 format(2x,'guerin: aterm,e1,e2',/,& 2x,3(1x,f12.5)) write(iout,180) pval 180 format(2x,'guerin: pval=p0*gc',/,& 2x,1(1x,f12.5)) endif ! ****************************** ! Back to Breon's paper (since guerin just discusses the probability term in ! in equation 9 ! Tilt angle ! equation 2 of Breon ! sum= (zx*zx) + (zy*zy) ! is this correct? sum= (su*su) + (sc*sc) tnbeta=sqrt(sum) betarad=atan(tnbeta) betadeg=betarad/convr csbeta=cos(betarad) cs4=csbeta**4.00 ! ****************************** ! equation 4 of Breon top=pi*fresnel*pval denom=4.0*csthetasun*csthetaview*cs4 rsp=top/denom ! ****************************** ! See equation 87 and 88 of Mischchenko and Travis ! for the shadow term (seems to be 1.00) call mtshadow(nopr,iout,& thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,& sig2,shadow) ! include the shadow term rsp=rsp*shadow ! ****************************** ! Gave result similar to gisscoxmunk when you calculated ratiom12m11 ! used below call hafermanmatrix(nopr,iout,& thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,& m11g,m12g,m33g,m44g,ratiom12m11,ratiom33m11) reflm11=rsp reflm12=rsp*ratiom12m11 reflm33=rsp*ratiom33m11 reflm44=rsp*ratiom33m11 ! ****************************** if (nopr .eq. 1) then write(iout,185) top,denom 185 format(/,2x,'guerin: top,denom',/,& 2x,1(1x,f12.5)) write(iout,190) rsp 190 format(2x,'guerin: rsp',/,& 2x,1(1x,f12.5)) write(iout,192) m11g,m12g,m33g,m44g 192 format(2x,'guerin: m11g,m12g,m33g,m44g',/,& 2x,4(1x,f12.5)) write(iout,194) ratiom12m11,ratiom33m11 194 format(2x,'guerin: ratiom12m11,ratiom33m11',/,& 2x,6(1x,f12.5)) write(iout,195) reflm11,reflm12,reflm33,reflm44 195 format(2x,'guerin: reflm11,reflm12,reflm33,reflm44',/,& 2x,4(1x,f12.5)) endif ! ****************************** return end