1module slantProfileLatLon_mod
2 ! MODULE slantProfileLatLon_mod (prefix='slp' category='5. Observation operators')
3 !
4 !:Purpose: To calculate latitudes/longitudes on slant-path based on
5 ! ColumnData.
6 !
7 use midasMpi_mod
8 use earthConstants_mod
9 use mathPhysConstants_mod
10 use utilities_mod
11 use obsSpaceData_mod
12 use horizontalCoord_mod
13 use tovsNL_mod
14 use codtyp_mod
15 use getGridPosition_mod
16 use radialVelocity_mod
17
18 implicit none
19 save
20 private
21
22 ! public procedures
23 public :: slp_calcLatLonTovs, slp_calcLatLonRO, slp_calcLatLonRadar
24
25 ! private module variables and derived types
26 logical :: nmlAlreadyRead = .false.
27
28 ! namelist variables
29 real(4) :: toleranceHeightDiff ! threshold of height diff (in m) for convergence of slant path
30 integer :: maxNumIteration ! maximum number of iterations for determining slant path
31
32
33contains
34
35 subroutine slp_calcLatLonTovs(obsSpaceData, hco, headerIndex, height3D_T_r4, height3D_M_r4, &
36 latSlantLev_T, lonSlantLev_T, latSlantLev_M, lonSlantLev_M, latSlantLev_S, lonSlantLev_S)
37 !
38 !:Purpose: call the computation of lat/lon on the slant path for radiance
39 ! observations, iteratively. To replace the vertical columns with
40 ! line-of-sight slanted columns.
41 !
42 implicit none
43
44 ! Arguments:
45 type(struct_obs), intent(in) :: obsSpaceData
46 type(struct_hco), intent(in) :: hco
47 integer, intent(in) :: headerIndex
48 real(4), intent(in) :: height3D_T_r4(:,:,:)
49 real(4), intent(in) :: height3D_M_r4(:,:,:)
50 real(8), intent(out) :: latSlantLev_T(:)
51 real(8), intent(out) :: lonSlantLev_T(:)
52 real(8), intent(out) :: latSlantLev_M(:)
53 real(8), intent(out) :: lonSlantLev_M(:)
54 real(8), intent(out) :: latSlantLev_S
55 real(8), intent(out) :: lonSlantLev_S
56
57 ! Locals:
58 real(8) :: lat, lon, azimuthAngle
59 integer :: idatyp
60 integer :: ierr, fnom, fclos, nulnam
61 integer :: nlev_T, lev_T, nlev_M, lev_M
62
63 namelist /namSlantPath/ toleranceHeightDiff, maxNumIteration
64
65 if ( .not. nmlAlreadyRead ) then
66 nmlAlreadyRead = .true.
67
68 ! default values
69 toleranceHeightDiff = 10.0
70 maxNumIteration = 1
71
72 ! reading namelist variables
73 nulnam = 0
74 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
75 read(nulnam, nml = namSlantPath, iostat = ierr)
76 if (ierr /= 0) write(*,*) 'slp_calcLatLonTovs: namSlantPath is missing in the namelist. The default value will be taken.'
77 if (mmpi_myid == 0) write(*, nml = namSlantPath)
78 ierr = fclos(nulnam)
79 end if
80
81 nlev_M = size(height3D_M_r4,3)
82 nlev_T = size(height3D_T_r4,3)
83
84 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
85 lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
86 if (lon < 0.0d0 ) lon = lon + 2.0d0*MPC_PI_R8
87 if (lon >= 2.0d0*MPC_PI_R8) lon = lon - 2.0d0*MPC_PI_R8
88 azimuthAngle = tvs_getCorrectedSatelliteAzimuth(obsSpaceData, headerIndex)
89
90 !SSMIS special case
91 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
92 if ( idatyp == codtyp_get_codtyp('ssmis') ) then
93 azimuthAngle = obs_missingValue_R
94 end if
95
96 ! loop through all thermo levels
97 do lev_T = 1, nlev_T
98 call solveIterativelyForLatLon(height3D_T_r4(:,:,lev_T), latSlantLev_T(lev_T), lonSlantLev_T(lev_T))
99 end do
100
101 ! loop through all momentum levels
102 do lev_M = 1, nlev_M
103 call solveIterativelyForLatLon(height3D_M_r4(:,:,lev_M), latSlantLev_M(lev_M), lonSlantLev_M(lev_M))
104 end do
105
106 ! then surface level (geoid in the case of radiances)
107 latSlantLev_S = latSlantLev_T(nlev_T)
108 lonSlantLev_S = lonSlantLev_T(nlev_T)
109
110 contains
111
112 subroutine solveIterativelyForLatLon(heightField2D, latSlant, lonSlant)
113 implicit none
114
115 ! Arguments:
116 real(4), intent(in) :: heightField2D(:,:)
117 real(8), intent(out) :: latSlant
118 real(8), intent(out) :: lonSlant
119
120 ! Locals:
121 logical :: doIteration
122 integer :: numIteration
123 real(4) :: heightInterp_r4, heightIntersect_r4, heightDiff_r4
124
125 ! find the interpolated height
126 call heightBilinearInterp(lat, lon, hco, heightField2D, heightInterp_r4)
127
128 doIteration = .true.
129 numIteration = 0
130 while_doIteration: do while (doIteration)
131
132 numIteration = numIteration + 1
133
134 call findIntersectLatlon(obsSpaceData, headerIndex, heightInterp_r4, azimuthAngle, latSlant, lonSlant)
135
136 ! find the interpolated height
137 call heightBilinearInterp(latSlant, lonSlant, hco, heightField2D, heightIntersect_r4)
138
139 heightDiff_r4 = abs(heightInterp_r4-heightIntersect_r4)
140 if ( heightDiff_r4 < toleranceHeightDiff .or. &
141 numIteration >= maxNumIteration ) doIteration = .false.
142
143 end do while_doIteration
144
145 end subroutine solveIterativelyForLatLon
146
147
148 end subroutine slp_calcLatLonTovs
149
150
151 subroutine findIntersectLatlon(obsSpaceData, headerIndex, height_r4, azimuthAngle, latSlant, lonSlant)
152 !
153 !:Purpose: Computation of lat/lon of the intersection between model level
154 ! and the slant line-of-sight for radiance observations.
155 !
156 implicit none
157
158 ! Arguments:
159 type(struct_obs), intent(in) :: obsSpaceData
160 real(4), intent(in) :: height_r4
161 integer, intent(in) :: headerIndex
162 real(8), intent(in) :: azimuthAngle
163 real(8), intent(out) :: latSlant
164 real(8), intent(out) :: lonSlant
165
166 ! Locals:
167 real(8) :: lat, lon, geometricHeight
168 real(8) :: zenithAngle, zenithAngle_rad, azimuthAngle_rad, elevationAngle_rad, distAlongPath
169 real(8) :: obsCordGlb(3), slantPathCordGlb(3), unitx(3), unity(3), unitz(3), unitSatLoc(3), unitSatGlb(3)
170
171 ! read lat/lon/angles from obsSpaceData
172 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
173 lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
174 if (lon < 0.0d0 ) lon = lon + 2.0d0*MPC_PI_R8
175 if (lon >= 2.0d0*MPC_PI_R8) lon = lon - 2.0d0*MPC_PI_R8
176
177 if ( azimuthAngle /= obs_missingValue_R) then
178
179 zenithAngle = obs_headElem_r(obsSpaceData,OBS_SZA,headerIndex)
180
181 ! convert angles to radian unit
182 azimuthAngle_rad = azimuthAngle * MPC_RADIANS_PER_DEGREE_R8
183 zenithAngle_rad = zenithAngle * MPC_RADIANS_PER_DEGREE_R8
184 elevationAngle_rad = 0.5d0 * MPC_PI_R8 - zenithAngle_rad
185
186 obsCordGlb = ec_ra * (/ cos(lat)*cos(lon), cos(lat)*sin(lon), sin(lat) /)
187 unitz = (/ cos(lat)*cos(lon), cos(lat)*sin(lon) , sin(lat) /)
188 unitx = (/ -sin(lon) , cos(lon) , 0.d0 /)
189 unity = (/ -sin(lat)*cos(lon), -sin(lat)*sin(lon), cos(lat) /)
190
191 ! unit vector towards satellite in local coordinate
192 unitSatLoc = (/ cos(elevationAngle_rad)*sin(azimuthAngle_rad) , &
193 cos(elevationAngle_rad)*cos(azimuthAngle_rad) , &
194 sin(elevationAngle_rad) /)
195 ! unit vector towards satellite in global coordinate
196 unitSatGlb = unitSatLoc(1) * unitx + unitSatLoc(2) * unity + unitSatLoc(3) * unitz
197
198 ! Geometric altitude
199 geometricHeight = real(height_r4,8)
200
201 ! distance along line of sight
202 distAlongPath = geometricHeight / cos(zenithAngle_rad)
203
204 slantPathCordGlb(:) = obsCordGlb(:) + distAlongPath * unitSatGlb(:)
205
206 latSlant = atan(slantPathCordGlb(3)/sqrt(slantPathCordGlb(1)**2+slantPathCordGlb(2)**2))
207 lonSlant = atan2(slantPathCordGlb(2),slantPathCordGlb(1))
208 if (lonSlant < 0.0d0 ) lonSlant = lonSlant + 2.0d0*MPC_PI_R8
209 if (lonSlant >= 2.0d0*MPC_PI_R8) lonSlant = lonSlant - 2.0d0*MPC_PI_R8
210
211 else
212
213 latSlant = lat
214 lonSlant = lon
215
216 end if
217
218 end subroutine findIntersectLatlon
219
220
221 subroutine heightBilinearInterp(lat, lon, hco, height_r4, heightInterp_r4)
222 !
223 !:Purpose: To interpolate the 2D height field to the obs location
224 !
225 implicit none
226
227 ! Arguments:
228 real(8), intent(in) :: lat
229 real(8), intent(in) :: lon
230 type(struct_hco), intent(in) :: hco
231 real(4), intent(in) :: height_r4(:,:)
232 real(4), intent(out) :: heightInterp_r4
233
234 ! Locals:
235 integer :: ierr, niP1
236 integer :: latIndex, lonIndex, latIndex2, lonIndex2, lonIndexP1
237 integer :: subGridIndex, subGridForInterp, numSubGridsForInterp
238 integer :: ipoint, gridptCount
239 integer :: latIndexVec(8), lonIndexVec(8)
240 real(8) :: WeightVec(8)
241 real(8) :: dldx, dldy
242 real(8) :: weightsSum
243 real(4) :: lon_deg_r4, lat_deg_r4
244 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
245 logical :: latlonOutsideGrid
246
247 lat_deg_r4 = real(lat * MPC_DEGREES_PER_RADIAN_R8)
248 lon_deg_r4 = real(lon * MPC_DEGREES_PER_RADIAN_R8)
249 ierr = gpos_getPositionXY( hco%EZscintID, &
250 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
251 lat_deg_r4, lon_deg_r4, subGridIndex )
252
253 ! Allow for periodicity in Longitude for global Gaussian grid
254 if ( hco%grtyp == 'G' .or. (hco%grtyp == 'Z' .and. hco%global) ) then
255 niP1 = hco%ni + 1
256 else
257 niP1 = hco%ni
258 end if
259
260 ! Check if the interpolated lat/lon is outside the hco domain
261 latlonOutsideGrid = ( xpos_r4 < 1.0 .or. &
262 xpos_r4 > real(niP1) .or. &
263 ypos_r4 < 1.0 .or. &
264 ypos_r4 > real(hco%nj) )
265
266 if ( latlonOutsideGrid ) then
267 write(*,*) 'heightBilinearInterp: interpolated lat/lon outside the hco domain.'
268 write(*,*) ' position lon, lat = ', lon_deg_r4, lat_deg_r4
269 write(*,*) ' position x, y = ', xpos_r4, ypos_r4
270
271 ! if above or below domain
272 if( ypos_r4 < 1.0 ) ypos_r4 = 1.0
273 if( ypos_r4 > real(hco%nj) ) ypos_r4 = real(hco%nj)
274
275 ! if on the left or right longitude band, move it to the edge of this longitude band
276 if( xpos_r4 < 1.0 ) xpos_r4 = 1.0
277 if( xpos_r4 > real(niP1) ) xpos_r4 = real(niP1)
278 write(*,*) ' new position x, y = ', xpos_r4, ypos_r4
279
280 end if
281
282 ! Find the lower-left grid point next to the observation
283 if ( xpos_r4 == real(niP1) ) then
284 lonIndex = floor(xpos_r4) - 1
285 else
286 lonIndex = floor(xpos_r4)
287 end if
288 if ( xpos2_r4 == real(niP1) ) then
289 lonIndex2 = floor(xpos2_r4) - 1
290 else
291 lonIndex2 = floor(xpos2_r4)
292 end if
293
294 if ( ypos_r4 == real(hco%nj) ) then
295 latIndex = floor(ypos_r4) - 1
296 else
297 latIndex = floor(ypos_r4)
298 end if
299 if ( ypos2_r4 == real(hco%nj) ) then
300 latIndex2 = floor(ypos2_r4) - 1
301 else
302 latIndex2 = floor(ypos2_r4)
303 end if
304
305 if ( hco%grtyp == 'U' ) then
306 if ( ypos_r4 == real(hco%nj/2) ) then
307 latIndex = floor(ypos_r4) - 1
308 else
309 latIndex = floor(ypos_r4)
310 end if
311 if ( ypos2_r4 == real(hco%nj/2) ) then
312 latIndex2 = floor(ypos2_r4) - 1
313 else
314 latIndex2 = floor(ypos2_r4)
315 end if
316 end if
317
318 ! Handle periodicity in longitude
319 lonIndexP1 = lonIndex + 1
320 if ( lonIndexP1 == hco%ni + 1 ) lonIndexP1 = 1
321
322 ! Check if location is in between Yin and Yang (should not happen)
323 if ( hco%grtyp == 'U' ) then
324 if ( ypos_r4 > real(hco%nj/2) .and. &
325 ypos_r4 < real((hco%nj/2)+1) ) then
326 write(*,*) 'heightBilinearInterp: WARNING, obs position in between Yin and Yang!!!'
327 write(*,*) ' xpos, ypos = ', xpos_r4, ypos_r4
328 end if
329 if ( ypos2_r4 > real(hco%nj/2) .and. &
330 ypos2_r4 < real((hco%nj/2)+1) ) then
331 write(*,*) 'heightBilinearInterp: WARNING, obs position in between Yin and Yang!!!'
332 write(*,*) ' xpos2, ypos2 = ', xpos2_r4, ypos2_r4
333 end if
334 end if
335
336 if ( subGridIndex == 3 ) then
337 ! both subGrids involved in interpolation, so first treat subGrid 1
338 numSubGridsForInterp = 2
339 subGridIndex = 1
340 else
341 ! only 1 subGrid involved in interpolation
342 numSubGridsForInterp = 1
343 end if
344
345 gridptCount = 0
346
347 do subGridForInterp = 1, numSubGridsForInterp
348
349 ! Compute the 4 weights of the bilinear interpolation
350 if ( subGridForInterp == 1 ) then
351 ! when only 1 subGrid involved, subGridIndex can be 1 or 2
352 dldx = real(xpos_r4,8) - real(lonIndex,8)
353 dldy = real(ypos_r4,8) - real(latIndex,8)
354 else
355 ! when 2 subGrids, subGridIndex is set to 1 for 1st iteration, 2 for second
356 subGridIndex = 2
357 lonIndex = lonIndex2
358 latIndex = latIndex2
359 lonIndexP1 = lonIndex2 + 1
360 dldx = real(xpos2_r4,8) - real(lonIndex,8)
361 dldy = real(ypos2_r4,8) - real(latIndex,8)
362 end if
363
364 gridptCount = gridptCount + 1
365 latIndexVec(gridptCount) = latIndex
366 lonIndexVec(gridptCount) = lonIndex
367 WeightVec(gridptCount) = (1.d0-dldx) * (1.d0-dldy)
368
369 gridptCount = gridptCount + 1
370 latIndexVec(gridptCount) = latIndex
371 lonIndexVec(gridptCount) = lonIndexP1
372 WeightVec(gridptCount) = dldx * (1.d0-dldy)
373
374 gridptCount = gridptCount + 1
375 latIndexVec(gridptCount) = latIndex + 1
376 lonIndexVec(gridptCount) = lonIndex
377 WeightVec(gridptCount) = (1.d0-dldx) * dldy
378
379 gridptCount = gridptCount + 1
380 latIndexVec(gridptCount) = latIndex + 1
381 lonIndexVec(gridptCount) = lonIndexP1
382 WeightVec(gridptCount) = dldx * dldy
383
384 end do ! subGrid
385
386 weightsSum = sum(WeightVec(1:gridptCount))
387 if ( weightsSum > 0.d0 ) then
388
389 WeightVec(1:gridptCount) = WeightVec(1:gridptCount) / weightsSum
390
391 else
392
393 call utl_abort('heightBilinearInterp: weightsSum smaller than 0.')
394
395 end if
396
397 ! perform the bi-linear interpolation
398 heightInterp_r4 = 0.0
399 do ipoint = 1, gridptCount
400
401 if ( latIndexVec(ipoint) > hco%nj ) then
402 write(*,*) 'heightBilinearInterp: latIndexVec(ipoint) > hco%nj: latIndex=',latIndexVec(ipoint),', lonIndex=',lonIndexVec(ipoint)
403 write(*,*) 'lat_deg_r4=',lat_deg_r4,', lon_deg_r4=',lon_deg_r4
404 write(*,*) 'ypos_r4=',ypos_r4,', xpos_r4=',xpos_r4
405 end if
406 if ( latIndexVec(ipoint) < 1 ) then
407 write(*,*) 'heightBilinearInterp: latIndexVec(ipoint) < 1: latIndex=',latIndexVec(ipoint),', lonIndex=',lonIndexVec(ipoint)
408 write(*,*) 'lat_deg_r4=',lat_deg_r4,', lon_deg_r4=',lon_deg_r4
409 write(*,*) 'ypos_r4=',ypos_r4,', xpos_r4=',xpos_r4
410 end if
411
412 if ( lonIndexVec(ipoint) > hco%ni ) then
413 write(*,*) 'heightBilinearInterp: lonIndexVec(ipoint) > hco%ni: latIndex=',latIndexVec(ipoint),', lonIndex=',lonIndexVec(ipoint)
414 write(*,*) 'lat_deg_r4=',lat_deg_r4,', lon_deg_r4=',lon_deg_r4
415 write(*,*) 'ypos_r4=',ypos_r4,', xpos_r4=',xpos_r4
416 end if
417 if ( lonIndexVec(ipoint) < 1 ) then
418 write(*,*) 'heightBilinearInterp: lonIndexVec(ipoint) < 1: latIndex=',latIndexVec(ipoint),', lonIndex=',lonIndexVec(ipoint)
419 write(*,*) 'lat_deg_r4=',lat_deg_r4,', lon_deg_r4=',lon_deg_r4
420 write(*,*) 'ypos_r4=',ypos_r4,', xpos_r4=',xpos_r4
421 end if
422
423 latlonOutsideGrid = ( latIndexVec(ipoint) < 1 .or. &
424 latIndexVec(ipoint) > hco%nj .or. &
425 lonIndexVec(ipoint) < 1 .or. &
426 lonIndexVec(ipoint) > hco%ni )
427
428 if ( latlonOutsideGrid ) &
429 call utl_abort('heightBilinearInterp: lat/lon outside the domain.')
430
431 heightInterp_r4 = heightInterp_r4 + &
432 real(WeightVec(ipoint),4) * &
433 height_r4(lonIndexVec(ipoint), latIndexVec(ipoint))
434 end do
435
436 end subroutine heightBilinearInterp
437
438 subroutine slp_calcLatLonRO(obsSpaceData, hco, headerIndex, &
439 height3D_T_r4, height3D_M_r4, &
440 latSlantLev_T, lonSlantLev_T, latSlantLev_M, lonSlantLev_M, latSlantLev_S, lonSlantLev_S )
441 !
442 !:Purpose: call the computation of lat/lon on the slant path for GPSRO
443 ! observations, iteratively. To replace the vertical columns with
444 ! slanted columns.
445 !
446 implicit none
447
448 ! Arguments:
449 type(struct_obs), intent(inout) :: obsSpaceData
450 type(struct_hco), intent(in) :: hco
451 integer, intent(in) :: headerIndex
452 real(4), intent(in) :: height3D_T_r4(:,:,:)
453 real(4), intent(in) :: height3D_M_r4(:,:,:)
454 real(8), intent(out) :: latSlantLev_T(:)
455 real(8), intent(out) :: lonSlantLev_T(:)
456 real(8), intent(out) :: latSlantLev_M(:)
457 real(8), intent(out) :: lonSlantLev_M(:)
458 real(8), intent(out) :: latSlantLev_S
459 real(8), intent(out) :: lonSlantLev_S
460
461 ! Locals:
462 real(8) :: latr, lonr, height, rad, dH, hmin
463 integer :: bodyIndex, imin
464 real(8), allocatable :: Lat_Obs(:), Lon_Obs(:), Hgt_Obs(:)
465 real(4), allocatable :: H_M(:,:), H_T(:,:)
466 integer :: nObs, iObs
467 integer :: levIndex, nlev_M, nlev_T
468
469 nlev_M = size(height3D_M_r4,3)
470 nlev_T = size(height3D_T_r4,3)
471
472 ! Header, and Obs counting:
473 call obs_set_current_header_list(obsSpaceData,'RO')
474 call obs_set_current_body_list(obsSpaceData, headerIndex)
475 Rad = obs_headElem_r(obsSpaceData,OBS_TRAD,headerIndex)
476 nObs = 0
477 BODY: do
478 bodyIndex = obs_getBodyIndex(obsSpaceData)
479 if (bodyIndex < 0) exit BODY
480 nObs = nObs + 1
481 end do BODY
482
483 call obs_set_current_body_list(obsSpaceData, headerIndex)
484
485 allocate(Hgt_Obs(nObs), Lat_Obs(nObs), Lon_Obs(nObs))
486 allocate(H_M(nObs,nlev_M), H_T(nObs,nlev_T))
487
488 do iObs = 1,nObs
489 bodyIndex = obs_getBodyIndex(obsSpaceData)
490 if (bodyIndex < 0) exit
491
492 height = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
493 ! If the vertical coordinate is an impact parameter (6e6<himp<7e6), subtract radius:
494 if (6.e6 < height .and. height < 7.e6) height = height-rad
495 Hgt_Obs(iObs) = height
496 latr = obs_bodyElem_r(obsSpaceData,OBS_LATD,bodyIndex)
497 lonr = obs_bodyElem_r(obsSpaceData,OBS_LOND,bodyIndex)
498 if (lonr < 0.d0 ) lonr = lonr + 2.d0*MPC_PI_R8
499 if (lonr >= 2.d0*MPC_PI_R8) lonr = lonr - 2.d0*MPC_PI_R8
500 Lat_Obs(iObs) = latr
501 Lon_Obs(iObs) = lonr
502 do levIndex = 1, nlev_M
503 call heightBilinearInterp(latr, lonr, hco, height3D_M_r4(:,:,levIndex), &
504 H_M(iObs,levIndex))
505 end do
506 do levIndex = 1, nlev_T
507 call heightBilinearInterp(latr, lonr, hco, height3D_T_r4(:,:,levIndex), &
508 H_T(iObs,levIndex))
509 end do
510 end do
511
512 do levIndex = 1, nlev_M
513 hmin = 1.e30
514 imin = -1
515 do iObs = 1, nObs
516 dH = abs(H_M(iObs,levIndex)-Hgt_Obs(iObs))
517 if (dH < hmin) then
518 hmin = dH
519 imin = iObs
520 end if
521 end do
522 latSlantLev_M(levIndex) = Lat_Obs(imin)
523 lonSlantLev_M(levIndex) = Lon_Obs(imin)
524 end do
525
526 do levIndex = 1, nlev_T
527 hmin = 1.e30
528 imin = -1
529 do iObs = 1, nObs
530 dH = abs(H_T(iObs,levIndex)-Hgt_Obs(iObs))
531 if (dH < hmin) then
532 hmin = dH
533 imin = iObs
534 end if
535 end do
536 latSlantLev_T(levIndex) = Lat_Obs(imin)
537 lonSlantLev_T(levIndex) = Lon_Obs(imin)
538 end do
539 latSlantLev_S = latSlantLev_T(nlev_T)
540 lonSlantLev_S = lonSlantLev_T(nlev_T)
541
542 deallocate(H_T, H_M)
543 deallocate(Lon_Obs, Lat_Obs, Hgt_Obs)
544
545 end subroutine slp_calcLatLonRO
546
547 subroutine slp_calcLatLonRadar(obsSpaceData, hco, headerIndex, height3D_T_r4, height3D_M_r4, &
548 latSlantLev_T, lonSlantLev_T, latSlantLev_M, lonSlantLev_M, latSlantLev_S, lonSlantLev_S)
549 !
550 !:Purpose: call the computation of lat/lon on the slant path for radar
551 ! observations, iteratively. To replace the vertical columns with
552 ! radar beam columns .
553 !
554 implicit none
555
556 ! Arguments:
557 type(struct_obs), intent(in) :: obsSpaceData
558 type(struct_hco), intent(in) :: hco
559 integer, intent(in) :: headerIndex
560 real(4), intent(in) :: height3D_T_r4(:,:,:)
561 real(4), intent(in) :: height3D_M_r4(:,:,:)
562 real(8), intent(out) :: latSlantLev_T(:)
563 real(8), intent(out) :: lonSlantLev_T(:)
564 real(8), intent(out) :: latSlantLev_M(:)
565 real(8), intent(out) :: lonSlantLev_M(:)
566 real(8), intent(out) :: latSlantLev_S
567 real(8), intent(out) :: lonSlantLev_S
568
569 ! Locals:
570 real(8) :: antennaLat, antennaLon, latSlant, lonSlant, beamElevation, beamAzimuth
571 real(8) :: radarAltitude, beamRangeStart, beamRangeEnd
572 integer :: nlev_M,lev_M, nlev_T,lev_T
573
574 nlev_M = size(height3D_M_r4,3)
575 nlev_T = size(height3D_T_r4,3)
576
577 antennaLat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
578 antennaLon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
579 if (antennaLon < 0.0d0 ) antennaLon = antennaLon + 2.0d0*MPC_PI_R8
580 if (antennaLon >= 2.0d0*MPC_PI_R8) antennaLon = antennaLon - 2.0d0*MPC_PI_R8
581
582 beamElevation = obs_headElem_r(obsSpaceData, OBS_RELE, headerIndex)
583 beamAzimuth = obs_headElem_r(obsSpaceData, OBS_RZAM, headerIndex)
584 radarAltitude = obs_headElem_r(obsSpaceData, OBS_ALT, headerIndex)
585 beamRangeStart = obs_headElem_r(obsSpaceData, OBS_RANS, headerIndex)
586 beamRangeEnd = obs_headElem_r(obsSpaceData, OBS_RANE, headerIndex)
587
588 ! put the last antennaLat/antennaLon at the surface of thermo level
589 latSlantLev_T(nlev_T) = antennaLat
590 lonSlantLev_T(nlev_T) = antennaLon
591 ! Loop through rest of thermo level
592 do lev_T = 1, nlev_T-1
593 call findIntersectLatlonRadar(antennaLat, antennaLon, beamElevation, beamAzimuth, radarAltitude, & !in
594 beamRangeStart, beamRangeEnd, hco, height3D_T_r4(:,:,lev_T), & !in
595 latSlant, lonSlant) !out
596 latSlantLev_T(lev_T) = latSlant
597 lonSlantLev_T(lev_T) = lonSlant
598 end do
599
600 ! put the last antennaLat/beamlon at the surface of momentum levels
601 latSlantLev_M(nlev_M) = antennaLat
602 lonSlantLev_M(nlev_M) = antennaLon
603 ! loop through rest of momentum levels
604 do lev_M = 1, nlev_M-1
605 ! find the intersection
606 call findIntersectLatlonRadar(antennaLat, antennaLon , beamElevation, beamAzimuth, radarAltitude, & !in
607 beamRangeStart, beamRangeEnd, hco, height3D_M_r4(:,:,lev_M), & !in
608 latSlant, lonSlant) !out
609 latSlantLev_M(lev_M) = latSlant
610 lonSlantLev_M(lev_M) = lonSlant
611 end do
612
613 latSlantLev_S = antennaLat
614 lonSlantLev_S = antennaLon
615
616 end subroutine slp_calcLatLonRadar
617
618 subroutine findIntersectLatlonRadar(antennaLat, antennaLon, beamElevation, beamAzimuth, radarAltitude, & !in
619 beamRangeStart, beamRangeEnd, hco, field2d_height, & !in
620 latSlant, lonSlant) !out
621 !
622 !:Purpose: Computation of lat lon of the intersection
623 ! of radar beam with the levels model
624 !
625 !:NOTE: Bisection method
626 !
627 implicit none
628
629 ! Arguments:
630 type(struct_hco), intent(in) :: hco
631 real(8), intent(in) :: antennaLat
632 real(8), intent(in) :: antennaLon
633 real(8), intent(in) :: beamElevation
634 real(8), intent(in) :: beamAzimuth
635 real(8), intent(in) :: radarAltitude
636 real(8), intent(in) :: beamRangeStart
637 real(8), intent(in) :: beamRangeEnd
638 real(4), intent(in) :: field2d_height(hco%ni,hco%nj)
639 real(8), intent(out) :: latSlant
640 real(8), intent(out) :: lonSlant
641
642 ! Locals:
643 real(8) :: upper_bound, lower_bound, tolerance
644 real(8) :: beamHeight, mid, difference_heights, beamDistance
645 integer :: iteration,maximum_iteration
646 real(4) :: heightInterp_r4
647
648 lower_bound = beamRangeStart ! m max useable range for weather radars
649 upper_bound = beamRangeEnd ! m minimum range
650 tolerance = 1 ! tolerance of the bisection method
651 iteration = 0
652 maximum_iteration = 30
653 ! bisection method
654 mid=0.
655 do while ((abs((upper_bound - lower_bound)/2.)>tolerance).or.(iteration>maximum_iteration))
656 mid = (upper_bound + lower_bound)/2
657 call rdv_getlatlonHRfromRange(antennaLat, antennaLon, beamElevation, beamAzimuth, & !in
658 radarAltitude, mid, & !in
659 latSlant, lonSlant, beamHeight, beamDistance) !out
660 call heightBilinearInterp(latSlant, lonSlant, hco, field2d_height, heightInterp_r4)
661 difference_heights = beamHeight - heightInterp_r4
662 if (difference_heights>0.) then
663 upper_bound = mid
664 else
665 lower_bound = mid
666 end if
667 iteration = iteration+1
668 end do
669
670 if (lonSlant < 0.0d0 ) lonSlant = lonSlant + 2.0d0*MPC_PI_R8
671 if (lonSlant >= 2.0d0*MPC_PI_R8) lonSlant = lonSlant - 2.0d0*MPC_PI_R8
672
673 end subroutine findIntersectLatlonRadar
674
675end module slantProfileLatLon_mod
