1module tovsLin_mod
2 ! MODULE tovsLin_mod (prefix='tvslin' category='5. Observation operators')
3 !
4 !:Purpose: Derived types, public variables and procedures related to the
5 ! tangent-linear and adjoint versions of RTTOV
6 !
7 use rttovInterfaces_mod
8 use rttov_types, only : &
9 rttov_profile ,&
10 rttov_profile_cloud ,&
11 rttov_radiance ,&
12 rttov_transmission ,&
13 rttov_chanprof ,&
14 rttov_emissivity
15 use rttov_const, only : &
16 gas_unit_specconc ,&
17 sensor_id_mw ,&
18 surftype_sea ,&
19 errorStatus_success
20 use parkind1, only : jpim, jprb
21 use verticalCoord_mod
22 use tovsNL_mod
23 use utilities_mod
24 use MathPhysConstants_mod
25 use obsSpaceData_mod
26 use columnData_mod
27
28 implicit none
29 save
30 private
31
32 public :: tvslin_rttov_tl, tvslin_rttov_ad
33
34
35contains
36
37 !--------------------------------------------------------------------------
38 ! tvslin_rttov_tl
39 !--------------------------------------------------------------------------
40 subroutine tvslin_rttov_tl(columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData)
41 !
42 !:Purpose: Tangent linear of computation of radiance with rttov_tl
43 !
44 implicit none
45
46 ! Arguments:
47 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData structure
48 type(struct_columnData), intent(in) :: columnAnlInc ! column structure for pertubation profile
49 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev ! column structure for background profile
50
51 ! Locals:
52 type(struct_vco), pointer :: vco_anl
53 integer, allocatable :: sensorTovsIndexes(:)
54 integer, allocatable :: sensorHeaderIndexes(:)
55 integer :: allocStatus(5)
56 integer :: nobmax
57 integer :: sensorIndex, tovsIndex
58 integer :: ilowlvl_M,ilowlvl_T,profileCount,headerIndex,levelIndex,nlv_M,nlv_T
59 integer :: profileIndex
60 integer :: Vcode
61 character(len=4) :: ozoneVarName
62 logical, allocatable :: surfTypeIsWater(:)
63 real(8), pointer :: delTT(:), delHU(:), delP(:)
64 real(8), pointer :: delO3(:)
65 real(8), pointer :: delCLW(:)
66 real(8), pointer :: delCIW(:), delRF(:), delSF(:)
67 integer :: btCount
68 integer,external :: omp_get_num_threads
69 integer :: nthreads, max_nthreads
70 integer :: btIndex, bodyIndex
71 integer :: instrum
72 integer :: sensorType !sensor type(1=infrared; 2=microwave; 3=high resolution, 4=polarimetric)
73 integer :: errorStatus
74 integer,allocatable :: sensorBodyIndexes(:)
75 logical,allocatable :: lchannel_subset(:,:)
76 real(8), allocatable :: surfem1(:)
77 integer, allocatable :: frequencies(:)
78 type(rttov_emissivity), pointer :: emissivity_local(:)
79 type(rttov_emissivity), pointer :: emissivity_tl(:)
80 type(rttov_radiance) :: radiancedata_d ! radiances full structure buffer used in rttov calls
81 type(rttov_radiance) :: radiancedata_tl ! tl radiances full structure buffer used in rttov calls
82 type(rttov_transmission) :: transmission ! transmission
83 type(rttov_transmission) :: transmission_tl ! transmission tl
84 type(rttov_profile), pointer :: profilesdata_tl(:) ! tl profiles buffer used in rttov calls
85 type(rttov_profile_cloud), pointer :: cld_profiles_tl(:) !tl profiles buffer used in RttovScatt calls
86 type(rttov_chanprof), pointer :: chanprof(:)
87 logical, pointer :: calcemis(:)
88 logical :: runObsOperatorWithClw_tl
89 logical :: runObsOperatorWithHydrometeors_tl
90 integer :: asw
91 real(8) :: obsOMP
92 type (rttov_profile), pointer :: profiles(:)
93 type(rttov_profile_cloud), pointer :: cld_profiles(:)
94
95 if (tvs_nobtov == 0) return ! exit if there are not tovs data
96
97 write(*,*) 'tvslin_rttov_tl: Starting'
98
99 call tvs_getProfile(profiles, 'tlad', cld_profiles)
100
101 if (.not. tvs_useO3Climatology .and. .not. col_varExist(columnTrlOnAnlIncLev,'TO3') .and. .not. col_varExist(columnTrlOnAnlIncLev,'O3L') ) then
102 call utl_abort('tvslin_rttov_tl: if tvs_useO3Climatology is set to .true. the ozone variable must be included as an analysis variable in NAMSTATE.')
103 else if (.not.tvs_useO3Climatology) then
104 if (col_varExist(columnTrlOnAnlIncLev,'TO3')) then
105 ozoneVarName = 'TO3'
106 else
107 ozoneVarName = 'O3L'
108 end if
109 end if
110
111 ! 1. Set index for model's lowest level and model top
112
113 nlv_M = col_getNumLev(columnTrlOnAnlIncLev,'MM')
114 nlv_T = col_getNumLev(columnTrlOnAnlIncLev,'TH')
115
116 if ( col_getPressure(columnTrlOnAnlIncLev,1,1,'TH') < col_getPressure(columnTrlOnAnlIncLev,nlv_T,1,'TH') ) then
117 ilowlvl_M = nlv_M
118 ilowlvl_T = nlv_T
119 else
120 ilowlvl_M = 1
121 ilowlvl_T = 1
122 end if
123
124 vco_anl => col_getVco(columnTrlOnAnlIncLev)
125 Vcode = vco_anl%Vcode
126
127
128 ! 1. Get number of threads available and allocate memory for some variables
129 ! . ----------------------------------------------------------------------
130
131 !$omp parallel
132 max_nthreads = omp_get_num_threads()
133 !$omp end parallel
134
135 allocStatus(:) = 0
136 allocate ( sensorTovsIndexes(tvs_nobtov), stat = allocStatus(1) )
137 call utl_checkAllocationStatus(allocStatus(1:1), ' tvslin_rttov_tl sensorTovsIndexes')
138
139 ! 2. Computation of hx for tovs data only
140
141
142 ! Loop over all sensors specified by user
143
144 sensor_loop: do sensorIndex = 1, tvs_nsensors
145
146 runObsOperatorWithClw_tl = col_varExist(columnTrlOnAnlIncLev,'LWCR') .and. &
147 tvs_isInstrumUsingCLW(tvs_instruments(sensorIndex)) .and. &
148 tvs_mwInstrumUsingCLW_tl
149 runObsOperatorWithHydrometeors_tl = col_varExist(columnTrlOnAnlIncLev,'LWCR') .and. &
150 col_varExist(columnTrlOnAnlIncLev,'IWCR') .and. &
151 tvs_isInstrumUsingHydrometeors(tvs_instruments(sensorIndex)) .and. &
152 tvs_mwInstrumUsingHydrometeors_tl
153
154 sensorType = tvs_coefs(sensorIndex) % coef % id_sensor
155 instrum = tvs_coefs(sensorIndex) % coef % id_inst
156 ! loop over all obs.
157 profileCount = 0
158 do tovsIndex = 1, tvs_nobtov
159 ! Currently processed sensor?
160 if ( tvs_lsensor(tovsIndex) == sensorIndex ) then
161 profileCount = profileCount + 1
162 sensorTovsIndexes(profileCount) = tovsIndex
163 end if
164 end do
165
166 if (profileCount == 0) cycle sensor_loop
167 nobmax = sensorTovsIndexes(profileCount)
168 ! compute the number of calculated radiances for one call
169 btCount = tvs_countRadiances(sensorTovsIndexes(1:profileCount), obsSpaceData, &
170 assim_flag_val_opt=obs_assimilated)
171 if ( btCount == 0 ) cycle sensor_loop
172
173 allocate (sensorHeaderIndexes (profileCount), stat= allocStatus(1))
174 allocate (profilesdata_tl(profileCount), stat= allocStatus(2))
175 if (runObsOperatorWithClw_tl) write(*,*) 'tvslin_rttov_tl: using clw_data'
176 if (runObsOperatorWithHydrometeors_tl) write(*,*) 'tvslin_rttov_tl: using hydrometeor data'
177 allocate (surfTypeIsWater(profileCount),stat= allocStatus(3))
178 call utl_checkAllocationStatus(allocStatus, ' tvslin_rttov_tl')
179
180 sensorHeaderIndexes(:) = 0
181
182 surfTypeIsWater(:) = .false.
183
184 ! allocate profiledata_tl structures
185 asw = 1 ! 1 to allocate
186 call rttov_alloc_tl( &
187 allocStatus(1), &
188 asw, &
189 nprofiles=profileCount, &
190 nchanprof=btCount, &
191 nlevels=nlv_T, &
192 chanprof=chanprof, &
193 opts=tvs_opts(sensorIndex), &
194 profiles_tl=profilesdata_tl, &
195 coefs=tvs_coefs(sensorIndex), &
196 transmission=transmission, &
197 transmission_tl=transmission_tl, &
198 radiance=radiancedata_d, &
199 radiance_tl=radiancedata_tl, &
200 calcemis=calcemis, &
201 emissivity=emissivity_local, &
202 emissivity_tl=emissivity_tl, &
203 init=.true.)
204 if (runObsOperatorWithHydrometeors_tl) then
205 allocate(cld_profiles_tl(profileCount))
206 call rttov_alloc_scatt_prof ( allocStatus(2), &
207 profileCount, &
208 cld_profiles_tl, &
209 nlv_T, &
210 nhydro=5, &
211 nhydro_frac=1, &
212 asw=asw, &
213 init=.false., &
214 flux_conversion=[1,2,0,0,0])
215 end if
216
217 call utl_checkAllocationStatus(allocStatus(1:2), ' tvslin_rtttov_tl rttov_alloc_tl 1')
218
219 profileCount = 0
220
221 obs_loop: do tovsIndex = 1, nobmax
222 if (tvs_lsensor(tovsIndex) /= sensorIndex) cycle obs_loop
223 headerIndex = tvs_headerIndex(tovsIndex)
224 profileCount = profileCount + 1
225 surfTypeIsWater(profileCount) = ( tvs_ChangedStypValue(obsSpaceData,headerIndex) == surftype_sea )
226 sensorHeaderIndexes(profileCount) = headerIndex
227 end do obs_loop
228
229 do profileIndex = 1 , profileCount
230 profilesdata_tl(profileIndex) % gas_units = gas_unit_specconc ! all gas profiles should be provided in kg/kg
231 profilesdata_tl(profileIndex) % nlevels = nlv_T
232 profilesdata_tl(profileIndex) % nlayers = nlv_T - 1
233 if (tvs_coefs(sensorIndex)%coef%nozone > 0) then
234 if (tvs_useO3Climatology) then
235 profilesdata_tl(profileIndex) % o3(:) = 0.0d0
236 else
237 delO3 => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),trim(ozoneVarName))
238 profilesdata_tl(profileIndex) % o3(1:nlv_T) = delO3(1:nlv_T) * 1.0d-9 ! Assumes model ozone in ug/kg
239 profilesdata_tl(profileIndex) % s2m % o = col_getElem(columnAnlInc,ilowlvl_T,sensorHeaderIndexes(profileIndex),trim(ozoneVarName)) * 1.0d-9 ! Assumes model ozone in ug/kg
240 end if
241 end if
242
243 ! using the zero CLW value for land FOV
244 if (runObsOperatorWithClw_tl) then
245 if (surfTypeIsWater(profileIndex)) then
246 delCLW => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'LWCR')
247 profilesdata_tl(profileIndex) % clw(1:nlv_T) = delCLW(:)
248 else
249 profilesdata_tl(profileIndex) % clw(1:nlv_T) = 0.d0
250 end if
251 end if
252
253 if (runObsOperatorWithHydrometeors_tl) then
254 if (surfTypeIsWater(profileIndex)) then
255 ! rain flux
256 if (col_varExist(columnAnlInc,'RF')) then
257 delRF => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'RF')
258 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,1) = delRF(:)
259 else
260 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,1) = 0.0d0
261 end if
262
263 ! snow flux
264 if (col_varExist(columnAnlInc,'SF')) then
265 delSF => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'SF')
266 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,2) = delSF(:)
267 else
268 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,2) = 0.0d0
269 end if
270
271 ! graupel
272 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,3) = 0.d0 ! no information for graupel
273
274 ! cloud liquid water content
275 delCLW => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'LWCR')
276 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,4) = delCLW(:)
277
278 ! cloud ice water content
279 delCIW => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'IWCR')
280 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,5) = delCIW(:)
281 else
282 cld_profiles_tl(profileIndex) % hydro(1:nlv_T,1:5) = 0.d0
283 end if ! surfTypeIsWater
284
285 cld_profiles_tl(profileIndex) % hydro_frac(1:nlv_T,1) = 0.d0 ! no perturbation on cloud fraction as it is a binary variable (or or 1.0) in this implementation
286 end if ! runObsOperatorWithHydrometeors_tl
287
288 profilesdata_tl(profileIndex) % ctp = 0.0d0
289 profilesdata_tl(profileIndex) % cfraction = 0.0d0
290 profilesdata_tl(profileIndex) % zenangle = 0.0d0
291 profilesdata_tl(profileIndex) % azangle = 0.0d0
292 profilesdata_tl(profileIndex) % skin % surftype = 0
293 profilesdata_tl(profileIndex) % skin % t = col_getElem(columnAnlInc,1,sensorHeaderIndexes(profileIndex),'TG')
294 profilesdata_tl(profileIndex) % skin % fastem(:)= 0.0d0
295 profilesdata_tl(profileIndex) % skin % salinity = 0.0d0
296 profilesdata_tl(profileIndex) % s2m % t = col_getElem(columnAnlInc,ilowlvl_T,sensorHeaderIndexes(profileIndex),'TT')
297 profilesdata_tl(profileIndex) % s2m % q = 0.d0
298
299 profilesdata_tl(profileIndex) % s2m % p = col_getElem(columnAnlInc,1,sensorHeaderIndexes(profileIndex),'P0')*MPC_MBAR_PER_PA_R8
300 profilesdata_tl(profileIndex) % s2m % u = col_getElem(columnAnlInc,ilowlvl_M,sensorHeaderIndexes(profileIndex),'UU')
301 profilesdata_tl(profileIndex) % s2m % v = col_getElem(columnAnlInc,ilowlvl_M,sensorHeaderIndexes(profileIndex),'VV')
302
303 delP => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'P_T')
304 profilesdata_tl(profileIndex) % p(1:nlv_T) = delP(:) * MPC_MBAR_PER_PA_R8
305 delTT => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'TT')
306 profilesdata_tl(profileIndex) % t(1:nlv_T) = delTT(:)
307 delHU => col_getColumn(columnAnlInc,sensorHeaderIndexes(profileIndex),'HU')
308 profilesdata_tl(profileIndex) % q(1:nlv_T) = delHU(:)
309 if (runObsOperatorWithHydrometeors_tl) then
310 cld_profiles_tl(profileIndex) % ph (1) = 0.d0
311 cld_profiles_tl(profileIndex) % cfrac = 0.d0
312 do levelIndex = 1, nlv_T - 1
313 cld_profiles_tl(profileIndex) % ph (levelIndex+1) = 0.5d0 * (profilesdata_tl(profileIndex) % p(levelIndex) + profilesdata_tl(profileIndex) % p(levelIndex+1) )
314 end do
315 cld_profiles_tl(profileIndex) % ph (nlv_T+1) = profilesdata_tl(profileIndex) % s2m % p
316 end if
317 end do
318
319 deallocate (sensorHeaderIndexes, stat= allocStatus(1) )
320 deallocate (surfTypeIsWater,stat= allocStatus(2))
321 call utl_checkAllocationStatus(allocStatus, ' tvslin_rttov_tl', .false.)
322
323 ! set nthreads to actual number of threads which will be used.
324
325 nthreads = min(max_nthreads, profileCount)
326
327 ! 2.2 Prepare all input variables required by rttov.
328
329 allocate ( surfem1(btCount) ,stat=allocStatus(1))
330 allocate ( sensorBodyIndexes(btCount) ,stat=allocStatus(2))
331 if (runObsOperatorWithHydrometeors_tl) then
332 allocate (frequencies(btCount), stat=allocStatus(3))
333 end if
334 call utl_checkAllocationStatus(allocStatus(1:3), ' tvslin_rtttov_tl')
335
336 ! get Hyperspecral IR emissivities
337 if ( tvs_isInstrumHyperSpectral(instrum) ) call tvs_getHIREmissivities(sensorTovsIndexes(1:profileCount), &
338 obsSpaceData, surfem1)
339 allocate( lchannel_subset(profileCount,tvs_nchan(sensorIndex)) )
340 call tvs_getChanprof(sensorTovsIndexes(1:profileCount), obsSpaceData, chanprof, &
341 iptobs_cma_opt=sensorBodyIndexes, lchannel_subset_opt = lchannel_subset)
342 if (runObsOperatorWithHydrometeors_tl) then
343 call rttov_scatt_setupindex ( &
344 errorStatus, &
345 profileCount, & ! number of profiles
346 tvs_nchan(sensorIndex), & ! number of channels
347 tvs_coefs(sensorIndex), & ! coef structure read in from rttov coef file
348 tvs_coef_scatt(sensorIndex), & ! coef structure read in from rttov coef file
349 btcount, & ! number of calculated channels
350 chanprof, & ! channels and profile numbers
351 frequencies, & ! array, frequency number for each channel
352 lchannel_subset ) ! OPTIONAL array of logical flags to indicate a subset of channels
353 if (errorStatus /= errorStatus_success) then
354 write(*,*) 'tvslin_rttov_tl: fatal error in rttov_scatt_setupindex ', errorStatus
355 call utl_abort('tvslin_rttov_tl')
356 end if
357 end if
358 deallocate( lchannel_subset )
359
360 call tvs_getOtherEmissivities(chanprof, sensorTovsIndexes, sensorType, instrum, surfem1, calcemis)
361
362 if (sensorType == sensor_id_mw) then
363 call tvs_getMWemissivityFromAtlas(surfem1(1:btcount), emissivity_local, sensorIndex, chanprof, sensorTovsIndexes(1:profileCount))
364 else
365 emissivity_local(:)%emis_in = surfem1(:)
366 end if
367
368 ! 2.3 Compute tl radiance with rttov_tl
369
370 errorStatus = errorStatus_success
371 emissivity_tl(:)%emis_in = 0.0d0
372
373 if (runObsOperatorWithHydrometeors_tl) then
374 call rttov_scatt_tl( &
375 errorStatus, & ! out
376 tvs_opts_scatt(sensorIndex), & ! in
377 nlv_T, & ! in
378 chanprof, & ! in
379 frequencies, & ! in
380 profiles(sensorTovsIndexes(1:profileCount)), & ! in
381 cld_profiles(sensorTovsIndexes(1:profileCount)),& ! in
382 tvs_coefs(sensorIndex), & ! in
383 tvs_coef_scatt(sensorIndex), & ! in
384 calcemis, & ! in
385 emissivity_local, & ! inout
386 profilesdata_tl, & ! in
387 cld_profiles_tl, & ! in
388 emissivity_tl, & ! inout
389 radiancedata_d, & ! inout
390 radiancedata_tl) ! inout
391 else
392 call rttov_parallel_tl( &
393 errorStatus, & ! out
394 chanprof, & ! in
395 tvs_opts(sensorIndex), & ! in
396 profiles(sensorTovsIndexes(1:profileCount)), & ! in
397 profilesdata_tl, & ! inout
398 tvs_coefs(sensorIndex), & ! in
399 transmission, & ! inout
400 transmission_tl, & ! inout
401 radiancedata_d, & ! inout
402 radiancedata_tl, & ! inout
403 calcemis=calcemis, & ! in
404 emissivity=emissivity_local, & ! in
405 emissivity_tl=emissivity_tl, & ! inout
406 nthreads=nthreads ) ! in
407 end if
408
409 if (errorStatus /= errorStatus_success) then
410 write(*,*) 'Error in rttov_parallel_tl', errorStatus
411 write(*,*) 'temperature profile=',profiles(sensorTovsIndexes(1)) % t(:)
412 write(*,*) 'temperature increment profile=',profilesdata_tl(1) % t(:)
413 call utl_abort('tvslin_rttov_tl')
414 end if
415
416 ! 2.4 Store hx in obsSpaceData,OBS_WORK
417
418 do btIndex = 1, btCount
419
420 bodyIndex = sensorBodyIndexes(btIndex)
421 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex, &
422 radiancedata_tl % bt(btIndex) )
423 if ( tvs_debug ) then
424 obsOMP = obs_bodyElem_r(obsSpaceData,OBS_OMP,bodyIndex)
425 write(*,'(a,i4,2f8.2)') ' ichn,sim,obs= ', &
426 chanprof(btIndex)%chan, radiancedata_tl % bt(btIndex), obsOMP
427 end if
428
429 end do
430
431 ! de-allocate memory
432
433 asw = 0 ! 0 to deallocate
434 if (runObsOperatorWithHydrometeors_tl) then
435 call rttov_alloc_scatt_prof (allocStatus(1), &
436 profileCount, &
437 cld_profiles_tl, &
438 nlv_T, &
439 nhydro=5, &
440 nhydro_frac=1, &
441 asw=asw, &
442 flux_conversion=[1,2,0,0,0])
443 deallocate(cld_profiles_tl)
444 end if
445 call rttov_alloc_tl( &
446 allocStatus(2), &
447 asw, &
448 nprofiles=profileCount, &
449 nchanprof=btCount, &
450 nlevels=nlv_T, &
451 chanprof=chanprof, &
452 opts=tvs_opts(sensorIndex), &
453 profiles_tl=profilesdata_tl, &
454 coefs=tvs_coefs(sensorIndex), &
455 transmission=transmission, &
456 transmission_tl=transmission_tl, &
457 radiance=radiancedata_d, &
458 radiance_tl=radiancedata_tl, &
459 calcemis=calcemis, &
460 emissivity=emissivity_local, &
461 emissivity_tl=emissivity_tl )
462
463 deallocate ( surfem1, stat=allocStatus(3) )
464 if (allocated(frequencies)) deallocate(frequencies, stat=allocStatus(4))
465 deallocate ( sensorBodyIndexes,stat=allocStatus(5) )
466 call utl_checkAllocationStatus(allocStatus(1:5), ' tvslin_rtttov_tl', .false.)
467
468 end do sensor_loop
469
470 deallocate ( sensorTovsIndexes )
471 nullify( profiles )
472 write(*,*) 'tvslin_rttov_tl: Finished'
473
474 end subroutine tvslin_rttov_tl
475
476 !--------------------------------------------------------------------------
477 ! tvslin_rttov_ad
478 !--------------------------------------------------------------------------
479 subroutine tvslin_rttov_ad( columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData )
480 !
481 !:Purpose: Adjoint of computation of radiance with rttov_ad
482 !
483
484 implicit none
485
486 ! Arguments:
487 type(struct_columnData), intent(inout) :: columnAnlInc
488 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
489 type(struct_obs), intent(in) :: obsSpaceData
490
491 ! Locals:
492 type(struct_vco), pointer :: vco_anl
493 integer, allocatable :: sensorTovsIndexes(:)
494 integer, allocatable :: sensorHeaderIndexes(:)
495 integer :: allocStatus(17)
496 integer :: omp_get_num_threads, nthreads
497 integer :: nobmax
498 integer :: sensorIndex, tovsIndex
499 integer :: ilowlvl_T,ilowlvl_M,profileCount,headerIndex,nlv_M,nlv_T
500 integer :: profileIndex, levelIndex
501 integer :: Vcode
502 real(8), allocatable :: tt_ad(:,:)
503 real(8), allocatable :: hu_ad(:,:)
504 real(8), allocatable :: pressure_ad(:,:)
505 real(8), allocatable :: ozone_ad(:,:)
506 character(len=4) :: ozoneVarName
507 real(8), allocatable :: clw_ad(:,:)
508 real(8), allocatable :: ciw_ad(:,:), rf_ad(:,:), sf_ad(:,:)
509 logical, allocatable :: surfTypeIsWater(:), lchannel_subset(:,:)
510 real(8), pointer :: uu_column(:),vv_column(:),tt_column(:),hu_column(:),ps_column(:)
511 real(8), pointer :: tg_column(:),p_column(:),o3_column(:),clw_column(:)
512 real(8), pointer :: ciw_column(:), rf_column(:),sf_column(:)
513 integer :: btCount
514 integer :: max_nthreads
515 integer :: instrum
516 integer :: btIndex, bodyIndex
517 integer :: sensorType ! sensor type (1=infrared; 2=microwave; 3=high resolution, 4=polarimetric)
518 integer, allocatable :: sensorBodyIndexes(:)
519 integer :: errorStatus
520 real(8), allocatable :: surfem1(:)
521 integer, allocatable :: frequencies(:)
522 type(rttov_emissivity), pointer :: emissivity_local(:)
523 type(rttov_emissivity), pointer :: emissivity_ad(:)
524 type(rttov_transmission) :: transmission,transmission_ad
525 type(rttov_radiance) :: radiancedata_ad, radiancedata_d
526 type(rttov_profile), pointer :: profilesdata_ad(:) ! ad profiles buffer used in rttov calls
527 type(rttov_profile), pointer :: profiles(:)
528 type(rttov_profile_cloud), pointer :: cld_profiles(:)
529 type(rttov_profile_cloud), pointer :: cld_profiles_ad(:)
530 type(rttov_chanprof), pointer :: chanprof(:)
531 integer :: asw
532 logical, pointer :: calcemis (:)
533 logical :: runObsOperatorWithClw_ad
534 logical :: runObsOperatorWithHydrometeors_ad
535
536 if (tvs_nobtov == 0) return ! exit if there are not tovs data
537 write(*,*) 'tvslin_rttov_ad: Starting'
538
539 call tvs_getProfile(profiles, 'tlad', cld_profiles)
540
541 if (.not. tvs_useO3Climatology .and. .not. col_varExist(columnTrlOnAnlIncLev,'TO3') .and. .not. col_varExist(columnTrlOnAnlIncLev,'O3L') ) then
542 call utl_abort('tvslin_rttov_ad: if tvs_useO3Climatology is set to .true. the ozone variable must be included as an analysis variable in NAMSTATE.')
543 else if (.not.tvs_useO3Climatology) then
544 if (col_varExist(columnTrlOnAnlIncLev,'TO3')) then
545 ozoneVarName = 'TO3'
546 else
547 ozoneVarName = 'O3L'
548 end if
549 end if
550
551 ! 1. Set index for model's lowest level and model top
552
553 nlv_M = col_getNumLev(columnTrlOnAnlIncLev,'MM')
554 nlv_T = col_getNumLev(columnTrlOnAnlIncLev,'TH')
555
556 if ( col_getPressure(columnTrlOnAnlIncLev,1,1,'TH') < col_getPressure(columnTrlOnAnlIncLev,nlv_T,1,'TH') ) then
557 ilowlvl_M = nlv_M
558 ilowlvl_T = nlv_T
559 else
560 ilowlvl_M = 1
561 ilowlvl_T = 1
562 end if
563
564 vco_anl => col_getVco(columnTrlOnAnlIncLev)
565 Vcode = vco_anl%Vcode
566
567 ! 1. Get number of threads available and allocate memory for some variables
568
569 !$omp parallel
570 max_nthreads = omp_get_num_threads()
571 !$omp end parallel
572
573 allocate ( sensorTovsIndexes(tvs_nobtov) )
574
575 ! 2. Computation of adjoint hx for tovs data only
576
577 ! Loop over all sensors specified by user
578
579 sensor_loop:do sensorIndex = 1, tvs_nsensors
580
581 runObsOperatorWithClw_ad = col_varExist(columnTrlOnAnlIncLev,'LWCR') .and. &
582 tvs_isInstrumUsingCLW(tvs_instruments(sensorIndex)) .and. &
583 tvs_mwInstrumUsingCLW_tl
584 runObsOperatorWithHydrometeors_ad = col_varExist(columnTrlOnAnlIncLev,'LWCR') .and. &
585 col_varExist(columnTrlOnAnlIncLev,'IWCR') .and. &
586 tvs_isInstrumUsingHydrometeors(tvs_instruments(sensorIndex)) .and. &
587 tvs_mwInstrumUsingHydrometeors_tl
588
589 sensorType = tvs_coefs(sensorIndex) % coef% id_sensor
590 instrum = tvs_coefs(sensorIndex) % coef% id_inst
591
592 profileCount = 0
593 do tovsIndex = 1, tvs_nobtov
594 ! Currently processed sensor?
595 if ( tvs_lsensor(tovsIndex) == sensorIndex ) then
596 profileCount = profileCount + 1
597 sensorTovsIndexes(profileCount) = tovsIndex
598 end if
599 end do
600
601 if (profileCount == 0) cycle sensor_loop
602
603 nobmax = sensorTovsIndexes(profileCount)
604
605 ! compute the number of radiances/tbs to be calculated
606 btCount = tvs_countRadiances(sensorTovsIndexes(1:profileCount), obsSpaceData)
607
608 if (btCount == 0) cycle sensor_loop
609
610 allocStatus(:) = 0
611 allocate (sensorHeaderIndexes(profileCount), stat= allocStatus(1))
612 allocate (tt_ad (nlv_T,profileCount), stat= allocStatus(2))
613 allocate (hu_ad (nlv_T,profileCount), stat= allocStatus(3))
614 allocate (pressure_ad (nlv_T,profileCount), stat= allocStatus(4))
615 if (.not. tvs_useO3Climatology) then
616 if (tvs_coefs(sensorIndex) %coef %nozone > 0) then
617 allocate (ozone_ad(nlv_T,profileCount), stat= allocStatus(5))
618 end if
619 end if
620 if (runObsOperatorWithClw_ad .or. runObsOperatorWithHydrometeors_ad) then
621 allocate (clw_ad(nlv_T,profileCount), stat= allocStatus(6))
622 end if
623 allocate (surfTypeIsWater(profileCount),stat= allocStatus(7))
624 surfTypeIsWater(:) = .false.
625 if (runObsOperatorWithHydrometeors_ad) then
626 allocate (ciw_ad(nlv_T,profileCount), stat= allocStatus(8))
627 allocate (rf_ad(nlv_T,profileCount), stat= allocStatus(9))
628 allocate (sf_ad(nlv_T,profileCount), stat= allocStatus(10))
629 end if
630
631 call utl_checkAllocationStatus(allocStatus, ' tvslin_rttov_ad')
632
633 ! loop over all obs.
634 profileCount = 0
635
636 ! loop over all obs.
637 obs_loop: do tovsIndex = 1, nobmax
638 if (tvs_lsensor(tovsIndex)/=sensorIndex) cycle obs_loop
639 headerIndex = tvs_headerIndex(tovsIndex)
640 profileCount = profileCount + 1
641 sensorHeaderIndexes(profileCount) = headerIndex
642 end do obs_loop
643
644 ! 2.1 Calculate the actual number of threads which will be used.
645
646 nthreads = min(max_nthreads, profileCount )
647
648 ! 2.2 Prepare all input variables required by rttov_ad.
649
650 asw = 1 ! 1 for allocation, 0 for deallocation
651
652 call rttov_alloc_ad( &
653 allocStatus(1), &
654 asw, &
655 profileCount, &
656 btCount, &
657 nlv_T, &
658 chanprof, &
659 opts=tvs_opts(sensorIndex), &
660 profiles_ad=profilesdata_ad, &
661 coefs=tvs_coefs(sensorIndex), &
662 transmission= transmission, &
663 transmission_ad= transmission_ad,&
664 radiance=radiancedata_d, &
665 radiance_ad=radiancedata_ad, &
666 calcemis=calcemis, &
667 emissivity=emissivity_local, &
668 emissivity_ad=emissivity_ad, &
669 init=.true.)
670
671 allocate (surfem1(btCount), stat=allocStatus(2))
672 if (runObsOperatorWithHydrometeors_ad) allocate ( frequencies(btCount), stat=allocStatus(3))
673 allocate (sensorBodyIndexes(btCount), stat=allocStatus(4))
674 if (runObsOperatorWithHydrometeors_ad) then
675 allocate(cld_profiles_ad(profileCount), stat=allocStatus(5))
676 call rttov_alloc_scatt_prof (allocStatus(6), &
677 profileCount, &
678 cld_profiles_ad, &
679 nlv_T, &
680 nhydro=5, &
681 nhydro_frac=1, &
682 asw=asw, &
683 flux_conversion=[1,2,0,0,0])
684 end if
685 call utl_checkAllocationStatus(allocStatus(1:6), ' tvslin_rttov_ad')
686
687 ! get Hyperspectral IR emissivities
688
689 if ( tvs_isInstrumHyperSpectral(instrum) ) call tvs_getHIREmissivities(sensorTovsIndexes(1:profileCount), obsSpaceData, surfem1)
690
691 ! Build the list of channels/profiles indices
692 allocate( lchannel_subset(profileCount,tvs_nchan(sensorIndex)) )
693 call tvs_getChanprof(sensorTovsIndexes(1:profileCount), obsSpaceData, chanprof, &
694 iptobs_cma_opt = sensorBodyIndexes, lchannel_subset_opt = lchannel_subset)
695 if (runObsOperatorWithHydrometeors_ad) then
696 call rttov_scatt_setupindex ( &
697 errorStatus, &
698 profileCount, & ! number of profiles
699 tvs_nchan(sensorIndex), & ! number of channels
700 tvs_coefs(sensorIndex), & ! coef structure read in from rttov coef file
701 tvs_coef_scatt(sensorIndex), & ! coef structure read in from rttov coef file
702 btcount, & ! number of calculated channels
703 chanprof, & ! channels and profile numbers
704 frequencies, & ! array, frequency number for each channel
705 lchannel_subset ) ! OPTIONAL array of logical flags to indicate a subset of channels
706 if (errorStatus /= errorStatus_success) then
707 write(*,*) 'tvslin_rttov_ad: fatal error in rttov_scatt_setupindex ', errorStatus
708 call utl_abort('tvslin_rttov_ad')
709 end if
710 end if
711 deallocate( lchannel_subset )
712 ! get non Hyperspectral IR emissivities
713 call tvs_getOtherEmissivities(chanprof, sensorTovsIndexes, sensorType, instrum, surfem1, calcemis)
714
715 if (sensorType == sensor_id_mw) then
716 call tvs_getMWemissivityFromAtlas(surfem1(1:btcount), emissivity_local, sensorIndex, chanprof, sensorTovsIndexes(1:profileCount))
717 else
718 emissivity_local(:)%emis_in = surfem1(:)
719 end if
720
721 do btIndex = 1, btCount
722 bodyIndex = sensorBodyIndexes(btIndex)
723 radiancedata_ad % bt( btIndex ) = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
724 end do
725
726 ! 2.3 Compute ad radiance with rttov_ad
727
728 errorStatus = errorStatus_success
729 emissivity_ad(:) % emis_in = 0.0d0
730 emissivity_ad(:) % emis_out = 0.0d0
731
732 if (runObsOperatorWithHydrometeors_ad) then
733 call rttov_scatt_ad( &
734 errorStatus, &! out
735 tvs_opts_scatt(sensorIndex), &! in
736 nlv_T, &! in
737 chanprof, &! in
738 frequencies, &! in
739 profiles(sensorTovsIndexes(1:profileCount)), &! in
740 cld_profiles(sensorTovsIndexes(1:profileCount)),&! in
741 tvs_coefs(sensorIndex), &! in
742 tvs_coef_scatt(sensorIndex), &! in
743 calcemis, &! in
744 emissivity_local, &! inout
745 profilesdata_ad, &! inout
746 cld_profiles_ad, &! inout
747 emissivity_ad, &! inout
748 radiancedata_d, &! inout
749 radiancedata_ad) ! inout
750 else
751 call rttov_parallel_ad( &
752 errorstatus, &! out
753 chanprof, &! in
754 tvs_opts(sensorIndex), &! in
755 profiles(sensorTovsIndexes(1:profileCount)), &! in
756 profilesdata_ad, &! in
757 tvs_coefs(sensorIndex), &! in
758 transmission, &! inout
759 transmission_ad, &! inout
760 radiancedata_d, &! inout
761 radiancedata_ad, &! inout
762 calcemis=calcemis, &! in
763 emissivity=emissivity_local, &! inout
764 emissivity_ad=emissivity_ad, &! inout
765 nthreads = nthreads )
766 end if
767 if (errorStatus /= errorStatus_success) then
768 write(*,*) 'Error in rttov_parallel_ad', errorStatus
769 call utl_abort('tvslin_rttov_ad')
770 end if
771
772 ! 2.0 Store adjoints in columnData object
773 tt_ad(:,:) = 0.d0
774 hu_ad(:,:) = 0.d0
775 pressure_ad(:,:) = 0.d0
776 if (.not. tvs_useO3Climatology) then
777 if (tvs_coefs(sensorIndex) %coef %nozone > 0) ozone_ad(:,:) = 0.d0
778 endif
779 if (runObsOperatorWithClw_ad .or. runObsOperatorWithHydrometeors_ad) clw_ad(:,:) = 0.d0
780 if (runObsOperatorWithHydrometeors_ad) then
781 ciw_ad(:,:) = 0.d0
782 rf_ad(:,:) = 0.d0
783 sf_ad(:,:) = 0.d0
784 end if
785
786 do btIndex = 1, btCount
787 profileIndex = chanprof(btIndex)%prof
788 headerIndex = sensorHeaderIndexes(profileIndex)
789
790 ps_column => col_getColumn(columnAnlInc,headerIndex,'P0')
791 p_column => col_getColumn(columnAnlInc,headerIndex,'P_T')
792 tg_column => col_getColumn(columnAnlInc,headerIndex,'TG')
793 tt_column => col_getColumn(columnAnlInc,headerIndex,'TT')
794 hu_column => col_getColumn(columnAnlInc,headerIndex,'HU')
795 uu_column => col_getColumn(columnAnlInc,headerIndex,'UU')
796 vv_column => col_getColumn(columnAnlInc,headerIndex,'VV')
797
798 tt_ad(:,profileIndex) = profilesdata_ad(profileIndex) % t(:)
799 hu_ad(:,profileIndex) = profilesdata_ad(profileIndex) % q(:)
800 pressure_ad(:,profileIndex) = profilesdata_ad(profileIndex) % p(:)
801 tg_column(1) = profilesdata_ad(profileIndex) % skin % t
802 tt_column(ilowlvl_T) = profilesdata_ad(profileIndex) % s2m % t
803 ps_column(1) = profilesdata_ad(profileIndex) % s2m % p * MPC_MBAR_PER_PA_R8
804 hu_column(ilowlvl_T) = 0.d0
805 uu_column(ilowlvl_M) = profilesdata_ad(profileIndex) % s2m % u
806 vv_column(ilowlvl_M) = profilesdata_ad(profileIndex) % s2m % v
807
808 if (.not. tvs_useO3Climatology) then
809 if (tvs_coefs(sensorIndex) %coef %nozone > 0) then
810 ! This step is just to transfer the value for ilowlvl_T to the memory space defined by 'col_getColumn(...trim(ozoneVarName))
811 o3_column => col_getColumn(columnAnlInc,headerIndex,trim(ozoneVarName))
812 o3_column(ilowlvl_T) = profilesdata_ad(profileIndex) % s2m % o * 1.0d-9
813 ozone_ad(:,profileIndex) = profilesdata_ad(profileIndex) % o3(:)
814 end if
815 end if
816
817 if (runObsOperatorWithClw_ad) then
818 clw_ad(:,profileIndex) = profilesdata_ad(profileIndex) % clw(:)
819 end if
820
821 if (runObsOperatorWithHydrometeors_ad) then
822 rf_ad(:,profileIndex) = cld_profiles_ad(profileIndex) % hydro(:,1)
823 sf_ad(:,profileIndex) = cld_profiles_ad(profileIndex) % hydro(:,2)
824 clw_ad(:,profileIndex) = cld_profiles_ad(profileIndex) % hydro(:,4)
825 ciw_ad(:,profileIndex) = cld_profiles_ad(profileIndex) % hydro(:,5)
826 end if
827 end do
828
829 ! . 2.1 Store adjoints in columnData object
830 ! . -----------------------------------
831
832 do profileIndex = 1 , profileCount
833 ps_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex), 'P0')
834 p_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex), 'P_T')
835 tt_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex), 'TT')
836 hu_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex), 'HU')
837
838 do levelIndex = 1, nlv_T
839 p_column(levelIndex) = p_column(levelIndex) + pressure_ad (levelIndex,profileIndex) * MPC_MBAR_PER_PA_R8
840 tt_column(levelIndex) = tt_column(levelIndex) + tt_ad (levelIndex,profileIndex)
841 hu_column(levelIndex) = hu_column(levelIndex) + hu_ad (levelIndex,profileIndex)
842 end do
843 end do
844
845 if (.not. tvs_useO3Climatology) then
846 if (tvs_coefs(sensorIndex) %coef %nozone > 0) then
847 do profileIndex = 1 , profileCount
848 o3_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex),trim(ozoneVarName))
849 do levelIndex = 1, col_getNumLev(columnAnlInc,'TH')
850 o3_column(levelIndex) = o3_column(levelIndex) + ozone_ad(levelIndex,profileIndex) * 1.0d-9
851 end do
852 end do
853 end if
854 end if
855
856 if (runObsOperatorWithClw_ad) then
857 do profileIndex = 1 , profileCount
858 surfTypeIsWater(profileIndex) = ( tvs_ChangedStypValue(obsSpaceData,sensorHeaderIndexes(profileIndex)) == surftype_sea )
859 if (surfTypeIsWater(profileIndex)) then
860 clw_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex),'LWCR')
861 do levelIndex = 1, col_getNumLev(columnAnlInc,'TH')
862 clw_column(levelIndex) = clw_column(levelIndex) + &
863 clw_ad(levelIndex,profileIndex)
864 end do
865 end if
866 end do
867 end if
868
869 if (runObsOperatorWithHydrometeors_ad) then
870 do profileIndex = 1 , profileCount
871 surfTypeIsWater(profileIndex) = (tvs_ChangedStypValue(obsSpaceData,sensorHeaderIndexes(profileIndex)) == surftype_sea)
872 if (surfTypeIsWater(profileIndex)) then
873 ! rain flux
874 if (col_varExist(columnAnlInc,'RF')) then
875 rf_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex),'RF')
876 do levelIndex = 1, col_getNumLev(columnAnlInc,'TH')
877 rf_column(levelIndex) = rf_column(levelIndex) + rf_ad(levelIndex,profileIndex)
878 end do
879 end if
880
881 ! snow flux
882 if (col_varExist(columnAnlInc,'SF')) then
883 sf_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex),'SF')
884 do levelIndex = 1, col_getNumLev(columnAnlInc,'TH')
885 sf_column(levelIndex) = sf_column(levelIndex) + sf_ad(levelIndex,profileIndex)
886 end do
887 end if
888
889 ! cloud liquid/ice water content
890 clw_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex),'LWCR')
891 ciw_column => col_getColumn(columnAnlInc, sensorHeaderIndexes(profileIndex),'IWCR')
892 do levelIndex = 1, col_getNumLev(columnAnlInc,'TH')
893 clw_column(levelIndex) = clw_column(levelIndex) + clw_ad(levelIndex,profileIndex)
894 ciw_column(levelIndex) = ciw_column(levelIndex) + ciw_ad(levelIndex,profileIndex)
895 end do
896 end if ! surfTypeIsWater
897 end do ! profileIndex
898 end if ! runObsOperatorWithHydrometeors_ad
899
900 deallocate (sensorHeaderIndexes, stat= allocStatus(1) )
901 deallocate (tt_ad, stat= allocStatus(2) )
902 deallocate (hu_ad, stat= allocStatus(3) )
903 deallocate (pressure_ad, stat= allocStatus(4) )
904 if (.not. tvs_useO3Climatology) then
905 if (tvs_coefs(sensorIndex) %coef %nozone > 0) then
906 deallocate (ozone_ad, stat= allocStatus(5) )
907 end if
908 end if
909 if ( allocated(clw_ad) ) then
910 deallocate (clw_ad,stat=allocStatus(6))
911 end if
912 deallocate (surfTypeIsWater,stat=allocStatus(7))
913 if ( allocated(ciw_ad) ) then
914 deallocate (ciw_ad, stat= allocStatus(8))
915 deallocate (rf_ad, stat= allocStatus(9))
916 deallocate (sf_ad, stat= allocStatus(10))
917 end if
918
919 call utl_checkAllocationStatus(allocStatus, ' tvslin_fill_profiles_ad', .false.)
920
921 ! de-allocate memory
922
923 asw = 0 ! 0 to deallocate
924 if (runObsOperatorWithHydrometeors_ad) then
925 call rttov_alloc_scatt_prof (allocStatus(1), &
926 profileCount, &
927 cld_profiles_ad, &
928 nlv_T, &
929 nhydro=5, &
930 nhydro_frac=1, &
931 asw=asw, &
932 flux_conversion= [1,2,0,0,0])
933 deallocate(cld_profiles_ad, stat=allocStatus(2))
934 end if
935 call rttov_alloc_ad( &
936 allocStatus(3), &
937 asw, &
938 profileCount, &
939 btCount, &
940 nlv_T, &
941 chanprof, &
942 opts=tvs_opts(sensorIndex), &
943 profiles_ad=profilesdata_ad, &
944 coefs=tvs_coefs(sensorIndex), &
945 transmission= transmission, &
946 transmission_ad= transmission_ad,&
947 radiance=radiancedata_d, &
948 radiance_ad=radiancedata_ad, &
949 calcemis=calcemis, &
950 emissivity=emissivity_local, &
951 emissivity_ad=emissivity_ad )
952
953
954 deallocate ( surfem1, stat=allocStatus(4))
955 deallocate ( sensorBodyIndexes, stat=allocStatus(5))
956 if (allocated(frequencies)) deallocate ( frequencies, stat=allocStatus(6))
957 call utl_checkAllocationStatus(allocStatus(1:6), ' tvslin_rttov_ad', .false.)
958
959 end do sensor_loop
960
961 ! 3. Close up
962
963 deallocate ( sensorTovsIndexes )
964 nullify( profiles )
965 write(*,*) 'tvslin_rttov_ad: Finished'
966
967 end subroutine tvslin_rttov_ad
968
969end module tovsLin_mod
