da_transform_xtoy_crtm_adj.inc
References to this file elsewhere.
1 subroutine da_transform_xtoy_crtm_adj ( iv, jo_grad_y, jo_grad_x )
2
3 !---------------------------------------------------------------------------
4 ! PURPOSE: transform gradient from obs space to model grid space.
5 !
6 ! METHOD: jo_grad_x = H^T jo_grad_y = - H^T R^-1 ( d - H delta_x )
7 ! 1. input gradient in obs space and reference state of RTTOV
8 ! 2. call adjoint of RTM
9 ! 3. adjoint of interpolation from model grid to obs loc
10 !
11 ! HISTORY: 11/19/2006 - Creation Zhiquan Liu
12 !
13 !---------------------------------------------------------------------------
14
15 implicit none
16
17 type (x_type), intent(inout) :: jo_grad_x !
18 type (y_type), intent(in) :: jo_grad_y ! H' delta_x
19 type (iv_type), intent(in) :: iv ! O-B structure.
20
21 #ifdef CRTM
22
23 integer :: k ! Index dimension.
24 integer :: num_rad ! Number of radiance obs
25 integer :: inst, nchanl, n
26 real, allocatable :: q_ad(:,:)
27 real, allocatable :: t(:,:)
28 real, allocatable :: p(:)
29 ! type(infa_type), pointer :: info
30 ! integer :: i,j
31 ! real :: dx,dy,dxm,dym
32
33
34 ! CRTM local varaibles and types
35 integer :: wmo_sensor_id, Error_Status, Allocate_Status
36 type (CRTM_RTSolution_type ), allocatable :: RTSolution(:,:),RTSolution_AD(:,:)
37 type (CRTM_Atmosphere_type ) :: Atmosphere(1), Atmosphere_AD(1)
38 type (CRTM_Surface_type ) :: Surface(1), Surface_AD(1)
39 type (CRTM_GeometryInfo_type ) :: GeometryInfo(1)
40 !type(CRTM_Options_type ) :: Options
41
42 !---------------------------------------------------------
43
44 if ( iv%num_inst < 1 ) return
45
46 if (trace_use) call da_trace_entry("da_transform_xtoy_crtm_adj")
47
48 !----------------------------------------------------------------------------
49 ! CRTM allocation
50 !
51 ! Atmosphere structure
52 Atmosphere(1)%n_Layers=(kte-kts)+1 ! number of vertical levels
53 Atmosphere(1)%n_Absorbers=2
54 Atmosphere(1)%n_Clouds=0
55 Atmosphere(1)%n_Aerosols=0
56 if (crtm_cloud) Atmosphere(1)%n_Clouds=6
57
58 Error_Status = CRTM_Allocate_Atmosphere( Atmosphere(1)%n_Layers, &
59 Atmosphere(1)%n_Absorbers, &
60 Atmosphere(1)%n_Clouds, &
61 Atmosphere(1)%n_Aerosols, &
62 Atmosphere)
63 if ( Error_Status /= 0 ) then
64 call da_error(__FILE__,__LINE__, &
65 (/"Error in allocatting CRTM Atmosphere Structure"/))
66 end if
67 Atmosphere(1)%Absorber_ID(1)=H2O_ID
68 Atmosphere(1)%Absorber_ID(2)=O3_ID
69
70 if (crtm_cloud) then
71 Atmosphere(1)%Cloud(1)%Type=WATER_CLOUD
72 Atmosphere(1)%Cloud(2)%Type=ICE_CLOUD
73 Atmosphere(1)%Cloud(3)%Type=RAIN_CLOUD
74 Atmosphere(1)%Cloud(4)%Type=SNOW_CLOUD
75 Atmosphere(1)%Cloud(5)%Type=GRAUPEL_CLOUD
76 Atmosphere(1)%Cloud(6)%Type=HAIL_CLOUD
77 end if
78
79 !-------------------------------------------------------------------------------
80
81 do inst = 1, iv%num_inst ! loop for sensor
82 num_rad = iv%instid(inst)%info%n2 - iv%instid(inst)%info%n1 + 1
83 if ( num_rad < 1 ) cycle
84 allocate (q_ad(kts:kte,iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
85 allocate (t(kts:kte,iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
86 allocate (p(iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
87
88 ! info => iv%instid(inst)%info
89
90 ! CRTM channel information structure
91 ! Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(inst),ChannelInfo)
92 ! if ( Error_Status /= 0) then
93 ! call da_error(__FILE__,__LINE__, &
94 ! (/"Error in calling CRTM_Set_ChannelInfo"/))
95 ! end if
96 nchanl = ChannelInfo(inst)%n_channels
97
98 ! Allocate forward model solution RTSolution array to number of channels
99 allocate( RTSolution( ChannelInfo(inst)%n_Channels, 1 ) , &
100 RTSolution_AD( ChannelInfo(inst)%n_Channels, 1 ), &
101 STAT = Allocate_Status )
102 if ( Allocate_Status /= 0 ) then
103 call da_error(__FILE__,__LINE__, &
104 (/"Error in allocatting RTSolution"/))
105 END IF
106
107 ! CRTM Surface Structure
108 if (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsua') then
109 wmo_sensor_id=WMO_AMSUA
110 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
111 wmo_sensor_id=WMO_AMSUB
112 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
113 wmo_sensor_id=WMO_AMSRE
114 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
115 wmo_sensor_id=WMO_SSMI
116 else
117 wmo_sensor_id=INVALID_WMO_SENSOR_ID
118 end if
119
120 Error_Status = CRTM_Allocate_Surface( nchanl, & ! Input
121 Surface) ! Output
122 if ( Error_Status /= 0 ) then
123 call da_error(__FILE__,__LINE__, &
124 (/"Error in allocatting CRTM Surface Structure"/))
125 end if
126
127 ! CRTM Options structure
128 !Options%n_channels = nchanl
129 !Error_Status = CRTM_Allocate_Options( nchanl, & ! Input
130 ! Options) ! Output
131 !if ( Error_Status /= 0 ) then
132 ! call da_error(__FILE__,__LINE__, &
133 ! (/"Error in allocatting CRTM Options Structure"/))
134 !endif
135
136 !do n=1,num_rad
137 do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2
138
139 ! [1.0] Extract base state Atmosphere variables
140 Atmosphere(1)%level_pressure(0) = iv%instid(inst)%pf(0,n)
141 do k=1,Atmosphere(1)%n_layers
142 Atmosphere(1)%pressure(k) = iv%instid(inst)%pm(k,n)
143 Atmosphere(1)%level_pressure(k) = iv%instid(inst)%pf(k,n)
144 Atmosphere(1)%temperature(k) = iv%instid(inst)%tm(k,n)
145 Atmosphere(1)%absorber(k,1) = iv%instid(inst)%qm(k,n)
146 if (crtm_cloud) then
147 Atmosphere(1)%cloud(1)%water_content(k)=iv%instid(inst)%qcw(k,n)
148 Atmosphere(1)%cloud(2)%water_content(k)=iv%instid(inst)%qci(k,n)
149 Atmosphere(1)%cloud(3)%water_content(k)=iv%instid(inst)%qrn(k,n)
150 Atmosphere(1)%cloud(4)%water_content(k)=iv%instid(inst)%qsn(k,n)
151 Atmosphere(1)%cloud(5)%water_content(k)=iv%instid(inst)%qgr(k,n)
152 Atmosphere(1)%cloud(6)%water_content(k)=0.0
153 Atmosphere(1)%cloud(1)%effective_radius(k)=10.0
154 Atmosphere(1)%cloud(2)%effective_radius(k)=200.0
155 Atmosphere(1)%cloud(3)%effective_radius(k)=200.0
156 Atmosphere(1)%cloud(4)%effective_radius(k)=200.0
157 Atmosphere(1)%cloud(5)%effective_radius(k)=200.0
158 Atmosphere(1)%cloud(6)%effective_radius(k)=200.0
159 end if
160 end do
161
162 ! [1.1] User-supplied emissivity
163 ! Options%emissivity_switch = 1
164 ! Options%emissivity(1:Options%n_channels) = &
165 ! iv%instid(inst)%emiss(1:Options%n_channels,n)
166
167 ! [1.1] CRTM Surface parameter data
168 Surface(1)%Land_Coverage=iv%instid(inst)%land_coverage(n)
169 Surface(1)%Water_Coverage=iv%instid(inst)%water_coverage(n)
170 Surface(1)%Snow_Coverage=iv%instid(inst)%snow_coverage(n)
171 Surface(1)%Ice_Coverage=iv%instid(inst)%ice_coverage(n)
172
173 if (Surface(1)%Land_Coverage > 0.0) then
174 Surface(1)%Land_Type=GRASS_SOIL ! land type (User guide appendix 3)
175 Surface(1)%Land_Temperature=iv%instid(inst)%ts(n) ! K
176 Surface(1)%Soil_Moisture_Content=iv%instid(inst)%smois(n) !0.05 ! volumetric water content (g/cm**3)
177 !Surface(1)%Canopy_Water_Content=0.05 ! gravimetric water content
178 Surface(1)%Vegetation_Fraction=iv%instid(inst)%vegtyp(n)
179 Surface(1)%Soil_Temperature=iv%instid(inst)%tslb(n)
180 end if
181 if (Surface(1)%Water_Coverage > 0.0) then
182 !Surface(1)%Water_Type=SEA_WATER ! (Currently NOT used)
183 Surface(1)%Water_Temperature=iv%instid(inst)%ts(n) ! K
184 Surface(1)%Wind_Speed=sqrt((iv%instid(inst)%u10(n))**2+ &
185 (iv%instid(inst)%v10(n))**2) ! m/sec
186 !surface(1)%Wind_Direction=0.0 ! NOT used
187 Surface(1)%Salinity=33.0 ! ppmv
188 end if
189 if (Surface(1)%Snow_Coverage > 0.0) then
190 Surface(1)%Snow_Type=NEW_SNOW ! User guide appendix 3
191 Surface(1)%Snow_Temperature=iv%instid(inst)%ts(n) ! K
192 Surface(1)%Snow_Depth=iv%instid(inst)%snowh(n) ! mm
193 !Surface(1)%Snow_Density=0.2 ! g/cm**3
194 !Surface(1)%Snow_Grain_Size=2.0 ! mm
195 end if
196 if (Surface(1)%Ice_Coverage > 0.0) then
197 !Surface(1)%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
198 Surface(1)%Ice_Temperature=iv%instid(inst)%ts(n) ! K
199 Surface(1)%Ice_Thickness=10.0 ! mm
200 !Surface(1)%Ice_Density=0.9 ! g/cm**3
201 !Surface(1)%Ice_Roughness=0.0 ! NO Table offered, single example is ZERO
202 end if
203 Surface(1)%SensorData%n_channels = nchanl
204 Surface(1)%SensorData%Sensor_ID = wmo_sensor_id
205 Surface(1)%SensorData%Tb(1:nchanl) = iv%instid(inst)%tb_inv(1:nchanl,n) + &
206 iv%instid(inst)%tb_xb(1:nchanl,n)
207
208 ! -- Copy the adjoint atmosphere structure
209 Error_Status = CRTM_Assign_Atmosphere( Atmosphere, Atmosphere_AD )
210
211 if ( Error_Status /= 0 ) then
212 call da_error(__FILE__,__LINE__, &
213 (/"Error copying Atmosphere_AD structure"/))
214 END IF
215
216 ! -- Copy the adjoint surface structure
217 Error_Status = CRTM_Assign_Surface( Surface, Surface_AD )
218
219 if ( Error_Status /= 0 ) then
220 call da_error(__FILE__,__LINE__, &
221 (/"Error copying Surface_AD structure"/))
222 END IF
223
224 ! -- Zero the Adjoint outputs
225 ! Important: adjoint variables must be initialized
226 call CRTM_Zero_Atmosphere( Atmosphere_AD )
227 call CRTM_Zero_Surface( Surface_AD )
228
229 ! [1.2] CRTM GeometryInfo Structure
230 GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
231 GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
232 ! GeometryInfo(1)%Satellite_Height=830.0
233 ! GeometryInfo(1)%Sensor_Scan_Angle=
234 ! GeometryInfo(1)%Sensor_Zenith_Angle=
235 ! GeometryInfo(1)%Sensor_Scan_Angle=
236 ! GeometryInfo(1)%Source_Zenith_Angle=
237
238
239 ! [1.3] assign tb = R^-1 Re :
240
241 RTSolution_AD(:,1)%brightness_temperature = jo_grad_y%instid(inst)%tb(:,n)
242 RTSolution_AD(:,1)%radiance = 0.0 ! must assign zero, since each call of AD model will return non-zero value
243 !if (n <= 10) write(6,'(15f8.3)') RTSolution_AD(:,1)%brightness_temperature
244
245 ! [1.4] Call CRTM_AD model
246
247 call da_crtm_ad (1, nchanl, 1, Atmosphere, &
248 Surface, &
249 RTSolution_AD,&
250 GeometryInfo, &
251 ChannelInfo(inst:inst), &
252 Atmosphere_AD,&
253 Surface_AD, &
254 RTSolution) !, &
255 !Options = Options)
256
257 !Error_Status = CRTM_Adjoint(Atmosphere, &
258 ! Surface, &
259 ! RTSolution_AD,&
260 ! GeometryInfo, &
261 ! ChannelInfo, &
262 ! Atmosphere_AD,&
263 ! Surface_AD, &
264 ! RTSolution) !, &
265 ! !Options = Options)
266 !if ( Error_Status /= 0 ) then
267 ! call da_error(__FILE__,__LINE__, &
268 ! (/"Error in calling CRTM_Adjoint"/))
269 !endif
270
271 ! [1.5] Scale transformation and fill zero for no-control variable
272
273 Atmosphere_AD(1)%Level_Pressure(Atmosphere(1)%n_layers) = &
274 Atmosphere_AD(1)%Level_Pressure(Atmosphere(1)%n_layers)*0.01
275 Atmosphere_AD(1)%Absorber(:,1) = 1000.0*Atmosphere_AD(1)%Absorber(:,1) ! in g/kg
276
277 ! [1.7] Adjoint of Interpolate horizontally from ob to grid:
278
279 do k=kts,kte ! from bottom to top
280 if (atmosphere(1)%pressure(kte-k+1) < 75.0) then
281 q_ad(k,n) = 0.0
282 else
283 q_ad(k,n) = Atmosphere_AD(1)%Absorber(kte-k+1,1)
284 end if
285 t(k,n)=Atmosphere_AD(1)%Temperature(kte-k+1)
286 end do
287
288 p(n) = Atmosphere_AD(1)%Level_Pressure(atmosphere(1)%n_layers)
289
290 !if (n <=10 ) then
291 ! write(6,'(15f8.2)') RTSolution(:,1)%brightness_temperature
292 ! write(6,'(e15.5)') jo_grad_x% psfc(i,j)
293 ! do k=kts,kte
294 ! write(6,'(2e15.5)') jo_grad_x%t(i,j,k), jo_grad_x%q(i,j,k)
295 ! end do
296 !endif
297
298 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere_AD )
299 if ( Error_Status /= 0 ) then
300 call da_error(__FILE__,__LINE__, &
301 (/"Error in deallocatting CRTM Atmosphere_AD Structure"/))
302 end if
303
304 Error_Status = CRTM_Destroy_Surface(Surface_AD)
305 if ( Error_Status /= 0 ) then
306 call da_error(__FILE__,__LINE__, &
307 (/"Error in deallocatting CRTM Surface_AD Structure"/))
308 end if
309
310 end do ! end loop for pixels
311
312 call da_interp_lin_2d_adj_partial(jo_grad_x%t(:,:,kts:kte), iv%instid(inst)%info, kts,kte, t)
313 call da_interp_lin_2d_adj_partial(jo_grad_x%q(:,:,kts:kte), iv%instid(inst)%info, kts,kte, q_ad)
314 call da_interp_lin_2d_adj_partial(jo_grad_x%psfc, iv%instid(inst)%info, 1,1, p)
315
316 deallocate (q_ad)
317 deallocate (t)
318 deallocate (p)
319
320 !-------------------------------------------------------------------
321 ! [2.0] Deallocating CRTM structures
322 !-------------------------------------------------------------------
323 deallocate( RTSolution, RTSolution_AD, STAT = Allocate_Status )
324 if ( Allocate_Status /= 0 ) then
325 call da_error(__FILE__,__LINE__, &
326 (/"Error in deallocatting RTSolution"/))
327 END IF
328
329 Error_Status = CRTM_Destroy_Surface(Surface)
330 if ( Error_Status /= 0 ) then
331 call da_error(__FILE__,__LINE__, &
332 (/"Error in deallocatting CRTM Surface Structure"/))
333 end if
334 end do ! end loop for sensor
335
336 !-------------------------------------------------------------------
337 ! [3.0] Deallocating CRTM Atmosphere structures
338 !-------------------------------------------------------------------
339 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere )
340 if ( Error_Status /= 0 ) then
341 call da_error(__FILE__,__LINE__, &
342 (/"Error in deallocatting CRTM Atmosphere Structure"/))
343 end if
344
345 if (trace_use) call da_trace_exit("da_transform_xtoy_crtm_adj")
346 #else
347 call da_error(__FILE__,__LINE__, &
348 (/"Must compile with $CRTM option for radiances"/))
349 #endif
350
351 end subroutine da_transform_xtoy_crtm_adj
352