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 :: i, j, k ! Index dimension.
24 REAL :: dx, dxm ! Interpolation weights.
25 REAL :: dy, dym ! Interpolation weights.
26
27 INTEGER :: inst, num_rad, nchanl, n
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_TL(:)
33 TYPE (CRTM_ChannelInfo_type) :: ChannelInfo
34 TYPE( CRTM_Atmosphere_type ) :: Atmosphere, Atmosphere_TL
35 TYPE( CRTM_Surface_type ) :: Surface, Surface_TL
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
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%n_Layers=(xp%kte-xp%kts)+1 ! number of vertical levels
55 Atmosphere%n_Absorbers=2
56 Atmosphere%n_Clouds=0
57 Atmosphere%n_Aerosols=0
58 if (crtm_cloud) Atmosphere%n_Clouds=6
59
60 Error_Status = CRTM_Allocate_Atmosphere( Atmosphere%n_Layers, &
61 Atmosphere%n_Absorbers, &
62 Atmosphere%n_Clouds, &
63 Atmosphere%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%Absorber_ID(1)=H2O_ID
71 Atmosphere%Absorber_ID(2)=O3_ID
72
73 if (crtm_cloud) then
74 Atmosphere%Cloud(1)%Type=WATER_CLOUD
75 Atmosphere%Cloud(2)%Type=ICE_CLOUD
76 Atmosphere%Cloud(3)%Type=RAIN_CLOUD
77 Atmosphere%Cloud(4)%Type=SNOW_CLOUD
78 Atmosphere%Cloud(5)%Type=GRAUPEL_CLOUD
79 Atmosphere%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%n_channels
95
96 ! Allocate forward model solution RTSolution array to number of channels
97 ALLOCATE( RTSolution( ChannelInfo%n_Channels ), &
98 RTSolution_TL( ChannelInfo%n_Channels ), &
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 nchan_emis=4
108 wmo_sensor_id=WMO_AMSUA
109 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
110 nchan_emis=2
111 wmo_sensor_id=WMO_AMSUB
112 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
113 nchan_emis=12
114 wmo_sensor_id=WMO_AMSRE
115 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
116 nchan_emis=7
117 wmo_sensor_id=WMO_SSMI
118 else
119 nchan_emis=0
120 wmo_sensor_id=INVALID_WMO_SENSOR_ID
121 endif
122
123 Error_Status = CRTM_Allocate_Surface( nchan_emis, & ! Input
124 Surface) ! Output
125 if ( Error_Status /= 0 ) THEN
126 call da_error(__FILE__,__LINE__, &
127 (/"Error in allocatting CRTM Surface Structure"/))
128 endif
129
130 ! CRTM Options structure
131 !Options%n_channels = nchanl
132 !Error_Status = CRTM_Allocate_Options( nchanl, & ! Input
133 ! Options) ! Output
134 !if ( Error_Status /= 0 ) THEN
135 ! call da_error(__FILE__,__LINE__, &
136 ! (/"Error in allocatting CRTM Options Structure"/))
137 !endif
138
139 do n= 1, num_rad ! loop for pixel
140
141 ! [1.1] Get horizontal interpolation weights:
142
143 i = iv%instid(inst)%loc_i(n)
144 j = iv%instid(inst)%loc_j(n)
145 dx = iv%instid(inst)%loc_dx(n)
146 dy = iv%instid(inst)%loc_dy(n)
147 dxm = iv%instid(inst)%loc_dxm(n)
148 dym = iv%instid(inst)%loc_dym(n)
149
150 ! [1.3] Extract base state Atmosphere variables
151 atmosphere%level_pressure(0) = iv%instid(inst)%pf(0,n)
152 do k=1,atmosphere%n_layers
153 atmosphere%pressure(k) = iv%instid(inst)%pm(k,n)
154 atmosphere%level_pressure(k) = iv%instid(inst)%pf(k,n)
155 atmosphere%temperature(k) = iv%instid(inst)%tm(k,n)
156 atmosphere%absorber(k,1) = iv%instid(inst)%qm(k,n)
157 if (crtm_cloud) then
158 atmosphere%cloud(1)%water_content(k)=iv%instid(inst)%qcw(k,n)
159 atmosphere%cloud(2)%water_content(k)=iv%instid(inst)%qci(k,n)
160 atmosphere%cloud(3)%water_content(k)=iv%instid(inst)%qrn(k,n)
161 atmosphere%cloud(4)%water_content(k)=iv%instid(inst)%qsn(k,n)
162 atmosphere%cloud(5)%water_content(k)=iv%instid(inst)%qgr(k,n)
163 atmosphere%cloud(6)%water_content(k)=0.
164 atmosphere%cloud(1)%effective_radius(k)=10.
165 atmosphere%cloud(2)%effective_radius(k)=200.
166 atmosphere%cloud(3)%effective_radius(k)=200.
167 atmosphere%cloud(4)%effective_radius(k)=200.
168 atmosphere%cloud(5)%effective_radius(k)=200.
169 atmosphere%cloud(6)%effective_radius(k)=200.
170 end if
171 enddo
172
173 ! [1.4] User-supplied emissivity
174 !Options%emissivity_switch = 1
175 !Options%emissivity(1:Options%n_channels) = &
176 ! iv%instid(inst)%emiss(1:Options%n_channels,n)
177
178 ! [1.4] CRTM Surface parameter data
179 Surface%Land_Coverage=iv%instid(inst)%land_coverage(n)
180 Surface%Water_Coverage=iv%instid(inst)%water_coverage(n)
181 Surface%Snow_Coverage=iv%instid(inst)%snow_coverage(n)
182 Surface%Ice_Coverage=iv%instid(inst)%ice_coverage(n)
183
184 if (Surface%Land_Coverage > 0.) then
185 Surface%Land_Type=GRASS_SOIL ! land type (User guide appendix 3)
186 Surface%Land_Temperature=iv%instid(inst)%ts(n) ! K
187 Surface%Soil_Moisture_Content=iv%instid(inst)%smois(n) !0.05 ! volumetric water content (g/cm**3)
188 !Surface%Canopy_Water_Content=0.05 ! gravimetric water content
189 Surface%Vegetation_Fraction=iv%instid(inst)%vegtyp(n)
190 Surface%Soil_Temperature=iv%instid(inst)%tslb(n)
191 endif
192 if (Surface%Water_Coverage > 0.) then
193 !Surface%Water_Type=SEA_WATER ! (Currently NOT used)
194 Surface%Water_Temperature=iv%instid(inst)%ts(n) ! K
195 Surface%Wind_Speed=sqrt((iv%instid(inst)%u10(n))**2+ &
196 (iv%instid(inst)%v10(n))**2) ! m/sec
197 !surface%Wind_Direction=0.0 ! NOT used
198 Surface%Salinity=33. ! ppmv
199 endif
200 if (Surface%Snow_Coverage > 0.) then
201 Surface%Snow_Type=NEW_SNOW ! User guide appendix 3
202 Surface%Snow_Temperature=iv%instid(inst)%ts(n) ! K
203 Surface%Snow_Depth=iv%instid(inst)%snowh(n) ! mm
204 !Surface%Snow_Density=0.2 ! g/cm**3
205 !Surface%Snow_Grain_Size=2.0 ! mm
206 endif
207 if (Surface%Ice_Coverage > 0.) then
208 !Surface%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
209 Surface%Ice_Temperature=iv%instid(inst)%ts(n) ! K
210 Surface%Ice_Thickness=10. ! mm
211 !Surface%Ice_Density=0.9 ! g/cm**3
212 !Surface%Ice_Roughness=0. ! NO Table offered, single example is ZERO
213 endif
214 Surface%SensorData%n_channels = nchan_emis
215 Surface%SensorData%Sensor_ID = wmo_sensor_id
216 Surface%SensorData%Tb(1:nchan_emis) = iv%instid(inst)%tb_inv(1:nchan_emis,n) + &
217 iv%instid(inst)%tb_xb(1:nchan_emis,n)
218
219 ! -- Copy the TL atmosphere structure
220 Error_Status = CRTM_Assign_Atmosphere( Atmosphere, Atmosphere_TL )
221
222 IF ( Error_Status /= 0 ) THEN
223 call da_error(__FILE__,__LINE__, &
224 (/"Error copying Atmosphere_TL structure"/))
225 END IF
226
227 ! -- Copy the TL surface structure
228 Error_Status = CRTM_Assign_Surface( Surface, Surface_TL )
229
230 IF ( Error_Status /= 0 ) THEN
231 call da_error(__FILE__,__LINE__, &
232 (/"Error copying Surface_TL structure"/))
233 END IF
234
235 ! -- Zero the TL outputs
236 ! Important: adjoint variables must be initialized
237 CALL CRTM_Zero_Atmosphere( Atmosphere_TL )
238 CALL CRTM_Zero_Surface( Surface_TL )
239
240 ! [1.2] Interpolate horizontally Atmoshere_TL variables to ob:
241 do k=xp%kts,xp%kte ! from bottem to top
242 call da_interp_lin_2d( xa%t(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
243 i, j, dx, dy, dxm, dym, & ! temperature (in K)
244 Atmosphere_TL%Temperature(xp%kte-k+1) )
245 call da_interp_lin_2d( xa%q(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
246 i, j, dx, dy, dxm, dym, & ! specific humidity (in kg/kg)
247 Atmosphere_TL%Absorber(xp%kte-k+1,1) )
248 if (atmosphere%pressure(xp%kte-k+1) < 75.) Atmosphere_TL%Absorber(xp%kte-k+1,1) = 0.
249 end do
250 call da_interp_lin_2d( xa%psfc, xp%ims, xp%ime, xp%jms, xp%jme, &
251 i, j, dx, dy, dxm, dym, & ! in Pa
252 Atmosphere_TL%Level_Pressure(Atmosphere_TL%n_Layers) )
253
254 Atmosphere_TL%Level_Pressure(Atmosphere%n_Layers)= &
255 0.01*Atmosphere_TL%Level_Pressure(Atmosphere%n_Layers)
256 Atmosphere_TL%Absorber(:,1) = 1000.*Atmosphere_TL%Absorber(:,1) ! in g/kg
257
258 ! [1.5] CRTM GeometryInfo Structure
259 GeometryInfo%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
260 GeometryInfo%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
261 ! GeometryInfo%Satellite_Height=830.
262 ! GeometryInfo%Sensor_Scan_Angle=
263 ! GeometryInfo%Sensor_Zenith_Angle=
264 ! GeometryInfo%Sensor_Scan_Angle=
265 ! GeometryInfo%Source_Zenith_Angle=
266
267 ! [1.6] Call CRTM_TL model
268
269 call da_crtm_tl (nchanl, Atmosphere, &
270 Surface, &
271 Atmosphere_TL,&
272 Surface_TL, &
273 GeometryInfo, &
274 ChannelInfo, &
275 RTSolution, &
276 RTSolution_TL)
277
278 !Error_Status = CRTM_Tangent_Linear(Atmosphere, &
279 ! Surface, &
280 ! Atmosphere_TL,&
281 ! Surface_TL, &
282 ! GeometryInfo, &
283 ! ChannelInfo, &
284 ! RTSolution, &
285 ! RTSolution_TL) !,&
286 ! !Options = Options)
287 !if ( Error_Status /= 0 ) THEN
288 ! call da_error(__FILE__,__LINE__, &
289 ! (/"Error in calling CRTM_Tangent_Linear"/))
290 !endif
291
292 !-------------------------------------------------------------------
293 ! [1.7] assign Hdx :
294 !-------------------------------------------------------------------
295 y%instid(inst)%tb(:,n) = RTSolution_TL(:)%brightness_temperature
296 !if (n <= 10 ) write(6,'(15f8.3)') RTSolution_TL(:)%brightness_temperature
297
298 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere_TL )
299 if ( Error_Status /= 0 ) THEN
300 call da_error(__FILE__,__LINE__, &
301 (/"Error in deallocatting CRTM Atmosphere_TL Structure"/))
302 endif
303
304 Error_Status = CRTM_Destroy_Surface(Surface_TL)
305 if ( Error_Status /= 0 ) THEN
306 call da_error(__FILE__,__LINE__, &
307 (/"Error in deallocatting CRTM Surface_TL Structure"/))
308 endif
309
310 end do ! end loop for pixels
311 !-------------------------------------------------------------------
312 ! [2.0] Deallocating CRTM structures
313 !-------------------------------------------------------------------
314 deallocate( rtsolution, rtsolution_tl, STAT = Allocate_Status )
315 IF ( Allocate_Status /= 0 ) THEN
316 call da_error(__FILE__,__LINE__, &
317 (/"Error in deallocatting RTSolution"/))
318 END IF
319
320 Error_Status = CRTM_Destroy_Surface(Surface)
321 if ( Error_Status /= 0 ) THEN
322 call da_error(__FILE__,__LINE__, &
323 (/"Error in deallocatting CRTM Surface Structure"/))
324 endif
325
326 end do ! end loop for sensor
327
328 !-------------------------------------------------------------------
329 ! [3.0] Deallocating CRTM Atmosphere structures
330 !-------------------------------------------------------------------
331 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere )
332 if ( Error_Status /= 0 ) THEN
333 call da_error(__FILE__,__LINE__, &
334 (/"Error in deallocatting CRTM Atmosphere Structure"/))
335 endif
336
337 if (trace_use) call da_trace_exit("da_transform_xtoy_crtm")
338 #endif
339
340 end subroutine da_transform_xtoy_crtm
341