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 RTTOV
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 ! For Microwave sensors
216 !-----------------------------
217 if (coef%id_sensor == 2) then
218 !
219 ! 1.0 over sea
220 !
221 ! print *, 'mw_emis_sea=', mw_emis_sea
222 if ( profiles(1) % skin % surftype == 1) then
223 if (mw_emis_sea == 0) then
224 calcemis (1:nchannels) = .true.
225 emissivity(1:nchannels) = -1. ! RTTOV fastem-2
226 else ! Weng's model
227 calcemis (1:nchannels) = .false.
228 do i = 1, nchanl ! loop for channels
229 ichannel = polarisations(i,1) ! position of first polar at chanl i
230 pol_id = coef%fastem_polar(i) + 1
231 if (polarisations(i,3) == 2) then ! number of polar at chanl i
232 emissivity(ichannel) = emisv(i)
233 emissivity(ichannel+1) = emish(i)
234 else if (polarisations(i,3) == 4) then
235 emissivity(ichannel) = emisv(i)
236 emissivity(ichannel+1) = emish(i)
237 emissivity(ichannel+2) = 0.
238 emissivity(ichannel+3) = 0.
239 else ! == 1 V or H polar
240 if (pol_id == 4) emissivity(ichannel) = emisv(i)
241 if (pol_id == 5) emissivity(ichannel) = emish(i)
242 end if
243 end do
244 end if
245
246 ! 1.0 over sea-ice/snow/land
247 else
248 calcemis (1:nchannels) = .false.
249 do i = 1, nchanl ! loop for channels
250 ichannel = polarisations(i,1) ! position of first polar at chanl i
251 pol_id = coef%fastem_polar(i) + 1
252 if (polarisations(i,3) == 2) then ! number of polar at chanl i
253 emissivity(ichannel) = emisv(i)
254 emissivity(ichannel+1) = emish(i)
255 else if (polarisations(i,3) == 4) then
256 emissivity(ichannel) = emisv(i)
257 emissivity(ichannel+1) = emish(i)
258 emissivity(ichannel+2) = 0.
259 emissivity(ichannel+3) = 0.
260 else ! == 1 V or H polar
261 if (pol_id == 4) emissivity(ichannel) = emisv(i)
262 if (pol_id == 5) emissivity(ichannel) = emish(i)
263 end if
264 end do
265 end if
266 end if
267
268 !-----------------------------------
269 ! calling RTTOV forward model
270 !----------------------------------
271
272 #ifdef RTTOV
273 call rttov_direct( &
274 & errorstatus, &! out
275 & nfrequencies, &! in
276 & nchannels, &! in
277 & nbtout, &! in
278 & nprofiles, &! in
279 & channels, &! in
280 & polarisations, &! in
281 & lprofiles, &! in
282 & profiles, &! in
283 & coef, &! in
284 & addcloud, &! in
285 & calcemis, &! in
286 & emissivity, &! inout
287 & transmission, &! inout
288 & radiance ) ! inout
289 #endif
290
291 ! rttov87 generates warnings we want to ignore
292 if ( any(errorstatus(:) == errorstatus_fatal )) then
293 WRITE (UNIT=stderr,FMT=*) 'rttov_direct error code = ', errorstatus(:)
294 WRITE (UNIT=stderr,FMT=*) 'nfrequencies = ', nfrequencies
295 WRITE (UNIT=stderr,FMT=*) 'nchannels = ', nchannels
296 WRITE (UNIT=stderr,FMT=*) 'nbtout = ', nbtout
297 WRITE (UNIT=stderr,FMT=*) 'nprofiles = ', nprofiles
298 WRITE (UNIT=stderr,FMT=*) 'channels = ', channels
299 WRITE (UNIT=stderr,FMT=*) 'polarisations= ', polarisations
300 WRITE (UNIT=stderr,FMT=*) 'lprofiles = ', lprofiles
301 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%t = ', profiles(1)%s2m%t
302 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%q = ', profiles(1)%s2m%q
303 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%o = ', profiles(1)%s2m%o
304 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%p = ', profiles(1)%s2m%p
305 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%u = ', profiles(1)%s2m%u
306 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%v = ', profiles(1)%s2m%v
307 WRITE (UNIT=stderr,FMT=*) 'profiles%skin%surftype = ', profiles(1)%skin%surftype
308 WRITE (UNIT=stderr,FMT=*) 'profiles%skin%t = ', profiles(1)%skin%t
309 WRITE (UNIT=stderr,FMT=*) 'profiles%skin%fastem = ', profiles(1)%skin%fastem
310 WRITE (UNIT=stderr,FMT=*) 'profiles%zenangle = ', profiles(1)%zenangle
311 WRITE (UNIT=stderr,FMT=*) 'profiles%azangle = ', profiles(1)%azangle
312 WRITE (UNIT=stderr,FMT=*) 'profiles%p = ', profiles(1)%p
313 WRITE (UNIT=stderr,FMT=*) 'profiles%t = ', profiles(1)%t
314 WRITE (UNIT=stderr,FMT=*) 'profiles%q = ', profiles(1)%q
315 do i=coef%nlevels,1,-1
316 write(UNIT=stderr,FMT='(i4,3f12.2)') i, profiles(1)%p(i),profiles(1)%t(i),profiles(1)%q(i)
317 end do
318 WRITE (UNIT=stderr,FMT=*) 'addcloud = ', addcloud
319 WRITE (UNIT=stderr,FMT=*) 'calcemis = ', calcemis
320 WRITE (UNIT=stderr,FMT=*) 'emissivity = ', emissivity
321 WRITE (UNIT=stderr,FMT=*) 'radiance = ', radiance%out_clear
322 call da_warning(__FILE__,__LINE__,(/"Problem in rttov_direct"/))
323 endif
324
325 nc = nbtout / nprofiles
326 do n = 1, nprofiles
327 tb(1:nc) = radiance % out_clear((n-1)*nc+1:n*nc)
328 end do
329 ! calcemis_out (:) = calcemis(:)
330 emissivity_out(:) = emissivity(:)
331
332 deallocate ( lprofiles )
333 deallocate ( channels )
334 deallocate ( polarisations )
335 deallocate ( emissivity )
336 deallocate ( calcemis )
337 deallocate ( surfem )
338 do n = 1, nprofiles
339 deallocate ( profiles(n)%p )
340 deallocate ( profiles(n)%t )
341 deallocate ( profiles(n)%q )
342 deallocate ( profiles(n)%o3 )
343 deallocate ( profiles(n)%co2 )
344 deallocate ( profiles(n)%clw )
345 end do
346
347 ! deallocate transmittance structure
348 Deallocate( transmission % tau_surf ,stat= alloc_status(6))
349 Deallocate( transmission % tau_layer ,stat= alloc_status(7))
350 Deallocate( transmission % od_singlelayer,stat= alloc_status(8))
351
352 ! deallocate radiance results arrays with number of channels
353 Deallocate( radiance % clear ,stat=alloc_status(9))
354 Deallocate( radiance % cloudy ,stat=alloc_status(10))
355 Deallocate( radiance % total ,stat=alloc_status(11))
356 Deallocate( radiance % bt ,stat=alloc_status(12))
357 Deallocate( radiance % bt_clear ,stat=alloc_status(13))
358 Deallocate( radiance % upclear ,stat=alloc_status(14))
359 Deallocate( radiance % dnclear ,stat=alloc_status(15))
360 Deallocate( radiance % reflclear,stat=alloc_status(16))
361 Deallocate( radiance % overcast ,stat=alloc_status(17))
362 Deallocate( radiance % downcld ,stat=alloc_status(18))
363 Deallocate( radiance % out ,stat= alloc_status(19))
364 Deallocate( radiance % out_clear ,stat= alloc_status(20))
365 Deallocate( radiance % total_out ,stat= alloc_status(21))
366 Deallocate( radiance % clear_out ,stat= alloc_status(22))
367
368 If( Any(alloc_status /= 0) ) Then
369 call da_error(__FILE__,__LINE__, &
370 (/"mem deallocation error"/))
371 End If
372
373 call da_trace_exit("da_rttov_direct")
374
375 end subroutine da_rttov_direct
376 #endif