var/da/da_radiance/gsi_emiss.inc

   1 subroutine gsi_emiss(inst,knchpf,kprof,kchan,knchan,indexn,zasat,zlsat,&
   2                 polar, isflg,pemsv,pemsh,pems5, ts5,soiltype5,soilt5,soilm5, &
   3                 vegtype5,vegf5,snow5,itype,nsdata,nchan,  &
   4                 btm,amsua,amsub,ssmi,u10,v10)
   5 !                .      .    .                                       .
   6 ! subprogram:    emiss       compute emissivity for IR and microwave
   7 !   prgmmr: treadon          org: np23                date: 1998-02-20
   8 !
   9 ! abstract: compute surface emissivity for IR and microwave channels.
  10 !
  11 ! program history log:
  12 !   1998-02-20  treadon - gather all emissivity calculations from
  13 !                         setuprad and move into this subroutine
  14 !   2004-07-23  weng,yan,okamoto - incorporate MW land and snow/ice emissivity
  15 !                         models for AMSU-A/B and SSM/I
  16 !   2004-08-02  treadon - add only to module use, add intent in/out
  17 !   2004-11-22  derber - modify for openMP
  18 !   2005-01-20  okamoto- add preprocessing for ocean MW emissivity jacobian
  19 !   2005-03-05  derber- add adjoint of surface emissivity to this routine
  20 !   2005-09-25  Zhiquan Liu- adopted for wrf-var
  21 !
  22 !
  23 !   input argument list:
  24 !     inst     - wrf-var internal instrument number
  25 !     knchpf   - total number of profiles (obs) times number of
  26 !                channels for each profile.
  27 !     kprof    - profile number array
  28 !     kchan    - channel number array
  29 !     kochan   - old channel number array
  30 !     knchan   - number of channels for each profile
  31 !     indexn   - index pointing from channel,profile to location in array
  32 !     zasat    - local satellite zenith angle in radians
  33 !     zlsat    - satellite look angle in radians
  34 !     isflg    - surface flag
  35 !                0 sea
  36 !                1 sea ice
  37 !                2 land
  38 !                3 snow
  39 !                4 mixed predominately sea
  40 !                5 mixed predominately sea ice
  41 !                6 mixed predominately land
  42 !                7 mixed predominately snow
  43 !     ts5      - skin temperature
  44 !     soiltype5- soil type
  45 !     soilt5   - soil temperature
  46 !     soilm5   - soil moisture (volumatic ratio 0.0 ~ 1.0)
  47 !     vegtype5 - vegetation type
  48 !     vegf5    - vegetation fraction
  49 !     snow5    - snow depth   [mm]
  50 !     itype    - ir/microwave instrument type
  51 !     nsdata   - number of profiles
  52 !     nchan    - maximum number of channels for this satellite/instrument
  53 !     kidsat   - satellite id
  54 !     btm      - observation tb
  55 !     amsua    - logical true if amsua is processed
  56 !     amsub    - logical true if amsub is processed
  57 !     ssmi     - logical true if ssmi  is processed
  58 !     u10      - 10m u wind
  59 !     v10      - 10m v wind
  60 !     f10      - factor reducing lowest sigma level to 10m wind
  61 !
  62 !   output argument list:
  63 !     emissav  - same as pems5 but for diagnostic array
  64 !     uuk      - adjoint of lowest sigma level u wind(microwave over ocean only)
  65 !     vvk      - adjoint of lowest sigma level v wind(microwave over ocean only)
  66 !     pems5    - surface emissivity at obs location
  67 !
  68 !     NOTE:  pems5 is passed through include file "prfvark3.h"
  69 !
  70 !   other important variables
  71 !    polar: channel polarization (0=vertical, 1=horizontal, or
  72 !                                  0 to 1=mix of V and H)
  73 !
  74 !
  75 ! ...............................................................
  76 !             land snow ice sea
  77 ! landem()     o   o    x   x  : for all MW but f<80GHz,lower accuracy over snow
  78 ! snwem_amsu() x   o    x   x  : for AMSU-A/B
  79 ! iceem_amsu() x   x    o   x  : for AMSU-A/B
  80 ! emiss_ssmi() x   o    o   x  : for SSM/I
  81 !
  82 ! if(snow)
  83 !   if(amsua .or. amsub) call snwem_amsu()
  84 !   else if(ssmi)        call emiss_ssmi()
  85 !   else                 call landem()
  86 ! if(land)               call landem()
  87 ! if(ice)
  88 !   if(amsua .or. amsub) call iceem_amsu()
  89 !   else if(ssmi)        call emiss_ssmi()
  90 !   else                 emissivity=0.92
  91 ! if(ocean) calculate
  92 ! ..............................................................
  93 !
  94 ! attributes:
  95 !   language: f90
  96 !   machine:  IBM sp
  97 !
  98 !$$$
  99 !  use kinds, only: r_kind,r_single,i_kind
 100 !  use error_handler, only: failure  ! change to local parameter
 101 !  use irsse_model, only: forward_irsse
 102 !  use radinfo, only: polar,jppf,newchn_ir
 103 !  use spectral_coefficients, only: sc
 104 !  use constants, only: zero,one,rad2deg,two,pi
 105   implicit none
 106
 107   integer, parameter  :: FAILURE = 3, jppf = 1
 108
 109 ! Declare passed variables.
 110   integer(i_kind),intent(in):: inst, knchpf,nsdata,nchan,itype
 111   integer(i_kind),dimension(nchan*nsdata),intent(in):: kprof,kchan
 112   integer(i_kind),dimension(nchan,nsdata),intent(in):: indexn
 113   integer(i_kind),dimension(nsdata),intent(in):: isflg,knchan
 114   real(r_kind),dimension(nchan,jppf), intent(in) :: btm
 115   real(r_kind),dimension(nsdata),intent(in):: ts5,snow5,soiltype5,&
 116        soilt5,soilm5,vegtype5,vegf5
 117   real(r_kind),dimension(nsdata),intent(in):: zasat,zlsat,u10,v10
 118   real , dimension(nchan),intent(in)       :: polar
 119 !  real(r_single),dimension(nchan,nsdata),intent(out):: emissav
 120 !  real(r_kind),dimension(nchan,nsdata),intent(out):: emissav
 121   real(r_kind),dimension(nsdata*nchan),intent(out):: pemsv, pemsh, pems5
 122 !  real(r_kind),dimension(59):: emc
 123
 124 ! Declare local variables
 125   integer(i_kind) kcho,n,kch,nn,nnp,i
 126   integer(i_kind) error_status
 127   integer(i_kind) quiet
 128   integer(i_kind),dimension(nchan)::indx
 129
 130   real(r_kind) zch4,xcorr2v,evertr,ehorzr,xcorr2h,ffoam,zcv2,zcv3
 131   real(r_kind) xcorr1,zcv1,zcv4,zch1,zch2,zcv5,zcv6,tau2,degre
 132   real(r_kind) wind,ehorz,evert,sec,sec2,freqghz2,dtde
 133   real(r_kind) u10mps2,usec,tccub,tau1,tc,tcsq,term2,freqghz
 134   real(r_kind) term1,u10mps,ps2,pc2,pcc,pss,rvertsi,rverts,rvertsr
 135   real(r_kind) rverts5,rhorzs5,xcorr15,ffoam5,evertr5,ehorzr5
 136   real(r_kind) perm_real,perm_imag,rhorzsr,zch5,zch6,zch3,rhorzsi
 137   real(r_kind) rhorzs,perm_imag2,einf,fen,del2,del1,fen2,perm_real2
 138   real(r_kind) perm_imag1,perm_real1,den1,den2
 139   real(r_kind),dimension(1):: emissir
 140
 141   complex(r_kind) perm1,perm2,rvth,rhth,xperm
 142
 143 ! -------- ice/snow-MW em  ---------
 144   integer(i_kind)     :: surf_type
 145 !  real(r_kind),dimension(nchan,jppf):: btm
 146   real(r_kind):: tbasnw(4),tbbsnw(2),tbbmi(nchan)
 147   logical      :: amsuab
 148   logical      :: amsua,amsub,ssmi
 149   real(r_kind):: ice_depth
 150
 151   integer     :: ipolar(nchan)
 152 !  real        :: polar(nchan)
 153
 154 !
 155 !  Start emiss here
 156 !
 157 !  uuk=zero
 158 !  vvk=zero
 159
 160 !  ////////////  IR emissivity //////////////////
 161 #ifdef RTTOV
 162
 163 ! Compute or set the emissivity for IR channels.
 164   if(itype == 0)then
 165
 166 !$omp parallel do  schedule(dynamic,1) private(n) &
 167 !$omp private(nn,nnp,wind,degre,kch,i,indx,error_status)
 168      do n = 1,nsdata
 169         nn   = indexn(1,n)
 170         nnp  = nn+knchan(n)-1
 171         if (isflg(n)==0 .or. isflg(n) == 4) then          ! sea
 172
 173 !      use RTTOV ISEM calculated values over sea .
 174           pemsv(nn:nnp) = zero
 175           pemsh(nn:nnp) = zero
 176
 177 !       Use fixed IR emissivity values over land, snow, and ice.
 178
 179         else if(isflg(n) == 1 .or. isflg(n) == 5) then
 180           pemsv(nn:nnp) = 0.98_r_kind                     ! sea ice
 181           pemsh(nn:nnp) = 0.98_r_kind
 182         else if(isflg(n) == 2 .or. isflg(n) == 6) then
 183           pemsv(nn:nnp) = 0.97_r_kind                     ! land
 184           pemsh(nn:nnp) = 0.97_r_kind
 185         else
 186           pemsv(nn:nnp) = one                             ! snow
 187           pemsh(nn:nnp) = one
 188         end if
 189 !        do i=1,knchan(n)
 190 !           emissav(i,n) = pemsv(nn+i-1)
 191 !        end do
 192      end do
 193
 194
 195 !  ////////////  MW emissivity //////////////////
 196   else if(itype == 1)then
 197
 198 ! Set emissivity for microwave channels
 199 !
 200      surf_type=0
 201      amsuab = amsua .or. amsub
 202
 203 !$omp parallel do  schedule(dynamic,1) private(nn) &
 204 !$omp private(n,kch,kcho,tbasnw,tbbsnw,tbbmi,evert,ehorz,surf_type,ice_depth) &
 205 !$omp private(freqghz,u10mps,pcc,pss,ps2,pc2,freqghz2,u10mps2,sec) &
 206 !$omp private(sec2,usec,tc,tcsq,tccub,tau1,tau2,del1,del2,einf,fen,fen2,den1) &
 207 !$omp private(den2,perm_real1,perm_real2,perm_imag1,perm_imag2,perm_real) &
 208 !$omp private(perm_imag,xperm,perm1,perm2,rhth,rvth,rvertsr,rvertsi) &
 209 !$omp private(rverts,rhorzsr,rhorzsi,rhorzs,xcorr1,zcv1,zcv2) &
 210 !$omp private(zcv3,zcv4,zcv5,zcv6,zch1,zch2,zch3,zch4,zch5,zch6) &
 211 !$omp private(xcorr2v,xcorr2h,ffoam,evertr,ehorzr) &
 212 !$omp private(rverts5,rhorzs5,xcorr15,ffoam5,evertr5,ehorzr5,dtde)
 213      do nn = 1,knchpf
 214         n    = kprof(nn)
 215         kch  = kchan(nn)
 216 !        ipolar(kch) = coefs(inst)%fastem_polar(kch)
 217 !        if (ipolar(kch) == 3) polar(kch) = 0.0   ! vertical polar
 218 !        if (ipolar(kch) == 4) polar(kch) = 1.0   ! horizontal polar
 219 !        kcho = kochan(nn)
 220
 221         pemsv(nn) = zero  !default value
 222         pemsh(nn) = zero
 223         pems5(nn) = zero
 224
 225 !        freqghz = sc%frequency(kch)
 226         freqghz = coefs(inst)%frequency_ghz(kch)  ! GHz
 227
 228         if((isflg(n)==3 .or. isflg(n) == 7).and.snow5(n)>0.1_r_kind)then
 229
 230 !      ----- snow MW -------
 231
 232 !          land/snow points
 233
 234            if(amsuab) then
 235               if(amsua) then
 236                  tbasnw(1:3) = btm(1:3,n)
 237                  tbasnw(4)   = btm(15,n)
 238                  tbbsnw(1:2) = -999.9_r_kind
 239               else if(amsub) then
 240                  tbasnw(1:4) = -999.9_r_kind
 241                  tbbsnw(1:2) = btm(1:2,n)
 242               end if
 243               call snwem_amsu(zasat(n),freqghz, &
 244                    snow5(n),ts5(n),tbasnw,tbbsnw,ehorz,evert )
 245
 246 !               call ehv2pem( ehorz,evert,zlsat(n),polar(kch), pems5(nn) )
 247               pemsv(nn) = evert !because esv=esh
 248               pemsh(nn) = ehorz
 249               pems5(nn) = evert
 250            else if(ssmi) then
 251               surf_type=3
 252               tbbmi(1:knchan(n)) = btm(1:knchan(n),n)
 253               call emiss_ssmi(   &
 254                    surf_type,zasat(n),freqghz, &
 255                    soilm5(n),vegf5(n),vegtype5(n),soiltype5(n), &
 256                    soilt5(n),ts5(n),snow5(n),tbbmi,ehorz,evert )
 257               pems5(nn) = (one-polar(kch))*evert + polar(kch)*ehorz
 258               pemsv(nn) = evert !because esv=esh
 259               pemsh(nn) = ehorz
 260            else
 261               if(freqghz<80._r_kind)then  !snow & f<80GHz
 262 !               currently landem only works for frequencies < 80 Ghz
 263                if (use_landem) then
 264                  call landem(zasat(n),freqghz,  &
 265                       soilm5(n),vegf5(n),vegtype5(n),soiltype5(n),soilt5(n), &
 266                       ts5(n),snow5(n),ehorz,evert)
 267                  call ehv2pem( ehorz,evert,zlsat(n),polar(kch), pems5(nn) )
 268                  pemsv(nn) = evert !because esv=esh
 269                  pemsh(nn) = ehorz
 270                else
 271                  pemsv(nn) = 0.90_r_kind
 272                  pemsh(nn) = 0.90_r_kind
 273                  pems5(nn) = 0.90_r_kind
 274                end if
 275               else  !snow & f>=80GHz
 276                  pemsv(nn) = 0.90_r_kind
 277                  pemsh(nn) = 0.90_r_kind
 278                  pems5(nn) = 0.90_r_kind
 279               end if
 280            end if  !snow & amsuab/ssmi/othermw
 281
 282
 283         else if(isflg(n) == 1 .or. isflg(n) == 5)then
 284 !      ----- sea ice MW  -------
 285
 286            if(amsuab) then
 287               ice_depth=zero
 288               if(amsua) then
 289                  tbasnw(1:3) = btm(1:3,n)
 290                  tbasnw(4)   = btm(15,n)
 291                  tbbsnw(1:2) = -999.9_r_kind
 292               else if(amsub) then
 293                  tbasnw(1:4) = -999.9_r_kind
 294                  tbbsnw(1:2) = btm(1:2,n)
 295               end if
 296               call iceem_amsu(  &
 297                    zasat(n),freqghz,ice_depth,ts5(n),&
 298                    tbasnw,tbbsnw,ehorz,evert )
 299 !                call ehv2pem( ehorz,evert,zlsat(n),polar(kch), pems5(nn) )
 300               pemsv(nn) = evert  !because esv=esh
 301               pemsh(nn) = ehorz
 302               pems5(nn) = evert
 303            else if(ssmi) then
 304               surf_type=2
 305               tbbmi(1:knchan(n)) = btm(1:knchan(n),n)
 306               call emiss_ssmi(   &
 307                    surf_type,zasat(n),freqghz, &
 308                    soilm5(n),vegf5(n),vegtype5(n),soiltype5(n), &
 309                    soilt5(n),ts5(n),snow5(n),tbbmi,ehorz,evert )
 310               pems5(nn) = (one-polar(kch))*evert + polar(kch)*ehorz
 311               pemsv(nn) = evert  !because esv=esh
 312               pemsh(nn) = ehorz
 313            else
 314               pemsv(nn) = 0.92_r_kind
 315               pemsh(nn) = 0.92_r_kind
 316               pems5(nn) = 0.92_r_kind
 317            end if
 318         else if(isflg(n) == 0 .or. isflg(n) == 4)then
 319 !      ----- sea (ice-free) MW  -------
 320
 321 !          Open ocean points
 322              call seaem(freqghz,zasat(n),zlsat(n),ts5(n), &
 323                         u10(n),v10(n),ehorz,evert)
 324
 325              pemsv(nn) = evert
 326              pemsh(nn) = ehorz
 327              call ehv2pem( ehorz,evert,zlsat(n),polar(kch), pems5(nn) )
 328
 329         else
 330 !      ----- land (snow-free or snow thin ) MW -------
 331            if(freqghz<80._r_kind)then  !land & f<80GHz
 332 !             currently landem only works for frequencies < 80 Ghz
 333              if (use_landem) then
 334               call landem(zasat(n),freqghz,  &
 335                    soilm5(n),vegf5(n),vegtype5(n),soiltype5(n),soilt5(n), &
 336                    ts5(n),snow5(n),ehorz,evert  )
 337               call ehv2pem( ehorz,evert,zlsat(n),polar(kch), pems5(nn) )
 338               pemsv(nn) = evert
 339               pemsh(nn) = ehorz
 340              else
 341               pemsv(nn) = 0.95_r_kind
 342               pemsh(nn) = 0.95_r_kind
 343               pems5(nn) = 0.95_r_kind
 344              end if
 345            else  !land & f>=80GHz
 346               pemsv(nn) = 0.95_r_kind
 347               pemsh(nn) = 0.95_r_kind
 348               pems5(nn) = 0.95_r_kind
 349            end if  !land & if f><80
 350
 351         end if
 352
 353
 354 !        Load emissivity into array for radiative transfer model
 355 !       (pems5) and diagnostic output file (emissav).
 356
 357         pemsv(nn)       = max(zero,min(pemsv(nn),one))
 358         pemsh(nn)       = max(zero,min(pemsh(nn),one))
 359         pems5(nn)       = max(zero,min(pems5(nn),one))
 360
 361      end do
 362   else
 363     WRITE(UNIT=message(1), FMT='(A,I5)') &
 364       "ILLEGAL surface emissivity type",itype
 365     CALL da_error(__FILE__,__LINE__,message(1:1))
 366   end if   !IR or MW
 367
 368
 369 ! End of routine.
 370   return
 371
 372 #endif
 373
 374
 375 contains
 376
 377 #include "ehv2pem.inc"
 378 #include "adm_ehv2pem.inc"
 379
 380 end subroutine gsi_emiss