da_transform_xtoy_crtmk_adj.inc
References to this file elsewhere.
1 subroutine da_transform_xtoy_crtmk_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/16/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 :: l, 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, allocatable :: q_k(:),t_k(:)
28 REAL :: ps_k
29
30 #if defined(CRTM)
31 ! CRTM local varaibles and types
32 INTEGER :: wmo_sensor_id, nchan_emis, Error_Status, Allocate_Status
33 TYPE( CRTM_RTSolution_type ), ALLOCATABLE :: RTSolution(:),RTSolution_K(:)
34 TYPE (CRTM_ChannelInfo_type) :: ChannelInfo
35 TYPE( CRTM_Atmosphere_type ) :: Atmosphere
36 TYPE( CRTM_Surface_type ) :: Surface
37 TYPE( CRTM_Atmosphere_type ), ALLOCATABLE :: Atmosphere_K(:)
38 TYPE( CRTM_Surface_type ), ALLOCATABLE :: Surface_K(:)
39 TYPE( CRTM_GeometryInfo_type ) :: GeometryInfo
40 !TYPE( CRTM_Options_type ) :: Options
41 #endif
42 !---------------------------------------------------------
43
44 #if !defined(CRTM)
45 call da_error(__FILE__,__LINE__, &
46 (/"Must compile with $CRTM option for radiances"/))
47 #else
48 IF ( iv%num_inst < 1 ) return
49
50 if (trace_use) call da_trace_entry("da_transform_xtoy_crtmk_adj")
51
52 !-------------------------------------------------------------------------------
53
54 do inst = 1, iv%num_inst ! loop for sensor
55 if ( iv%instid(inst)%num_rad < 1 ) cycle
56 num_rad = iv%instid(inst)%num_rad
57
58 ! CRTM channel information structure
59 Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(inst),ChannelInfo)
60 if ( Error_Status /= 0) then
61 call da_error(__FILE__,__LINE__, &
62 (/"Error in calling CRTM_Set_ChannelInfo"/))
63 endif
64 nchanl = ChannelInfo%n_channels
65
66 ! Allocate forward model solution RTSolution array to number of channels
67 ALLOCATE (RTSolution(ChannelInfo%n_Channels))
68 ALLOCATE (RTSolution_K(ChannelInfo%n_Channels))
69 ALLOCATE (Atmosphere_K(ChannelInfo%n_Channels))
70 ALLOCATE (Surface_K(ChannelInfo%n_Channels))
71 ALLOCATE (t_k(xp%kte-xp%kts+1))
72 ALLOCATE (q_k(xp%kte-xp%kts+1))
73
74 !----------------------------------------------------------------------------
75 ! CRTM Atmosphere structure
76 Atmosphere%n_Layers=(xp%kte-xp%kts)+1 ! number of vertical levels
77 Atmosphere%n_Absorbers=2
78 Atmosphere%n_Clouds=0
79 Atmosphere%n_Aerosols=0
80
81 Error_Status = CRTM_Allocate_Atmosphere( Atmosphere%n_Layers, &
82 Atmosphere%n_Absorbers, &
83 Atmosphere%n_Clouds, &
84 Atmosphere%n_Aerosols, &
85 Atmosphere)
86 if ( Error_Status /= 0 ) THEN
87 call da_error(__FILE__,__LINE__, &
88 (/"Error in allocatting CRTM Atmosphere Structure"/))
89 endif
90 Atmosphere%Absorber_ID(1)=H2O_ID
91 Atmosphere%Absorber_ID(2)=O3_ID
92
93 ! CRTM Surface Structure
94 if (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsua') then
95 nchan_emis=4
96 wmo_sensor_id=WMO_AMSUA
97 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
98 nchan_emis=2
99 wmo_sensor_id=WMO_AMSUB
100 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
101 nchan_emis=12
102 wmo_sensor_id=WMO_AMSRE
103 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
104 nchan_emis=7
105 wmo_sensor_id=WMO_SSMI
106 else
107 nchan_emis=0
108 wmo_sensor_id=INVALID_WMO_SENSOR_ID
109 endif
110
111 Error_Status = CRTM_Allocate_Surface( nchan_emis, & ! Input
112 Surface) ! Output
113 if ( Error_Status /= 0 ) THEN
114 call da_error(__FILE__,__LINE__, &
115 (/"Error in allocatting CRTM Surface Structure"/))
116 endif
117
118 do n=1,num_rad
119
120 ! [1.0] Extract base state Atmosphere variables
121 atmosphere%level_pressure(0) = iv%instid(inst)%pf(0,n)
122 do k=1,atmosphere%n_layers
123 atmosphere%pressure(k) = iv%instid(inst)%pm(k,n)
124 atmosphere%level_pressure(k) = iv%instid(inst)%pf(k,n)
125 atmosphere%temperature(k) = iv%instid(inst)%tm(k,n)
126 atmosphere%absorber(k,1) = iv%instid(inst)%qm(k,n)
127 enddo
128
129 ! [1.1] Extract base Surface parameter data
130 Surface%Land_Coverage=iv%instid(inst)%land_coverage(n)
131 Surface%Water_Coverage=iv%instid(inst)%water_coverage(n)
132 Surface%Snow_Coverage=iv%instid(inst)%snow_coverage(n)
133 Surface%Ice_Coverage=iv%instid(inst)%ice_coverage(n)
134
135 if (Surface%Land_Coverage > 0.) then
136 Surface%Land_Type=GRASS_SOIL ! land type (User guide appendix 3)
137 Surface%Land_Temperature=iv%instid(inst)%ts(n) ! K
138 Surface%Soil_Moisture_Content=iv%instid(inst)%smois(n) !0.05 ! volumetric water content (g/cm**3)
139 !Surface%Canopy_Water_Content=0.05 ! gravimetric water content
140 Surface%Vegetation_Fraction=iv%instid(inst)%vegtyp(n)
141 Surface%Soil_Temperature=iv%instid(inst)%tslb(n)
142 endif
143 if (Surface%Water_Coverage > 0.) then
144 !Surface%Water_Type=SEA_WATER ! (Currently NOT used)
145 Surface%Water_Temperature=iv%instid(inst)%ts(n) ! K
146 Surface%Wind_Speed=sqrt((iv%instid(inst)%u10(n))**2+ &
147 (iv%instid(inst)%v10(n))**2) ! m/sec
148 !surface%Wind_Direction=0.0 ! NOT used
149 Surface%Salinity=33. ! ppmv
150 endif
151 if (Surface%Snow_Coverage > 0.) then
152 Surface%Snow_Type=NEW_SNOW ! User guide appendix 3
153 Surface%Snow_Temperature=iv%instid(inst)%ts(n) ! K
154 Surface%Snow_Depth=iv%instid(inst)%snowh(n) ! mm
155 !Surface%Snow_Density=0.2 ! g/cm**3
156 !Surface%Snow_Grain_Size=2.0 ! mm
157 endif
158 if (Surface%Ice_Coverage > 0.) then
159 !Surface%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
160 Surface%Ice_Temperature=iv%instid(inst)%ts(n) ! K
161 Surface%Ice_Thickness=10. ! mm
162 !Surface%Ice_Density=0.9 ! g/cm**3
163 !Surface%Ice_Roughness=0. ! NO Table offered, single example is ZERO
164 endif
165 Surface%SensorData%n_channels = nchan_emis
166 Surface%SensorData%Sensor_ID = wmo_sensor_id
167 Surface%SensorData%Tb(1:nchan_emis) = iv%instid(inst)%tb_inv(1:nchan_emis,n) + &
168 iv%instid(inst)%tb_xb(1:nchan_emis,n)
169
170 do l = 1, ChannelInfo%n_Channels
171 ! -- Copy the adjoint atmosphere structure
172 Error_Status = CRTM_Assign_Atmosphere( Atmosphere, Atmosphere_K(l) )
173
174 IF ( Error_Status /= 0 ) THEN
175 call da_error(__FILE__,__LINE__, &
176 (/"Error copying Atmosphere_K structure"/))
177 END IF
178
179 ! -- Copy the adjoint surface structure
180 Error_Status = CRTM_Assign_Surface( Surface, Surface_K(l) )
181
182 IF ( Error_Status /= 0 ) THEN
183 call da_error(__FILE__,__LINE__, &
184 (/"Error copying Surface_K structure"/))
185 END IF
186 end do
187
188 ! -- Zero the Adjoint outputs
189 ! Important: adjoint variables must be initialized
190 CALL CRTM_Zero_Atmosphere( Atmosphere_K(:) )
191 CALL CRTM_Zero_Surface( Surface_K(:) )
192
193 ! [1.2] CRTM GeometryInfo Structure
194 GeometryInfo%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
195 GeometryInfo%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
196
197 ! [1.3] assign tb = R^-1 Re :
198
199 RTSolution_K(:)%brightness_temperature = jo_grad_y%instid(inst)%tb(:,n)
200 RTSolution_K(:)%radiance = 0.
201
202 ! [1.4] Call CRTM_K model
203
204 Error_Status = CRTM_K_Matrix(Atmosphere, &
205 Surface, &
206 RTSolution_K,&
207 GeometryInfo, &
208 ChannelInfo, &
209 Atmosphere_K,&
210 Surface_K, &
211 RTSolution) !, &
212 !Options = Options)
213 if ( Error_Status /= 0 ) THEN
214 call da_error(__FILE__,__LINE__, &
215 (/"Error in calling CRTM_Adjoint"/))
216 endif
217
218 ! [1.5] Scale transformation and fill zero for no-control variable
219 ps_k = 0.
220 do l=1, ChannelInfo%n_Channels
221 ps_k = ps_k + Atmosphere_K(l)%Level_Pressure(atmosphere%n_layers)
222 end do
223 ps_k = ps_k*0.01
224
225 do k=1,atmosphere%n_layers
226 t_k(k) = 0.
227 q_k(k) = 0.
228 do l=1, ChannelInfo%n_Channels
229 t_k(k) = t_k(k) + Atmosphere_K(l)%Temperature(k)
230 q_k(k) = q_k(k) + Atmosphere_K(l)%Absorber(k,1)
231 end do
232 q_k(k) = 1000.*q_k(k)
233 end do
234
235 !-----------------------------------------------------
236 ! [1.6] Get horizontal interpolation weights:
237 !-----------------------------------------------------
238
239 i = iv%instid(inst)%loc_i(n)
240 j = iv%instid(inst)%loc_j(n)
241 dx = iv%instid(inst)%loc_dx(n)
242 dy = iv%instid(inst)%loc_dy(n)
243 dxm = iv%instid(inst)%loc_dxm(n)
244 dym = iv%instid(inst)%loc_dym(n)
245
246 ! [1.7] Adjoint of Interpolate horizontally from ob to grid:
247
248 do k=xp%kts,xp%kte ! from bottem to top
249 if (atmosphere%pressure(xp%kte-k+1) < 75.) q_k(xp%kte-k+1) = 0.
250 call da_interp_lin_2d_adj( jo_grad_x%t(:,:,k), xp%ims, xp%ime, xp%jms, &
251 xp%jme, i, j, dx, dy, dxm, dym, &
252 t_k(xp%kte-k+1))
253
254 call da_interp_lin_2d_adj( jo_grad_x%q(:,:,k), xp%ims, xp%ime, xp%jms, &
255 xp%jme, i, j, dx, dy, dxm, dym, &
256 q_k(xp%kte-k+1) )
257 enddo
258
259 call da_interp_lin_2d_adj(jo_grad_x% psfc, xp%ims, xp%ime, xp%jms, &
260 xp%jme, i, j, dx, dy, dxm, dym, &
261 ps_k )
262
263 !if (n <=10 ) then
264 ! write(6,'(15f8.2)') rtsolution(:)%brightness_temperature
265 ! write(6,'(e15.5)') jo_grad_x% psfc(i,j)
266 ! do k=xp%kts,xp%kte
267 ! write(6,'(2e15.5)') jo_grad_x%t(i,j,k), jo_grad_x%q(i,j,k)
268 ! enddo
269 !endif
270
271 end do ! end loop for pixels
272
273 !-------------------------------------------------------------------
274 ! [2.0] Deallocating CRTM structures
275 !-------------------------------------------------------------------
276 Error_Status = CRTM_Destroy_Surface(Surface)
277 if ( Error_Status /= 0 ) THEN
278 call da_error(__FILE__,__LINE__, &
279 (/"Error in deallocatting CRTM Surface Structure"/))
280 endif
281
282 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere )
283 if ( Error_Status /= 0 ) THEN
284 call da_error(__FILE__,__LINE__, &
285 (/"Error in deallocatting CRTM Atmosphere Structure"/))
286 endif
287
288 Error_Status = CRTM_Destroy_Surface(Surface_K)
289 if ( Error_Status /= 0 ) THEN
290 call da_error(__FILE__,__LINE__, &
291 (/"Error in deallocatting CRTM Surface_K Structure"/))
292 endif
293
294 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere_K )
295 if ( Error_Status /= 0 ) THEN
296 call da_error(__FILE__,__LINE__, &
297 (/"Error in deallocatting CRTM Atmosphere_K Structure"/))
298 endif
299
300 deallocate( RTSolution, RTSolution_K, Atmosphere_K, Surface_K, t_k, q_k, &
301 STAT = Allocate_Status )
302 if ( Allocate_Status /= 0 ) THEN
303 call da_error(__FILE__,__LINE__, &
304 (/"Error in deallocatting CRTM Surface Structure"/))
305 endif
306
307 end do ! end loop for sensor
308
309 if (trace_use) call da_trace_exit("da_transform_xtoy_crtmk_adj")
310 #endif
311
312 end subroutine da_transform_xtoy_crtmk_adj
313