metgrid/src/rotate_winds_module.F

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