gsi_emiss.inc
References to this file elsewhere.
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