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