da_rttov_ad.inc
References to this file elsewhere.
1 #ifdef RTTOV
2 subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
3 aux_vars, con_vars_ad, tb )
4
5 !---------------------------------------------------------------------------
6 ! PURPOSE: interface to the adjoint subroutine of RTTOV
7 !---------------------------------------------------------------------------
8
9 implicit none
10
11 #include "rttov_ad.interface"
12
13 integer , intent (in) :: inst, nchanl, nprofiles
14 type (con_vars_type), intent (in) :: con_vars (nprofiles)
15 type (con_vars_type), intent (inout) :: con_vars_ad (nprofiles)
16 type (aux_vars_type), intent (in) :: aux_vars (nprofiles)
17 real , intent (in) :: tb(nchanl,nprofiles)
18
19 ! local variables
20 integer :: n, nc
21 Integer :: alloc_status(140)
22
23 ! RTTOV input parameters
24 integer :: nfrequencies, nchannels, nbtout
25 integer :: nchan(nprofiles)
26 integer , pointer :: lprofiles(:)
27 type(rttov_coef) :: coef
28 type(profile_type) :: profiles(nprofiles), profiles_ad(nprofiles)
29 logical :: addcloud, switchrad
30 real, allocatable :: surfem(:)
31 integer , pointer :: channels (:), polarisations(:,:)
32 logical, allocatable :: calcemis (:)
33
34 ! RTTOV out parameters
35 integer :: errorstatus(nprofiles)
36
37 ! RTTOV inout parameters
38 real , pointer :: emissivity (:), emissivity_ad (:)
39 type (radiance_type) :: radiance, radiance_ad
40 type (transmission_type) :: transmission, transmission_ad
41
42
43 call da_trace_entry("da_rttov_ad")
44
45 nchan (:) = nchanl
46 coef = coefs(inst)
47 addcloud = .false.
48 switchrad = .TRUE.
49 alloc_status (:) = 0.0
50
51 do n = 1, nprofiles
52 profiles(n) % nlevels = con_vars(n) % nlevels
53 allocate ( profiles(n)%p(profiles(n) % nlevels), stat=alloc_status(1) )
54 allocate ( profiles(n)%t(profiles(n) % nlevels), stat=alloc_status(2) )
55 allocate ( profiles(n)%q(profiles(n) % nlevels), stat=alloc_status(3) )
56 allocate ( profiles(n)%o3(profiles(n) % nlevels), stat=alloc_status(4) )
57 allocate ( profiles(n)%co2(profiles(n) % nlevels), stat=alloc_status(5) )
58 allocate ( profiles(n)%clw(profiles(n) % nlevels), stat=alloc_status(6) )
59 if ( Any(alloc_status /= 0) ) then
60 WRITE(UNIT=message(1),FMT='(A,I5)') &
61 "ad mem allocation error to profile",n
62 call da_error(__FILE__,__LINE__,message(1:1))
63 End if
64
65 profiles(n) % ozone_data = .false.
66 profiles(n) % co2_data = .false.
67 profiles(n) % clw_data = .false.
68
69 profiles(n) % p(:) = coef%ref_prfl_p(:)
70 profiles(n) % t(:) = con_vars(n)%t(:)
71 profiles(n) % q(:) = con_vars(n)%q(:)
72 profiles(n) % o3(:) = 0.0 !con_vars(n)%o3(:)
73 profiles(n) % co2(:) = 0.0 !con_vars(n)%co2(:)
74 profiles(n) % clw(:) = 0.0 !con_vars(n)%clw(:)
75
76 profiles(n) % skin % surftype = aux_vars (n) % surftype
77 profiles(n) % skin % t = aux_vars (n) % surft
78 profiles(n) % skin % fastem (:) = 0.0 ! aux_vars (n) % fastem (:)
79
80 profiles(n) % s2m % t = aux_vars (n) % t2m
81 profiles(n) % s2m % q = aux_vars (n) % q2m
82 profiles(n) % s2m % o = 0.0 !aux_vars (n) % o3
83 profiles(n) % s2m % p = con_vars (n) % ps
84 profiles(n) % s2m % u = aux_vars (n) % u10
85 profiles(n) % s2m % v = aux_vars (n) % v10
86
87 profiles(n) % zenangle = aux_vars (n) % satzen
88 profiles(n) % azangle = aux_vars (n) % satazi
89
90 profiles(n) % ctp = 500.0
91 profiles(n) % cfraction = 0.0
92
93 profiles_ad(n) % nlevels = con_vars(n) % nlevels
94 allocate ( profiles_ad(n)%p(profiles_ad(n) % nlevels), stat=alloc_status(1) )
95 allocate ( profiles_ad(n)%t(profiles_ad(n) % nlevels), stat=alloc_status(2) )
96 allocate ( profiles_ad(n)%q(profiles_ad(n) % nlevels), stat=alloc_status(3) )
97 allocate ( profiles_ad(n)%o3(profiles_ad(n) % nlevels), stat=alloc_status(4) )
98 allocate ( profiles_ad(n)%co2(profiles_ad(n) % nlevels), stat=alloc_status(5) )
99 allocate ( profiles_ad(n)%clw(profiles_ad(n) % nlevels), stat=alloc_status(6) )
100 if (Any(alloc_status /= 0) ) then
101 WRITE(UNIT=message(1),FMT='(A,I5)') &
102 "mem allocation error to profiles_ad",n
103 call da_error(__FILE__,__LINE__,message(1:1))
104 End if
105
106 profiles_ad(n) % ozone_Data = .False. ! no meaning
107 profiles_ad(n) % co2_Data = .False. ! no meaning
108 profiles_ad(n) % clw_Data = .False. ! no meaning
109 profiles_ad(n) % zenangle = -1 ! no meaning
110 profiles_ad(n) % p(:) = 0.0 ! no AD on pressure levels
111 profiles_ad(n) % t(:) = 0.0 ! temperarure
112 profiles_ad(n) % o3(:) = 0.0 ! O3 ppmv
113 profiles_ad(n) % clw(:) = 0.0 ! clw
114 profiles_ad(n) % q(:) = 0.0 ! WV
115 profiles_ad(n) % s2m % t = 0.0 ! temperarure
116 profiles_ad(n) % s2m % q = 0 ! WV
117 profiles_ad(n) % s2m % p = 0.0 ! pressure
118 profiles_ad(n) % s2m % u = 0.0 ! wind components
119 profiles_ad(n) % s2m % v = 0.0 ! wind components
120 profiles_ad(n) % skin % surftype = -1 ! no meaning
121 profiles_ad(n) % skin % t = 0.0 ! on temperarure
122 profiles_ad(n) % skin % fastem = 0.0 ! Fastem
123 profiles_ad(n) % ctp = 0.0 ! cloud top pressure
124 profiles_ad(n) % cfraction = 0.0 ! cloud fraction
125
126 end do
127
128 #ifdef RTTOV
129 call rttov_setupchan(nprofiles, nchan, coef, & ! in
130 nfrequencies, nchannels, nbtout ) ! out
131 #endif
132
133 allocate (lprofiles(nfrequencies), stat = alloc_status(31) )
134 allocate (channels (nfrequencies), stat = alloc_status(32) )
135 allocate (polarisations(nchannels, 3), stat = alloc_status(33) )
136 allocate (emissivity( nchannels ), stat = alloc_status(34) )
137 allocate (emissivity_ad( nchannels ), stat = alloc_status(134) )
138 allocate (calcemis( nchannels ), stat = alloc_status(35) )
139 allocate (surfem( nchannels ), stat = alloc_status(36) )
140
141 ! allocate transmittance structure
142 allocate (transmission % tau_surf ( nchannels ) ,stat= alloc_status(8))
143 allocate (transmission % tau_layer ( coef % nlevels, nchannels ) ,stat= alloc_status(9))
144 allocate (transmission % od_singlelayer( coef % nlevels, nchannels ),stat= alloc_status(10))
145
146 allocate (transmission_ad % tau_surf ( nchannels ) ,stat= alloc_status(108))
147 allocate (transmission_ad % tau_layer ( coef % nlevels, nchannels ) ,stat= alloc_status(109))
148 allocate (transmission_ad % od_singlelayer( coef % nlevels, nchannels ),stat= alloc_status(110))
149
150
151 ! allocate radiance results arrays with number of channels
152 allocate (radiance % clear ( nchannels ) ,stat= alloc_status(11))
153 allocate (radiance % cloudy ( nchannels ) ,stat= alloc_status(12))
154 allocate (radiance % total ( nchannels ) ,stat= alloc_status(13))
155 allocate (radiance % bt ( nchannels ) ,stat= alloc_status(14))
156 allocate (radiance % bt_clear ( nchannels ) ,stat= alloc_status(15))
157 allocate (radiance % upclear ( nchannels ) ,stat= alloc_status(16))
158 allocate (radiance % dnclear ( nchannels ) ,stat= alloc_status(17))
159 allocate (radiance % reflclear( nchannels ) ,stat= alloc_status(18))
160 allocate (radiance % overcast ( coef % nlevels, nchannels ) ,stat= alloc_status(19))
161 ! allocate the cloudy radiances with full size even
162 ! if not used
163 allocate (radiance % downcld ( coef % nlevels, nchannels ) ,stat= alloc_status(20))
164
165 allocate (radiance % out ( nbtout ) ,stat= alloc_status(121))
166 allocate (radiance % out_clear( nbtout ) ,stat= alloc_status(122))
167 allocate (radiance % total_out( nbtout ) ,stat= alloc_status(123))
168 allocate (radiance % clear_out( nbtout ) ,stat= alloc_status(124))
169
170 ! allocate radiance results arrays with number of channels
171 allocate (radiance_ad % clear ( nchannels ) ,stat= alloc_status(111))
172 allocate (radiance_ad % cloudy ( nchannels ) ,stat= alloc_status(112))
173 allocate (radiance_ad % total ( nchannels ) ,stat= alloc_status(113))
174 allocate (radiance_ad % bt ( nchannels ) ,stat= alloc_status(114))
175 allocate (radiance_ad % bt_clear ( nchannels ) ,stat= alloc_status(115))
176 allocate (radiance_ad % upclear ( nchannels ) ,stat= alloc_status(116))
177 allocate (radiance_ad % dnclear ( nchannels ) ,stat= alloc_status(117))
178 allocate (radiance_ad % reflclear( nchannels ) ,stat= alloc_status(118))
179 allocate (radiance_ad % overcast ( coef % nlevels, nchannels ) ,stat= alloc_status(119))
180 ! allocate the cloudy radiances with full size even
181 ! if not used
182 allocate (radiance_ad % downcld ( coef % nlevels, nchannels ) ,stat= alloc_status(120))
183
184 allocate (radiance_ad % out ( nbtout ) ,stat= alloc_status(121))
185 allocate (radiance_ad % out_clear( nbtout ) ,stat= alloc_status(122))
186 allocate (radiance_ad % total_out( nbtout ) ,stat= alloc_status(123))
187 allocate (radiance_ad % clear_out( nbtout ) ,stat= alloc_status(124))
188
189 if ( Any(alloc_status /= 0) ) then
190 call da_error(__FILE__,__LINE__, &
191 (/"mem allocation error prior to rttov_ad"/))
192 End if
193
194 ! initialise all radiance increments to 0
195 radiance_ad % clear(:) = 0.0
196 radiance_ad % clear_out(:) = 0.0
197 radiance_ad % cloudy(:) = 0.0
198 radiance_ad % total(:) = 0.0
199 radiance_ad % total_out(:) = 0.0
200 radiance_ad % bt(:) = 0.0
201 radiance_ad % bt_clear(:) = 0.0
202 radiance_ad % out(:) = 0.0
203 radiance_ad % out_clear(:) = 0.0
204 radiance_ad % upclear(:) = 0.0
205 radiance_ad % reflclear(:) = 0.0
206 radiance_ad % overcast(:,:) = 0.0
207 radiance_ad % downcld(:,:) = 0.0
208
209 transmission_ad % tau_surf(:) = 0.0
210 transmission_ad % tau_layer(:,:)= 0.0
211 transmission_ad % od_singlelayer(:,:) = 0.0
212
213 nc = nbtout / nprofiles
214 do n = 1, nprofiles
215 radiance_ad % out((n-1)*nc+1:n*nc) = tb(1:nc,n)
216 end do
217
218 surfem (:) = 0.0
219 #ifdef RTTOV
220 call rttov_setupindex(nchan, nprofiles, nfrequencies, & ! in
221 nchannels, nbtout, coef, surfem, & ! in
222 lprofiles, channels, polarisations, & ! out
223 emissivity ) ! out
224 #endif
225
226 nc = nchannels/nprofiles
227
228 if (coef%id_sensor == 1) then ! infrared sensor
229 calcemis (1:nchannels) = .true.
230 emissivity (1:nchannels) = 0.0
231 emissivity_ad (1:nchannels) = 0.0
232 else if (coef%id_sensor == 2) then ! microwave sensor
233 do n = 1, nprofiles
234 if ( profiles(n) % skin % surftype == 1) then ! sea
235 calcemis ((n-1)*nc+1:n*nc) = .true.
236 emissivity ((n-1)*nc+1:n*nc) = 0.0
237 emissivity_ad ((n-1)*nc+1:n*nc) = 0.0
238 else ! 0:land ; 2:sea-ice
239 calcemis ((n-1)*nc+1:n*nc) = .false.
240 emissivity ((n-1)*nc+1:n*nc) = 0.9
241 emissivity_ad ((n-1)*nc+1:n*nc) = 0.0
242 end if
243 end do
244 end if
245
246 #ifdef RTTOV
247 call rttov_ad( &
248 & errorstatus, &! out
249 & nfrequencies, &! in
250 & nchannels, &! in
251 & nbtout, &! in
252 & nprofiles, &! in
253 & channels, &! in
254 & polarisations, &! in
255 & lprofiles, &! in
256 & profiles, &! in
257 & coef, &! in
258 & addcloud, &! in
259 & switchrad, &! in
260 & calcemis, &! in
261 & emissivity, &! inout direct model
262 & profiles_ad, &! inout adjoint
263 & emissivity_ad, &! inout adjoint
264 & transmission, &! inout direct model
265 & transmission_ad, &! inout adjoint input
266 & radiance, &! inout direct model (input due to pointers alloc)
267 & radiance_ad ) ! inout adjoint input (output if converstion Bt -> rad)
268 #endif
269
270 ! rttov87 generates warnings we want to ignore
271 if (any(errorstatus(:) == errorstatus_fatal)) then
272 write (message(1),*) 'rttov_direct error code = ', errorstatus(:)
273 write (message(2),*) 'nfrequencies = ', nfrequencies
274 write (message(3),*) 'nchannels = ', nchannels
275 write (message(4),*) 'nbtout = ', nbtout
276 write (message(5),*) 'nprofiles = ', nprofiles
277 write (message(6),*) 'channels = ', channels
278 write (message(7),*) 'polarisations = ', polarisations
279 write (message(8),*) 'lprofiles = ', lprofiles
280 write (message(9),*) 'addcloud = ', addcloud
281 write (message(10),*) 'switchrad = ', switchrad
282 write (message(11),*) 'calcemis = ', calcemis
283 write (message(12),*) 'profiles%s2m = ', profiles(1)%s2m
284 write (message(13),*) 'profiles%skin = ', profiles(1)%skin
285 write (message(14),*) 'profiles%zenangle = ', profiles(1)%zenangle
286 write (message(15),*) 'profiles%azangle = ', profiles(1)%azangle
287 write (message(16),*) 'profiles%p = ', profiles(1)%p
288 write (message(17),*) 'profiles%t = ', profiles(1)%t
289 write (message(18),*) 'profiles%q = ', profiles(1)%q
290 write (message(19),*) 'profiles%o3 = ', profiles(1)%o3
291 write (message(20),*) 'profiles%co2 = ', profiles(1)%co2
292 write (message(21),*) 'profiles%clw = ', profiles(1)%clw
293 write (message(22),*) 'emissivity = ', emissivity
294 write (message(23),*) 'radiance = ', radiance%out_clear
295 write (message(24),*) 'profiles_ad%s2m = ', profiles_ad(1)%s2m
296 write (message(25),*) 'profiles_ad%skin = ', profiles_ad(1)%skin
297 write (message(26),*) 'profiles_ad%zenangle = ', profiles_ad(1)%zenangle
298 write (message(27),*) 'profiles_ad%azangle = ', profiles_ad(1)%azangle
299 write (message(28),*) 'profiles_ad%p = ', profiles_ad(1)%p
300 write (message(29),*) 'profiles_ad%t = ', profiles_ad(1)%t
301 write (message(30),*) 'profiles_ad%q = ', profiles_ad(1)%q
302 write (message(31),*) 'profiles_ad%o3 = ', profiles_ad(1)%o3
303 write (message(32),*) 'profiles_ad%co2 = ', profiles_ad(1)%co2
304 write (message(33),*) 'profiles_ad%clw = ', profiles_ad(1)%clw
305 write (message(34),*) 'emissivity_ad = ', emissivity_ad
306 write (message(35),*) 'radiance_ad = ', radiance_ad%out_clear
307 call da_warning(__FILE__,__LINE__,message(1:35))
308 end if
309
310 do n = 1, nprofiles
311 con_vars_ad(n)%t(:) = profiles_ad(n) % t(:)
312 con_vars_ad(n)%q(:) = profiles_ad(n) % q(:)
313 con_vars_ad(n)%ps = profiles_ad(n) % s2m % p
314 end do
315
316 deallocate (lprofiles )
317 deallocate (channels )
318 deallocate (polarisations )
319 deallocate (emissivity )
320 deallocate (emissivity_ad )
321 deallocate (calcemis )
322 deallocate (surfem)
323
324 do n = 1, nprofiles
325 deallocate (profiles(n)%p )
326 deallocate (profiles(n)%t )
327 deallocate (profiles(n)%q )
328 deallocate (profiles(n)%o3 )
329 deallocate (profiles(n)%co2 )
330 deallocate (profiles(n)%clw )
331
332 deallocate (profiles_ad(n)%p )
333 deallocate (profiles_ad(n)%t )
334 deallocate (profiles_ad(n)%q )
335 deallocate (profiles_ad(n)%o3 )
336 deallocate (profiles_ad(n)%co2 )
337 deallocate (profiles_ad(n)%clw )
338 end do
339
340 ! deallocate transmittance structure
341 deallocate (transmission % tau_surf ,stat= alloc_status(6))
342 deallocate (transmission % tau_layer ,stat= alloc_status(7))
343 deallocate (transmission % od_singlelayer,stat= alloc_status(8))
344
345 ! deallocate adjoint transmittance structure
346 deallocate (transmission_ad % tau_surf ,stat= alloc_status(106))
347 deallocate (transmission_ad % tau_layer ,stat= alloc_status(107))
348 deallocate (transmission_ad % od_singlelayer,stat= alloc_status(108))
349
350 ! deallocate radiance results arrays with number of channels
351 deallocate (radiance % clear ,stat=alloc_status(9))
352 deallocate (radiance % cloudy ,stat=alloc_status(10))
353 deallocate (radiance % total ,stat=alloc_status(11))
354 deallocate (radiance % bt ,stat=alloc_status(12))
355 deallocate (radiance % bt_clear ,stat=alloc_status(13))
356 deallocate (radiance % upclear ,stat=alloc_status(14))
357 deallocate (radiance % dnclear ,stat=alloc_status(15))
358 deallocate (radiance % reflclear,stat=alloc_status(16))
359 deallocate (radiance % overcast ,stat=alloc_status(17))
360 deallocate (radiance % downcld ,stat=alloc_status(18))
361 deallocate (radiance % out ,stat= alloc_status(19))
362 deallocate (radiance % out_clear ,stat= alloc_status(20))
363 deallocate (radiance % total_out ,stat= alloc_status(21))
364 deallocate (radiance % clear_out ,stat= alloc_status(22))
365
366 deallocate (radiance_ad % clear ,stat=alloc_status(109))
367 deallocate (radiance_ad % cloudy ,stat=alloc_status(110))
368 deallocate (radiance_ad % total ,stat=alloc_status(111))
369 deallocate (radiance_ad % bt ,stat=alloc_status(112))
370 deallocate (radiance_ad % bt_clear ,stat=alloc_status(113))
371 deallocate (radiance_ad % upclear ,stat=alloc_status(114))
372 deallocate (radiance_ad % dnclear ,stat=alloc_status(115))
373 deallocate (radiance_ad % reflclear,stat=alloc_status(116))
374 deallocate (radiance_ad % overcast ,stat=alloc_status(117))
375 deallocate (radiance_ad % downcld ,stat=alloc_status(118))
376 deallocate (radiance_ad % out ,stat= alloc_status(119))
377 deallocate (radiance_ad % out_clear ,stat= alloc_status(120))
378 deallocate (radiance_ad % total_out ,stat= alloc_status(121))
379 deallocate (radiance_ad % clear_out ,stat= alloc_status(122))
380
381 if (Any(alloc_status /= 0) ) then
382 call da_error(__FILE__,__LINE__, &
383 (/"mem deallocation error"/))
384 End if
385
386 call da_trace_exit("da_rttov_ad")
387
388 end subroutine da_rttov_ad
389 #endif
390
391