WPS/metgrid/src/rotate_winds_module.f90

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