!
!-------------------------------------------------------------------------------------------------------------
!M+
! NAME:
! CRTM_FASTEM1
!
! PURPOSE:
! This module computes ocean emissivity and its jacobian over water. The code is adopted
! from RTTOV Fastem version 1.
!
! Method:
! FASTEM-1 English and Hewison 1998.
! http://www.metoffice.com/research/interproj/nwpsaf/rtm/evalfastems.pdf
!
! History:
!
! CATEGORY:
! CRTM : Surface : MW OPEN OCEAN EM
!
! LANGUAGE:
! Fortran-95
!
! CREATION HISTORY:
!
! 1998-02-20 treadon - gather all emissivity calculations from
! setuprad and move into this subroutine
! 2004-07-23 weng,yan,okamoto - incorporate MW land and snow/ice emissivity
! models for AMSU-A/B and SSM/I
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-11-22 derber - modify for openMP
! 2005-01-20 okamoto- add preprocessing for ocean MW emissivity jacobian
! 2005-03-05 derber- add adjoint of surface emissivity to this routine
! 2005-07-15 derber- include mhs
! 2005-09-28 derber - modify to handle percent surface type and mixed conditions
!
! 2005-12-15 Modified for CRTM
! by: Quanhua Liu, QSS Group Inc., Quanhua.Liu@noaa.gov
! Yong Han, NOAA/NESDIS; Yong.Han@noaa.gov
! Paul van Delst, CIMSS/SSEC; paul.vandelst@ssec.wisc.edu
!
!------------------------------------------------------------------------------------------------------------
MODULE CRTM_Fastem1
! ----------
! Module use
! ----------
USE Type_Kinds
! -- CRTM modules
USE CRTM_Parameters
! -----------------------
! Disable implicit typing
! -----------------------
IMPLICIT NONE
PRIVATE
PUBLIC :: Fastem1
! Explanation for Fastem1_Coef :
! emc(59): Emissivity model data
! Permittivity model data (Lamkaouchi model)
! [1-3]: Temperature polynomial coefficients for Tau1 - Lamkaouchi (1996)
! [4-7]: Temperature polynomial coefficients for Tau2 - Lamkaouchi (1996)
! [8-11]: Temperature polynomial coefficients for Del1 - Lamkaouchi (1996)
! [12-15]: Temperature polynomial coefficients for Del2 - Lamkaouchi (1996)
! [16-17]: Temperature polynomial coefficients for static permittivity - Lamkaouchi (1996)
! [18-19]: Temperature polynomial coefficients for infinite freq. permittivity - Lamkaouchi (1996)
! Pi is stored for good measure
! [20]: Stored value of Pi
! Bragg scattering correction coefficients
! [21]: Scaling factor for small scale correction - see English (1997)
! Foam model coefficients for Monahan model
! [22]: First coefficient in Monahan foam model (neutral stability) - see English (1997)
! [23]: Second coefficient in Monahan foam model (neutral stability) - see English (1997)
! Alternative permittivity model (Liebe)
! [30]: a1 in Liebe's dielectric model - see Liebe (1989)
! [31]: b1 in Liebe's dielectric model - see Liebe (1989)
! [32]: c1 in Liebe's dielectric model - see Liebe (1989)
! [33]: c2 in Liebe's dielectric model - see Liebe (1989)
! [34]: d1 in Liebe's dielectric model - see Liebe (1989)
! [35]: d2 in Liebe's dielectric model - see Liebe (1989)
! [36]: d3 in Liebe's dielectric model - see Liebe (1989)
! [37]: e1 in Liebe's dielectric model - see Liebe (1989)
! [38]: e2 in Liebe's dielectric model - see Liebe (1989)
! Version 2 of large scale correction which *DOES��* take account of
! hemispherical scattering.
! 1.) Vertical polarisation mode
! [24]: Term a00 in vertical pol of large scale correction model
! [25]: Term a01 in vertical pol mode of large scale correction model
! [26]: Term a02 in vertical pol mode of large scale correction model
! [27]: Term a10 in vertical pol mode of large scale correction model
! [28]: Term a11 in vertical pol mode of large scale correction model
! [29]: Term a12 in vertical pol mode of large scale correction model
! [30]: Term a20 in vertical pol mode of large scale correction model
! [31]: Term a21 in vertical pol mode of large scale correction model
! [32]: Term a22 in vertical pol mode of large scale correction model
! [33]: Term a30 in vertical pol mode of large scale correction model
! [34]: Term a31 in vertical pol mode of large scale correction model
! [35]: Term a32 in vertical pol mode of large scale correction model
! [36]: Term a40 in vertical pol mode of large scale correction model
! [37]: Term a41 in vertical pol mode of large scale correction model
! [38]: Term a42 in vertical pol mode of large scale correction model
! [39]: Term a50 in vertical pol mode of large scale correction model
! [40]: Term a51 in vertical pol mode of large scale correction model
! [41]: Term a52 in vertical pol mode of large scale correction model
! 2. ) Horizontal polarisation mode
! [42]: Term a00 in horizontal pol mode of large scale correction model
! [43]: Term a01 in horizontal pol mode of large scale correction model
! [44]: Term a02 in horizontal pol mode of large scale correction model
! [45]: Term a10 in horizontal pol mode of large scale correction model
! [46]: Term a11 in horizontal pol mode of large scale correction model
! [47]: Term a12 in horizontal pol mode of large scale correction model
! [48]: Term a20 in horizontal pol mode of large scale correction model
! [49]: Term a21 in horizontal pol mode of large scale correction model
! [50]: Term a22 in horizontal pol mode of large scale correction model
! [51]: Term a30 in horizontal pol mode of large scale correction model
! [52]: Term a31 in horizontal pol mode of large scale correction model
! [53]: Term a32 in horizontal pol mode of large scale correction model
! [54]: Term a40 in horizontal pol mode of large scale correction model
! [55]: Term a41 in horizontal pol mode of large scale correction model
! [56]: Term a42 in horizontal pol mode of large scale correction model
! [57]: Term a50 in horizontal pol mode of large scale correction model
! [58]: Term a51 in horizontal pol mode of large scale correction model
! [59]: Term a52 in horizontal pol mode of large scale correction model
REAL( fp ), PARAMETER :: emc(59) = Reshape( &
& (/0.175350E+02_fp, -.617670E+00_fp, .894800E-02_fp, .318420E+01_fp,&
0.191890E-01_fp, -.108730E-01_fp, .258180E-03_fp, .683960E+02_fp,&
-.406430E+00_fp, .228320E-01_fp, -.530610E-03_fp, .476290E+01_fp,&
0.154100E+00_fp, -.337170E-01_fp, .844280E-03_fp, .782870E+02_fp,&
-.434630E-02_fp, .531250E+01_fp, -.114770E-01_fp, .314159E+01_fp,&
-.100000E+01_fp, .195000E-04_fp, .255000E+01_fp, -.637182E+01_fp,&
0.253918E-01_fp, .357569E-04_fp, .942928E+01_fp, -.332839E-01_fp,&
-.647724E-04_fp, -.329282E+01_fp, .965450E-02_fp, .281588E-04_fp,&
0.252676E+00_fp, .343867E-02_fp, -.156362E-04_fp, -.156669E-03_fp,&
0.139485E-04_fp, -.407633E-07_fp, -.141316E+00_fp, -.356556E-02_fp,&
0.142869E-04_fp, -.240701E+01_fp, -.563888E-01_fp, .325227E-03_fp,&
0.296005E+01_fp, .704675E-01_fp, -.426440E-03_fp, -.751252E+00_fp,&
-.191934E-01_fp, .125937E-03_fp, -.288253E+00_fp, -.102655E-02_fp,&
0.226701E-05_fp, -.119072E-02_fp, -.263165E-04_fp, .114597E-06_fp,&
0.406300E+00_fp, .200031E-02_fp, -.781635E-05_fp/), (/59/) )
CONTAINS
!-------------------------------------------------------------------------------------------------------------
!M+
! NAME:
! FASTEM1
!
! PURPOSE:
! Subroutine to compute ocean emissivity and its derivative to wind speed. This code is adopted from
! RTTOV Fastem version 1.
!
! CATEGORY:
! CRTM : Surface : MW OPEN OCEAN EM
!
! LANGUAGE:
! Fortran-95
!
! CALLING SEQUENCE:
! CALL NESDIS_OCeanEM
!
! INPUT ARGUMENTS:
!
! Frequency Frequency User defines
! This is the "I" dimension
! UNITS: GHz
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
!
! Sat_Zenith_Angle The angle values in degree
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: Degrees
! TYPE: REAL( fp )
! DIMENSION: Rank-1, (I)
!
! SST Ocean surface temperature
! UNITS: Kelvin, K
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
! Wind_Speed Ocean surface wind speed
! UNITS: m/s
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
!
! OUTPUT ARGUMENTS:
!
! Emissivity: The surface emissivity at vertical and horizontal polarizations.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: N/A
! TYPE: REAL( fp )
! DIMENSION: ONE
!
! dEH_dWindSpeed: The surface horizontally polarized emissivity derivative to wind speed.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: N/A
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
! dEV_dWindSpeed: The surface vertically polarized emissivity derivative to wind speed.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: N/A
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
!
! CALLS:
! None
!
! SIDE EFFECTS:
! None.
!
! RESTRICTIONS:
! None.
!
!M-
!------------------------------------------------------------------------------------------------------------
SUBROUTINE Fastem1(Frequency, & ! INPUT
Sat_Zenith_Angle, & ! INPUT
SST, & ! INPUT
Wind_Speed, & ! INPUT
Emissivity, & ! OUTPUT
dEH_dWindSpeed, & ! OUTPUT)
dEV_dWindSpeed) ! OUTPUT)
! ---------------------------------------------------------------------------------------------------
!
REAL( fp ), INTENT( IN ) :: Frequency, Sat_Zenith_Angle
REAL( fp ), INTENT( IN ) :: SST, Wind_Speed
REAL( fp ), INTENT( IN OUT ) :: Emissivity(:), dEH_dWindSpeed, dEV_dWindSpeed
!
! Declare passed variables.
! Declare local variables
INTEGER :: kcho,n,kch,nn,nnp,i
INTEGER :: error_status
real(fp) zch4,xcorr2v,evertr,ehorzr,xcorr2h,ffoam,zcv2,zcv3
real(fp) xcorr1,zcv1,zcv4,zch1,zch2,zcv5,zcv6,tau2,degre
real(fp) wind,ehorz,evert,sec,sec2,freqghz2,dtde
real(fp) u10mps2,usec,tccub,tau1,tc,tcsq,term2,freqghz
real(fp) term1,u10mps,ps2,pc2,pcc,pss,rvertsi,rverts,rvertsr
real(fp) rverts5,rhorzs5,xcorr15,ffoam5,evertr5,ehorzr5
real(fp) perm_real,perm_imag,rhorzsr,zch5,zch6,zch3,rhorzsi
real(fp) rhorzs,perm_imag2,einf,fen,del2,del1,fen2,perm_real2
real(fp) perm_imag1,perm_real1,den1,den2,emisst
real(fp) ffoamv,ffoamh,xcorr1h,xcorr1v
complex(fp) perm1,perm2,rvth,rhth,xperm
!
! First set constants. Then perform the calculation.
wind = Wind_Speed
u10mps = wind
pcc=cos(Sat_Zenith_Angle * DEGREES_TO_RADIANS)
pss=sin(Sat_Zenith_Angle * DEGREES_TO_RADIANS)
ps2=pss*pss
pc2=pcc*pcc
freqghz = Frequency
freqghz2=freqghz*freqghz
u10mps2=u10mps*u10mps
sec=one/pcc
sec2=sec*sec
usec=u10mps*sec
! calculate piom (ellison et al.) xperm
! to calculate xperm of saline water based on piom model.
! convert from kelvin to centigrate and define quadratic and
! cubic functions for later polynomials
tc=SST-273.15_fp
tcsq=tc*tc
tccub=tcsq*tc
! define two relaxation frequencies, tau1 and tau2
tau1=emc(1)+emc(2)*tc+emc(3)*tcsq
tau2=emc(4)+emc(5)*tc+emc(6)*tcsq+emc(7)*tccub
! static xperm estatic=del1+del2+einf
del1=emc(8)+emc(9)*tc+emc(10)*tcsq+emc(11)*tccub
del2=emc(12)+emc(13)*tc+emc(14)*tcsq+emc(15)*tccub
einf=emc(18)+emc(19)*tc
! calculate xperm using double-debye formula
fen=two*pi*freqghz*0.001_fp
fen2=fen**two
den1=one+fen2*tau1*tau1
den2=one+fen2*tau2*tau2
perm_real1=del1/den1
perm_real2=del2/den2
perm_imag1=del1*fen*tau1/den1
perm_imag2=del2*fen*tau2/den2
perm_real=perm_real1+perm_real2+einf
perm_imag=perm_imag1+perm_imag2
xperm=cmplx(perm_real,perm_imag)
! calculate complex fresnel reflection coefficients
! to calculate vertical and horizontal polarised reflectivities
! given xperm at local incidencence angle for all channels
! and profiles
perm1 = sqrt(xperm - cmplx(ps2,zero,fp))
perm2 = xperm*pcc
rhth = (pcc - perm1)/(pcc + perm1)
rvth = (perm2 - perm1)/(perm2 + perm1)
rvertsr=real(rvth,fp)
! rvertsi=dimag(rvth)
rvertsi=aimag(rvth)
rverts=rvertsr*rvertsr+rvertsi*rvertsi
rhorzsr=real(rhth,fp)
! rhorzsi=dimag(rhth)
rhorzsi=aimag(rhth)
rhorzs=rhorzsr*rhorzsr+rhorzsi*rhorzsi
! calculate small scale xcorr to reflection coefficients
xcorr1=exp(emc(21)*u10mps*pc2/freqghz2)
! calculate large scale geometric correction
! to calculate a correction to the fresnel reflection coefficients
! allowing for the presence of large scale roughness
! jc: six coefficients (constant, u, u^2, sec, sec^2, u*sec)
zcv1=emc(24)+emc(25)*freqghz+emc(26)*freqghz2
zcv2=(emc(27)+emc(28)*freqghz+emc(29)*freqghz2)*sec
zcv3=(emc(30)+emc(31)*freqghz+emc(32)*freqghz2)*sec2
zcv4=(emc(33)+emc(34)*freqghz+emc(35)*freqghz2)*u10mps
zcv5=(emc(36)+emc(37)*freqghz+emc(38)*freqghz2)*u10mps2
zcv6=(emc(39)+emc(40)*freqghz+emc(41)*freqghz2)*usec
zch1=emc(42)+emc(43)*freqghz+emc(44)*freqghz2
zch2=(emc(45)+emc(46)*freqghz+emc(47)*freqghz2)*sec
zch3=(emc(48)+emc(49)*freqghz+emc(50)*freqghz2)*sec2
zch4=(emc(51)+emc(52)*freqghz+emc(53)*freqghz2)*u10mps
zch5=(emc(54)+emc(55)*freqghz+emc(56)*freqghz2)*u10mps2
zch6=(emc(57)+emc(58)*freqghz+emc(59)*freqghz2)*usec
! calculate correction for this polarisation
xcorr2v=.01_fp*(zcv1+zcv2+zcv3+zcv4+zcv5+zcv6)
xcorr2h=.01_fp*(zch1+zch2+zch3+zch4+zch5+zch6)
evertr=one-rverts*xcorr1+xcorr2v
ehorzr=one-rhorzs*xcorr1+xcorr2h
! write(769,'(6f12.6)') xcorr1,xcorr2v,xcorr2h,usec,u10mps,evertr
! calculate foam emissivity correction
ffoam=emc(22)*(u10mps**emc(23))
evert=evertr - ffoam*evertr+ ffoam
ehorz=ehorzr - ffoam*ehorzr + ffoam
! write(769,'(6f12.6)') ffoam,evertr,ehorzr,evert,ehorz
rverts5 = rverts
rhorzs5 = rhorzs
xcorr15 = xcorr1
ffoam5 = ffoam
evertr5 = evert
ehorzr5 = ehorz
Emissivity(1) = evertr5
Emissivity(2) = ehorzr5
! write(769,'(5E15.6)') SST,Wind_Speed,Sat_Zenith_Angle,Emissivity(1),Emissivity(2)
!
! Begin K matrix calculation
! Combine horizontal and vertical polarizations.
ehorz = ONE
evert = ONE
! calculate corrected emissivity from corrected refectivity
ehorzr=ehorz - ffoam5*ehorz
ffoamh =-ehorz*ehorzr5 + ehorz
evertr=evert - ffoam5*evert
ffoamv =-evert*evertr5 + evert
! calculate corrected emissivity from corrected refectivity
rhorzs = -ehorzr*xcorr15
xcorr1h = -rhorzs5*ehorzr
xcorr2h = ehorzr
rverts = -evertr*xcorr15
xcorr1v = -rverts5*evertr
xcorr2v = evertr
! calculate foam emissivity correction
! calculate correction for this polarisation
zch4=.01_fp*xcorr2h
zch5=.01_fp*xcorr2h
zch6=.01_fp*xcorr2h
zcv4=.01_fp*xcorr2v
zcv5=.01_fp*xcorr2v
zcv6=.01_fp*xcorr2v
! calculate large scale geometric correction
! to calculate a correction to the fresnel reflection coefficients
! allowing for the presence of large scale roughness
dEH_dWindSpeed = emc(23)*ffoam5/wind*ffoamh + &
zch4*(emc(51)+emc(52)*freqghz+emc(53)*freqghz2) + &
xcorr1h*emc(21)*pc2/freqghz2*xcorr15 + &
zch6*(emc(57)+emc(58)*freqghz+emc(59)*freqghz2)*sec + &
zch5*(emc(54)+emc(55)*freqghz+emc(56)*freqghz2)*two*wind
dEV_dWindSpeed = emc(23)*ffoam5/wind*ffoamv + &
zcv4*(emc(33)+emc(34)*freqghz+emc(35)*freqghz2) + &
xcorr1v*emc(21)*pc2/freqghz2*xcorr15 + &
zcv6*(emc(39)+emc(40)*freqghz+emc(41)*freqghz2)*sec + &
zcv5*(emc(36)+emc(37)*freqghz+emc(38)*freqghz2)*two*wind
END SUBROUTINE Fastem1
END MODULE CRTM_Fastem1