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                   diff = idir * (xlon_u(i,j) - proj_stack(current_nest_number)%stdlon)
 165                   if (diff > 180.) then
 166                      diff = diff - 360.
 167                   else if (diff < -180.) then
 168                      diff = diff + 360.
 169                   end if
 170
 171                   ! Calculate the rotation angle, alpha, in radians
 172                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
 173                      alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi
 174                   else if ( (proj_stack(SOURCE_PROJ)%code == PROJ_CASSINI) .and. ( idir == 1 ) ) then
 175                         call latlon_to_ij ( proj_stack(SOURCE_PROJ) , &
 176                                             xlat_u(i,j), xlon_u(i,j), x_u, y_u )
 177                         call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
 178                                             x_u, y_u + 0.1 , xlat_u_p1 , xlon_u_p1 )
 179                         call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
 180                                             x_u, y_u - 0.1 , xlat_u_m1 , xlon_u_m1 )
 181                         diff_lon = xlon_u_p1-xlon_u_m1
 182                         if (diff_lon > 180.) then
 183                            diff_lon = diff_lon - 360.
 184                         else if (diff_lon < -180.) then
 185                            diff_lon = diff_lon + 360.
 186                         end if
 187                         diff_lat = xlat_u_p1-xlat_u_m1
 188                         alpha =-atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
 189                                                         diff_lat*rad_per_deg    &
 190                                      )
 191                   else if ( (proj_stack(current_nest_number)%code == PROJ_CASSINI) .and. ( idir == -1 ) ) then
 192                         if (j == ue2) then
 193                            diff = xlon_u(i,j)-xlon_u(i,j-1)
 194                            if (diff > 180.) then
 195                               diff = diff - 360.
 196                            else if (diff < -180.) then
 197                               diff = diff + 360.
 198                            end if
 199                            alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &
 200                                                            (xlat_u(i,j)-xlat_u(i,j-1))*rad_per_deg    &
 201                                         )
 202                         else if (j == us2) then
 203                            diff = xlon_u(i,j+1)-xlon_u(i,j)
 204                            if (diff > 180.) then
 205                               diff = diff - 360.
 206                            else if (diff < -180.) then
 207                               diff = diff + 360.
 208                            end if
 209                            alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &
 210                                                            (xlat_u(i,j+1)-xlat_u(i,j))*rad_per_deg    &
 211                                         )
 212                         else
 213                            diff = xlon_u(i,j+1)-xlon_u(i,j-1)
 214                            if (diff > 180.) then
 215                               diff = diff - 360.
 216                            else if (diff < -180.) then
 217                               diff = diff + 360.
 218                            end if
 219                            alpha = atan2(   (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg),   &
 220                                                            (xlat_u(i,j+1)-xlat_u(i,j-1))*rad_per_deg    &
 221                                         )
 222                         end if
 223                   else
 224                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
 225                   end if
 226
 227                   v_weight = 0.
 228
 229                   ! On C grid, take U_ij, and get V value at the same lat/lon
 230                   !   by averaging the four surrounding V points
 231                   if (bitarray_test(u_mask, i-us1+1, j-us2+1)) then
 232                      u_map = u(i,j)
 233                      if (i == us1) then
 234                         if (j == ue2 .and. do_last_row_u) then
 235                            v_weight = v_mult(i,j)
 236                            v_map = v(i,j)*v_mult(i,j)
 237                         else
 238                            v_weight = v_mult(i,j) + v_mult(i,j+1)
 239                            v_map = v(i,j)*v_mult(i,j) + v(i,j+1)*v_mult(i,j+1)
 240                         end if
 241                      else if (i == ue1 .and. do_last_col_u) then
 242                         if (j == ue2 .and. do_last_row_u) then
 243                            v_weight = v_mult(i-1,j)
 244                            v_map = v(i-1,j)
 245                         else
 246                            v_weight = v_mult(i-1,j) + v_mult(i-1,j+1)
 247                            v_map = v(i-1,j)*v_mult(i-1,j) + v(i-1,j+1)*v_mult(i-1,j+1)
 248                         end if
 249                      else if (j == ue2 .and. do_last_row_u) then
 250                         v_weight = v_mult(i-1,j-1) + v_mult(i,j-1)
 251                         v_map = v(i-1,j-1)*v_mult(i-1,j-1) + v(i,j-1)*v_mult(i,j-1)
 252                      else
 253                         v_weight = v_mult(i-1,j) + v_mult(i-1,j+1) + v_mult(i,j) + v_mult(i,j+1)
 254                         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)
 255                      end if
 256                      if (v_weight > 0.) then
 257                         u_new(i,j) = cos(alpha)*u_map + sin(alpha)*v_map/v_weight
 258                      else
 259                         u_new(i,j) = u(i,j)
 260                      end if
 261                   else
 262                      u_new(i,j) = u(i,j)
 263                   end if
 264
 265                end do
 266             end do
 267
 268             ! Rotate V field
 269             do j=vs2,ve2
 270                do i=vs1,ve1
 271
 272                   diff = idir * (xlon_v(i,j) - proj_stack(current_nest_number)%stdlon)
 273                   if (diff > 180.) then
 274                      diff = diff - 360.
 275                   else if (diff < -180.) then
 276                      diff = diff + 360.
 277                   end if
 278
 279                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
 280                      alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi
 281                   else if ( (proj_stack(SOURCE_PROJ)%code == PROJ_CASSINI) .and. ( idir == 1 ) ) then
 282                         call latlon_to_ij ( proj_stack(SOURCE_PROJ) , &
 283                                             xlat_v(i,j), xlon_v(i,j), x_v, y_v )
 284                         call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
 285                                             x_v, y_v + 0.1 , xlat_v_p1 , xlon_v_p1 )
 286                         call ij_to_latlon ( proj_stack(SOURCE_PROJ) , &
 287                                             x_v, y_v - 0.1 , xlat_v_m1 , xlon_v_m1 )
 288                         diff_lon = xlon_v_p1-xlon_v_m1
 289                         if (diff_lon > 180.) then
 290                            diff_lon = diff_lon - 360.
 291                         else if (diff_lon < -180.) then
 292                            diff_lon = diff_lon + 360.
 293                         end if
 294                         diff_lat = xlat_v_p1-xlat_v_m1
 295                         alpha =-atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff_lon*rad_per_deg),   &
 296                                                         diff_lat*rad_per_deg    &
 297                                      )
 298                   else if ( (proj_stack(current_nest_number)%code == PROJ_CASSINI) .and. ( idir == -1 ) ) then
 299                         if (j == ve2) then
 300                            diff = xlon_v(i,j)-xlon_v(i,j-1)
 301                            if (diff > 180.) then
 302                               diff = diff - 360.
 303                            else if (diff < -180.) then
 304                               diff = diff + 360.
 305                            end if
 306                            alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
 307                                                            (xlat_v(i,j)-xlat_v(i,j-1))*rad_per_deg    &
 308                                         )
 309                         else if (j == vs2) then
 310                            diff = xlon_v(i,j+1)-xlon_v(i,j)
 311                            if (diff > 180.) then
 312                               diff = diff - 360.
 313                            else if (diff < -180.) then
 314                               diff = diff + 360.
 315                            end if
 316                            alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
 317                                                            (xlat_v(i,j+1)-xlat_v(i,j))*rad_per_deg    &
 318                                         )
 319                         else
 320                            diff = xlon_v(i,j+1)-xlon_v(i,j-1)
 321                            if (diff > 180.) then
 322                               diff = diff - 360.
 323                            else if (diff < -180.) then
 324                               diff = diff + 360.
 325                            end if
 326                            alpha = atan2(   (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg),   &
 327                                                            (xlat_v(i,j+1)-xlat_v(i,j-1))*rad_per_deg    &
 328                                         )
 329                         end if
 330                   else
 331                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
 332                   end if
 333
 334                   u_weight = 0.
 335
 336                   if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
 337                      v_map = v(i,j)
 338                      if (j == vs2) then
 339                         if (i == ve1 .and. do_last_col_v) then
 340                            u_weight = u_mult(i,j)
 341                            u_map = u(i,j)*u_mult(i,j)
 342                         else
 343                            u_weight = u_mult(i,j) + u_mult(i+1,j)
 344                            u_map = u(i,j)*u_mult(i,j) + u(i+1,j)*u_mult(i+1,j)
 345                         end if
 346                      else if (j == ve2 .and. do_last_row_v) then
 347                         if (i == ve1 .and. do_last_col_v) then
 348                            u_weight = u_mult(i,j-1)
 349                            u_map = u(i,j-1)*u_mult(i,j-1)
 350                         else
 351                            u_weight = u_mult(i,j-1) + u_mult(i+1,j-1)
 352                            u_map = u(i,j-1)*u_mult(i,j-1) + u(i+1,j-1)*u_mult(i+1,j-1)
 353                         end if
 354                      else if (i == ve1 .and. do_last_col_v) then
 355                         u_weight = u_mult(i,j) + u_mult(i,j-1)
 356                         u_map = u(i,j)*u_mult(i,j) + u(i,j-1)*u_mult(i,j-1)
 357                      else
 358                         u_weight = u_mult(i,j-1) + u_mult(i,j) + u_mult(i+1,j-1) + u_mult(i+1,j)
 359                         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)
 360                      end if
 361                      if (u_weight > 0.) then
 362                         v_new(i,j) = -sin(alpha)*u_map/u_weight + cos(alpha)*v_map
 363                      else
 364                         v_new(i,j) = v(i,j)
 365                      end if
 366                   else
 367                      v_new(i,j) = v(i,j)
 368                   end if
 369
 370                end do
 371             end do
 372
 373             ! Copy rotated winds back into argument arrays
 374             u = u_new
 375             v = v_new
 376
 377             deallocate(u_new)
 378             deallocate(v_new)
 379             deallocate(u_mult)
 380             deallocate(v_mult)
 381          end if
 382
 383       else
 384          call mprintf(.true.,ERROR,'In metmap_xform(), uninitialized proj_info structure.')
 385       end if
 386
 387    end subroutine metmap_xform
 388
 389
 390    !
 391    ! NMM Wind Rotation Code
 392    !
 393
 394    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 395    ! Name: map_to_met_nmm                                                         !
 396    !                                                                              !
 397    ! Purpose: Rotate grid-relative winds to Earth-relative winds                  !
 398    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 399    subroutine map_to_met_nmm(u, u_mask, v, v_mask, &
 400                              vs1, vs2, ve1, ve2, &
 401                              xlat_v, xlon_v)
 402
 403       implicit none
 404
 405       ! Arguments
 406       integer, intent(in) :: vs1, vs2, ve1, ve2
 407       real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
 408       type (bitarray), intent(in) :: u_mask, v_mask
 409
 410       orig_selected_projection = iget_selected_domain()
 411       call select_domain(SOURCE_PROJ)
 412       call metmap_xform_nmm(u, u_mask, v, v_mask, &
 413                             vs1, vs2, ve1, ve2, &
 414                             xlat_v, xlon_v, 1)
 415       call select_domain(orig_selected_projection)
 416
 417    end subroutine map_to_met_nmm
 418
 419
 420    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 421    ! Name: met_to_map_nmm                                                         !
 422    !                                                                              !
 423    ! Purpose: Rotate Earth-relative winds to grid-relative winds                  !
 424    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 425    subroutine met_to_map_nmm(u, u_mask, v, v_mask, &
 426                              vs1, vs2, ve1, ve2, &
 427                              xlat_v, xlon_v)
 428
 429       implicit none
 430
 431       ! Arguments
 432       integer, intent(in) :: vs1, vs2, ve1, ve2
 433       real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
 434       type (bitarray), intent(in) :: u_mask, v_mask
 435
 436       orig_selected_projection = iget_selected_domain()
 437       call select_domain(1)
 438       call metmap_xform_nmm(u, u_mask, v, v_mask, &
 439                             vs1, vs2, ve1, ve2, &
 440                             xlat_v, xlon_v, -1)
 441       call select_domain(orig_selected_projection)
 442
 443    end subroutine met_to_map_nmm
 444
 445
 446    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 447    ! Name: metmap_xform_nmm                                                       !
 448    !                                                                              !
 449    ! Purpose: Do the actual work of rotating winds for E grid.                    !
 450    !          If idir= 1, rotate grid-relative winds to Earth-relative winds      !
 451    !          If idir=-1, rotate Earth-relative winds to grid-relative winds      !
 452    !                                                                              !
 453    ! ASSUMPTIONS: 1) MEMORY ORDER IS XY.                                          !
 454    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 455    subroutine metmap_xform_nmm(u, u_mask, v, v_mask, &
 456                                vs1, vs2, ve1, ve2, &
 457                                xlat_v, xlon_v, idir)
 458
 459       implicit none
 460
 461       ! Arguments
 462       integer, intent(in) :: vs1, vs2, ve1, ve2, idir
 463       real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
 464       type (bitarray), intent(in) :: u_mask, v_mask
 465
 466       ! Local variables
 467       integer :: i, j
 468       real :: u_map, v_map, diff, alpha
 469       real :: phi0, lmbd0, big_denominator, relm, rlat_v,rlon_v, clontemp
 470       real :: sin_phi0, cos_phi0,  cos_alpha, sin_alpha
 471       real, pointer, dimension(:,:) :: u_new, v_new
 472
 473
 474       ! If the proj_info structure has not been initialized, we don't have
 475       !   information about the projection and standard longitude.
 476       if (proj_stack(current_nest_number)%init) then
 477
 478          if (proj_stack(current_nest_number)%code == PROJ_ROTLL) then
 479
 480             call mprintf((idir ==  1),LOGFILE,'Rotating map winds to earth winds.')
 481             call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
 482
 483             ! Create arrays to hold rotated winds
 484             allocate(u_new(vs1:ve1, vs2:ve2))
 485             allocate(v_new(vs1:ve1, vs2:ve2))
 486
 487             phi0  = proj_stack(current_nest_number)%lat1 * rad_per_deg
 488
 489             clontemp= proj_stack(current_nest_number)%lon1
 490
 491             if (clontemp .lt. 0.) then
 492                lmbd0 = (clontemp + 360.) * rad_per_deg
 493             else
 494                lmbd0 = (clontemp) * rad_per_deg
 495             endif
 496
 497             sin_phi0  = sin(phi0)
 498             cos_phi0  = cos(phi0)
 499
 500             do j=vs2,ve2
 501                do i=vs1,ve1
 502
 503                   ! Calculate the sine and cosine of rotation angle
 504                   rlat_v = xlat_v(i,j) * rad_per_deg
 505                   rlon_v = xlon_v(i,j) * rad_per_deg
 506                   relm = rlon_v - lmbd0
 507                   big_denominator = cos(asin( &
 508                                        cos_phi0 * sin(rlat_v) - &
 509                                        sin_phi0 * cos(rlat_v) * cos(relm) &
 510                                         )   )
 511
 512                   sin_alpha = sin_phi0 * sin(relm)  /  &
 513                                          big_denominator
 514
 515                   cos_alpha = (cos_phi0 * cos(rlat_v) + &
 516                                sin_phi0 * sin(rlat_v) * cos(relm))  /  &
 517                                             big_denominator
 518
 519                   ! Rotate U field
 520                   if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
 521                      u_map = u(i,j)
 522                      if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
 523                         v_map = v(i,j)
 524                      else
 525                         v_map = 0.
 526                      end if
 527
 528                      u_new(i,j) = cos_alpha*u_map + idir*sin_alpha*v_map
 529                   else
 530                      u_new(i,j) = u(i,j)
 531                   end if
 532
 533                   ! Rotate V field
 534                   if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
 535                      v_map = v(i,j)
 536                      if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
 537                         u_map = u(i,j)
 538                      else
 539                         u_map = 0.
 540                      end if
 541
 542                      v_new(i,j) = -idir*sin_alpha*u_map + cos_alpha*v_map
 543                   else
 544                      v_new(i,j) = v(i,j)
 545                   end if
 546
 547                end do
 548             end do
 549
 550             ! Copy rotated winds back into argument arrays
 551             u = u_new
 552             v = v_new
 553
 554             deallocate(u_new)
 555             deallocate(v_new)
 556
 557          ! Only rotate winds for Lambert conformal, polar stereographic, or Cassini
 558          else if ((proj_stack(current_nest_number)%code == PROJ_LC) .or. &
 559                   (proj_stack(current_nest_number)%code == PROJ_PS) .or. &
 560                   (proj_stack(current_nest_number)%code == PROJ_CASSINI)) then
 561
 562             call mprintf((idir ==  1),LOGFILE,'Rotating map winds to earth winds.')
 563             call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
 564
 565             ! Create arrays to hold rotated winds
 566             allocate(u_new(vs1:ve1, vs2:ve2))
 567             allocate(v_new(vs1:ve1, vs2:ve2))
 568
 569             do j=vs2,ve2
 570                do i=vs1,ve1
 571
 572                   diff = idir * (xlon_v(i,j) - proj_stack(current_nest_number)%stdlon)
 573                   if (diff > 180.) then
 574                      diff = diff - 360.
 575                   else if (diff < -180.) then
 576                      diff = diff + 360.
 577                   end if
 578
 579                   ! Calculate the rotation angle, alpha, in radians
 580                   if (proj_stack(current_nest_number)%code == PROJ_LC) then
 581                      alpha = diff * proj_stack(current_nest_number)%cone * &
 582                              rad_per_deg * proj_stack(current_nest_number)%hemi
 583                   else
 584                      alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
 585                   end if
 586
 587                   ! Rotate U field
 588                   if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
 589                      u_map = u(i,j)
 590                      if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
 591                         v_map = v(i,j)
 592                      else
 593                         v_map = 0.
 594                      end if
 595
 596                      u_new(i,j) = cos(alpha)*u_map + idir*sin(alpha)*v_map
 597                   else
 598                      u_new(i,j) = u(i,j)
 599                   end if
 600
 601                   ! Rotate V field
 602                   if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
 603                      v_map = v(i,j)
 604                      if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
 605                         u_map = u(i,j)
 606                      else
 607                         u_map = 0.
 608                      end if
 609
 610                      v_new(i,j) = -idir*sin(alpha)*u_map + cos(alpha)*v_map
 611                   else
 612                      v_new(i,j) = v(i,j)
 613                   end if
 614
 615                end do
 616             end do
 617
 618             ! Copy rotated winds back into argument arrays
 619             u = u_new
 620             v = v_new
 621
 622             deallocate(u_new)
 623             deallocate(v_new)
 624
 625          end if
 626
 627       else
 628          call mprintf(.true.,ERROR,'In metmap_xform_nmm(), uninitialized proj_info structure.')
 629       end if
 630
 631    end subroutine metmap_xform_nmm
 632
 633 end module rotate_winds_module