diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index bc3062084..5703bcc76 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -1037,7 +1037,9 @@ module fv_arrays_mod ! to .true. in the next city release if desired !integer, pointer :: test_case - !real, pointer :: alpha + !real, pointer :: alpha + + logical :: compute_coords_locally = .false. end type fv_flags_type diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 74b4eb651..84d54c0c2 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -321,6 +321,8 @@ module fv_control_mod deglat_start, deglat_stop real(kind=R_GRID), pointer :: deglat + logical, pointer :: compute_coords_locally + logical, pointer :: nested, twowaynest integer, pointer :: parent_tile, refinement, nestbctype, nestupdate, nsponge, ioffset, joffset real, pointer :: s_weight, update_blend @@ -676,7 +678,8 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split) nested, twowaynest, parent_grid_num, parent_tile, nudge_qv, & refinement, nestbctype, nestupdate, nsponge, s_weight, & ioffset, joffset, check_negative, nudge_ic, halo_update_type, gfs_phil, agrid_vel_rst, & - do_uni_zfull, adj_mass_vmr + do_uni_zfull, adj_mass_vmr, & + compute_coords_locally namelist /test_case_nml/test_case, bubble_do, alpha, nsolitons, soliton_Umax, soliton_size @@ -1329,6 +1332,7 @@ subroutine setup_pointers(Atm) layout => Atm%layout io_layout => Atm%io_layout + compute_coords_locally => Atm%flagstruct%compute_coords_locally end subroutine setup_pointers diff --git a/model/fv_grid_utils.F90 b/model/fv_grid_utils.F90 index 0100eefc5..be853b98d 100644 --- a/model/fv_grid_utils.F90 +++ b/model/fv_grid_utils.F90 @@ -111,7 +111,7 @@ module fv_grid_utils_mod public f_p public ptop_min, big_number !CLEANUP: OK to keep since they are constants? public cos_angle - public latlon2xyz, gnomonic_grids, & + public latlon2xyz, gnomonic_grids, gnomonic_grids_local, & global_mx, unit_vect_latlon, & cubed_to_latlon, c2l_ord2, g_sum, global_qsum, great_circle_dist, & v_prod, get_unit_vect2, project_sphere_v @@ -3373,4 +3373,168 @@ subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng) end subroutine init_mq #endif + subroutine gnomonic_grids_local(grid_type, im, start, lon, lat) + integer, intent(in) :: grid_type + integer, intent(in) :: im + integer, intent(in) :: start(:) + real(R_GRID), intent(out):: lon(:,:) + real(R_GRID), intent(out):: lat(:,:) + + integer i, j + + select case (grid_type) + case (0) + call get_gnomonic_ed_coords(im, start, lon, lat) + case (1) + call get_gnomonic_dist_coords(im, start, lon, lat) + case (2) + call get_gnomonic_angl_coords(im, start, lon, lat) + end select + + if (grid_type < 3) then + lon = lon - pi + end if + + end subroutine gnomonic_grids_local + + subroutine get_gnomonic_ed_coords(im, start, lambda, theta) + integer, intent(in) :: im + integer, intent(in) :: start(:) + real(R_GRID), intent(out) :: lambda(start(1):, start(2):) + real(R_GRID), intent(out) :: theta(start(1):, start(2):) + + + ! Local: + real(R_GRID) :: pp(3,lbound(lambda,1):ubound(lambda,1), lbound(lambda,2):ubound(lambda,2)) + real(R_GRID) :: rsq3, alpha, beta, dely + real(R_GRID) :: b(1:im+1) + integer i, j + + rsq3 = 1./sqrt(3.) + alpha = asin( rsq3 ) + dely = alpha / im + + ! Compute distances along cube edge (segment) at equispaced + ! angles. The formula is same for i and for j, but for local + ! region these may or may not overlap. Worst case we + ! do this twice for a "diagonal subdomain" + + do i = lbound(lambda,1), ubound(lambda,1) + beta = -dely * (im + 2 - 2*i) + b(i) = tan(beta)*cos(alpha) + end do + + do j = lbound(lambda,2), ubound(lambda,2) + beta = -dely * (im + 2 - 2*j) + b(j) = tan(beta)*cos(alpha) + end do + + ! Use the edge position to construct an array of + ! cartesian points in the local domain. + do j = lbound(lambda,2), ubound(lambda,2) + do i = lbound(lambda,1), ubound(lambda,1) + pp(:,i,j) = [-rsq3, -b(i), b(j)] + end do + end do + + call cart_to_latlon_new(pp, lambda, theta) + + end subroutine get_gnomonic_ed_coords + + subroutine get_gnomonic_angl_coords(im, start, lambda, theta) + integer, intent(in) :: im + integer, intent(in) :: start(:) + real(R_GRID), intent(out) :: lambda(start(1):, start(2):) + real(R_GRID), intent(out) :: theta(start(1):, start(2):) + + ! Local + real(R_GRID) p(3,im+1,im+1) + real(R_GRID) rsq3 + real(R_GRID) dp + integer i,j + + dp = (pi/4) /real(im,R_GRID) + rsq3 = 1./sqrt(3.) + do j = lbound(lambda,2), ubound(lambda,2) + do i = lbound(lambda,1), ubound(lambda,1) + p(1,i,j) =-rsq3 ! constant + p(2,i,j) =-rsq3*tan(dp * (im + 2 - 2*i)) + p(2,i,j) =+rsq3*tan(dp * (im + 2 - 2*j)) + enddo + enddo + + call cart_to_latlon_new(p, lambda, theta) + + end subroutine get_gnomonic_angl_coords + + + ! It is important to use symmetric expressions for coordinates here. + ! This avoids the need for a subsequent forced symmetrization which in turn + ! thwarts fully local computation of coordinates. + subroutine get_gnomonic_dist_coords(im, start, lambda, theta) + integer, intent(in) :: im + integer, intent(in) :: start(:) + real(R_GRID), intent(out) :: lambda(start(1):, start(2):) + real(R_GRID), intent(out) :: theta(start(1):, start(2):) + + ! Local + real(R_GRID) rsq3, xf + real(R_GRID) dx + real(R_GRID) p(3,lbound(lambda,1):ubound(lambda,1),lbound(lambda,2):ubound(lambda,2)) + integer i, j + + ! Face-2 + + rsq3 = 1./sqrt(3.) + xf = -rsq3 + dx = rsq3/im + + do j = lbound(lambda,2), ubound(lambda,2) + do i = lbound(lambda,1), ubound(lambda,1) + p(1,i,j) = xf + p(2,i,j) = +dx * (im + 2 - 2*i) + p(3,i,j) = -dx * (im + 2 - 2*j) + enddo + enddo + call cart_to_latlon_new(p, lambda, theta) + + end subroutine get_gnomonic_dist_coords + + subroutine cart_to_latlon_new(q, xs, ys) + ! vector version of cart_to_latlon1 + real(R_GRID), intent(inout) :: q(:,:,:) + real(R_GRID), intent(inout) :: xs(:,:), ys(:,:) + + ! local + real(R_GRID), parameter:: esl=1.e-10 + real(R_GRID) :: p(3) + real(R_GRID) :: dist, lat, lon + integer i, j, k + + + do j = 1, size(q,3) + do i = 1, size(q,2) + p = q(:,i,j) + + dist = sqrt(p(1)**2 + p(2)**2 + p(3)**2) + p = p/dist + + if ( (abs(p(1))+abs(p(2))) < esl ) then + lon = 0. + else + lon = atan2( p(2), p(1) ) ! range [-pi,pi] + endif + + if ( lon < 0.) lon = 2.*pi + lon + lat = asin(p(3)) + + xs(i,j) = lon + ys(i,j) = lat + ! q Normalized: + q(:,i,j) = p + enddo + end do + + end subroutine cart_to_latlon_new + end module fv_grid_utils_mod diff --git a/tools/fv_grid_tools.F90 b/tools/fv_grid_tools.F90 index 4a869d521..40550921d 100644 --- a/tools/fv_grid_tools.F90 +++ b/tools/fv_grid_tools.F90 @@ -117,7 +117,7 @@ module fv_grid_tools_mod use constants_mod, only: grav, omega, pi=>pi_8, cnst_radius=>radius use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_grid_bounds_type, R_GRID - use fv_grid_utils_mod, only: gnomonic_grids, great_circle_dist, & + use fv_grid_utils_mod, only: gnomonic_grids, gnomonic_grids_local, great_circle_dist, & mid_pt_sphere, spherical_angle, & cell_center2, get_area, inner_prod, fill_ghost, & direct_transform, dist2side_latlon, & @@ -172,7 +172,7 @@ module fv_grid_tools_mod END INTERFACE public :: todeg, missing, init_grid, spherical_to_cartesian - public :: mirror_grid, get_unit_vector + public :: mirror_grid, mirror_grid_local, get_unit_vector contains @@ -565,6 +565,7 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, integer, pointer, dimension(:,:,:) :: iinta, jinta, iintb, jintb real(kind=R_GRID), pointer, dimension(:,:,:,:) :: grid_global + real(kind=R_GRID), allocatable, dimension(:,:,:) :: grid_local integer, pointer :: npx_g, npy_g, ntiles_g, tile logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner @@ -574,6 +575,11 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, integer :: is, ie, js, je integer :: isd, ied, jsd, jed + + logical :: local_algorithm + + local_algorithm = atm%flagstruct%compute_coords_locally + is = Atm%bd%is ie = Atm%bd%ie js = Atm%bd%js @@ -610,9 +616,11 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, e2 => Atm%gridstruct%e2 if (Atm%neststruct%nested .or. ANY(Atm%neststruct%child_grids)) then - grid_global => Atm%grid_global - else if( trim(grid_file) .NE. 'INPUT/grid_spec.nc') then - allocate(grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)) + grid_global => Atm%grid_global + else if( trim(grid_file) .ne. 'INPUT/grid_spec.nc') then + if (.not. local_algorithm) then + allocate(grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)) + end if endif iinta => Atm%gridstruct%iinta @@ -675,79 +683,104 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, call read_grid(Atm, grid_file, ndims, nregions, ng) else - if (Atm%flagstruct%grid_type>=0) call gnomonic_grids(Atm%flagstruct%grid_type, npx-1, xs, ys) + if (.not. local_algorithm) then + if (Atm%flagstruct%grid_type>=0) call gnomonic_grids(Atm%flagstruct%grid_type, npx-1, xs, ys) - if (is_master()) then + if (is_master()) then - if (Atm%flagstruct%grid_type>=0) then - do j=1,npy - do i=1,npx - grid_global(i,j,1,1) = xs(i,j) - grid_global(i,j,2,1) = ys(i,j) - enddo - enddo -! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] - call mirror_grid(grid_global, ng, npx, npy, 2, 6) - do n=1,nregions - do j=1,npy - do i=1,npx -!--------------------------------- -! Shift the corner away from Japan -!--------------------------------- -!--------------------- This will result in the corner close to east coast of China ------------------ - if ( .not.Atm%flagstruct%do_schmidt .and. (Atm%flagstruct%shift_fac)>1.E-4 ) & - grid_global(i,j,1,n) = grid_global(i,j,1,n) - pi/Atm%flagstruct%shift_fac -!---------------------------------------------------------------------------------------------------- - if ( grid_global(i,j,1,n) < 0. ) & - grid_global(i,j,1,n) = grid_global(i,j,1,n) + 2.*pi - if (ABS(grid_global(i,j,1,1)) < 1.d-10) grid_global(i,j,1,1) = 0.0 - if (ABS(grid_global(i,j,2,1)) < 1.d-10) grid_global(i,j,2,1) = 0.0 + if (Atm%flagstruct%grid_type>=0) then + do j=1,npy + do i=1,npx + grid_global(i,j,1,1) = xs(i,j) + grid_global(i,j,2,1) = ys(i,j) enddo enddo + ! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] + call mirror_grid(grid_global, ng, npx, npy, 2, 6) + do n=1,nregions + do j=1,npy + do i=1,npx + !--------------------------------- + ! Shift the corner away from Japan + !--------------------------------- + !--------------------- This will result in the corner close to east coast of China ------------------ + if ( .not.Atm%flagstruct%do_schmidt .and. (Atm%flagstruct%shift_fac)>1.E-4 ) & + grid_global(i,j,1,n) = grid_global(i,j,1,n) - pi/Atm%flagstruct%shift_fac + !---------------------------------------------------------------------------------------------------- + if ( grid_global(i,j,1,n) < 0. ) & + grid_global(i,j,1,n) = grid_global(i,j,1,n) + 2.*pi + if (abs(grid_global(i,j,1,1)) < 1.d-10) grid_global(i,j,1,1) = 0.0 + if (abs(grid_global(i,j,2,1)) < 1.d-10) grid_global(i,j,2,1) = 0.0 + enddo + enddo + enddo + else + call mpp_error(FATAL, "fv_grid_tools: reading of ASCII grid files no longer supported") + endif + + grid_global( 1,1:npy,:,2)=grid_global(npx,1:npy,:,1) + grid_global( 1,1:npy,:,3)=grid_global(npx:1:-1,npy,:,1) + grid_global(1:npx,npy,:,5)=grid_global(1,npy:1:-1,:,1) + grid_global(1:npx,npy,:,6)=grid_global(1:npx,1,:,1) + + grid_global(1:npx, 1,:,3)=grid_global(1:npx,npy,:,2) + grid_global(1:npx, 1,:,4)=grid_global(npx,npy:1:-1,:,2) + grid_global(npx,1:npy,:,6)=grid_global(npx:1:-1,1,:,2) + + grid_global( 1,1:npy,:,4)=grid_global(npx,1:npy,:,3) + grid_global( 1,1:npy,:,5)=grid_global(npx:1:-1,npy,:,3) + + grid_global(npx,1:npy,:,3)=grid_global(1,1:npy,:,4) + grid_global(1:npx, 1,:,5)=grid_global(1:npx,npy,:,4) + grid_global(1:npx, 1,:,6)=grid_global(npx,npy:1:-1,:,4) + + grid_global( 1,1:npy,:,6)=grid_global(npx,1:npy,:,5) + + !------------------------ + ! Schmidt transformation: + !------------------------ + if ( Atm%flagstruct%do_schmidt ) then + do n=1,nregions + call direct_transform(Atm%flagstruct%stretch_fac, 1, npx, 1, npy, & + Atm%flagstruct%target_lon, Atm%flagstruct%target_lat, & + n, grid_global(1:npx,1:npy,1,n), grid_global(1:npx,1:npy,2,n)) + enddo + endif + endif ! is_master + + call mpp_broadcast(grid_global, size(grid_global), mpp_root_pe()) + !--- copy grid to compute domain + do n=1,ndims + do j=js,je+1 + do i=is,ie+1 + grid(i,j,n) = grid_global(i,j,n,tile) + enddo enddo - else - call mpp_error(FATAL, "fv_grid_tools: reading of ASCII grid files no longer supported") - endif - - grid_global( 1,1:npy,:,2)=grid_global(npx,1:npy,:,1) - grid_global( 1,1:npy,:,3)=grid_global(npx:1:-1,npy,:,1) - grid_global(1:npx,npy,:,5)=grid_global(1,npy:1:-1,:,1) - grid_global(1:npx,npy,:,6)=grid_global(1:npx,1,:,1) - - grid_global(1:npx, 1,:,3)=grid_global(1:npx,npy,:,2) - grid_global(1:npx, 1,:,4)=grid_global(npx,npy:1:-1,:,2) - grid_global(npx,1:npy,:,6)=grid_global(npx:1:-1,1,:,2) - - grid_global( 1,1:npy,:,4)=grid_global(npx,1:npy,:,3) - grid_global( 1,1:npy,:,5)=grid_global(npx:1:-1,npy,:,3) - - grid_global(npx,1:npy,:,3)=grid_global(1,1:npy,:,4) - grid_global(1:npx, 1,:,5)=grid_global(1:npx,npy,:,4) - grid_global(1:npx, 1,:,6)=grid_global(npx,npy:1:-1,:,4) - - grid_global( 1,1:npy,:,6)=grid_global(npx,1:npy,:,5) - -!------------------------ -! Schmidt transformation: -!------------------------ - if ( Atm%flagstruct%do_schmidt ) then - do n=1,nregions - call direct_transform(Atm%flagstruct%stretch_fac, 1, npx, 1, npy, & - Atm%flagstruct%target_lon, Atm%flagstruct%target_lat, & - n, grid_global(1:npx,1:npy,1,n), grid_global(1:npx,1:npy,2,n)) - enddo + enddo + else ! local algorithm + allocate(grid_local(is:ie+1,js:je+1,2)) + call gnomonic_grids_local(Atm%flagstruct%grid_type, npx-1, [is,js], grid_local(:,:,1), grid_local(:,:,2)) + call mirror_grid_local(grid_local, tile) + !--------------------------------- + ! Shift the corner away from Japan + !--------------------------------- + !--------------------- This will result in the corner close to east coast of China ------------------ + if ( .not.Atm%flagstruct%do_schmidt .and. (Atm%flagstruct%shift_fac)>1.E-4 ) & + grid_local(:,:,1) = grid_local(:,:,1) - pi/Atm%flagstruct%shift_fac + !---------------------------------------------------------------------------------------------------- + where (grid_local(:,:,1) < 0.) + grid_local(:,:,1) = grid_local(:,:,1) + 2*pi + end where + where (abs(grid_local(:,:,:)) < 1.d-10) grid_local = 0 + do n=1,nregions + call direct_transform(Atm%flagstruct%stretch_fac, is, ie+1, js, je+1, & + Atm%flagstruct%target_lon, Atm%flagstruct%target_lat, & + n, grid_local(:,:,1), grid_local(:,:,2)) + enddo + grid(is:ie+1,js:je+1,:) = grid_local(:,:,:) + deallocate(grid_local) endif - endif - call mpp_broadcast(grid_global, size(grid_global), mpp_root_pe()) -!--- copy grid to compute domain - do n=1,ndims - do j=js,je+1 - do i=is,ie+1 - grid(i,j,n) = grid_global(i,j,n,tile) - enddo - enddo - enddo - endif + end if ! ! SJL: For phys/exchange grid, etc ! @@ -1115,9 +1148,11 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, endif!if gridtype > 3 if (Atm%neststruct%nested .or. ANY(Atm%neststruct%child_grids)) then - nullify(grid_global) - else if( trim(grid_file) .NE. 'INPUT/grid_spec.nc') then - deallocate(grid_global) + nullify(grid_global) + else if( trim(grid_file) .ne. 'INPUT/grid_spec.nc') then + if (.not. local_algorithm) then + deallocate(grid_global) + end if endif nullify(agrid) @@ -2457,5 +2492,160 @@ subroutine normalize_vect(np, e) end subroutine normalize_vect + !----------- + ! New version of coordinate calculations + ! + ! These only require allocation of arrays that span local domain and + ! are far more accurate due to careful construction of symmetric + ! loops. + + subroutine mirror_grid_local(local_tile, tileno) + real(R_GRID) , intent(INOUT) :: local_tile(:,:,:) + integer, intent(IN) :: tileno + + integer :: i,j,n,n1,n2,nreg, npx, npy + real(R_GRID) :: x1,y1,z1, x2,y2,z2, ang, sa, ca + + if (tileno == 1) then + ! no op + else + do j = 1, size(local_tile,2) + do i = 1, size(local_tile,1) + + x1 = local_tile(i,j,1) + y1 = local_tile(i,j,2) + z1 = radius + + select case (tileno) + case (2) + ang = -90. + sa = -1 + ca = 0 + call rot_3d_new( 3, x1, y1, z1, sa, ca, x2, y2, z2, 1) ! rotate about the z-axis + + case (3) + ang = -90. + sa = -1 + ca = 0 + call rot_3d_new( 3, x1, y1, z1, sa, ca, x2, y2, z2, 1) ! rotate about the z-axis + ang = 90. + sa = +1 + ca = 0 + call rot_3d_new( 1, x2, y2, z2, sa, ca, x1, y1, z1, 1) ! rotate about the x-axis + x2=x1 + y2=y1 + z2=z1 + + case (4) + ang = -180. + sa = 0 + ca = -1 + call rot_3d_new( 3, x1, y1, z1, sa, ca, x2, y2, z2, 1) ! rotate about the z-axis + ang = 90. + sa = 1 + ca = 0 + call rot_3d_new( 1, x2, y2, z2, sa, ca, x1, y1, z1, 1) ! rotate about the x-axis + x2=x1 + y2=y1 + z2=z1 + + case (5) + ang = 90. + sa = 1 + ca = 0 + call rot_3d_new( 3, x1, y1, z1, sa, ca, x2, y2, z2, 1) ! rotate about the z-axis + ang = 90. + sa = 1 + ca = 0 + call rot_3d_new( 2, x2, y2, z2, sa, ca, x1, y1, z1, 1) ! rotate about the y-axis + x2=x1 + y2=y1 + z2=z1 + + case (6) + ang = 90. + sa = 1 + ca = 0 + call rot_3d_new( 2, x1, y1, z1, sa, ca, x2, y2, z2, 1) ! rotate about the y-axis + ang = 0. + sa = 0 + ca = 1 + call rot_3d_new( 3, x2, y2, z2, sa, ca, x1, y1, z1, 1) ! rotate about the z-axis + x2=x1 + y2=y1 + z2=z1 + + end select + + local_tile(i,j,1) = x2 + local_tile(i,j,2) = y2 + + enddo + enddo + endif + end subroutine mirror_grid_local + + + !------------------------------------------------------------------------------- + ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! + ! + ! rot_3d_new :: rotate points on a sphere in xyz coords given sin + ! and cos of angle. Only works for {0,-1,+1} - 90 + ! degree rotations. This approach guarantees + ! symmetry by avoiding roundoff associated with + ! inexact values of pi (and thus cos(pi/2)). + + ! + subroutine rot_3d_new(axis, x1in, y1in, z1in, sa, ca, x2out, y2out, z2out, convert) + + integer, intent(IN) :: axis ! axis of rotation 1=x, 2=y, 3=z + real(R_GRID) , intent(IN) :: x1in, y1in, z1in + real(R_GRID) , intent(IN) :: sa, ca ! sin and cos of angle to rotate in radians + real(R_GRID) , intent(OUT) :: x2out, y2out, z2out + integer, intent(IN), optional :: convert ! if present convert input point + ! from spherical to cartesian, rotate, + ! and convert back + + real(R_GRID) :: x1,y1,z1, x2,y2,z2 + + if ( present(convert) ) then + call spherical_to_cartesian(x1in, y1in, z1in, x1, y1, z1) + else + x1=x1in + y1=y1in + z1=z1in + endif + + + SELECT CASE(axis) + + CASE(1) + x2 = x1 + y2 = ca*y1 + sa*z1 + z2 = -sa*y1 + ca*z1 + CASE(2) + x2 = ca*x1 - sa*z1 + y2 = y1 + z2 = sa*x1 + ca*z1 + CASE(3) + x2 = ca*x1 + sa*y1 + y2 = -sa*x1 + ca*y1 + z2 = z1 + CASE DEFAULT + write(*,*) "Invalid axis: must be 1 for X, 2 for Y, 3 for Z." + + END SELECT + + if ( present(convert) ) then + call cartesian_to_spherical(x2, y2, z2, x2out, y2out, z2out) + else + x2out=x2 + y2out=y2 + z2out=z2 + endif + + end subroutine rot_3d_new + + end module fv_grid_tools_mod