1module windRotation_mod
2 ! MODULE windRotation_mod (prefix='uvr' category='4. Data Object transformations')
3 !
4 !:Purpose: To transform winds FROM the rotated spherical coordinate system
5 ! TO the non-rotated spherical coordinate system. Also includes
6 ! tangent-linear and adjoint versions of the transformation.
7 !
8 use mathPhysConstants_mod
9 use horizontalCoord_mod
10 use utilities_mod
11 use midasMpi_mod
12 implicit none
13 save
14 private
15
16 ! Public derived type definition
17 public :: struct_uvr
18
19 ! Public Subroutines
20 public :: uvr_setup, uvr_rotateWind_nl, uvr_rotateWind_tl, uvr_rotateWind_ad, uvr_rotateLatLon
21
22 integer, parameter :: msize = 3
23 integer, parameter :: maxNumSubGrid = 2
24
25 type :: struct_uvr
26 real(8) :: grd_rot_8 (msize,msize,maxNumSubGrid) ! Real Earth to Rotated Earth
27 real(8) :: grd_rotinv_8(msize,msize,maxNumSubGrid) ! Rotated Earth to Real Earth
28 logical :: initialized = .false.
29 end type struct_uvr
30
31 ! Minimum value for latitude used in TL/AD wind rotation routines
32 real(8), parameter :: lat_minValue = 1D-10
33
34 contains
35
36 subroutine uvr_setup( uvr, hco_in )
37 !
38 !:Purpose: Setup the information for wind rotation
39 !
40 implicit none
41
42 ! Arguments:
43 type(struct_uvr), pointer, intent(inout) :: uvr ! Wind rotation object
44 type(struct_hco), intent(in) :: hco_in ! Horizontal grid object
45
46 ! Locals:
47 logical, save :: firstCall = .true.
48
49 if ( associated(uvr) ) then
50 if ( uvr%initialized ) then
51 if ( mmpi_myid == 0 ) write(*,*) 'uvr_setup: already initialized, returning'
52 return
53 else
54 call utl_abort('uvr_setup: the supplied non-null uvr pointer is not initialized!')
55 end if
56 end if
57
58 !
59 !- Compute the rotation matrices (grd_rot_8 and grd_rotinv_8)
60 !
61 if ( mmpi_myid == 0 ) then
62 write(*,*)
63 write(*,*) 'uvr_setup: Starting... for grid type = ', hco_in%grtyp
64 end if
65
66 allocate(uvr)
67
68 if ( mmpi_myid == 0 .and. firstCall ) then
69 write(*,*) hco_in % xlon1
70 write(*,*) hco_in % xlat1
71 write(*,*) hco_in % xlon2
72 write(*,*) hco_in % xlat2
73 end if
74
75 call sugrdpar(uvr, 1, hco_in % xlon1, hco_in % xlat1, hco_in % xlon2, hco_in % xlat2 )
76
77 if ( hco_in%grtyp == 'U' ) then
78
79 if ( mmpi_myid == 0 .and. firstCall ) then
80 write(*,*) 'uvr_setup: doing setup for YAN grid'
81 write(*,*) hco_in % xlon1_yan
82 write(*,*) hco_in % xlat1_yan
83 write(*,*) hco_in % xlon2_yan
84 write(*,*) hco_in % xlat2_yan
85 end if
86
87 call sugrdpar(uvr, 2, hco_in % xlon1_yan, hco_in % xlat1_yan, hco_in % xlon2_yan, hco_in % xlat2_yan )
88
89 end if
90
91 uvr%initialized = .true.
92
93 if ( mmpi_myid == 0 ) then
94 write(*,*)
95 write(*,*) 'uvr_setup: Done!'
96 end if
97
98 firstCall = .false.
99
100 end subroutine uvr_setup
101
102 subroutine sugrdpar( uvr, subGridIndex, grd_xlon1, grd_xlat1, grd_xlon2, grd_xlat2 )
103 !
104 !:Purpose: Compute the rotation matrix (r_8) that allows transformation
105 ! from the non-rotated to the rotated spherical coordinate system.
106 !
107 implicit none
108
109 ! Arguments:
110 type(struct_uvr), pointer, intent(inout) :: uvr ! Wind rotation object
111 integer, intent(in) :: subGridIndex ! Horizontal subGrid index = 2
112 real(8), intent(in) :: grd_xlon1 ! Horizontal grid xlon1_yan
113 real(8), intent(in) :: grd_xlat1 ! Horizontal grid xlat1_yan
114 real(8), intent(in) :: grd_xlon2 ! Horizontal grid xlon2_yan
115 real(8), intent(in) :: grd_xlat2 ! Horizontal grid xlat2_yan
116
117 ! Locals:
118 integer :: j1, j2
119 real(8) :: zxlon1_8,zxlat1_8,zxlon2_8,zxlat2_8
120 real(8) :: a_8, b_8, c_8, d_8, xyz1_8(msize), xyz2_8(msize)
121 real(8) :: zunit(msize,msize)
122 logical, save :: firstCall = .true.
123
124 zxlon1_8 = grd_xlon1
125 zxlat1_8 = grd_xlat1
126 zxlon2_8 = grd_xlon2
127 zxlat2_8 = grd_xlat2
128
129 if ( zxlon1_8 == zxlon2_8 .and. zxlat1_8 == zxlat2_8 ) then
130 !
131 !- 1. Non rotated grid (for test case only)
132 !
133 if ( mmpi_myid == 0 ) then
134 write(*,*)
135 write(*,*) 'uvr_sugrdpar: Warning: This grid is not rotated !!!'
136 end if
137 uvr%grd_rot_8(1,1,subGridIndex) = 1.0d0
138 uvr%grd_rot_8(1,2,subGridIndex) = 0.0d0
139 uvr%grd_rot_8(1,3,subGridIndex) = 0.0d0
140 uvr%grd_rot_8(2,1,subGridIndex) = 0.0d0
141 uvr%grd_rot_8(2,2,subGridIndex) = 1.0d0
142 uvr%grd_rot_8(2,3,subGridIndex) = 0.0d0
143 uvr%grd_rot_8(3,1,subGridIndex) = 0.0d0
144 uvr%grd_rot_8(3,2,subGridIndex) = 0.0d0
145 uvr%grd_rot_8(3,3,subGridIndex) = 1.0d0
146
147 else
148 !
149 !- 2. Rotated grid
150 !
151 call vllacar( xyz1_8, zxlon1_8, zxlat1_8 )
152 call vllacar( xyz2_8, zxlon2_8, zxlat2_8 )
153
154 !- 2.1 Compute a = cos(alpha) & b = sin(alpha)
155 a_8 = (xyz1_8(1)*xyz2_8(1)) + (xyz1_8(2)*xyz2_8(2)) &
156 + (xyz1_8(3)*xyz2_8(3))
157
158 b_8 = sqrt (((xyz1_8(2)*xyz2_8(3)) - (xyz2_8(2)*xyz1_8(3)))**2 &
159 + ((xyz2_8(1)*xyz1_8(3)) - (xyz1_8(1)*xyz2_8(3)))**2 &
160 + ((xyz1_8(1)*xyz2_8(2)) - (xyz2_8(1)*xyz1_8(2)))**2)
161
162 !- 2.2 Compute c = norm(-r1) & d = norm(r4)
163 c_8 = sqrt ( xyz1_8(1)**2 + xyz1_8(2)**2 + xyz1_8(3)**2 )
164
165 d_8 = sqrt ( ( ( (a_8*xyz1_8(1)) - xyz2_8(1) ) / b_8 )**2 + &
166 ( ( ( a_8*xyz1_8(2)) - xyz2_8(2) ) / b_8 )**2 + &
167 ( ( ( a_8*xyz1_8(3)) - xyz2_8(3) ) / b_8 )**2 )
168
169 !- 2.3 Compute the forward rotation matrix
170 uvr%grd_rot_8(1,1,subGridIndex) = -xyz1_8(1) / c_8
171 uvr%grd_rot_8(1,2,subGridIndex) = -xyz1_8(2) / c_8
172 uvr%grd_rot_8(1,3,subGridIndex) = -xyz1_8(3) / c_8
173 uvr%grd_rot_8(2,1,subGridIndex) = ( ((a_8*xyz1_8(1)) - xyz2_8(1)) / b_8 ) / d_8
174 uvr%grd_rot_8(2,2,subGridIndex) = ( ((a_8*xyz1_8(2)) - xyz2_8(2)) / b_8 ) / d_8
175 uvr%grd_rot_8(2,3,subGridIndex) = ( ((a_8*xyz1_8(3)) - xyz2_8(3)) / b_8 ) / d_8
176 uvr%grd_rot_8(3,1,subGridIndex) = ( (xyz1_8(2)*xyz2_8(3)) - (xyz2_8(2)*xyz1_8(3))) / b_8
177 uvr%grd_rot_8(3,2,subGridIndex) = ( (xyz2_8(1)*xyz1_8(3)) - (xyz1_8(1)*xyz2_8(3))) / b_8
178 uvr%grd_rot_8(3,3,subGridIndex) = ( (xyz1_8(1)*xyz2_8(2)) - (xyz2_8(1)*xyz1_8(2))) / b_8
179
180 end if
181
182 do j1 = 1, msize
183 do j2 = 1, msize
184 if ( mmpi_myid == 0 .and. firstCall ) then
185 write(*,*) 'sugrdpar: grd_rot_8(j1,j2) =',j1,j2,uvr%grd_rot_8(j1,j2,subGridIndex)
186 end if
187 end do
188 end do
189
190 !
191 !- 3. Compute the inverse matrix (transpose is the inverse)
192 !
193 do j1 = 1, msize
194 do j2 = 1, msize
195 uvr%grd_rotinv_8(j1,j2,subGridIndex) = uvr%grd_rot_8(j2,j1,subGridIndex)
196 end do
197 end do
198
199 zunit(:,:) = 0.d0
200 call mxma8x(zunit,uvr%grd_rotinv_8(:,:,subGridIndex),uvr%grd_rot_8(:,:,subGridIndex),msize,msize,msize)
201 do j1 = 1, msize
202 do j2 = 1, msize
203 if ( mmpi_myid == 0 .and. firstCall ) then
204 write(*,*) 'sugrdpar: unit = ', j1, j2, zunit(j1,j2)
205 end if
206 end do
207 end do
208
209 firstCall = .false.
210
211 end subroutine sugrdpar
212
213 subroutine vllacar( F_xyz_8, F_lon, F_lat )
214 !
215 !:Purpose: Compute parameters of rotated grid
216 !
217 implicit none
218
219 ! Arguments:
220 real(8), intent(out) :: F_xyz_8(msize) ! output
221 real(8), intent(in) :: F_lon ! Input in degrees
222 real(8), intent(in) :: F_lat ! Input in degrees
223
224 F_xyz_8(1) = cos(MPC_RADIANS_PER_DEGREE_R8*F_lat) * cos(MPC_RADIANS_PER_DEGREE_R8*F_lon)
225 F_xyz_8(2) = cos(MPC_RADIANS_PER_DEGREE_R8*F_lat) * sin(MPC_RADIANS_PER_DEGREE_R8*F_lon)
226 F_xyz_8(3) = sin(MPC_RADIANS_PER_DEGREE_R8*F_lat)
227
228 end subroutine vllacar
229
230 subroutine mxma8x( pmat3, pmat1, pmat2, kdimi1, kdimj1, kdimj2 )
231 !
232 !:Purpose: Compute a product of two matrices.
233 !
234 implicit none
235
236 ! Arguments:
237 real(8), intent(out) :: pmat3(kdimi1,kdimj2) ! output
238 real(8), intent(in) :: pmat1(kdimi1,kdimj1) ! input matrix one
239 real(8), intent(in) :: pmat2(kdimj1,kdimj2) ! input matrix two
240 integer, intent(in) :: kdimi1 ! first dimension of the first matrix
241 integer, intent(in) :: kdimj1 ! second dimension of the first matrix
242 integer, intent(in) :: kdimj2 ! second dimension of the second matrix
243
244 ! Locals:
245 integer :: ji1,jj2,jj
246
247 pmat3(:,:) = 0.d0
248
249 do jj2 = 1, kdimj2
250 do jj = 1, kdimj1
251 do ji1 = 1, kdimi1
252 pmat3(ji1,jj2) = pmat3(ji1,jj2) + pmat1(ji1,jj) * pmat2(jj,jj2)
253 end do
254 end do
255 end do
256
257 end subroutine mxma8x
258
259 subroutine uvr_rotateWind_nl( uvr, subGridIndex, uwind, vwind, Lat, Lon, LatRot, LonRot, mode )
260 !
261 !:Purpose: Go from tangential wind components from one sphere to another
262 ! (same origin!). Original ezsint version used for computing innovation.
263 !
264 implicit none
265
266 ! Arguments:
267 type(struct_uvr), pointer, intent(in) :: uvr ! Wind rotation object
268 integer, intent(in) :: subGridIndex ! Current subgrid index
269 real(8), intent(inout) :: uwind ! interpUU
270 real(8), intent(inout) :: vwind ! interpVV
271 real(8), intent(in) :: Lat ! Latitude in radians
272 real(8), intent(in) :: Lon ! Longitude in radians
273 real(8), intent(in) :: LatRot ! Rotated latitude in radians
274 real(8), intent(in) :: LonRot ! Rotated longitude in radians
275 character(*), intent(in) :: mode ! ToMetWind or ToRotWind
276
277 ! Locals:
278 integer :: index1, index2
279 real(8) :: coslatr, sinlatr, coslonr, sinlonr, coslat, sinlat, coslon, sinlon, ezCoeff_C, ezCoeff_D
280 real(8) :: xyz(msize), uvcart(msize)
281
282 if ( .not. uvr%initialized ) then
283 write(*,*)
284 call utl_abort('uvr_rotateWind_nl: WindRotation module is not initialize')
285 endif
286
287 coslatr = cos(LatRot)
288 sinlatr = sin(LatRot)
289 coslonr = cos(LonRot)
290 sinlonr = sin(LonRot)
291 coslat = cos(Lat)
292 sinlat = sin(Lat)
293 coslon = cos(Lon)
294 sinlon = sin(Lon)
295
296 if ( trim(mode) == 'ToMetWind' ) then
297
298 xyz(1) = -uwind*sinlonr - vwind*coslonr*sinlatr
299 xyz(2) = uwind*coslonr - vwind*sinlonr*sinlatr
300 xyz(3) = vwind*coslatr
301
302 uvcart(:) = 0.0d0
303 do index2 = 1, msize
304 do index1 = 1, msize
305 uvcart(index1) = uvcart(index1) + &
306 uvr%grd_rotinv_8(index1,index2,subGridIndex)*xyz(index2)
307 end do
308 end do
309
310 uwind = uvcart(2)*coslon - uvcart(1)*sinlon
311 ezCoeff_C = uvcart(1)*coslon + uvcart(2)*sinlon
312 ezCoeff_D = sqrt( ezCoeff_C**2 + uvcart(3)**2 )
313 vwind = sign( ezCoeff_D, uvcart(3)*coslat - ezCoeff_C*sinlat )
314
315 else if ( trim(mode) == 'ToRotWind' ) then
316 write(*,*)
317 call utl_abort('uvr_rotateWind_nl: mode ToRotWind is not available yet')
318 else
319 write(*,*)
320 write(*,*) 'uvr_rotateWind_nl: Unknown transform name: ', trim(mode)
321 write(*,*) ' mode = ToMetWind or ToRotWind'
322 call utl_abort('uvr_rotateWind_nl')
323 end if
324
325 end subroutine uvr_rotateWind_nl
326
327 subroutine uvr_rotateWind_tl( uvr, subGridIndex, uwind, vwind, Lat_in, Lon_in, LatRot_in, LonRot_in, mode )
328 !
329 !:Purpose: Go from tangential wind components from one sphere to another
330 ! (same origin!). Fast version used by Variational analysis.
331 !
332 implicit none
333
334 ! Arguments:
335 type(struct_uvr), pointer, intent(in) :: uvr ! Wind rotation object
336 integer, intent(in) :: subGridIndex ! Current subgrid index
337 real(8), intent(inout) :: uwind ! interpUU
338 real(8), intent(inout) :: vwind ! interpVV
339 real(8), intent(in) :: Lat_in ! Latitude in radians
340 real(8), intent(in) :: Lon_in ! Longitude in radians
341 real(8), intent(in) :: LatRot_in ! Rotated latitude in radians
342 real(8), intent(in) :: LonRot_in ! Rotated longitude in radians
343 character(*), intent(in) :: mode ! ToMetWind or ToRotWind
344
345 ! Locals:
346 integer :: index1, index2
347 real(8) :: coslatr, sinlatr, coslonr, sinlonr, coslon, sinlon, rsinlat
348 real(8) :: lat, lon, latRot, lonRot
349 real(8) :: xyz(msize), uvcart(msize)
350
351 if ( .not. uvr%initialized ) then
352 write(*,*)
353 call utl_abort('uvr_rotateWind_tl: WindRotation module is not initialize')
354 endif
355
356 ! Special case for the equator
357 if (abs(lat_in) < lat_minValue) then
358 lat = lat_minValue
359 lon = lon_in
360 call uvr_RotateLatLon( uvr, & ! INOUT
361 subGridIndex, & ! IN
362 latRot, lonRot, & ! OUT (radians)
363 lat, lon, & ! IN (radians)
364 'ToLatLonRot') ! IN
365 else
366 lat = lat_in
367 lon = lon_in
368 latRot = latRot_in
369 lonRot = lonRot_in
370 end if
371
372 coslatr = cos(LatRot)
373 sinlatr = sin(LatRot)
374 coslonr = cos(LonRot)
375 sinlonr = sin(LonRot)
376 coslon = cos(Lon)
377 sinlon = sin(Lon)
378 if ( Lat /= 0.0d0 ) then
379 rsinlat = 1.0d0/sin(Lat)
380 else
381 call utl_abort('uvr_rotateWind_tl: cannot be used for points on the equator')
382 end if
383
384 if ( trim(mode) == 'ToMetWind' ) then
385
386 xyz(1) = -uwind*sinlonr - vwind*coslonr*sinlatr
387 xyz(2) = uwind*coslonr - vwind*sinlonr*sinlatr
388 xyz(3) = vwind*coslatr
389
390 uvcart(:) = 0.0d0
391 do index2 = 1, msize
392 do index1 = 1, msize
393 uvcart(index1) = uvcart(index1) + &
394 uvr%grd_rotinv_8(index1,index2,subGridIndex)*xyz(index2)
395 end do
396 end do
397
398 uwind = -uvcart(1)*sinlon + uvcart(2)*coslon
399 vwind = -uvcart(1)*rsinlat*coslon - uvcart(2)*rsinlat*sinlon
400
401 else if ( trim(mode) == 'ToRotWind' ) then
402 write(*,*)
403 call utl_abort('uvr_rotateWind_tl: mode ToRotWind is not available yet')
404 else
405 write(*,*)
406 write(*,*) 'uvr_rotateWind_tl: Unknown transform name: ', trim(mode)
407 write(*,*) ' mode = ToMetWind or ToRotWind'
408 call utl_abort('uvr_rotateWind_tl')
409 end if
410
411 end subroutine uvr_rotateWind_tl
412
413 subroutine uvr_rotateWind_ad( uvr, subGridIndex, uwind, vwind, Lat_in, Lon_in, LatRot_in, LonRot_in, mode )
414 !
415 !:Purpose: Adjoint of : Go from tangential wind components from one sphere to another
416 ! (same origin!). Fast version used by Variational analysis.
417 !
418 implicit none
419
420 ! Arguments:
421 type(struct_uvr), pointer, intent(in) :: uvr ! Wind rotation object
422 integer, intent(in) :: subGridIndex ! Current subgrid index
423 real(8), intent(inout) :: uwind ! interpUU
424 real(8), intent(inout) :: vwind ! interpVV
425 real(8), intent(in) :: Lat_in ! Latitude in radians
426 real(8), intent(in) :: Lon_in ! Longitude in radians
427 real(8), intent(in) :: LatRot_in ! Rotated latitude in radians
428 real(8), intent(in) :: LonRot_in ! Rotated longitude in radians
429 character(*), intent(in) :: mode ! ToMetWind or ToRotWind
430
431 ! Locals:
432 integer :: index1, index2
433 real(8) :: coslatr, sinlatr, coslonr, sinlonr, coslon, sinlon, rsinlat
434 real(8) :: lat, lon, latRot, lonRot
435 real(8) :: xyz(msize), uvcart(msize)
436
437 if ( .not. uvr%initialized ) then
438 write(*,*)
439 call utl_abort('uvr_rotateWind_ad: WindRotation module is not initialize')
440 endif
441
442 ! Special case for the equator
443 if (abs(lat_in) < lat_minValue) then
444 lat = lat_minValue
445 lon = lon_in
446 call uvr_RotateLatLon( uvr, & ! INOUT
447 subGridIndex, & ! IN
448 latRot, lonRot, & ! OUT (radians)
449 lat, lon, & ! IN (radians)
450 'ToLatLonRot') ! IN
451 else
452 lat = lat_in
453 lon = lon_in
454 latRot = latRot_in
455 lonRot = lonRot_in
456 end if
457
458 coslatr = cos(LatRot)
459 sinlatr = sin(LatRot)
460 coslonr = cos(LonRot)
461 sinlonr = sin(LonRot)
462 coslon = cos(Lon)
463 sinlon = sin(Lon)
464 if ( Lat /= 0.0d0 ) then
465 rsinlat = 1.0d0/sin(Lat)
466 else
467 call utl_abort('uvr_rotateWind_ad: cannot be used for points on the equator')
468 end if
469
470 if ( trim(mode) == 'ToMetWind' ) then
471
472 uvcart(1) = -uwind*sinlon - vwind*rsinlat*coslon
473 uvcart(2) = uwind*coslon - vwind*rsinlat*sinlon
474 uvcart(3) = 0.0d0
475
476 xyz(:) = 0.0d0
477 do index2 = 1, msize
478 do index1 = 1, msize
479 xyz(index1) = xyz(index1) + uvr%grd_rot_8(index1,index2,subGridIndex)*uvcart(index2)
480 end do
481 end do
482
483 uwind = -xyz(1)*sinlonr + xyz(2)*coslonr
484 vwind = -xyz(1)*coslonr*sinlatr - xyz(2)*sinlonr*sinlatr + xyz(3)*coslatr
485
486 else if ( trim(mode) == 'ToRotWind' ) then
487 write(*,*)
488 call utl_abort('uvr_rotateWind_ad: mode ToRotWind is not available yet')
489 else
490 write(*,*)
491 write(*,*) 'uvr_rotateWind_ad: Unknown transform name: ', trim(mode)
492 write(*,*) ' mode = ToMetWind or ToRotWind'
493 call utl_abort('uvr_rotateWind_ad')
494 end if
495
496 end subroutine uvr_rotateWind_ad
497
498 subroutine uvr_rotateLatLon( uvr, subGridIndex, LatOut, LonOut, LatIn, LonIn, mode )
499 !
500 !:Purpose: Go from (lat,lon) of one Cartesian frame to (lat,lon) of another
501 ! Cartesian frame given the rotation matrix.
502 !
503 implicit none
504
505 ! Arguments:
506 type(struct_uvr), pointer, intent(in) :: uvr ! Wind rotation object
507 integer, intent(in) :: subGridIndex ! Current subgrid index
508 real(8), intent(in) :: LatIn ! Input latitude in radians
509 real(8), intent(in) :: LonIn ! Input longitude in radians
510 real(8), intent(out) :: LatOut ! Output latitude in radians
511 real(8), intent(out) :: LonOut ! Output longitude in radians
512 character(*), intent(in) :: mode ! ToLatLonRot or ToLatLon
513
514 ! Locals:
515 real(8) :: CartIn(msize),CartOut(msize)
516 real(8) :: rLon,rLat
517
518 if ( .not. uvr%initialized ) then
519 write(*,*)
520 call utl_abort('uvr_rotateLatLon: WindRotation module is not initialize')
521 endif
522
523 rLon = LonIn * MPC_DEGREES_PER_RADIAN_R8 ! To degress
524 rLat = LatIn * MPC_DEGREES_PER_RADIAN_R8 ! To degrees
525
526 call vllacar( CartIn, & ! OUT
527 rLon, rLat ) ! IN
528
529 if ( trim(mode) == 'ToLatLonRot' ) then
530 call mxv( CartOut, & ! OUT
531 uvr%grd_rot_8(:,:,subGridIndex), CartIn, & ! IN
532 msize, msize) ! IN
533 else if ( trim(mode) == 'ToLatLon' ) then
534 call mxv( CartOut, & ! OUT
535 uvr%grd_rotinv_8(:,:,subGridIndex), CartIn, & ! IN
536 msize, msize) ! IN
537 else
538 write(*,*)
539 write(*,*) 'uvr_rotateLatLon: Unknown transform name: ', trim(mode)
540 write(*,*) ' mode = ToLatLonRot or ToLatLon'
541 call utl_abort('uvr_rotateLatLon')
542 end if
543
544 call carall( LonOut, LatOut, & ! OUT
545 CartOut ) ! IN
546
547 LonOut = LonOut * MPC_RADIANS_PER_DEGREE_R8 ! To radians
548 LatOut = LatOut * MPC_RADIANS_PER_DEGREE_R8 ! To radians
549
550 end subroutine uvr_rotateLatLon
551
552
553 subroutine carall( plon, plat, pcart )
554 !
555 !:Purpose: Returns (lat,lon) (degrees) of an input Cartesian position vector
556 ! on the unit sphere.
557 !
558 implicit none
559
560 ! Arguments:
561 real(8), intent(out) :: plat ! output latitude
562 real(8), intent(out) :: plon ! output longitude
563 real(8), intent(in) :: pcart(msize) ! input Cartesian vector
564
565 plat = asin(pcart(3))
566 plat = plat * MPC_DEGREES_PER_RADIAN_R8
567
568 if ( pcart(1) == 0.d0 ) then
569 if ( pcart(2) == 0.d0 ) then ! point is located at the pole
570 plon = 0.0d0 ! can be any longitude... set it to zero simply.
571 else if ( pcart(2) > 0.0d0 ) then
572 plon = 90.0d0
573 else if ( pcart(2) < 0.0d0 ) then
574 plon = 270.0d0
575 end if
576 else
577 plon = atan2(pcart(2),pcart(1))
578 if (plon < 0.0d0) plon = plon + 2.d0 * MPC_PI_R8
579 plon = plon * MPC_DEGREES_PER_RADIAN_R8
580 end if
581
582 end subroutine carall
583
584 subroutine mxv( pvec2, pmat, pvec1, kdimi, kdimj )
585 !
586 !:Purpose: Compute a product : matrix times vector.
587 !
588 implicit none
589
590 ! Arguments:
591 real(8), intent(out) :: pvec2(kdimi) ! output vector
592 real(8), intent(in) :: pmat(kdimi,kdimj) ! input matrix
593 real(8), intent(in) :: pvec1(kdimj) ! input vector
594 integer, intent(in) :: kdimi ! first dimension
595 integer, intent(in) :: kdimj ! second dimension
596
597 ! Locals:
598 integer :: ji,jj
599
600 pvec2(:) = 0.0d0
601 do jj = 1,kdimj
602 do ji = 1,kdimi
603 pvec2(ji) = pvec2(ji) + pmat(ji,jj) * pvec1(jj)
604 end do
605 end do
606
607 end subroutine mxv
608
609end module windRotation_mod
