CRTM_Fastem1.f90

Go to the documentation of this file.

!
!-------------------------------------------------------------------------------------------------------------
!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
  real(fp) zch4,xcorr2v,evertr,ehorzr,xcorr2h,ffoam,zcv2,zcv3
  real(fp) xcorr1,zcv1,zcv4,zch1,zch2,zcv5,zcv6,tau2
  real(fp) wind,ehorz,evert,sec,sec2,freqghz2
  real(fp) u10mps2,usec,tccub,tau1,tc,tcsq,freqghz
  real(fp) 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
  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,fp)

!          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