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, i, ichannel, pol_id
31 Integer :: alloc_status(40)
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, allocatable :: surfem(:)
41 integer , pointer :: channels (:), polarisations(:,:)
42 logical, allocatable :: 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.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.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.0
118 profiles(n) % cfraction = 0.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.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.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 if (profiles(1) % skin % surftype == 1) then
222 if (mw_emis_sea == 0) then
223 calcemis (1:nchannels) = .true.
224 emissivity(1:nchannels) = -1.0 ! RTTOV fastem-2
225 else ! Weng's model
226 calcemis (1:nchannels) = .false.
227 do i = 1, nchanl ! loop for channels
228 ichannel = polarisations(i,1) ! position of first polar at chanl i
229 pol_id = coef%fastem_polar(i) + 1
230 if (polarisations(i,3) == 2) then ! number of polar at chanl i
231 emissivity(ichannel) = emisv(i)
232 emissivity(ichannel+1) = emish(i)
233 else if (polarisations(i,3) == 4) then
234 emissivity(ichannel) = emisv(i)
235 emissivity(ichannel+1) = emish(i)
236 emissivity(ichannel+2) = 0.0
237 emissivity(ichannel+3) = 0.0
238 else ! == 1 V or H polar
239 if (pol_id == 4) emissivity(ichannel) = emisv(i)
240 if (pol_id == 5) emissivity(ichannel) = emish(i)
241 end if
242 end do
243 end if
244
245 ! 1.0 over sea-ice/snow/land
246 else
247 calcemis (1:nchannels) = .false.
248 do i = 1, nchanl ! loop for channels
249 ichannel = polarisations(i,1) ! position of first polar at chanl i
250 pol_id = coef%fastem_polar(i) + 1
251 if (polarisations(i,3) == 2) then ! number of polar at chanl i
252 emissivity(ichannel) = emisv(i)
253 emissivity(ichannel+1) = emish(i)
254 else if (polarisations(i,3) == 4) then
255 emissivity(ichannel) = emisv(i)
256 emissivity(ichannel+1) = emish(i)
257 emissivity(ichannel+2) = 0.0
258 emissivity(ichannel+3) = 0.0
259 else ! == 1 V or H polar
260 if (pol_id == 4) emissivity(ichannel) = emisv(i)
261 if (pol_id == 5) emissivity(ichannel) = emish(i)
262 end if
263 end do
264 end if
265 end if
266
267 !-----------------------------------
268 ! calling RTTOV forward model
269 !----------------------------------
270
271 #ifdef RTTOV
272 call rttov_direct(&
273 & errorstatus, &! out
274 & nfrequencies, &! in
275 & nchannels, &! in
276 & nbtout, &! in
277 & nprofiles, &! in
278 & channels, &! in
279 & polarisations, &! in
280 & lprofiles, &! in
281 & profiles, &! in
282 & coef, &! in
283 & addcloud, &! in
284 & calcemis, &! in
285 & emissivity, &! inout
286 & transmission, &! inout
287 & radiance) ! inout
288 #endif
289
290 ! rttov87 generates warnings we want to ignore
291 if (any(errorstatus(:) == errorstatus_fatal)) then
292 WRITE (UNIT=stderr,FMT=*) 'rttov_direct error code = ', errorstatus(:)
293 WRITE (UNIT=stderr,FMT=*) 'nfrequencies = ', nfrequencies
294 WRITE (UNIT=stderr,FMT=*) 'nchannels = ', nchannels
295 WRITE (UNIT=stderr,FMT=*) 'nbtout = ', nbtout
296 WRITE (UNIT=stderr,FMT=*) 'nprofiles = ', nprofiles
297 WRITE (UNIT=stderr,FMT=*) 'channels = ', channels
298 WRITE (UNIT=stderr,FMT=*) 'polarisations= ', polarisations
299 WRITE (UNIT=stderr,FMT=*) 'lprofiles = ', lprofiles
300 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%t = ', profiles(1)%s2m%t
301 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%q = ', profiles(1)%s2m%q
302 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%o = ', profiles(1)%s2m%o
303 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%p = ', profiles(1)%s2m%p
304 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%u = ', profiles(1)%s2m%u
305 WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%v = ', profiles(1)%s2m%v
306 WRITE (UNIT=stderr,FMT=*) 'profiles%skin%surftype = ', profiles(1)%skin%surftype
307 WRITE (UNIT=stderr,FMT=*) 'profiles%skin%t = ', profiles(1)%skin%t
308 WRITE (UNIT=stderr,FMT=*) 'profiles%skin%fastem = ', profiles(1)%skin%fastem
309 WRITE (UNIT=stderr,FMT=*) 'profiles%zenangle = ', profiles(1)%zenangle
310 WRITE (UNIT=stderr,FMT=*) 'profiles%azangle = ', profiles(1)%azangle
311 WRITE (UNIT=stderr,FMT=*) 'profiles%p = ', profiles(1)%p
312 WRITE (UNIT=stderr,FMT=*) 'profiles%t = ', profiles(1)%t
313 WRITE (UNIT=stderr,FMT=*) 'profiles%q = ', profiles(1)%q
314 do i=coef%nlevels,1,-1
315 write(UNIT=stderr,FMT='(i4,3f12.2)') i, profiles(1)%p(i),profiles(1)%t(i),profiles(1)%q(i)
316 end do
317 WRITE (UNIT=stderr,FMT=*) 'addcloud = ', addcloud
318 WRITE (UNIT=stderr,FMT=*) 'calcemis = ', calcemis
319 WRITE (UNIT=stderr,FMT=*) 'emissivity = ', emissivity
320 WRITE (UNIT=stderr,FMT=*) 'radiance = ', radiance%out_clear
321 call da_warning(__FILE__,__LINE__,(/"Problem in rttov_direct"/))
322 end if
323
324 nc = nbtout / nprofiles
325 do n = 1, nprofiles
326 tb(1:nc) = radiance % out_clear((n-1)*nc+1:n*nc)
327 end do
328 ! calcemis_out (:) = calcemis(:)
329 emissivity_out(:) = emissivity(:)
330
331 deallocate (lprofiles)
332 deallocate (channels)
333 deallocate (polarisations)
334 deallocate (emissivity)
335 deallocate (calcemis)
336 deallocate (surfem)
337 do n = 1, nprofiles
338 deallocate (profiles(n)%p)
339 deallocate (profiles(n)%t)
340 deallocate (profiles(n)%q)
341 deallocate (profiles(n)%o3)
342 deallocate (profiles(n)%co2)
343 deallocate (profiles(n)%clw)
344 end do
345
346 ! deallocate transmittance structure
347 deallocate (transmission % tau_surf ,stat= alloc_status(6))
348 deallocate (transmission % tau_layer ,stat= alloc_status(7))
349 deallocate (transmission % od_singlelayer,stat= alloc_status(8))
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 if (any(alloc_status /= 0)) then
368 call da_error(__FILE__,__LINE__, &
369 (/"mem deallocation error"/))
370 end if
371
372 call da_trace_exit("da_rttov_direct")
373
374 end subroutine da_rttov_direct
375 #endif