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