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