PageRenderTime 99ms CodeModel.GetById 17ms RepoModel.GetById 0ms app.codeStats 0ms

/WPS/metgrid/src/rotate_winds_module.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 595 lines | 394 code | 99 blank | 102 comment | 97 complexity | 35175581bc2258903df5d12db8ed1c82 MD5 | raw file
Possible License(s): AGPL-1.0
  1. module rotate_winds_module
  2. use bitarray_module
  3. use constants_module
  4. use llxy_module
  5. use misc_definitions_module
  6. use module_debug
  7. integer :: orig_selected_projection
  8. contains
  9. !
  10. ! ARW Wind Rotation Code
  11. !
  12. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  13. ! Name: map_to_met !
  14. ! !
  15. ! Purpose: Rotate grid-relative winds to Earth-relative winds !
  16. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  17. subroutine map_to_met(u, u_mask, v, v_mask, &
  18. us1, us2, ue1, ue2, &
  19. vs1, vs2, ve1, ve2, &
  20. xlon_u, xlon_v, xlat_u, xlat_v)
  21. implicit none
  22. ! Arguments
  23. integer, intent(in) :: us1, us2, ue1, ue2, vs1, vs2, ve1, ve2
  24. real, pointer, dimension(:,:) :: u, v, xlon_u, xlon_v, xlat_u, xlat_v
  25. type (bitarray), intent(in) :: u_mask, v_mask
  26. orig_selected_projection = iget_selected_domain()
  27. call select_domain(SOURCE_PROJ)
  28. call metmap_xform(u, u_mask, v, v_mask, &
  29. us1, us2, ue1, ue2, &
  30. vs1, vs2, ve1, ve2, &
  31. xlon_u, xlon_v, xlat_u, xlat_v, 1)
  32. call select_domain(orig_selected_projection)
  33. end subroutine map_to_met
  34. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  35. ! Name: met_to_map !
  36. ! !
  37. ! Purpose: Rotate Earth-relative winds to grid-relative winds !
  38. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  39. subroutine met_to_map(u, u_mask, v, v_mask, &
  40. us1, us2, ue1, ue2, &
  41. vs1, vs2, ve1, ve2, &
  42. xlon_u, xlon_v, xlat_u, xlat_v)
  43. implicit none
  44. ! Arguments
  45. integer, intent(in) :: us1, us2, ue1, ue2, vs1, vs2, ve1, ve2
  46. real, pointer, dimension(:,:) :: u, v, xlon_u, xlon_v, xlat_u, xlat_v
  47. type (bitarray), intent(in) :: u_mask, v_mask
  48. orig_selected_projection = iget_selected_domain()
  49. call select_domain(1)
  50. call metmap_xform(u, u_mask, v, v_mask, &
  51. us1, us2, ue1, ue2, &
  52. vs1, vs2, ve1, ve2, &
  53. xlon_u, xlon_v, xlat_u, xlat_v, -1)
  54. call select_domain(orig_selected_projection)
  55. end subroutine met_to_map
  56. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  57. ! Name: metmap_xform !
  58. ! !
  59. ! Purpose: Do the actual work of rotating winds for C grid. !
  60. ! If idir= 1, rotate grid-relative winds to Earth-relative winds !
  61. ! If idir=-1, rotate Earth-relative winds to grid-relative winds !
  62. ! !
  63. ! ASSUMPTIONS: 1) MEMORY ORDER IS XY. !
  64. ! 2) U ARRAY HAS ONE MORE COLUMN THAN THE V ARRAY, AND V ARRAY !
  65. ! HAS ONE MORE ROW THAN U ARRAY. !
  66. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  67. subroutine metmap_xform(u, u_mask, v, v_mask, &
  68. us1, us2, ue1, ue2, vs1, vs2, ve1, ve2, &
  69. xlon_u, xlon_v, xlat_u, xlat_v, idir)
  70. implicit none
  71. ! Arguments
  72. integer, intent(in) :: us1, us2, ue1, ue2, vs1, vs2, ve1, ve2, idir
  73. real, pointer, dimension(:,:) :: u, v, xlon_u, xlon_v, xlat_u, xlat_v
  74. type (bitarray), intent(in) :: u_mask, v_mask
  75. ! Local variables
  76. integer :: i, j
  77. real :: u_weight, v_weight
  78. real :: u_map, v_map, alpha, diff
  79. real, pointer, dimension(:,:) :: u_new, v_new, u_mult, v_mult
  80. logical :: do_last_col_u, do_last_row_u, do_last_col_v, do_last_row_v
  81. ! If the proj_info structure has not been initialized, we don't have
  82. ! information about the projection and standard longitude.
  83. if (proj_stack(current_nest_number)%init) then
  84. ! Only rotate winds for Lambert conformal, polar stereographic, or Cassini
  85. if ((proj_stack(current_nest_number)%code == PROJ_LC) .or. &
  86. (proj_stack(current_nest_number)%code == PROJ_PS) .or. &
  87. (proj_stack(current_nest_number)%code == PROJ_CASSINI)) then
  88. call mprintf((idir == 1),LOGFILE,'Rotating map winds to earth winds.')
  89. call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
  90. allocate(u_mult(us1:ue1,us2:ue2))
  91. allocate(v_mult(vs1:ve1,vs2:ve2))
  92. do j=us2,ue2
  93. do i=us1,ue1
  94. if (bitarray_test(u_mask, i-us1+1, j-us2+1)) then
  95. u_mult(i,j) = 1.
  96. else
  97. u_mult(i,j) = 0.
  98. end if
  99. end do
  100. end do
  101. do j=vs2,ve2
  102. do i=vs1,ve1
  103. if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
  104. v_mult(i,j) = 1.
  105. else
  106. v_mult(i,j) = 0.
  107. end if
  108. end do
  109. end do
  110. if (ue1-us1 == ve1-vs1) then
  111. do_last_col_u = .false.
  112. do_last_col_v = .true.
  113. else
  114. do_last_col_u = .true.
  115. do_last_col_v = .false.
  116. end if
  117. if (ue2-us2 == ve2-vs2) then
  118. do_last_row_u = .true.
  119. do_last_row_v = .false.
  120. else
  121. do_last_row_u = .false.
  122. do_last_row_v = .true.
  123. end if
  124. ! Create arrays to hold rotated winds
  125. allocate(u_new(us1:ue1, us2:ue2))
  126. allocate(v_new(vs1:ve1, vs2:ve2))
  127. ! Rotate U field
  128. do j=us2,ue2
  129. do i=us1,ue1
  130. diff = idir * (xlon_u(i,j) - proj_stack(current_nest_number)%stdlon)
  131. if (diff > 180.) then
  132. diff = diff - 360.
  133. else if (diff < -180.) then
  134. diff = diff + 360.
  135. end if
  136. ! Calculate the rotation angle, alpha, in radians
  137. if (proj_stack(current_nest_number)%code == PROJ_LC) then
  138. alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi
  139. else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
  140. if (j == ue2) then
  141. diff = xlon_u(i,j)-xlon_u(i,j-1)
  142. if (diff > 180.) then
  143. diff = diff - 360.
  144. else if (diff < -180.) then
  145. diff = diff + 360.
  146. end if
  147. alpha = atan2( (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg), &
  148. (xlat_u(i,j)-xlat_u(i,j-1))*rad_per_deg &
  149. )
  150. else if (j == us2) then
  151. diff = xlon_u(i,j+1)-xlon_u(i,j)
  152. if (diff > 180.) then
  153. diff = diff - 360.
  154. else if (diff < -180.) then
  155. diff = diff + 360.
  156. end if
  157. alpha = atan2( (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg), &
  158. (xlat_u(i,j+1)-xlat_u(i,j))*rad_per_deg &
  159. )
  160. else
  161. diff = xlon_u(i,j+1)-xlon_u(i,j-1)
  162. if (diff > 180.) then
  163. diff = diff - 360.
  164. else if (diff < -180.) then
  165. diff = diff + 360.
  166. end if
  167. alpha = atan2( (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg), &
  168. (xlat_u(i,j+1)-xlat_u(i,j-1))*rad_per_deg &
  169. )
  170. end if
  171. else
  172. alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
  173. end if
  174. v_weight = 0.
  175. ! On C grid, take U_ij, and get V value at the same lat/lon
  176. ! by averaging the four surrounding V points
  177. if (bitarray_test(u_mask, i-us1+1, j-us2+1)) then
  178. u_map = u(i,j)
  179. if (i == us1) then
  180. if (j == ue2 .and. do_last_row_u) then
  181. v_weight = v_mult(i,j)
  182. v_map = v(i,j)*v_mult(i,j)
  183. else
  184. v_weight = v_mult(i,j) + v_mult(i,j+1)
  185. v_map = v(i,j)*v_mult(i,j) + v(i,j+1)*v_mult(i,j+1)
  186. end if
  187. else if (i == ue1 .and. do_last_col_u) then
  188. if (j == ue2 .and. do_last_row_u) then
  189. v_weight = v_mult(i-1,j)
  190. v_map = v(i-1,j)
  191. else
  192. v_weight = v_mult(i-1,j) + v_mult(i-1,j+1)
  193. v_map = v(i-1,j)*v_mult(i-1,j) + v(i-1,j+1)*v_mult(i-1,j+1)
  194. end if
  195. else if (j == ue2 .and. do_last_row_u) then
  196. v_weight = v_mult(i-1,j-1) + v_mult(i,j-1)
  197. v_map = v(i-1,j-1)*v_mult(i-1,j-1) + v(i,j-1)*v_mult(i,j-1)
  198. else
  199. v_weight = v_mult(i-1,j) + v_mult(i-1,j+1) + v_mult(i,j) + v_mult(i,j+1)
  200. v_map = v(i-1,j)*v_mult(i-1,j) + v(i-1,j+1)*v_mult(i-1,j+1) + v(i,j)*v_mult(i,j) + v(i,j+1)*v_mult(i,j+1)
  201. end if
  202. if (v_weight > 0.) then
  203. u_new(i,j) = cos(alpha)*u_map + sin(alpha)*v_map/v_weight
  204. else
  205. u_new(i,j) = u(i,j)
  206. end if
  207. else
  208. u_new(i,j) = u(i,j)
  209. end if
  210. end do
  211. end do
  212. ! Rotate V field
  213. do j=vs2,ve2
  214. do i=vs1,ve1
  215. diff = idir * (xlon_v(i,j) - proj_stack(current_nest_number)%stdlon)
  216. if (diff > 180.) then
  217. diff = diff - 360.
  218. else if (diff < -180.) then
  219. diff = diff + 360.
  220. end if
  221. if (proj_stack(current_nest_number)%code == PROJ_LC) then
  222. alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi
  223. else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
  224. if (j == ve2) then
  225. diff = xlon_v(i,j)-xlon_v(i,j-1)
  226. if (diff > 180.) then
  227. diff = diff - 360.
  228. else if (diff < -180.) then
  229. diff = diff + 360.
  230. end if
  231. alpha = atan2( (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg), &
  232. (xlat_v(i,j)-xlat_v(i,j-1))*rad_per_deg &
  233. )
  234. else if (j == vs2) then
  235. diff = xlon_v(i,j+1)-xlon_v(i,j)
  236. if (diff > 180.) then
  237. diff = diff - 360.
  238. else if (diff < -180.) then
  239. diff = diff + 360.
  240. end if
  241. alpha = atan2( (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg), &
  242. (xlat_v(i,j+1)-xlat_v(i,j))*rad_per_deg &
  243. )
  244. else
  245. diff = xlon_v(i,j+1)-xlon_v(i,j-1)
  246. if (diff > 180.) then
  247. diff = diff - 360.
  248. else if (diff < -180.) then
  249. diff = diff + 360.
  250. end if
  251. alpha = atan2( (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg), &
  252. (xlat_v(i,j+1)-xlat_v(i,j-1))*rad_per_deg &
  253. )
  254. end if
  255. else
  256. alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
  257. end if
  258. u_weight = 0.
  259. if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
  260. v_map = v(i,j)
  261. if (j == vs2) then
  262. if (i == ve1 .and. do_last_col_v) then
  263. u_weight = u_mult(i,j)
  264. u_map = u(i,j)*u_mult(i,j)
  265. else
  266. u_weight = u_mult(i,j) + u_mult(i+1,j)
  267. u_map = u(i,j)*u_mult(i,j) + u(i+1,j)*u_mult(i+1,j)
  268. end if
  269. else if (j == ve2 .and. do_last_row_v) then
  270. if (i == ve1 .and. do_last_col_v) then
  271. u_weight = u_mult(i,j-1)
  272. u_map = u(i,j-1)*u_mult(i,j-1)
  273. else
  274. u_weight = u_mult(i,j-1) + u_mult(i+1,j-1)
  275. u_map = u(i,j-1)*u_mult(i,j-1) + u(i+1,j-1)*u_mult(i+1,j-1)
  276. end if
  277. else if (i == ve1 .and. do_last_col_v) then
  278. u_weight = u_mult(i,j) + u_mult(i,j-1)
  279. u_map = u(i,j)*u_mult(i,j) + u(i,j-1)*u_mult(i,j-1)
  280. else
  281. u_weight = u_mult(i,j-1) + u_mult(i,j) + u_mult(i+1,j-1) + u_mult(i+1,j)
  282. u_map = u(i,j-1)*u_mult(i,j-1) + u(i,j)*u_mult(i,j) + u(i+1,j-1)*u_mult(i+1,j-1) + u(i+1,j)*u_mult(i+1,j)
  283. end if
  284. if (u_weight > 0.) then
  285. v_new(i,j) = -sin(alpha)*u_map/u_weight + cos(alpha)*v_map
  286. else
  287. v_new(i,j) = v(i,j)
  288. end if
  289. else
  290. v_new(i,j) = v(i,j)
  291. end if
  292. end do
  293. end do
  294. ! Copy rotated winds back into argument arrays
  295. u = u_new
  296. v = v_new
  297. deallocate(u_new)
  298. deallocate(v_new)
  299. deallocate(u_mult)
  300. deallocate(v_mult)
  301. end if
  302. else
  303. call mprintf(.true.,ERROR,'In metmap_xform(), uninitialized proj_info structure.')
  304. end if
  305. end subroutine metmap_xform
  306. !
  307. ! NMM Wind Rotation Code
  308. !
  309. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  310. ! Name: map_to_met_nmm !
  311. ! !
  312. ! Purpose: Rotate grid-relative winds to Earth-relative winds !
  313. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  314. subroutine map_to_met_nmm(u, u_mask, v, v_mask, &
  315. vs1, vs2, ve1, ve2, &
  316. xlat_v, xlon_v)
  317. implicit none
  318. ! Arguments
  319. integer, intent(in) :: vs1, vs2, ve1, ve2
  320. real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
  321. type (bitarray), intent(in) :: u_mask, v_mask
  322. orig_selected_projection = iget_selected_domain()
  323. call select_domain(SOURCE_PROJ)
  324. call metmap_xform_nmm(u, u_mask, v, v_mask, &
  325. vs1, vs2, ve1, ve2, &
  326. xlat_v, xlon_v, 1)
  327. call select_domain(orig_selected_projection)
  328. end subroutine map_to_met_nmm
  329. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  330. ! Name: met_to_map_nmm !
  331. ! !
  332. ! Purpose: Rotate Earth-relative winds to grid-relative winds !
  333. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  334. subroutine met_to_map_nmm(u, u_mask, v, v_mask, &
  335. vs1, vs2, ve1, ve2, &
  336. xlat_v, xlon_v)
  337. implicit none
  338. ! Arguments
  339. integer, intent(in) :: vs1, vs2, ve1, ve2
  340. real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
  341. type (bitarray), intent(in) :: u_mask, v_mask
  342. orig_selected_projection = iget_selected_domain()
  343. call select_domain(1)
  344. call metmap_xform_nmm(u, u_mask, v, v_mask, &
  345. vs1, vs2, ve1, ve2, &
  346. xlat_v, xlon_v, -1)
  347. call select_domain(orig_selected_projection)
  348. end subroutine met_to_map_nmm
  349. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  350. ! Name: metmap_xform_nmm !
  351. ! !
  352. ! Purpose: Do the actual work of rotating winds for E grid. !
  353. ! If idir= 1, rotate grid-relative winds to Earth-relative winds !
  354. ! If idir=-1, rotate Earth-relative winds to grid-relative winds !
  355. ! !
  356. ! ASSUMPTIONS: 1) MEMORY ORDER IS XY. !
  357. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  358. subroutine metmap_xform_nmm(u, u_mask, v, v_mask, &
  359. vs1, vs2, ve1, ve2, &
  360. xlat_v, xlon_v, idir)
  361. implicit none
  362. ! Arguments
  363. integer, intent(in) :: vs1, vs2, ve1, ve2, idir
  364. real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
  365. type (bitarray), intent(in) :: u_mask, v_mask
  366. ! Local variables
  367. integer :: i, j
  368. real :: u_map, v_map, diff, alpha
  369. real :: phi0, lmbd0, big_denominator, relm, rlat_v,rlon_v, clontemp
  370. real :: sin_phi0, cos_phi0, cos_alpha, sin_alpha
  371. real, pointer, dimension(:,:) :: u_new, v_new
  372. ! If the proj_info structure has not been initialized, we don't have
  373. ! information about the projection and standard longitude.
  374. if (proj_stack(current_nest_number)%init) then
  375. if (proj_stack(current_nest_number)%code == PROJ_ROTLL) then
  376. call mprintf((idir == 1),LOGFILE,'Rotating map winds to earth winds.')
  377. call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
  378. ! Create arrays to hold rotated winds
  379. allocate(u_new(vs1:ve1, vs2:ve2))
  380. allocate(v_new(vs1:ve1, vs2:ve2))
  381. phi0 = proj_stack(current_nest_number)%lat1 * rad_per_deg
  382. clontemp= proj_stack(current_nest_number)%lon1
  383. if (clontemp .lt. 0.) then
  384. lmbd0 = (clontemp + 360.) * rad_per_deg
  385. else
  386. lmbd0 = (clontemp) * rad_per_deg
  387. endif
  388. sin_phi0 = sin(phi0)
  389. cos_phi0 = cos(phi0)
  390. do j=vs2,ve2
  391. do i=vs1,ve1
  392. ! Calculate the sine and cosine of rotation angle
  393. rlat_v = xlat_v(i,j) * rad_per_deg
  394. rlon_v = xlon_v(i,j) * rad_per_deg
  395. relm = rlon_v - lmbd0
  396. big_denominator = cos(asin( &
  397. cos_phi0 * sin(rlat_v) - &
  398. sin_phi0 * cos(rlat_v) * cos(relm) &
  399. ) )
  400. sin_alpha = sin_phi0 * sin(relm) / &
  401. big_denominator
  402. cos_alpha = (cos_phi0 * cos(rlat_v) + &
  403. sin_phi0 * sin(rlat_v) * cos(relm)) / &
  404. big_denominator
  405. ! Rotate U field
  406. if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
  407. u_map = u(i,j)
  408. if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
  409. v_map = v(i,j)
  410. else
  411. v_map = 0.
  412. end if
  413. u_new(i,j) = cos_alpha*u_map + idir*sin_alpha*v_map
  414. else
  415. u_new(i,j) = u(i,j)
  416. end if
  417. ! Rotate V field
  418. if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
  419. v_map = v(i,j)
  420. if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
  421. u_map = u(i,j)
  422. else
  423. u_map = 0.
  424. end if
  425. v_new(i,j) = -idir*sin_alpha*u_map + cos_alpha*v_map
  426. else
  427. v_new(i,j) = v(i,j)
  428. end if
  429. end do
  430. end do
  431. ! Copy rotated winds back into argument arrays
  432. u = u_new
  433. v = v_new
  434. deallocate(u_new)
  435. deallocate(v_new)
  436. ! Only rotate winds for Lambert conformal, polar stereographic, or Cassini
  437. else if ((proj_stack(current_nest_number)%code == PROJ_LC) .or. &
  438. (proj_stack(current_nest_number)%code == PROJ_PS) .or. &
  439. (proj_stack(current_nest_number)%code == PROJ_CASSINI)) then
  440. call mprintf((idir == 1),LOGFILE,'Rotating map winds to earth winds.')
  441. call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
  442. ! Create arrays to hold rotated winds
  443. allocate(u_new(vs1:ve1, vs2:ve2))
  444. allocate(v_new(vs1:ve1, vs2:ve2))
  445. do j=vs2,ve2
  446. do i=vs1,ve1
  447. diff = idir * (xlon_v(i,j) - proj_stack(current_nest_number)%stdlon)
  448. if (diff > 180.) then
  449. diff = diff - 360.
  450. else if (diff < -180.) then
  451. diff = diff + 360.
  452. end if
  453. ! Calculate the rotation angle, alpha, in radians
  454. if (proj_stack(current_nest_number)%code == PROJ_LC) then
  455. alpha = diff * proj_stack(current_nest_number)%cone * &
  456. rad_per_deg * proj_stack(current_nest_number)%hemi
  457. else
  458. alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
  459. end if
  460. ! Rotate U field
  461. if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
  462. u_map = u(i,j)
  463. if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
  464. v_map = v(i,j)
  465. else
  466. v_map = 0.
  467. end if
  468. u_new(i,j) = cos(alpha)*u_map + idir*sin(alpha)*v_map
  469. else
  470. u_new(i,j) = u(i,j)
  471. end if
  472. ! Rotate V field
  473. if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
  474. v_map = v(i,j)
  475. if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
  476. u_map = u(i,j)
  477. else
  478. u_map = 0.
  479. end if
  480. v_new(i,j) = -idir*sin(alpha)*u_map + cos(alpha)*v_map
  481. else
  482. v_new(i,j) = v(i,j)
  483. end if
  484. end do
  485. end do
  486. ! Copy rotated winds back into argument arrays
  487. u = u_new
  488. v = v_new
  489. deallocate(u_new)
  490. deallocate(v_new)
  491. end if
  492. else
  493. call mprintf(.true.,ERROR,'In metmap_xform_nmm(), uninitialized proj_info structure.')
  494. end if
  495. end subroutine metmap_xform_nmm
  496. end module rotate_winds_module