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, nchan_emis, Error_Status, Allocate_Status
32 TYPE( CRTM_RTSolution_type ), ALLOCATABLE :: RTSolution(:),RTSolution_AD(:)
33 TYPE (CRTM_ChannelInfo_type) :: ChannelInfo
34 TYPE( CRTM_Atmosphere_type ) :: Atmosphere, Atmosphere_AD
35 TYPE( CRTM_Surface_type ) :: Surface, Surface_AD
36 TYPE( CRTM_GeometryInfo_type ) :: GeometryInfo
37 !TYPE( CRTM_Options_type ) :: Options
38 #endif
39 !---------------------------------------------------------
40
41 #if !defined(CRTM)
42 call da_error(__FILE__,__LINE__, &
43 (/"Must compile with $CRTM option for radiances"/))
44 #else
45 IF ( iv%num_inst < 1 ) return
46
47 if (trace_use) call da_trace_entry("da_transform_xtoy_crtm_adj")
48
49 !----------------------------------------------------------------------------
50 ! CRTM allocation
51 !
52 ! Atmosphere structure
53 Atmosphere%n_Layers=(xp%kte-xp%kts)+1 ! number of vertical levels
54 Atmosphere%n_Absorbers=2
55 Atmosphere%n_Clouds=0
56 Atmosphere%n_Aerosols=0
57 if (crtm_cloud) Atmosphere%n_Clouds=6
58
59 Error_Status = CRTM_Allocate_Atmosphere( Atmosphere%n_Layers, &
60 Atmosphere%n_Absorbers, &
61 Atmosphere%n_Clouds, &
62 Atmosphere%n_Aerosols, &
63 Atmosphere)
64 if ( Error_Status /= 0 ) THEN
65 call da_error(__FILE__,__LINE__, &
66 (/"Error in allocatting CRTM Atmosphere Structure"/))
67 endif
68 Atmosphere%Absorber_ID(1)=H2O_ID
69 Atmosphere%Absorber_ID(2)=O3_ID
70
71 if (crtm_cloud) then
72 Atmosphere%Cloud(1)%Type=WATER_CLOUD
73 Atmosphere%Cloud(2)%Type=ICE_CLOUD
74 Atmosphere%Cloud(3)%Type=RAIN_CLOUD
75 Atmosphere%Cloud(4)%Type=SNOW_CLOUD
76 Atmosphere%Cloud(5)%Type=GRAUPEL_CLOUD
77 Atmosphere%Cloud(6)%Type=HAIL_CLOUD
78 end if
79
80 !-------------------------------------------------------------------------------
81
82 do inst = 1, iv%num_inst ! loop for sensor
83 if ( iv%instid(inst)%num_rad < 1 ) cycle
84 num_rad = iv%instid(inst)%num_rad
85
86 ! CRTM channel information structure
87 Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(inst),ChannelInfo)
88 if ( Error_Status /= 0) then
89 call da_error(__FILE__,__LINE__, &
90 (/"Error in calling CRTM_Set_ChannelInfo"/))
91 endif
92 nchanl = ChannelInfo%n_channels
93
94 ! Allocate forward model solution RTSolution array to number of channels
95 ALLOCATE( RTSolution( ChannelInfo%n_Channels ) , &
96 RTSolution_AD( ChannelInfo%n_Channels ), &
97 STAT = Allocate_Status )
98 IF ( Allocate_Status /= 0 ) THEN
99 call da_error(__FILE__,__LINE__, &
100 (/"Error in allocatting RTSolution"/))
101 END IF
102
103 ! CRTM Surface Structure
104 if (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsua') then
105 nchan_emis=4
106 wmo_sensor_id=WMO_AMSUA
107 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
108 nchan_emis=2
109 wmo_sensor_id=WMO_AMSUB
110 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
111 nchan_emis=12
112 wmo_sensor_id=WMO_AMSRE
113 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
114 nchan_emis=7
115 wmo_sensor_id=WMO_SSMI
116 else
117 nchan_emis=0
118 wmo_sensor_id=INVALID_WMO_SENSOR_ID
119 endif
120
121 Error_Status = CRTM_Allocate_Surface( nchan_emis, & ! Input
122 Surface) ! Output
123 if ( Error_Status /= 0 ) THEN
124 call da_error(__FILE__,__LINE__, &
125 (/"Error in allocatting CRTM Surface Structure"/))
126 endif
127
128 ! CRTM Options structure
129 !Options%n_channels = nchanl
130 !Error_Status = CRTM_Allocate_Options( nchanl, & ! Input
131 ! Options) ! Output
132 !if ( Error_Status /= 0 ) THEN
133 ! call da_error(__FILE__,__LINE__, &
134 ! (/"Error in allocatting CRTM Options Structure"/))
135 !endif
136
137 do n=1,num_rad
138
139 ! [1.0] Extract base state Atmosphere variables
140 atmosphere%level_pressure(0) = iv%instid(inst)%pf(0,n)
141 do k=1,atmosphere%n_layers
142 atmosphere%pressure(k) = iv%instid(inst)%pm(k,n)
143 atmosphere%level_pressure(k) = iv%instid(inst)%pf(k,n)
144 atmosphere%temperature(k) = iv%instid(inst)%tm(k,n)
145 atmosphere%absorber(k,1) = iv%instid(inst)%qm(k,n)
146 if (crtm_cloud) then
147 atmosphere%cloud(1)%water_content(k)=iv%instid(inst)%qcw(k,n)
148 atmosphere%cloud(2)%water_content(k)=iv%instid(inst)%qci(k,n)
149 atmosphere%cloud(3)%water_content(k)=iv%instid(inst)%qrn(k,n)
150 atmosphere%cloud(4)%water_content(k)=iv%instid(inst)%qsn(k,n)
151 atmosphere%cloud(5)%water_content(k)=iv%instid(inst)%qgr(k,n)
152 atmosphere%cloud(6)%water_content(k)=0.
153 atmosphere%cloud(1)%effective_radius(k)=10.
154 atmosphere%cloud(2)%effective_radius(k)=200.
155 atmosphere%cloud(3)%effective_radius(k)=200.
156 atmosphere%cloud(4)%effective_radius(k)=200.
157 atmosphere%cloud(5)%effective_radius(k)=200.
158 atmosphere%cloud(6)%effective_radius(k)=200.
159 end if
160 enddo
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%Land_Coverage=iv%instid(inst)%land_coverage(n)
169 Surface%Water_Coverage=iv%instid(inst)%water_coverage(n)
170 Surface%Snow_Coverage=iv%instid(inst)%snow_coverage(n)
171 Surface%Ice_Coverage=iv%instid(inst)%ice_coverage(n)
172
173 if (Surface%Land_Coverage > 0.) then
174 Surface%Land_Type=GRASS_SOIL ! land type (User guide appendix 3)
175 Surface%Land_Temperature=iv%instid(inst)%ts(n) ! K
176 Surface%Soil_Moisture_Content=iv%instid(inst)%smois(n) !0.05 ! volumetric water content (g/cm**3)
177 !Surface%Canopy_Water_Content=0.05 ! gravimetric water content
178 Surface%Vegetation_Fraction=iv%instid(inst)%vegtyp(n)
179 Surface%Soil_Temperature=iv%instid(inst)%tslb(n)
180 endif
181 if (Surface%Water_Coverage > 0.) then
182 !Surface%Water_Type=SEA_WATER ! (Currently NOT used)
183 Surface%Water_Temperature=iv%instid(inst)%ts(n) ! K
184 Surface%Wind_Speed=sqrt((iv%instid(inst)%u10(n))**2+ &
185 (iv%instid(inst)%v10(n))**2) ! m/sec
186 !surface%Wind_Direction=0.0 ! NOT used
187 Surface%Salinity=33. ! ppmv
188 endif
189 if (Surface%Snow_Coverage > 0.) then
190 Surface%Snow_Type=NEW_SNOW ! User guide appendix 3
191 Surface%Snow_Temperature=iv%instid(inst)%ts(n) ! K
192 Surface%Snow_Depth=iv%instid(inst)%snowh(n) ! mm
193 !Surface%Snow_Density=0.2 ! g/cm**3
194 !Surface%Snow_Grain_Size=2.0 ! mm
195 endif
196 if (Surface%Ice_Coverage > 0.) then
197 !Surface%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
198 Surface%Ice_Temperature=iv%instid(inst)%ts(n) ! K
199 Surface%Ice_Thickness=10. ! mm
200 !Surface%Ice_Density=0.9 ! g/cm**3
201 !Surface%Ice_Roughness=0. ! NO Table offered, single example is ZERO
202 endif
203 Surface%SensorData%n_channels = nchan_emis
204 Surface%SensorData%Sensor_ID = wmo_sensor_id
205 Surface%SensorData%Tb(1:nchan_emis) = iv%instid(inst)%tb_inv(1:nchan_emis,n) + &
206 iv%instid(inst)%tb_xb(1:nchan_emis,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%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
231 GeometryInfo%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
232 ! GeometryInfo%Satellite_Height=830.
233 ! GeometryInfo%Sensor_Scan_Angle=
234 ! GeometryInfo%Sensor_Zenith_Angle=
235 ! GeometryInfo%Sensor_Scan_Angle=
236 ! GeometryInfo%Source_Zenith_Angle=
237
238
239 ! [1.3] assign tb = R^-1 Re :
240
241 RTSolution_AD(:)%brightness_temperature = jo_grad_y%instid(inst)%tb(:,n)
242 RTSolution_AD(:)%radiance = 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(:)%brightness_temperature
244
245 ! [1.4] Call CRTM_AD model
246
247 call da_crtm_ad (nchanl,Atmosphere, &
248 Surface, &
249 RTSolution_AD,&
250 GeometryInfo, &
251 ChannelInfo, &
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%Level_Pressure(atmosphere%n_layers) = &
274 Atmosphere_AD%Level_Pressure(atmosphere%n_layers)*0.01
275 Atmosphere_AD%Absorber(:,1) = 1000.*Atmosphere_AD%Absorber(:,1) ! in g/kg
276
277 !-----------------------------------------------------
278 ! [1.6] Get horizontal interpolation weights:
279 !-----------------------------------------------------
280
281 i = iv%instid(inst)%loc_i(n)
282 j = iv%instid(inst)%loc_j(n)
283 dx = iv%instid(inst)%loc_dx(n)
284 dy = iv%instid(inst)%loc_dy(n)
285 dxm = iv%instid(inst)%loc_dxm(n)
286 dym = iv%instid(inst)%loc_dym(n)
287
288 ! [1.7] Adjoint of Interpolate horizontally from ob to grid:
289
290 do k=xp%kts,xp%kte ! from bottem to top
291 q_ad = Atmosphere_AD%Absorber(xp%kte-k+1,1)
292 if (atmosphere%pressure(xp%kte-k+1) < 75.) q_ad = 0.
293 call da_interp_lin_2d_adj( jo_grad_x%t(:,:,k), xp%ims, xp%ime, xp%jms, &
294 xp%jme, i, j, dx, dy, dxm, dym, &
295 Atmosphere_AD%Temperature(xp%kte-k+1))
296
297 call da_interp_lin_2d_adj( jo_grad_x%q(:,:,k), xp%ims, xp%ime, xp%jms, &
298 xp%jme, i, j, dx, dy, dxm, dym, &
299 q_ad )
300 enddo
301
302 call da_interp_lin_2d_adj(jo_grad_x% psfc, xp%ims, xp%ime, xp%jms, &
303 xp%jme, i, j, dx, dy, dxm, dym, &
304 Atmosphere_AD%Level_Pressure(atmosphere%n_layers) )
305
306 !if (n <=10 ) then
307 ! write(6,'(15f8.2)') rtsolution(:)%brightness_temperature
308 ! write(6,'(e15.5)') jo_grad_x% psfc(i,j)
309 ! do k=xp%kts,xp%kte
310 ! write(6,'(2e15.5)') jo_grad_x%t(i,j,k), jo_grad_x%q(i,j,k)
311 ! enddo
312 !endif
313
314 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere_AD )
315 if ( Error_Status /= 0 ) THEN
316 call da_error(__FILE__,__LINE__, &
317 (/"Error in deallocatting CRTM Atmosphere_AD Structure"/))
318 endif
319
320 Error_Status = CRTM_Destroy_Surface(Surface_AD)
321 if ( Error_Status /= 0 ) THEN
322 call da_error(__FILE__,__LINE__, &
323 (/"Error in deallocatting CRTM Surface_AD Structure"/))
324 endif
325
326 end do ! end loop for pixels
327
328 !-------------------------------------------------------------------
329 ! [2.0] Deallocating CRTM structures
330 !-------------------------------------------------------------------
331 deallocate( rtsolution, rtsolution_ad, STAT = Allocate_Status )
332 IF ( Allocate_Status /= 0 ) THEN
333 call da_error(__FILE__,__LINE__, &
334 (/"Error in deallocatting RTSolution"/))
335 END IF
336
337 Error_Status = CRTM_Destroy_Surface(Surface)
338 if ( Error_Status /= 0 ) THEN
339 call da_error(__FILE__,__LINE__, &
340 (/"Error in deallocatting CRTM Surface Structure"/))
341 endif
342
343 end do ! end loop for sensor
344
345 !-------------------------------------------------------------------
346 ! [3.0] Deallocating CRTM Atmosphere structures
347 !-------------------------------------------------------------------
348 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere )
349 if ( Error_Status /= 0 ) THEN
350 call da_error(__FILE__,__LINE__, &
351 (/"Error in deallocatting CRTM Atmosphere Structure"/))
352 endif
353
354 if (trace_use) call da_trace_exit("da_transform_xtoy_crtm_adj")
355 #endif
356
357 end subroutine da_transform_xtoy_crtm_adj
358