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