var/da/da_radiance/seaem.inc

   1 subroutine seaem(freqghz,zasat,zlsat,ts5,u10,v10,ehorz,evert)
   2 !
   3 !------------------------------------------------------------------------------
   4 !  PURPOSE: Calculate microwave surface emissivity over sea.
   5 !
   6 !  METHOD:  adopted from GSI code emiss
   7 !       out V and H polarized emissivitys for used in RTTOV
   8 !
   9 !  HISTORY: 03/10/2005 - Creation            Zhiquan Liu
  10 !
  11 !------------------------------------------------------------------------------
  12 !   input argument list:
  13 !     freqghz  - microwave freqency
  14 !     zasat    - local satellite zenith angle in radians
  15 !     zlsat    - satellite look angle in radians (not used)
  16 !     ts5      - skin temperature
  17 !     u10      - 10m u wind
  18 !     v10      - 10m v wind
  19 !
  20 !   output argument list:
  21 !     ehorz    - horizontal polarization emissivity
  22 !     evert    - vertical polarization emissivity
  23 !
  24 ! ...............................................................
  25
  26   implicit none
  27
  28 ! Declare passed variables.
  29   real(r_kind), intent(in):: freqghz
  30   real(r_kind), intent(in):: ts5
  31   real(r_kind), intent(in):: zasat,zlsat,u10,v10
  32   real(r_kind),intent(out):: ehorz, evert
  33   real(r_kind),dimension(59):: emc
  34
  35 ! Declare local variables
  36   integer(i_kind) kcho,n,kch,nn,nnp,i
  37   integer(i_kind) error_status
  38   integer(i_kind) quiet
  39 !  integer(i_kind),dimension(nchan)::indx
  40
  41   real(r_kind) zch4,xcorr2v,evertr,ehorzr,xcorr2h,ffoam,zcv2,zcv3
  42   real(r_kind) xcorr1,zcv1,zcv4,zch1,zch2,zcv5,zcv6,tau2,degre
  43   real(r_kind) wind,sec,sec2,freqghz2,dtde
  44   real(r_kind) u10mps2,usec,tccub,tau1,tc,tcsq,term2
  45   real(r_kind) term1,u10mps,ps2,pc2,pcc,pss,rvertsi,rverts,rvertsr
  46   real(r_kind) rverts5,rhorzs5,xcorr15,ffoam5,evertr5,ehorzr5
  47   real(r_kind) perm_real,perm_imag,rhorzsr,zch5,zch6,zch3,rhorzsi
  48   real(r_kind) rhorzs,perm_imag2,einf,fen,del2,del1,fen2,perm_real2
  49   real(r_kind) perm_imag1,perm_real1,den1,den2
  50
  51   complex(r_kind) perm1,perm2,rvth,rhth,xperm
  52
  53 !  integer     :: ipolar(nchan)
  54 !  real        :: polar(nchan)
  55
  56 !
  57 !  Start emiss here
  58 !
  59 !  uuk=zero
  60 !  vvk=zero
  61
  62 ! Explanation for emc :
  63 ! emc(59): Emissivity model data
  64 ! Permittivity model data (Lamkaouchi model)
  65 !   [1-3]: Temperature polynomial coefficients for Tau1 - Lamkaouchi (1996)
  66 !   [4-7]: Temperature polynomial coefficients for Tau2 - Lamkaouchi (1996)
  67 !  [8-11]: Temperature polynomial coefficients for Del1 - Lamkaouchi (1996)
  68 ! [12-15]: Temperature polynomial coefficients for Del2 - Lamkaouchi (1996)
  69 ! [16-17]: Temperature polynomial coefficients for static permittivity - Lamkaouchi (1996)
  70 ! [18-19]: Temperature polynomial coefficients for infinite freq. permittivity - Lamkaouchi (1996)
  71 ! Pi is stored for good measure
  72 !    [20]: Stored value of Pi
  73 ! Bragg scattering correction coefficients
  74 !    [21]: Scaling factor for small scale correction - see English (1997)
  75 ! Foam model coefficients for Monahan model
  76 !    [22]: First coefficient in Monahan foam model (neutral stability)  - see English (1997)
  77 !    [23]: Second coefficient in Monahan foam model (neutral stability) - see English (1997)
  78 ! Alternative permittivity model (Liebe)
  79 !    [30]: a1 in Liebe's dielectric model - see Liebe (1989)
  80 !    [31]: b1 in Liebe's dielectric model - see Liebe (1989)
  81 !    [32]: c1 in Liebe's dielectric model - see Liebe (1989)
  82 !    [33]: c2 in Liebe's dielectric model - see Liebe (1989)
  83 !    [34]: d1 in Liebe's dielectric model - see Liebe (1989)
  84 !    [35]: d2 in Liebe's dielectric model - see Liebe (1989)
  85 !    [36]: d3 in Liebe's dielectric model - see Liebe (1989)
  86 !    [37]: e1 in Liebe's dielectric model - see Liebe (1989)
  87 !    [38]: e2 in Liebe's dielectric model - see Liebe (1989)
  88 ! Version 2 of large scale correction which *DOES�»* take account of
  89 ! hemispherical scattering.
  90 ! 1.) Vertical polarisation mode
  91 !    [24]: Term a00 in vertical pol of large scale correction model
  92 !    [25]: Term a01 in vertical pol mode of large scale correction model
  93 !    [26]: Term a02 in vertical pol mode of large scale correction model
  94 !    [27]: Term a10 in vertical pol mode of large scale correction model
  95 !    [28]: Term a11 in vertical pol mode of large scale correction model
  96 !    [29]: Term a12 in vertical pol mode of large scale correction model
  97 !    [30]: Term a20 in vertical pol mode of large scale correction model
  98 !    [31]: Term a21 in vertical pol mode of large scale correction model
  99 !    [32]: Term a22 in vertical pol mode of large scale correction model
 100 !    [33]: Term a30 in vertical pol mode of large scale correction model
 101 !    [34]: Term a31 in vertical pol mode of large scale correction model
 102 !    [35]: Term a32 in vertical pol mode of large scale correction model
 103 !    [36]: Term a40 in vertical pol mode of large scale correction model
 104 !    [37]: Term a41 in vertical pol mode of large scale correction model
 105 !    [38]: Term a42 in vertical pol mode of large scale correction model
 106 !    [39]: Term a50 in vertical pol mode of large scale correction model
 107 !    [40]: Term a51 in vertical pol mode of large scale correction model
 108 !    [41]: Term a52 in vertical pol mode of large scale correction model
 109 ! 2. ) Horizontal polarisation mode
 110 !    [42]: Term a00 in horizontal pol mode of large scale correction model
 111 !    [43]: Term a01 in horizontal pol mode of large scale correction model
 112 !    [44]: Term a02 in horizontal pol mode of large scale correction model
 113 !    [45]: Term a10 in horizontal pol mode of large scale correction model
 114 !    [46]: Term a11 in horizontal pol mode of large scale correction model
 115 !    [47]: Term a12 in horizontal pol mode of large scale correction model
 116 !    [48]: Term a20 in horizontal pol mode of large scale correction model
 117 !    [49]: Term a21 in horizontal pol mode of large scale correction model
 118 !    [50]: Term a22 in horizontal pol mode of large scale correction model
 119 !    [51]: Term a30 in horizontal pol mode of large scale correction model
 120 !    [52]: Term a31 in horizontal pol mode of large scale correction model
 121 !    [53]: Term a32 in horizontal pol mode of large scale correction model
 122 !    [54]: Term a40 in horizontal pol mode of large scale correction model
 123 !    [55]: Term a41 in horizontal pol mode of large scale correction model
 124 !    [56]: Term a42 in horizontal pol mode of large scale correction model
 125 !    [57]: Term a50 in horizontal pol mode of large scale correction model
 126 !    [58]: Term a51 in horizontal pol mode of large scale correction model
 127 !    [59]: Term a52 in horizontal pol mode of large scale correction model
 128
 129    emc = (/&
 130        0.175350E+02_r_kind, -.617670E+00_r_kind,  .894800E-02_r_kind,  .318420E+01_r_kind,&
 131        0.191890E-01_r_kind, -.108730E-01_r_kind,  .258180E-03_r_kind,  .683960E+02_r_kind,&
 132        -.406430E+00_r_kind,  .228320E-01_r_kind, -.530610E-03_r_kind,  .476290E+01_r_kind,&
 133        0.154100E+00_r_kind, -.337170E-01_r_kind,  .844280E-03_r_kind,  .782870E+02_r_kind,&
 134        -.434630E-02_r_kind,  .531250E+01_r_kind, -.114770E-01_r_kind,  .314159E+01_r_kind,&
 135        -.100000E+01_r_kind,  .195000E-04_r_kind,  .255000E+01_r_kind, -.637182E+01_r_kind,&
 136        0.253918E-01_r_kind,  .357569E-04_r_kind,  .942928E+01_r_kind, -.332839E-01_r_kind,&
 137        -.647724E-04_r_kind, -.329282E+01_r_kind,  .965450E-02_r_kind,  .281588E-04_r_kind,&
 138        0.252676E+00_r_kind,  .343867E-02_r_kind, -.156362E-04_r_kind, -.156669E-03_r_kind,&
 139        0.139485E-04_r_kind, -.407633E-07_r_kind, -.141316E+00_r_kind, -.356556E-02_r_kind,&
 140        0.142869E-04_r_kind, -.240701E+01_r_kind, -.563888E-01_r_kind,  .325227E-03_r_kind,&
 141        0.296005E+01_r_kind,  .704675E-01_r_kind, -.426440E-03_r_kind, -.751252E+00_r_kind,&
 142        -.191934E-01_r_kind,  .125937E-03_r_kind, -.288253E+00_r_kind, -.102655E-02_r_kind,&
 143        0.226701E-05_r_kind, -.119072E-02_r_kind, -.263165E-04_r_kind,  .114597E-06_r_kind,&
 144        0.406300E+00_r_kind,  .200031E-02_r_kind, -.781635E-05_r_kind/)
 145
 146 !      ----- sea (ice-free) MW  -------
 147
 148 !          Open ocean points
 149
 150 !          First set constants.  Then perform the calculation.
 151 !           wind  = f10(n)*sqrt(uu5(n)*uu5(n)+vv5(n)*vv5(n))  ! wind speed in m/s
 152            wind  = sqrt(u10*u10+v10*v10)  ! 10m wind speed in m/s
 153            u10mps  = wind
 154            pcc=cos(zasat)
 155            pss=sin(abs(zasat))
 156 !          pcl2=cos(zlsat)**2
 157 !          psl2=sin(zlsat)**2
 158            ps2=pss*pss
 159            pc2=pcc*pcc
 160            freqghz2=freqghz*freqghz
 161            u10mps2=u10mps*u10mps
 162            sec=one/pcc
 163            sec2=sec*sec
 164            usec=u10mps*sec
 165
 166 !          calculate piom (ellison et al.) xperm
 167 !          to calculate xperm of saline water based on piom model.
 168 !          convert from kelvin to centigrate and define quadratic and
 169 !          cubic functions for later polynomials
 170            tc=ts5-t_kelvin
 171            tcsq=tc*tc
 172            tccub=tcsq*tc
 173
 174 !          define two relaxation frequencies, tau1 and tau2
 175            tau1=emc(1)+emc(2)*tc+emc(3)*tcsq
 176            tau2=emc(4)+emc(5)*tc+emc(6)*tcsq+emc(7)*tccub
 177
 178 !          static xperm estatic=del1+del2+einf
 179            del1=emc(8)+emc(9)*tc+emc(10)*tcsq+emc(11)*tccub
 180            del2=emc(12)+emc(13)*tc+emc(14)*tcsq+emc(15)*tccub
 181            einf=emc(18)+emc(19)*tc
 182
 183 !          calculate xperm using double-debye formula
 184            fen=two*pi*freqghz*0.001_r_kind
 185            fen2=fen**two
 186            den1=one+fen2*tau1*tau1
 187            den2=one+fen2*tau2*tau2
 188            perm_real1=del1/den1
 189            perm_real2=del2/den2
 190            perm_imag1=del1*fen*tau1/den1
 191            perm_imag2=del2*fen*tau2/den2
 192            perm_real=perm_real1+perm_real2+einf
 193            perm_imag=perm_imag1+perm_imag2
 194            xperm=dcmplx(perm_real,perm_imag)
 195
 196 !          calculate complex fresnel reflection coefficients
 197 !          to calculate vertical and horizontal polarised reflectivities
 198 !          given xperm at local incidencence angle for all channels
 199 !          and profiles
 200            perm1 = cdsqrt(xperm - dcmplx(ps2,zero))
 201            perm2  = xperm*pcc
 202            rhth = (pcc - perm1)/(pcc + perm1)
 203            rvth = (perm2 - perm1)/(perm2 + perm1)
 204            rvertsr=real(rvth)
 205            rvertsi=imag(rvth)
 206            rverts=rvertsr*rvertsr+rvertsi*rvertsi
 207            rhorzsr=real(rhth)
 208            rhorzsi=imag(rhth)
 209            rhorzs=rhorzsr*rhorzsr+rhorzsi*rhorzsi
 210
 211 !          calculate small scale xcorr to reflection coefficients
 212            xcorr1=exp(emc(21)*u10mps*pc2/freqghz2)
 213
 214 !          calculate large scale geometric correction
 215 !          to calculate a correction to the fresnel reflection coefficients
 216 !          allowing for the presence of large scale roughness
 217
 218 !          jc: six coefficients (constant, u, u^2, sec, sec^2, u*sec)
 219            zcv1=emc(24)+emc(25)*freqghz+emc(26)*freqghz2
 220            zcv2=(emc(27)+emc(28)*freqghz+emc(29)*freqghz2)*sec
 221            zcv3=(emc(30)+emc(31)*freqghz+emc(32)*freqghz2)*sec2
 222            zcv4=(emc(33)+emc(34)*freqghz+emc(35)*freqghz2)*u10mps
 223            zcv5=(emc(36)+emc(37)*freqghz+emc(38)*freqghz2)*u10mps2
 224            zcv6=(emc(39)+emc(40)*freqghz+emc(41)*freqghz2)*usec
 225            zch1=emc(42)+emc(43)*freqghz+emc(44)*freqghz2
 226            zch2=(emc(45)+emc(46)*freqghz+emc(47)*freqghz2)*sec
 227            zch3=(emc(48)+emc(49)*freqghz+emc(50)*freqghz2)*sec2
 228            zch4=(emc(51)+emc(52)*freqghz+emc(53)*freqghz2)*u10mps
 229            zch5=(emc(54)+emc(55)*freqghz+emc(56)*freqghz2)*u10mps2
 230            zch6=(emc(57)+emc(58)*freqghz+emc(59)*freqghz2)*usec
 231
 232 !          calculate correction for this polarisation
 233            xcorr2v=.01_r_kind*(zcv1+zcv2+zcv3+zcv4+zcv5+zcv6)
 234            xcorr2h=.01_r_kind*(zch1+zch2+zch3+zch4+zch5+zch6)
 235
 236            evertr=one-rverts*xcorr1+xcorr2v
 237            ehorzr=one-rhorzs*xcorr1+xcorr2h
 238
 239 !          calculate foam emissivity correction
 240            ffoam=emc(22)*(u10mps**emc(23))
 241            evert=evertr - ffoam*evertr+ ffoam
 242            ehorz=ehorzr - ffoam*ehorzr + ffoam
 243
 244 !           rverts5 = rverts
 245 !           rhorzs5 = rhorzs
 246 !           xcorr15 = xcorr1
 247 !           ffoam5 = ffoam
 248 !           evertr5 = evertr
 249 !           ehorzr5 = ehorzr
 250
 251
 252 !          Combine horizontal and vertical polarizations.
 253 !           call ehv2pem( ehorz,evert,zlsat(n),polar(kch), pems5(nn) )
 254
 255 !          Begin K matrix calculation
 256
 257 !          Combine horizontal and vertical polarizations.
 258 !           dtde=one
 259 !           call adm_ehv2pem(zlsat(n),polar(kch),dtde, ehorz,evert )
 260
 261 !          calculate corrected emissivity from corrected refectivity
 262 !           ehorzr=ehorz - ffoam5*ehorz
 263 !           ffoam =-ehorz*ehorzr5 + ehorz
 264 !           evertr=evert - ffoam5*evert
 265 !           ffoam =ffoam-evert*evertr5 + evert
 266
 267 !          calculate corrected emissivity from corrected refectivity
 268 !           rhorzs = -ehorzr*xcorr15
 269 !           xcorr1 = -rhorzs5*ehorzr
 270 !           xcorr2h = ehorzr
 271 !           rverts = -evertr*xcorr15
 272 !           xcorr1 = xcorr1 - rverts5*evertr
 273 !           xcorr2v = evertr
 274
 275 !          calculate foam emissivity correction
 276 !          calculate correction for this polarisation
 277 !           zch4=.01_r_kind*xcorr2h
 278 !           zch5=.01_r_kind*xcorr2h
 279 !           zch6=.01_r_kind*xcorr2h
 280 !           zcv4=.01_r_kind*xcorr2v
 281 !           zcv5=.01_r_kind*xcorr2v
 282 !           zcv6=.01_r_kind*xcorr2v
 283
 284 !          calculate large scale geometric correction
 285 !          to calculate a correction to the fresnel reflection coefficients
 286 !          allowing for the presence of large scale roughness
 287 !
 288 !          jc: six coefficients (constant, u, u^2, sec, sec^2, u*sec)
 289 !           u10mps = emc(23)*ffoam5/wind*ffoam +                       &
 290 !                    zch4*(emc(51)+emc(52)*freqghz+emc(53)*freqghz2) + &
 291 !                    zcv4*(emc(33)+emc(34)*freqghz+emc(35)*freqghz2) + &
 292 !                    xcorr1*emc(21)*pc2/freqghz2*xcorr15
 293
 294 !           usec=zch6*(emc(57)+emc(58)*freqghz+emc(59)*freqghz2) + &
 295 !                zcv6*(emc(39)+emc(40)*freqghz+emc(41)*freqghz2)
 296 !           u10mps2=zch5*(emc(54)+emc(55)*freqghz+emc(56)*freqghz2) + &
 297 !                   zcv5*(emc(36)+emc(37)*freqghz+emc(38)*freqghz2)
 298 !          calculate small scale xcorr to reflection coefficients
 299 !          the following lines are commented out because a warning will
 300 !          be printed from dcalmkaouchi if freqghz<10.
 301
 302 !           u10mps  = u10mps+usec*sec+u10mps2*two*wind
 303 !           u10mps  = f10(n)*f10(n)/wind*u10mps
 304 !           uuk(nn) = uu5(n)*u10mps
 305 !           vvk(nn) = vv5(n)*u10mps
 306
 307
 308 !        Load emissivity into array for radiative transfer model
 309 !       (pems5) and diagnostic output file (emissav).
 310
 311 !        pems5(nn)       = max(zero,min(pems5(nn),one))
 312
 313 ! End of routine.
 314   return
 315
 316 end subroutine seaem