da_rttov_direct.inc
References to this file elsewhere.
1 #ifdef RTTOV
2 subroutine da_rttov_direct( inst, isflg, nchanl, nprofiles, &
3 ! nfrequencies, nchannels, nbtout , &
4 con_vars, aux_vars, &
5 ! tb, calcemis_out, emisv, emish, emissivity_out )
6 tb, emisv, emish, emissivity_out )
7
8 !---------------------------------------------------------------------------
9 ! PURPOSE: interface to the forward subroutine of RTTOV8_5
10 !---------------------------------------------------------------------------
11
12 implicit none
13
14 #include "rttov_direct.interface"
15
16 integer , intent (in) :: inst, isflg, nchanl, nprofiles
17 ! integer , intent (in) :: nfrequencies, nchannels, nbtout
18 type (con_vars_type), intent (in) :: con_vars (nprofiles)
19 type (aux_vars_type), intent (in) :: aux_vars (nprofiles)
20 ! real , intent (out) :: tb(nprofiles, nchanl)
21 ! real , intent (inout) :: tb(nchanl,nprofiles)
22 real , intent (inout) :: tb(nchanl)
23 ! logical , intent (out) :: calcemis_out(nprofiles*nchanl)
24 real , intent (in) :: emisv(nprofiles*nchanl)
25 real , intent (in) :: emish(nprofiles*nchanl)
26 real , intent (out) :: emissivity_out(:)
27
28 ! local variables
29 integer :: nfrequencies, nchannels, nbtout
30 integer :: n, nc, ios, ich, i, ichannel, pol_id
31 Integer :: alloc_status(40), emis_scheme
32
33 ! RTTOV input parameters
34 ! integer :: nfrequencies, nchannels, nbtout
35 integer :: nchan(nprofiles)
36 integer , pointer :: lprofiles(:)
37 type(rttov_coef) :: coef
38 type(profile_type) :: profiles(nprofiles)
39 logical :: addcloud
40 real , pointer :: surfem(:)
41 integer , pointer :: channels (:), polarisations(:,:)
42 logical , pointer :: calcemis (:)
43
44 ! RTTOV out parameters
45 integer :: errorstatus(nprofiles)
46
47 ! RTTOV inout parameters
48 real , pointer :: emissivity (:)
49 type (radiance_type) :: radiance
50 type (transmission_type) :: transmission
51
52 call da_trace_entry("da_rttov_direct")
53
54 nchan (:) = nchanl
55 coef = coefs(inst)
56 addcloud = .false.
57 alloc_status (:) = 0
58
59 do n = 1, nprofiles
60 profiles(n) % nlevels = con_vars(n) % nlevels
61 allocate ( profiles(n)%p(profiles(n) % nlevels), stat=alloc_status(1) )
62 allocate ( profiles(n)%t(profiles(n) % nlevels), stat=alloc_status(2) )
63 allocate ( profiles(n)%q(profiles(n) % nlevels), stat=alloc_status(3) )
64 allocate ( profiles(n)%o3(profiles(n) % nlevels), stat=alloc_status(4) )
65 allocate ( profiles(n)%co2(profiles(n) % nlevels), stat=alloc_status(5) )
66 allocate ( profiles(n)%clw(profiles(n) % nlevels), stat=alloc_status(6) )
67
68 If (Any(alloc_status /= 0) ) Then
69 WRITE(UNIT=message(1),FMT='(A,I5)') &
70 "mem allocation error to profile",n
71 call da_error(__FILE__,__LINE__,message(1:1))
72 End If
73
74 profiles(n) % ozone_data = .false.
75 profiles(n) % co2_data = .false.
76 profiles(n) % clw_data = .false.
77
78 profiles(n) % p(:) = coef%ref_prfl_p(:)
79 profiles(n) % t(:) = con_vars(n)%t(:)
80 profiles(n) % q(:) = con_vars(n)%q(:)
81 profiles(n) % o3(:) = 0.0 !con_vars(n)%o3(:)
82 profiles(n) % co2(:) = 0.0 !con_vars(n)%co2(:)
83 profiles(n) % clw(:) = 0.0 !con_vars(n)%clw(:)
84
85 if (isflg == 0 .or. isflg == 4) then ! sea
86 profiles(n) % skin % surftype = 1 ! aux_vars (n) % surftype
87 ! profiles(n) % skin % fastem (:) = 0.
88 else if (isflg == 1 .or. isflg == 5) then ! sea-ice with snow
89 profiles(n) % skin % surftype = 2
90 ! profiles(n) % skin % fastem (1) = 2.2
91 ! profiles(n) % skin % fastem (2) = 3.7
92 ! profiles(n) % skin % fastem (3) = 122.0
93 ! profiles(n) % skin % fastem (4) = 0.0
94 ! profiles(n) % skin % fastem (5) = 0.15
95 else
96 profiles(n) % skin % surftype = 0 ! land (Deep dry snow)
97 ! profiles(n) % skin % fastem (1) = 3.0
98 ! profiles(n) % skin % fastem (2) = 24.0
99 ! profiles(n) % skin % fastem (3) = 60.0
100 ! profiles(n) % skin % fastem (4) = 0.1
101 ! profiles(n) % skin % fastem (5) = 0.15
102 end if
103 ! profiles(n) % skin % surftype = aux_vars (n) % surftype
104 profiles(n) % skin % t = aux_vars (n) % surft
105 profiles(n) % skin % fastem (:) = 0. ! aux_vars (n) % fastem (:)
106
107 profiles(n) % s2m % t = aux_vars (n) % t2m
108 profiles(n) % s2m % q = aux_vars (n) % q2m
109 profiles(n) % s2m % o = 0.0 !aux_vars (n) % o3
110 profiles(n) % s2m % p = con_vars (n) % ps
111 profiles(n) % s2m % u = aux_vars (n) % u10
112 profiles(n) % s2m % v = aux_vars (n) % v10
113
114 profiles(n) % zenangle = aux_vars (n) % satzen
115 profiles(n) % azangle = aux_vars (n) % satazi
116
117 profiles(n) % ctp = 500.
118 profiles(n) % cfraction = 0.
119 end do
120
121 #ifdef RTTOV
122 call rttov_setupchan(nprofiles, nchan, coef, & ! in
123 nfrequencies, nchannels, nbtout ) ! out
124 #endif
125
126 Allocate ( lprofiles(nfrequencies), stat = alloc_status(31) )
127 Allocate ( channels (nfrequencies), stat = alloc_status(32) )
128 Allocate ( polarisations(nchannels, 3), stat = alloc_status(33) )
129 Allocate ( emissivity( nchannels ), stat = alloc_status(34) )
130 Allocate ( calcemis( nchannels ), stat = alloc_status(35) )
131 Allocate ( surfem( nchannels ), stat = alloc_status(36) )
132
133 emissivity = 0.0
134
135 ! allocate transmittance structure
136 Allocate( transmission % tau_surf ( nchannels ) ,stat= alloc_status(8))
137 Allocate( transmission % tau_layer ( coef % nlevels, nchannels ) ,stat= alloc_status(9))
138 Allocate( transmission % od_singlelayer( coef % nlevels, nchannels ),stat= alloc_status(10))
139
140 transmission % tau_surf = 0.0
141 transmission % tau_layer = 0.0
142 transmission % od_singlelayer = 0.0
143
144 ! allocate radiance results arrays with number of channels
145 Allocate( radiance % clear ( nchannels ) ,stat= alloc_status(11))
146 Allocate( radiance % cloudy ( nchannels ) ,stat= alloc_status(12))
147 Allocate( radiance % total ( nchannels ) ,stat= alloc_status(13))
148 Allocate( radiance % bt ( nchannels ) ,stat= alloc_status(14))
149 Allocate( radiance % bt_clear ( nchannels ) ,stat= alloc_status(15))
150 Allocate( radiance % upclear ( nchannels ) ,stat= alloc_status(16))
151 Allocate( radiance % dnclear ( nchannels ) ,stat= alloc_status(17))
152 Allocate( radiance % reflclear( nchannels ) ,stat= alloc_status(18))
153 Allocate( radiance % overcast ( coef % nlevels, nchannels ) ,stat= alloc_status(19))
154 ! allocate the cloudy radiances with full size even
155 ! if not used
156 Allocate( radiance % downcld ( coef % nlevels, nchannels ) ,stat= alloc_status(20))
157
158 Allocate( radiance % out ( nbtout ) ,stat= alloc_status(21))
159 Allocate( radiance % out_clear( nbtout ) ,stat= alloc_status(22))
160 Allocate( radiance % total_out( nbtout ) ,stat= alloc_status(23))
161 Allocate( radiance % clear_out( nbtout ) ,stat= alloc_status(24))
162
163 radiance % clear = 0.0
164 radiance % cloudy = 0.0
165 radiance % total = 0.0
166 radiance % bt = 0.0
167 radiance % bt_clear = 0.0
168 radiance % upclear = 0.0
169 radiance % dnclear = 0.0
170 radiance % reflclear = 0.0
171 radiance % overcast = 0.0
172 radiance % downcld = 0.0
173 radiance % out = 0.0
174 radiance % out_clear = 0.0
175 radiance % total_out = 0.0
176 radiance % clear_out = 0.0
177
178 If (Any(alloc_status /= 0)) Then
179 call da_error(__FILE__,__LINE__, &
180 (/"mem allocation error prior to rttov_direct"/))
181 End If
182
183 surfem (:) = 0.
184 #ifdef RTTOV
185 call rttov_setupindex(nchan, nprofiles, nfrequencies, & ! in
186 nchannels, nbtout, coef, surfem, & ! in
187 lprofiles, channels, polarisations, & ! out
188 emissivity ) ! out
189 #endif
190
191 ! surface emissivity scheme
192 !----------------------------------------------------------
193
194 ! For Infrared sensors
195 !-----------------------------
196 if (coef%id_sensor == 1 .or. coef%id_sensor == 3) then
197 if (profiles(1) % skin % surftype == 1) then
198 calcemis (1:nchannels) = .true. ! using ISSEM over sea
199 emissivity (1:nchannels) = 0.
200 else if (profiles(1) % skin % surftype == 2) then
201 calcemis (1:nchannels) = .false.
202 emissivity (1:nchannels) = 0.98 ! over sea-ice
203 else
204 if (isflg == 2 .or. isflg == 6) then
205 calcemis (1:nchannels) = .false.
206 emissivity (1:nchannels) = 0.97 ! land without snow
207 end if
208 if (isflg == 3 .or. isflg == 7) then
209 calcemis (1:nchannels) = .false.
210 emissivity (1:nchannels) = 1.0 ! land with snow
211 end if
212 end if
213 end if
214
215 !write (0,*) __FILE__,__LINE__,"sum(emissivity)",sum(emissivity)
216 !write (0,*) __FILE__,__LINE__,"profiles(1) % skin % surftype",profiles(1) % skin % surftype
217 !write (0,*) __FILE__,__LINE__,"mw_emis_sea,nchanl",mw_emis_sea,nchanl
218 !write (0,*) __FILE__,__LINE__,"emisv",emisv
219 !write (0,*) __FILE__,__LINE__,"emish",emish
220
221 ! For Microwave sensors
222 !-----------------------------
223 if (coef%id_sensor == 2) then
224 !
225 ! 1.0 over sea
226 !
227 ! print *, 'mw_emis_sea=', mw_emis_sea
228 if ( profiles(1) % skin % surftype == 1) then
229 if (mw_emis_sea == 0) then
230 calcemis (1:nchannels) = .true.
231 emissivity(1:nchannels) = -1. ! RTTOV fastem-2
232 else ! Weng's model
233 calcemis (1:nchannels) = .false.
234 do i = 1, nchanl ! loop for channels
235 ichannel = polarisations(i,1) ! position of first polar at chanl i
236 pol_id = coef%fastem_polar(i) + 1
237 if (polarisations(i,3) == 2) then ! number of polar at chanl i
238 emissivity(ichannel) = emisv(i)
239 emissivity(ichannel+1) = emish(i)
240 else if (polarisations(i,3) == 4) then
241 emissivity(ichannel) = emisv(i)
242 emissivity(ichannel+1) = emish(i)
243 emissivity(ichannel+2) = 0.
244 emissivity(ichannel+3) = 0.
245 else ! == 1 V or H polar
246 if (pol_id == 4) emissivity(ichannel) = emisv(i)
247 if (pol_id == 5) emissivity(ichannel) = emish(i)
248 end if
249 end do
250 end if
251
252 ! 1.0 over sea-ice/snow/land
253 else
254 calcemis (1:nchannels) = .false.
255 do i = 1, nchanl ! loop for channels
256 ichannel = polarisations(i,1) ! position of first polar at chanl i
257 pol_id = coef%fastem_polar(i) + 1
258 if (polarisations(i,3) == 2) then ! number of polar at chanl i
259 emissivity(ichannel) = emisv(i)
260 emissivity(ichannel+1) = emish(i)
261 else if (polarisations(i,3) == 4) then
262 emissivity(ichannel) = emisv(i)
263 emissivity(ichannel+1) = emish(i)
264 emissivity(ichannel+2) = 0.
265 emissivity(ichannel+3) = 0.
266 else ! == 1 V or H polar
267 if (pol_id == 4) emissivity(ichannel) = emisv(i)
268 if (pol_id == 5) emissivity(ichannel) = emish(i)
269 end if
270 end do
271 end if
272 end if
273
274 ! write (0,*) "sum(emissivity)",sum(emissivity)
275 ! write (0,*) "emissivity)",emissivity
276 ! call da_error(__FILE__,__LINE__,(/"bugger"/))
277
278 !-----------------------------------
279 ! calling RTTOV forward model
280 !----------------------------------
281 goto 37
282 write (0,*) __FILE__,__LINE__
283 write (0,*) __FILE__,__LINE__,"sum(channels)",sum(channels)
284 write (0,*) __FILE__,__LINE__,"sum(polarisations(:,1))",sum(polarisations(:,1))
285 write (0,*) __FILE__,__LINE__,"sum(polarisations(:,2))",sum(polarisations(:,2))
286 write (0,*) __FILE__,__LINE__,"sum(polarisations(:,3))",sum(polarisations(:,3))
287 !write (0,*) __FILE__,__LINE__,"sum(polarisations((1:nc,1))",sum(polarisations(1:nc,1))
288 !write (0,*) __FILE__,__LINE__,"sum(polarisations((1:nc,2))",sum(polarisations(1:nc,2))
289 !write (0,*) __FILE__,__LINE__,"sum(polarisations((1:nc,3))",sum(polarisations(1:nc,3))
290 write (0,*) __FILE__,__LINE__,"sum(emissivity)",sum(emissivity)
291 write (0,*) __FILE__,__LINE__,"sum(transmission % tau_surf)",sum(transmission % tau_surf)
292 write (0,*) __FILE__,__LINE__,"sum(transmission % tau_layer)",sum(transmission % tau_layer)
293 write (0,*) __FILE__,__LINE__,"sum(transmission % od_singlelayer)",sum(transmission %od_singlelayer)
294 write (0,*) __FILE__,__LINE__,"sum(radiance % clear)",sum(radiance % clear)
295 write (0,*) __FILE__,__LINE__,"sum(radiance % cloudy)",sum(radiance % cloudy)
296 write (0,*) __FILE__,__LINE__,"sum(radiance % total)",sum(radiance % total)
297 write (0,*) __FILE__,__LINE__,"sum(radiance % bt)",sum(radiance % bt)
298 write (0,*) __FILE__,__LINE__,"sum(radiance % bt_clear)",sum(radiance % bt_clear)
299 write (0,*) __FILE__,__LINE__,"sum(radiance % upclear)",sum(radiance % upclear)
300 write (0,*) __FILE__,__LINE__,"sum(radiance % dnclear)",sum(radiance % dnclear)
301 write (0,*) __FILE__,__LINE__,"sum(radiance % reflclear)",sum(radiance % reflclear)
302 write (0,*) __FILE__,__LINE__,"sum(radiance % overcast)",sum(radiance % overcast)
303 write (0,*) __FILE__,__LINE__,"sum(radiance % downcld)",sum(radiance % downcld)
304 write (0,*) __FILE__,__LINE__,"sum(radiance % out)",sum(radiance % out)
305 write (0,*) __FILE__,__LINE__,"sum(radiance % out_clear)",sum(radiance % out_clear)
306 write (0,*) __FILE__,__LINE__,"sum(radiance % total_out)",sum(radiance % total_out)
307 write (0,*) __FILE__,__LINE__,"sum(radiance % clear_out)",sum(radiance % clear_out)
308 write (0,*) __FILE__,__LINE__,"sum(profiles(1) % p(:))", sum(profiles(1) % p(:))
309 write (0,*) __FILE__,__LINE__,"sum(profiles(1) % t(:)) ", sum(profiles(1) % t(:))
310 write (0,*) __FILE__,__LINE__,"sum(profiles(1) % q(:)) ", sum(profiles(1) % q(:))
311 write (0,*) __FILE__,__LINE__,"sum(profiles(1) % o3(:))", sum(profiles(1) % o3(:))
312 write (0,*) __FILE__,__LINE__,"sum(profiles(1) % co2(:))", sum(profiles(1) % co2(:))
313 write (0,*) __FILE__,__LINE__,"sum(profiles(1) % clw(:))", sum(profiles(1) % clw(:))
314 write (0,*) __FILE__,__LINE__,"profiles(1) % skin % surftype",profiles(1) % skin % surftype
315 write (0,*) __FILE__,__LINE__,"profiles(1) % skin % t",profiles(1) % skin % t
316 write (0,*) __FILE__,__LINE__,"profiles(1) % skin % fastem)", sum(profiles(1) % skin % fastem)
317 write (0,*) __FILE__,__LINE__,"profiles(1) % s2m % t)",profiles(1) % s2m % t
318 write (0,*) __FILE__,__LINE__,"profiles(1) % s2m % q)",profiles(1) % s2m % q
319 write (0,*) __FILE__,__LINE__,"profiles(1) % s2m % o)",profiles(1) % s2m % o
320 write (0,*) __FILE__,__LINE__,"profiles(1) % s2m % p)",profiles(1) % s2m % p
321 write (0,*) __FILE__,__LINE__,"profiles(1) % s2m % u)",profiles(1) % s2m % u
322 write (0,*) __FILE__,__LINE__,"profiles(1) % s2m % v)",profiles(1) % s2m % v
323 write (0,*) __FILE__,__LINE__,"profiles(1) % zenangle)",profiles(1) % zenangle
324 write (0,*) __FILE__,__LINE__,"profiles(1) % azangle)",profiles(1) % azangle
325 write (0,*) __FILE__,__LINE__,"profiles(1) % ctp)",profiles(1) % ctp
326 write (0,*) __FILE__,__LINE__,"profiles(1) % cfraction)",profiles(1) % cfraction
327 37 continue
328 #ifdef RTTOV
329 call rttov_direct( &
330 & errorstatus, &! out
331 & nfrequencies, &! in
332 & nchannels, &! in
333 & nbtout, &! in
334 & nprofiles, &! in
335 & channels, &! in
336 & polarisations, &! in
337 & lprofiles, &! in
338 & profiles, &! in
339 & coef, &! in
340 & addcloud, &! in
341 & calcemis, &! in
342 & emissivity, &! inout
343 & transmission, &! inout
344 & radiance ) ! inout
345 #endif
346
347 if ( any(errorstatus(:) /= 0 )) then
348 !WRITE (UNIT=stderr,FMT=*) 'rttov_direct error code = ', errorstatus(:)
349 !WRITE (UNIT=stderr,FMT=*) 'nfrequencies = ', nfrequencies
350 !WRITE (UNIT=stderr,FMT=*) 'nchannels = ', nchannels
351 !WRITE (UNIT=stderr,FMT=*) 'nbtout = ', nbtout
352 !WRITE (UNIT=stderr,FMT=*) 'nprofiles = ', nprofiles
353 !WRITE (UNIT=stderr,FMT=*) 'channels = ', channels
354 !WRITE (UNIT=stderr,FMT=*) 'polarisations= ', polarisations
355 !WRITE (UNIT=stderr,FMT=*) 'lprofiles = ', lprofiles
356 !WRITE (UNIT=stderr,FMT=*) 'profiles%s2m = ', profiles(1)%s2m
357 !WRITE (UNIT=stderr,FMT=*) 'profiles%skin = ', profiles(1)%skin
358 write(UNIT=stderr,FMT='(6f12.2)') &
359 profiles(1)%s2m%p,profiles(1)%skin%t,profiles(1)%s2m%t, &
360 profiles(1)%s2m%q,profiles(1)%s2m%u,profiles(1)%s2m%v
361 !WRITE (UNIT=stderr,FMT=*) 'profiles%zenangle = ', profiles(1)%zenangle
362 !WRITE (UNIT=stderr,FMT=*) 'profiles%azangle = ', profiles(1)%azangle
363 !WRITE (UNIT=stderr,FMT=*) 'profiles%p = ', profiles(1)%p
364 !WRITE (UNIT=stderr,FMT=*) 'profiles%t = ', profiles(1)%t
365 !WRITE (UNIT=stderr,FMT=*) 'profiles%q = ', profiles(1)%q
366 !do i=coef%nlevels,1,-1
367 ! write(UNIT=stderr,FMT='(i4,3f12.2)') i, profiles(1)%p(i),profiles(1)%t(i),profiles(1)%q(i)
368 !end do
369 !WRITE (UNIT=stderr,FMT=*) 'addcloud = ', addcloud
370 !WRITE (UNIT=stderr,FMT=*) 'calcemis = ', calcemis
371 !WRITE (UNIT=stderr,FMT=*) 'emissivity = ', emissivity
372 !WRITE (UNIT=stderr,FMT=*) 'radiance = ', radiance%out_clear
373 !call da_error(__FILE__,__LINE__,(/"Erk"/))
374 endif
375
376 nc = nbtout / nprofiles
377 do n = 1, nprofiles
378 tb(1:nc) = radiance % out_clear((n-1)*nc+1:n*nc)
379 end do
380 ! calcemis_out (:) = calcemis(:)
381 emissivity_out(:) = emissivity(:)
382
383 deallocate ( lprofiles )
384 deallocate ( channels )
385 deallocate ( polarisations )
386 deallocate ( emissivity )
387 deallocate ( calcemis )
388 deallocate ( surfem )
389 do n = 1, nprofiles
390 deallocate ( profiles(n)%p )
391 deallocate ( profiles(n)%t )
392 deallocate ( profiles(n)%q )
393 deallocate ( profiles(n)%o3 )
394 deallocate ( profiles(n)%co2 )
395 deallocate ( profiles(n)%clw )
396 end do
397
398 ! deallocate transmittance structure
399 Deallocate( transmission % tau_surf ,stat= alloc_status(6))
400 Deallocate( transmission % tau_layer ,stat= alloc_status(7))
401 Deallocate( transmission % od_singlelayer,stat= alloc_status(8))
402
403 ! deallocate radiance results arrays with number of channels
404 Deallocate( radiance % clear ,stat=alloc_status(9))
405 Deallocate( radiance % cloudy ,stat=alloc_status(10))
406 Deallocate( radiance % total ,stat=alloc_status(11))
407 Deallocate( radiance % bt ,stat=alloc_status(12))
408 Deallocate( radiance % bt_clear ,stat=alloc_status(13))
409 Deallocate( radiance % upclear ,stat=alloc_status(14))
410 Deallocate( radiance % dnclear ,stat=alloc_status(15))
411 Deallocate( radiance % reflclear,stat=alloc_status(16))
412 Deallocate( radiance % overcast ,stat=alloc_status(17))
413 Deallocate( radiance % downcld ,stat=alloc_status(18))
414 Deallocate( radiance % out ,stat= alloc_status(19))
415 Deallocate( radiance % out_clear ,stat= alloc_status(20))
416 Deallocate( radiance % total_out ,stat= alloc_status(21))
417 Deallocate( radiance % clear_out ,stat= alloc_status(22))
418
419 If( Any(alloc_status /= 0) ) Then
420 call da_error(__FILE__,__LINE__, &
421 (/"mem deallocation error"/))
422 End If
423
424 call da_trace_exit("da_rttov_direct")
425
426 end subroutine da_rttov_direct
427 #endif