var/da/da_radiance/ossmem.inc

   1 subroutine ossmem(ntype_index, theta,frequency,ts,tv,th,em_vector)
   2 !
   3 !$$$  subprogram documentation block
   4 !                .      .    .                                       .
   5 ! subprogram: iceem_amsua  noaa/nesdis SSM/I emissivity model over snow/ice
   6 !
   7 !   prgmmr: Banghua Yan                 org: nesdis              date: 2004-02-12
   8 !
   9 ! abstract: Simulate microwave emissivity over ocean, sea ice and snow
  10 !           using SSM/I  measurements and surface temperature
  11 !
  12 ! program history log:
  13 !
  14 !      01/2004   : Implement the algorithm for snow/ice emissivity to F90 code by Banghua Yan
  15 !      02/2004   : Modify the code for SSI subsystem                           by Banghua Yan
  16 !      07/2004   : Modify the code for GSI subsystem                           by Kozo Okamoto
  17 !
  18 ! input argument list:
  19 !
  20 !      ntype_index  : surface type
  21 !                 1 : ocean
  22 !                 2 : sea ice
  23 !                 3 : snow over land
  24 !      theta    : local zenith angle in radian  (not used here)
  25 !      frequency: frequency in GHz
  26 !      ts       :  scattering layer temperature (K)
  27 !      tv[1] ~ tv[4]: brightness temperature at four SSM/I vertical polarization
  28 !          tv[1] : 19.35  GHz
  29 !          tv[2] : 22.235 GHz
  30 !          tv[3] : 37     GHz
  31 !          tv[4] : 85     GHz
  32 !
  33 !      th[1] ~ th[3]: brightness temperature at three SSM/I horizontal polarization
  34 !          th[1] : 19.35  GHz
  35 !          th[2] : 37     GHz
  36 !          th[3] : 85     GHz
  37 ! output argument list:
  38 !
  39 !      em_vector     :  emissivity at two polarizations
  40 !           em_vector[1] = eh
  41 !           em_vector[2] = ev
  42 !
  43 ! remarks:
  44 !
  45 ! attributes:
  46 !   language: f90
  47 !   machine:  ibm rs/6000 sp
  48 !
  49 !$$$
  50
  51 !  use kinds, only: r_kind,i_kind
  52 !  use constants, only: one
  53   implicit none
  54
  55   integer(i_kind),parameter :: ntype = 3, nv=4, nh=3,ncoev=5,ncoeh=4
  56   integer(i_kind) ntype_index,ich,k,lp,nch
  57   real(r_kind) theta,frequency,ts,tv(*),th(*),em_vector(*)
  58   real(r_kind) ev(nv),eh(nh),freq_v(nv),freq_h(nh),ev_22
  59   real(r_kind) coe_v(ntype,nv,ncoev),coe_h(ntype,nh,ncoeh),pe , ev_cor,eh_cor
  60   logical data_invalid
  61
  62   data freq_v/19.35_r_kind, 22.235_r_kind, 37.0_r_kind, 85.0_r_kind/
  63   data freq_h/19.35_r_kind, 37.0_r_kind, 85.0_r_kind/
  64
  65 ! ocean
  66   data (coe_v(1,1,k),k=1,5)/ &
  67        -5.650765e-002_r_kind,  2.796378e-003_r_kind,  4.603629e-004_r_kind ,  &
  68        1.163488e-003_r_kind, -6.402050e-004_r_kind/
  69   data (coe_v(1,2,k),k=1,5)/ &
  70        -7.773900e-002_r_kind, -1.390087e-003_r_kind,  4.374652e-003_r_kind , &
  71        1.893395e-003_r_kind, -1.053837e-003_r_kind/
  72   data (coe_v(1,3,k),k=1,5)/ &
  73        -1.774548e-001_r_kind, -1.280647e-003_r_kind,  7.487299e-004_r_kind , &
  74        5.565533e-003_r_kind, -7.623489e-004_r_kind/
  75   data (coe_v(1,4,k),k=1,5)/ &
  76        -2.941845e-001_r_kind, -1.522888e-003_r_kind,  6.942110e-004_r_kind , &
  77        1.798103e-003_r_kind,  3.735965e-003_r_kind/
  78   data (coe_h(1,1,k),k=1,4)/ &
  79        -7.468897e-002_r_kind,  3.759362e-003_r_kind,  1.529964e-004_r_kind , &
  80        -4.894926e-005_r_kind/
  81   data (coe_h(1,2,k),k=1,4)/ &
  82        -1.989357e-001_r_kind,  1.534271e-004_r_kind,  4.117263e-003_r_kind , &
  83        1.201543e-004_r_kind/
  84   data (coe_h(1,3,k),k=1,4)/ &
  85        -3.180339e-001_r_kind, -4.772533e-005_r_kind,  1.822194e-004_r_kind , &
  86        4.691318e-003_r_kind/
  87
  88 ! ice
  89   data (coe_v(2,1,k),k=1,5)/ -8.722723e-002_r_kind,  1.064573e-002_r_kind, &
  90        -5.333843e-003_r_kind, -1.394910e-003_r_kind,  4.007640e-004_r_kind/
  91   data (coe_v(2,2,k),k=1,5)/-1.373924e-001_r_kind,  6.580569e-003_r_kind, &
  92        -9.991220e-004_r_kind, -1.476022e-003_r_kind,  4.131816e-004_r_kind/
  93   data (coe_v(2,3,k),k=1,5)/ -2.329867e-001_r_kind,  6.419856e-003_r_kind, &
  94        -5.260987e-003_r_kind, 3.342582e-003_r_kind,  4.139272e-004_r_kind/
  95   data (coe_v(2,4,k),k=1,5)/ -3.528638e-001_r_kind,  6.342649e-003_r_kind, &
  96        -5.002575e-003_r_kind, -1.469298e-003_r_kind,  5.529711e-003_r_kind/
  97   data (coe_h(2,1,k),k=1,4)/ &
  98        -1.338736e-001_r_kind,  6.229798e-003_r_kind, -2.169491e-003_r_kind,  &
  99        5.706367e-004_r_kind/
 100   data (coe_h(2,2,k),k=1,4)/ &
 101        -2.747500e-001_r_kind,  2.041477e-003_r_kind,  2.581898e-003_r_kind,  &
 102        5.924890e-004_r_kind/
 103   data (coe_h(2,3,k),k=1,4)/ &
 104        -3.889575e-001_r_kind,  2.188889e-003_r_kind, -2.253243e-003_r_kind,  &
 105        5.750499e-003_r_kind/
 106
 107 !snow
 108   data (coe_v(3,1,k),k=1,5)/  1.109066e-001_r_kind,  5.449409e-003_r_kind,  &
 109        1.835799e-004_r_kind, -1.765248e-003_r_kind, -2.996101e-004_r_kind/
 110   data (coe_v(3,2,k),k=1,5)/ 9.356505e-002_r_kind,  1.320617e-003_r_kind,  &
 111        4.449195e-003_r_kind, -1.786960e-003_r_kind, -3.479687e-004_r_kind/
 112   data (coe_v(3,3,k),k=1,5)/ 6.387097e-002_r_kind,  1.252447e-003_r_kind,  &
 113        1.998846e-004_r_kind, 2.680219e-003_r_kind, -3.740141e-004_r_kind/
 114   data (coe_v(3,4,k),k=1,5)/ 4.150689e-002_r_kind,  1.420274e-003_r_kind,  &
 115        1.223339e-004_r_kind, -1.948946e-003_r_kind,  4.248289e-003_r_kind/
 116
 117   data (coe_h(3,1,k),k=1,4)/ &
 118        8.503807e-002_r_kind,  5.357374e-003_r_kind, -1.361660e-003_r_kind, &
 119        -3.319696e-004_r_kind/
 120   data (coe_h(3,2,k),k=1,4)/ &
 121        4.200333e-002_r_kind,  1.278894e-003_r_kind,  2.963129e-003_r_kind, &
 122        -4.087036e-004_r_kind/
 123   data (coe_h(3,3,k),k=1,4)/ &
 124        2.082461e-002_r_kind,  1.438480e-003_r_kind, -1.723992e-003_r_kind,  &
 125        4.194914e-003_r_kind/
 126
 127   save  coe_v,coe_h
 128
 129
 130 !*** Quality check
 131   if(ntype_index < 1) ntype_index = 1
 132   if(ntype_index > 3) ntype_index = 3
 133
 134
 135 ! Initialize
 136   if(ntype_index == 1) then
 137      em_vector(1) = 0.4_r_kind
 138      em_vector(2) = 0.6_r_kind
 139   else if(ntype_index == 2) then
 140      em_vector(1) = 0.7_r_kind
 141      em_vector(2) = 0.8_r_kind
 142   else if(ntype_index == 3) then
 143      em_vector(1) = 0.75_r_kind
 144      em_vector(2) = 0.8_r_kind
 145   end if
 146
 147 ! Data status check
 148   data_invalid = .False.
 149   if ( (ts <= 140.0_r_kind) .or. (ts >= 330.0_r_kind) ) data_invalid = .True.
 150   do ich = 1, nv
 151      if ( (tv(ich) .le. 50.0_r_kind) .or. (tv(ich) .ge. 330.0_r_kind) )  then
 152         data_invalid = .True.
 153         exit
 154      end if
 155   end do
 156   do ich = 1, nh
 157      if ( (th(ich) <= 50.0_r_kind) .or. (th(ich) >= 330.0_r_kind) )  then
 158         data_invalid = .True.
 159         exit
 160      end if
 161   end do
 162   if (data_invalid) RETURN
 163
 164
 165 !*** Get intial emissivity for each frequency
 166
 167 ! v components
 168   do ich=1,nv
 169      ev(ich) =  coe_v(ntype_index,ich,1) + coe_v(ntype_index,ich,2)*tv(1)  &
 170           + coe_v(ntype_index,ich,3)*tv(2) + coe_v(ntype_index,ich,4)*tv(3)  &
 171           + coe_v(ntype_index,ich,5)*tv(4)
 172   end do
 173
 174 ! h components
 175   do ich=1,nh
 176      eh(ich) =  coe_h(ntype_index,ich,1)
 177      do lp =2,4
 178         eh(ich) =  eh(ich) + coe_h(ntype_index,ich,lp)*th(lp-1)
 179      end do
 180   end do
 181
 182
 183 ! *** Emissivity bias value
 184
 185   if (ntype_index == 1) then         ! ocean
 186      pe= 0.001_r_kind + 3.885776e-003_r_kind*(tv(1) - th(1)) +  &
 187           3.727114e-005_r_kind*(tv(3) - th(2)) -  &
 188           1.141903e-004_r_kind*(tv(4) - th(3))
 189   else if (ntype_index == 2) then       ! seaice
 190      pe= 0.011_r_kind + 3.786080e-003_r_kind*(tv(1) - th(1)) -  &
 191           7.217788e-005_r_kind*(tv(3) - th(2)) +  &
 192           1.018791e-004_r_kind*(tv(4) - th(3))
 193   else            ! snow
 194      pe=     -0.002_r_kind + 4.470142e-003_r_kind*(tv(1) - th(1)) -  &
 195           1.991876e-004_r_kind*(tv(3) - th(2)) -  &
 196           1.704354e-005_r_kind*(tv(4) - th(3))
 197   end if
 198
 199   ev_cor = one - pe*(ts-tv(1))/(tv(1)-th(1))
 200   if (ev_cor >  one)         ev_cor = one
 201   if (ev_cor <= 0.2_r_kind) ev_cor = 0.2_r_kind
 202   eh_cor = ev_cor - pe
 203   ev_cor = ev(1) - ev_cor
 204   eh_cor = eh(1) - eh_cor
 205
 206 !*** Calculate emissivity
 207   do ich=1, nv
 208      ev(ich) = ev(ich) - ev_cor
 209      if(ich <= 3) eh(ich) = eh(ich) - eh_cor
 210   end do
 211
 212
 213 !*** Quality control at 22.235 GHz
 214   ev_22 = ev(1) + (ev(3)-ev(1))*(22.235_r_kind-19.35_r_kind)/(37.0_r_kind-19.35_r_kind)
 215 !/\ type
 216   if( (ev(2) .gt. ev(1)) .and. (ev(2) .gt. ev(3)) ) ev(2) = ev_22
 217 !\/ type
 218   if( (ev(2) .lt. ev(1)) .and. (ev(2) .lt. ev(3)) ) ev(2) = ev_22
 219
 220
 221 !*** Interpolate emissivity at a certain frequency
 222
 223 ! v-component
 224   nch = 4
 225   do ich=1,nv
 226      if(frequency <= freq_v(ich)) then
 227         nch = ich
 228         exit
 229      end if
 230   end do
 231
 232   if (nch == 1) then
 233      em_vector(2) = ev(1)
 234   else
 235      if (nch == 4) then
 236         em_vector(2) = ev(4)
 237      else
 238         em_vector(2) = ev(nch-1) + (ev(nch) - ev(nch-1))* &
 239              (frequency - freq_v(nch-1))/(freq_v(nch) - freq_v(nch-1))
 240      end if
 241   end if
 242
 243 ! h-component
 244   nch = 3
 245   do ich=1,nh
 246      if(frequency <= freq_h(ich)) then
 247         nch = ich
 248         exit
 249      end if
 250   end do
 251   if (nch == 1) then
 252      em_vector(1) = eh(1)
 253   else
 254      if (nch == 3) then
 255         em_vector(1) = eh(3)
 256      else
 257         em_vector(1) = eh(nch-1) + (eh(nch) - eh(nch-1))* &
 258              (frequency - freq_h(nch-1))/(freq_h(nch) - freq_h(nch-1))
 259      end if
 260
 261   end if
 262
 263 end subroutine ossmem