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