program winds ! Unit test for routine convert_winds, which ! converts u/v component winds to x/y/z component ! winds on a sphere. This test converts unit ! vectors at several latitude and longitudes. ! ! Author: George Gayno NCEP/EMC use esmf use model_grid, only : i_input, j_input, & input_grid, & latitude_input_grid, & longitude_input_grid use atm_input_data, only : lev_input, convert_winds_to_xyz, & xwind_input_grid, & ywind_input_grid, & zwind_input_grid, & u_input_grid, & v_input_grid implicit none integer, parameter :: IPTS=4 integer, parameter :: JPTS=3 real, parameter :: EPSILON=0.0001 integer :: clb(3), cub(3) integer :: ierr, localpet, npets, rc integer :: i, j, k real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) real(esmf_kind_r8), allocatable :: u_wind(:,:,:) real(esmf_kind_r8), allocatable :: v_wind(:,:,:) real(esmf_kind_r8), pointer :: xwindptr(:,:,:) real(esmf_kind_r8), pointer :: ywindptr(:,:,:) real(esmf_kind_r8), pointer :: zwindptr(:,:,:) real :: expected_x_component(IPTS,JPTS) real :: expected_y_component(IPTS,JPTS) real :: expected_z_component(IPTS,JPTS) type(esmf_vm) :: vm type(esmf_polekind_flag) :: polekindflag(2) ! These are the expected x/y/z components return from ! convert_winds for our test points. data expected_x_component/1.0, 0.0, -1.0, 0.0, & 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 1.0, 0.0 / data expected_y_component/0.0, 1.0, 0.0, -1.0, & 0.0, 0.0, 0.0, 0.0, & 1.0, -1.0, 0.0, -0.707106/ data expected_z_component/0.0, 0.0, 0.0, 0.0, & 1.0, 1.0, 1.0, 1.0, & 0.0, 0.0, 0.0, 0.707106/ print*,"Starting test of convert_winds." call mpi_init(ierr) call ESMF_Initialize(rc=ierr) call ESMF_VMGetGlobal(vm, rc=ierr) call ESMF_VMGet(vm, localPet=localpet, petCount=npets, rc=ierr) i_input = IPTS j_input = JPTS lev_input = 1 ! number of vertical levels polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE input_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & maxIndex=(/i_input,j_input/), & polekindflag=polekindflag, & periodicDim=1, & poleDim=2, & coordSys=ESMF_COORDSYS_SPH_DEG, & regDecomp=(/1,npets/), & indexflag=ESMF_INDEX_GLOBAL, rc=rc) latitude_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & name="input_grid_latitude", & rc=rc) longitude_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & name="input_grid_longitude", & rc=rc) xwind_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & ungriddedLBound=(/1/), & ungriddedUBound=(/lev_input/), rc=rc) ywind_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & ungriddedLBound=(/1/), & ungriddedUBound=(/lev_input/), rc=rc) zwind_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & ungriddedLBound=(/1/), & ungriddedUBound=(/lev_input/), rc=rc) u_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & ungriddedLBound=(/1/), & ungriddedUBound=(/lev_input/), rc=rc) v_input_grid = ESMF_FieldCreate(input_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & ungriddedLBound=(/1/), & ungriddedUBound=(/lev_input/), rc=rc) allocate(latitude(i_input,j_input)) allocate(longitude(i_input,j_input)) allocate(u_wind(i_input,j_input,lev_input)) allocate(v_wind(i_input,j_input,lev_input)) ! Test point 1 - A west wind at the Equator/Greenwich. latitude(1,1) = 0. longitude(1,1) = 0. u_wind(1,1,1) = 1.0 v_wind(1,1,1) = 0.0 ! Test point 2 - A west wind at the Equator/90 degrees East longitude. latitude(2,1) = 0. longitude(2,1) = 90. u_wind(2,1,1) = 1.0 v_wind(2,1,1) = 0.0 ! Test point 3 - A west wind at the Equator/Dateline. latitude(3,1) = 0. longitude(3,1) = 180. u_wind(3,1,1) = 1.0 v_wind(3,1,1) = 0.0 ! Test point 4 - A west wind at the Equator/90 degrees West longitude. latitude(4,1) = 0. longitude(4,1) = 270. u_wind(4,1,1) = 1.0 v_wind(4,1,1) = 0.0 ! Test point 5 - A south wind at the Equator/Greenwich. latitude(1,2) = 0. longitude(1,2) = 0. u_wind(1,2,1) = 0.0 v_wind(1,2,1) = 1.0 ! Test point 6 - A south wind at the Equator/90 degrees East longitude. latitude(2,2) = 0. longitude(2,2) = 90. u_wind(2,2,1) = 0.0 v_wind(2,2,1) = 1.0 ! Test point 7 - A south wind at the Equator/Dateline. latitude(3,2) = 0. longitude(3,2) = 180. u_wind(3,2,1) = 0.0 v_wind(3,2,1) = 1.0 ! Test point 8 - A south wind at the Equator/90 degrees West longitude. latitude(4,2) = 0. longitude(4,2) = 270. u_wind(4,2,1) = 0.0 v_wind(4,2,1) = 1.0 ! Test point 9 - A south wind at the North Pole/Greenwich. latitude(1,3) = 90. longitude(1,3) = 0. u_wind(1,3,1) = 0.0 v_wind(1,3,1) = 1.0 ! Test point 10 - A south wind at the South Pole/Greenwich. latitude(2,3) = -(90.0) longitude(2,3) = 0. u_wind(2,3,1) = 0.0 v_wind(2,3,1) = 1.0 ! Test point 11 - A west wind at the 45 degrees North latitude/Greenwich. latitude(3,3) = 45. longitude(3,3) = 0. u_wind(3,3,1) = 1.0 v_wind(3,3,1) = 0.0 ! Test point 12 - A south wind at 45 degrees North latitude/Dateline. latitude(4,3) = 45. longitude(4,3) = 180. u_wind(4,3,1) = 0.0 v_wind(4,3,1) = 1.0 call ESMF_FieldScatter(u_input_grid, u_wind, rootpet=0, rc=rc) call ESMF_FieldScatter(v_input_grid, v_wind, rootpet=0, rc=rc) call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc) call ESMF_FieldScatter(latitude_input_grid, latitude, rootpet=0, rc=rc) ! Call the routine to unit test. call convert_winds_to_xyz call ESMF_FieldGet(xwind_input_grid, & computationalLBound=clb, & computationalUBound=cub, & farrayPtr=xwindptr, rc=rc) call ESMF_FieldGet(ywind_input_grid, & computationalLBound=clb, & computationalUBound=cub, & farrayPtr=ywindptr, rc=rc) call ESMF_FieldGet(zwind_input_grid, & computationalLBound=clb, & computationalUBound=cub, & farrayPtr=zwindptr, rc=rc) print*,"Check results." do j = clb(2), cub(2) do i = clb(1), cub(1) do k = clb(3), cub(3) if (abs(xwindptr(i,j,k) - expected_x_component(i,j)) > EPSILON) stop 2 if (abs(ywindptr(i,j,k) - expected_y_component(i,j)) > EPSILON) stop 3 if (abs(zwindptr(i,j,k) - expected_z_component(i,j)) > EPSILON) stop 4 enddo enddo enddo print*,"OK" deallocate(u_wind, v_wind, latitude, longitude) call ESMF_finalize(endflag=ESMF_END_KEEPMPI) call mpi_finalize(rc) print*,"SUCCESS!" end program winds