/WPS/metgrid/src/rotate_winds_module.F
FORTRAN Legacy | 595 lines | 394 code | 99 blank | 102 comment | 97 complexity | 35175581bc2258903df5d12db8ed1c82 MD5 | raw file
Possible License(s): AGPL-1.0
- module rotate_winds_module
- use bitarray_module
- use constants_module
- use llxy_module
- use misc_definitions_module
- use module_debug
- integer :: orig_selected_projection
- contains
- !
- ! ARW Wind Rotation Code
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: map_to_met !
- ! !
- ! Purpose: Rotate grid-relative winds to Earth-relative winds !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine map_to_met(u, u_mask, v, v_mask, &
- us1, us2, ue1, ue2, &
- vs1, vs2, ve1, ve2, &
- xlon_u, xlon_v, xlat_u, xlat_v)
- implicit none
- ! Arguments
- integer, intent(in) :: us1, us2, ue1, ue2, vs1, vs2, ve1, ve2
- real, pointer, dimension(:,:) :: u, v, xlon_u, xlon_v, xlat_u, xlat_v
- type (bitarray), intent(in) :: u_mask, v_mask
- orig_selected_projection = iget_selected_domain()
- call select_domain(SOURCE_PROJ)
- call metmap_xform(u, u_mask, v, v_mask, &
- us1, us2, ue1, ue2, &
- vs1, vs2, ve1, ve2, &
- xlon_u, xlon_v, xlat_u, xlat_v, 1)
- call select_domain(orig_selected_projection)
- end subroutine map_to_met
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: met_to_map !
- ! !
- ! Purpose: Rotate Earth-relative winds to grid-relative winds !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine met_to_map(u, u_mask, v, v_mask, &
- us1, us2, ue1, ue2, &
- vs1, vs2, ve1, ve2, &
- xlon_u, xlon_v, xlat_u, xlat_v)
- implicit none
- ! Arguments
- integer, intent(in) :: us1, us2, ue1, ue2, vs1, vs2, ve1, ve2
- real, pointer, dimension(:,:) :: u, v, xlon_u, xlon_v, xlat_u, xlat_v
- type (bitarray), intent(in) :: u_mask, v_mask
- orig_selected_projection = iget_selected_domain()
- call select_domain(1)
- call metmap_xform(u, u_mask, v, v_mask, &
- us1, us2, ue1, ue2, &
- vs1, vs2, ve1, ve2, &
- xlon_u, xlon_v, xlat_u, xlat_v, -1)
- call select_domain(orig_selected_projection)
- end subroutine met_to_map
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: metmap_xform !
- ! !
- ! Purpose: Do the actual work of rotating winds for C grid. !
- ! If idir= 1, rotate grid-relative winds to Earth-relative winds !
- ! If idir=-1, rotate Earth-relative winds to grid-relative winds !
- ! !
- ! ASSUMPTIONS: 1) MEMORY ORDER IS XY. !
- ! 2) U ARRAY HAS ONE MORE COLUMN THAN THE V ARRAY, AND V ARRAY !
- ! HAS ONE MORE ROW THAN U ARRAY. !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine metmap_xform(u, u_mask, v, v_mask, &
- us1, us2, ue1, ue2, vs1, vs2, ve1, ve2, &
- xlon_u, xlon_v, xlat_u, xlat_v, idir)
- implicit none
- ! Arguments
- integer, intent(in) :: us1, us2, ue1, ue2, vs1, vs2, ve1, ve2, idir
- real, pointer, dimension(:,:) :: u, v, xlon_u, xlon_v, xlat_u, xlat_v
- type (bitarray), intent(in) :: u_mask, v_mask
- ! Local variables
- integer :: i, j
- real :: u_weight, v_weight
- real :: u_map, v_map, alpha, diff
- real, pointer, dimension(:,:) :: u_new, v_new, u_mult, v_mult
- logical :: do_last_col_u, do_last_row_u, do_last_col_v, do_last_row_v
- ! If the proj_info structure has not been initialized, we don't have
- ! information about the projection and standard longitude.
- if (proj_stack(current_nest_number)%init) then
- ! Only rotate winds for Lambert conformal, polar stereographic, or Cassini
- if ((proj_stack(current_nest_number)%code == PROJ_LC) .or. &
- (proj_stack(current_nest_number)%code == PROJ_PS) .or. &
- (proj_stack(current_nest_number)%code == PROJ_CASSINI)) then
- call mprintf((idir == 1),LOGFILE,'Rotating map winds to earth winds.')
- call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
- allocate(u_mult(us1:ue1,us2:ue2))
- allocate(v_mult(vs1:ve1,vs2:ve2))
- do j=us2,ue2
- do i=us1,ue1
- if (bitarray_test(u_mask, i-us1+1, j-us2+1)) then
- u_mult(i,j) = 1.
- else
- u_mult(i,j) = 0.
- end if
- end do
- end do
- do j=vs2,ve2
- do i=vs1,ve1
- if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
- v_mult(i,j) = 1.
- else
- v_mult(i,j) = 0.
- end if
- end do
- end do
- if (ue1-us1 == ve1-vs1) then
- do_last_col_u = .false.
- do_last_col_v = .true.
- else
- do_last_col_u = .true.
- do_last_col_v = .false.
- end if
- if (ue2-us2 == ve2-vs2) then
- do_last_row_u = .true.
- do_last_row_v = .false.
- else
- do_last_row_u = .false.
- do_last_row_v = .true.
- end if
- ! Create arrays to hold rotated winds
- allocate(u_new(us1:ue1, us2:ue2))
- allocate(v_new(vs1:ve1, vs2:ve2))
- ! Rotate U field
- do j=us2,ue2
- do i=us1,ue1
- diff = idir * (xlon_u(i,j) - proj_stack(current_nest_number)%stdlon)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- ! Calculate the rotation angle, alpha, in radians
- if (proj_stack(current_nest_number)%code == PROJ_LC) then
- alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi
- else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
- if (j == ue2) then
- diff = xlon_u(i,j)-xlon_u(i,j-1)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- alpha = atan2( (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg), &
- (xlat_u(i,j)-xlat_u(i,j-1))*rad_per_deg &
- )
- else if (j == us2) then
- diff = xlon_u(i,j+1)-xlon_u(i,j)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- alpha = atan2( (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg), &
- (xlat_u(i,j+1)-xlat_u(i,j))*rad_per_deg &
- )
- else
- diff = xlon_u(i,j+1)-xlon_u(i,j-1)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- alpha = atan2( (-cos(xlat_u(i,j)*rad_per_deg) * diff*rad_per_deg), &
- (xlat_u(i,j+1)-xlat_u(i,j-1))*rad_per_deg &
- )
- end if
- else
- alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
- end if
-
- v_weight = 0.
- ! On C grid, take U_ij, and get V value at the same lat/lon
- ! by averaging the four surrounding V points
- if (bitarray_test(u_mask, i-us1+1, j-us2+1)) then
- u_map = u(i,j)
- if (i == us1) then
- if (j == ue2 .and. do_last_row_u) then
- v_weight = v_mult(i,j)
- v_map = v(i,j)*v_mult(i,j)
- else
- v_weight = v_mult(i,j) + v_mult(i,j+1)
- v_map = v(i,j)*v_mult(i,j) + v(i,j+1)*v_mult(i,j+1)
- end if
- else if (i == ue1 .and. do_last_col_u) then
- if (j == ue2 .and. do_last_row_u) then
- v_weight = v_mult(i-1,j)
- v_map = v(i-1,j)
- else
- v_weight = v_mult(i-1,j) + v_mult(i-1,j+1)
- v_map = v(i-1,j)*v_mult(i-1,j) + v(i-1,j+1)*v_mult(i-1,j+1)
- end if
- else if (j == ue2 .and. do_last_row_u) then
- v_weight = v_mult(i-1,j-1) + v_mult(i,j-1)
- v_map = v(i-1,j-1)*v_mult(i-1,j-1) + v(i,j-1)*v_mult(i,j-1)
- else
- v_weight = v_mult(i-1,j) + v_mult(i-1,j+1) + v_mult(i,j) + v_mult(i,j+1)
- 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)
- end if
- if (v_weight > 0.) then
- u_new(i,j) = cos(alpha)*u_map + sin(alpha)*v_map/v_weight
- else
- u_new(i,j) = u(i,j)
- end if
- else
- u_new(i,j) = u(i,j)
- end if
- end do
- end do
- ! Rotate V field
- do j=vs2,ve2
- do i=vs1,ve1
- diff = idir * (xlon_v(i,j) - proj_stack(current_nest_number)%stdlon)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- if (proj_stack(current_nest_number)%code == PROJ_LC) then
- alpha = diff * proj_stack(current_nest_number)%cone * rad_per_deg * proj_stack(current_nest_number)%hemi
- else if (proj_stack(current_nest_number)%code == PROJ_CASSINI) then
- if (j == ve2) then
- diff = xlon_v(i,j)-xlon_v(i,j-1)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- alpha = atan2( (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg), &
- (xlat_v(i,j)-xlat_v(i,j-1))*rad_per_deg &
- )
- else if (j == vs2) then
- diff = xlon_v(i,j+1)-xlon_v(i,j)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- alpha = atan2( (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg), &
- (xlat_v(i,j+1)-xlat_v(i,j))*rad_per_deg &
- )
- else
- diff = xlon_v(i,j+1)-xlon_v(i,j-1)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- alpha = atan2( (-cos(xlat_v(i,j)*rad_per_deg) * diff*rad_per_deg), &
- (xlat_v(i,j+1)-xlat_v(i,j-1))*rad_per_deg &
- )
- end if
- else
- alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
- end if
- u_weight = 0.
- if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
- v_map = v(i,j)
- if (j == vs2) then
- if (i == ve1 .and. do_last_col_v) then
- u_weight = u_mult(i,j)
- u_map = u(i,j)*u_mult(i,j)
- else
- u_weight = u_mult(i,j) + u_mult(i+1,j)
- u_map = u(i,j)*u_mult(i,j) + u(i+1,j)*u_mult(i+1,j)
- end if
- else if (j == ve2 .and. do_last_row_v) then
- if (i == ve1 .and. do_last_col_v) then
- u_weight = u_mult(i,j-1)
- u_map = u(i,j-1)*u_mult(i,j-1)
- else
- u_weight = u_mult(i,j-1) + u_mult(i+1,j-1)
- u_map = u(i,j-1)*u_mult(i,j-1) + u(i+1,j-1)*u_mult(i+1,j-1)
- end if
- else if (i == ve1 .and. do_last_col_v) then
- u_weight = u_mult(i,j) + u_mult(i,j-1)
- u_map = u(i,j)*u_mult(i,j) + u(i,j-1)*u_mult(i,j-1)
- else
- u_weight = u_mult(i,j-1) + u_mult(i,j) + u_mult(i+1,j-1) + u_mult(i+1,j)
- 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)
- end if
- if (u_weight > 0.) then
- v_new(i,j) = -sin(alpha)*u_map/u_weight + cos(alpha)*v_map
- else
- v_new(i,j) = v(i,j)
- end if
- else
- v_new(i,j) = v(i,j)
- end if
- end do
- end do
- ! Copy rotated winds back into argument arrays
- u = u_new
- v = v_new
- deallocate(u_new)
- deallocate(v_new)
- deallocate(u_mult)
- deallocate(v_mult)
- end if
- else
- call mprintf(.true.,ERROR,'In metmap_xform(), uninitialized proj_info structure.')
- end if
-
- end subroutine metmap_xform
- !
- ! NMM Wind Rotation Code
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: map_to_met_nmm !
- ! !
- ! Purpose: Rotate grid-relative winds to Earth-relative winds !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine map_to_met_nmm(u, u_mask, v, v_mask, &
- vs1, vs2, ve1, ve2, &
- xlat_v, xlon_v)
- implicit none
- ! Arguments
- integer, intent(in) :: vs1, vs2, ve1, ve2
- real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
- type (bitarray), intent(in) :: u_mask, v_mask
- orig_selected_projection = iget_selected_domain()
- call select_domain(SOURCE_PROJ)
- call metmap_xform_nmm(u, u_mask, v, v_mask, &
- vs1, vs2, ve1, ve2, &
- xlat_v, xlon_v, 1)
- call select_domain(orig_selected_projection)
- end subroutine map_to_met_nmm
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: met_to_map_nmm !
- ! !
- ! Purpose: Rotate Earth-relative winds to grid-relative winds !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine met_to_map_nmm(u, u_mask, v, v_mask, &
- vs1, vs2, ve1, ve2, &
- xlat_v, xlon_v)
- implicit none
- ! Arguments
- integer, intent(in) :: vs1, vs2, ve1, ve2
- real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
- type (bitarray), intent(in) :: u_mask, v_mask
- orig_selected_projection = iget_selected_domain()
- call select_domain(1)
- call metmap_xform_nmm(u, u_mask, v, v_mask, &
- vs1, vs2, ve1, ve2, &
- xlat_v, xlon_v, -1)
- call select_domain(orig_selected_projection)
- end subroutine met_to_map_nmm
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: metmap_xform_nmm !
- ! !
- ! Purpose: Do the actual work of rotating winds for E grid. !
- ! If idir= 1, rotate grid-relative winds to Earth-relative winds !
- ! If idir=-1, rotate Earth-relative winds to grid-relative winds !
- ! !
- ! ASSUMPTIONS: 1) MEMORY ORDER IS XY. !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine metmap_xform_nmm(u, u_mask, v, v_mask, &
- vs1, vs2, ve1, ve2, &
- xlat_v, xlon_v, idir)
- implicit none
- ! Arguments
- integer, intent(in) :: vs1, vs2, ve1, ve2, idir
- real, pointer, dimension(:,:) :: u, v, xlat_v, xlon_v
- type (bitarray), intent(in) :: u_mask, v_mask
- ! Local variables
- integer :: i, j
- real :: u_map, v_map, diff, alpha
- real :: phi0, lmbd0, big_denominator, relm, rlat_v,rlon_v, clontemp
- real :: sin_phi0, cos_phi0, cos_alpha, sin_alpha
- real, pointer, dimension(:,:) :: u_new, v_new
- ! If the proj_info structure has not been initialized, we don't have
- ! information about the projection and standard longitude.
- if (proj_stack(current_nest_number)%init) then
- if (proj_stack(current_nest_number)%code == PROJ_ROTLL) then
- call mprintf((idir == 1),LOGFILE,'Rotating map winds to earth winds.')
- call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
-
- ! Create arrays to hold rotated winds
- allocate(u_new(vs1:ve1, vs2:ve2))
- allocate(v_new(vs1:ve1, vs2:ve2))
- phi0 = proj_stack(current_nest_number)%lat1 * rad_per_deg
- clontemp= proj_stack(current_nest_number)%lon1
- if (clontemp .lt. 0.) then
- lmbd0 = (clontemp + 360.) * rad_per_deg
- else
- lmbd0 = (clontemp) * rad_per_deg
- endif
- sin_phi0 = sin(phi0)
- cos_phi0 = cos(phi0)
- do j=vs2,ve2
- do i=vs1,ve1
- ! Calculate the sine and cosine of rotation angle
- rlat_v = xlat_v(i,j) * rad_per_deg
- rlon_v = xlon_v(i,j) * rad_per_deg
- relm = rlon_v - lmbd0
- big_denominator = cos(asin( &
- cos_phi0 * sin(rlat_v) - &
- sin_phi0 * cos(rlat_v) * cos(relm) &
- ) )
- sin_alpha = sin_phi0 * sin(relm) / &
- big_denominator
- cos_alpha = (cos_phi0 * cos(rlat_v) + &
- sin_phi0 * sin(rlat_v) * cos(relm)) / &
- big_denominator
-
- ! Rotate U field
- if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
- u_map = u(i,j)
- if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
- v_map = v(i,j)
- else
- v_map = 0.
- end if
-
- u_new(i,j) = cos_alpha*u_map + idir*sin_alpha*v_map
- else
- u_new(i,j) = u(i,j)
- end if
-
- ! Rotate V field
- if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
- v_map = v(i,j)
- if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
- u_map = u(i,j)
- else
- u_map = 0.
- end if
-
- v_new(i,j) = -idir*sin_alpha*u_map + cos_alpha*v_map
- else
- v_new(i,j) = v(i,j)
- end if
-
- end do
- end do
- ! Copy rotated winds back into argument arrays
- u = u_new
- v = v_new
-
- deallocate(u_new)
- deallocate(v_new)
-
- ! Only rotate winds for Lambert conformal, polar stereographic, or Cassini
- else if ((proj_stack(current_nest_number)%code == PROJ_LC) .or. &
- (proj_stack(current_nest_number)%code == PROJ_PS) .or. &
- (proj_stack(current_nest_number)%code == PROJ_CASSINI)) then
- call mprintf((idir == 1),LOGFILE,'Rotating map winds to earth winds.')
- call mprintf((idir == -1),LOGFILE,'Rotating earth winds to grid winds')
- ! Create arrays to hold rotated winds
- allocate(u_new(vs1:ve1, vs2:ve2))
- allocate(v_new(vs1:ve1, vs2:ve2))
- do j=vs2,ve2
- do i=vs1,ve1
- diff = idir * (xlon_v(i,j) - proj_stack(current_nest_number)%stdlon)
- if (diff > 180.) then
- diff = diff - 360.
- else if (diff < -180.) then
- diff = diff + 360.
- end if
- ! Calculate the rotation angle, alpha, in radians
- if (proj_stack(current_nest_number)%code == PROJ_LC) then
- alpha = diff * proj_stack(current_nest_number)%cone * &
- rad_per_deg * proj_stack(current_nest_number)%hemi
- else
- alpha = diff * rad_per_deg * proj_stack(current_nest_number)%hemi
- end if
- ! Rotate U field
- if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
- u_map = u(i,j)
- if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
- v_map = v(i,j)
- else
- v_map = 0.
- end if
-
- u_new(i,j) = cos(alpha)*u_map + idir*sin(alpha)*v_map
- else
- u_new(i,j) = u(i,j)
- end if
-
- ! Rotate V field
- if (bitarray_test(v_mask, i-vs1+1, j-vs2+1)) then
- v_map = v(i,j)
- if (bitarray_test(u_mask, i-vs1+1, j-vs2+1)) then
- u_map = u(i,j)
- else
- u_map = 0.
- end if
-
- v_new(i,j) = -idir*sin(alpha)*u_map + cos(alpha)*v_map
- else
- v_new(i,j) = v(i,j)
- end if
- end do
- end do
- ! Copy rotated winds back into argument arrays
- u = u_new
- v = v_new
-
- deallocate(u_new)
- deallocate(v_new)
- end if
- else
- call mprintf(.true.,ERROR,'In metmap_xform_nmm(), uninitialized proj_info structure.')
- end if
- end subroutine metmap_xform_nmm
- end module rotate_winds_module