! +-======-+ ! Copyright (c) 2003-2007 United States Government as represented by ! the Admistrator of the National Aeronautics and Space Administration. ! All Rights Reserved. ! ! THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, ! REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN ! COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS ! REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). ! THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN ! INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR ! REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, ! DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED ! HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE ! RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. ! ! Government Agency: National Aeronautics and Space Administration ! Government Agency Original Software Designation: GSC-15354-1 ! Government Agency Original Software Title: GEOS-5 GCM Modeling Software ! User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov ! Government Agency Point of Contact for Original Software: ! Dale Hithon, SRA Assistant, (301) 286-2691 ! ! +-======-+ ! $Id: CubeToLatLon.F90,v 1.1.2.38 2011-03-18 16:04:03 msuarez Exp $ #define SUCCESS 0 #define VERIFY_(I) if((I)/=0) then; if(present(rc)) rc=I; return; end if #define RETURN_(I) rc=I; return #define ASSERT_(I) if(.not.(I)) then; if(present(rc)) rc=1; return; end if Module CubeLatLonTransformMod implicit none private public T_CubeLatLonTransform public CubeLatLonIsCreated public CubeLatLonCreate public CubeLatLonDestroy public CubeToLatLon public LatLonToCube public CartesianToSpherical public SphericalToCartesian #define R8 8 type T_CubeLatLonTransform private real(R8),pointer :: weight(:,:,:),l2c(:,:,:) integer, pointer :: index (:,:,:) integer, pointer :: id1(:,:), id2(:,:), jdc(:,:) logical :: Created=.false. character(len=120) :: name integer :: npx, npy, nlon, nlat real(R8), pointer :: ee1(:,:,:) real(R8), pointer :: ee2(:,:,:) real(R8), pointer :: ff1(:,:,:) real(R8), pointer :: ff2(:,:,:) real(R8), pointer :: elon(:,:,:) real(R8), pointer :: elat(:,:,:) end type T_CubeLatLonTransform interface CubeToLatLon module procedure CubeToLatLonr8 module procedure CubeToLatLonr4 end interface interface LatLonToCube module procedure LatLonToCuber8 module procedure LatLonToCuber4 end interface integer, parameter :: ntiles=6 integer, parameter :: ndims=2 integer, parameter :: r8=R8 integer, parameter :: maxstring=120 real(R8), parameter :: PI=3.14159265358979323846 ! This EXTERNAL subroutine is in the fv directory ! and has real*8 interfaces interface subroutine GetWeights(npx, npy, nlat, nlon, & index, weight, id1, id2, jdc, l2c, & ee1, ee2, ff1, ff2) integer, intent(in ) :: npx, npy integer, intent(in ) :: nlon, nlat integer, intent(inout) :: index(3,nlon,nlat) real(R8), intent(inout) :: weight(4,nlon,nlat) integer, intent(inout) :: id1(npx,npy) integer, intent(inout) :: id2(npx,npy) integer, intent(inout) :: jdc(npx,npy) real(R8), intent(inout) :: l2c(4,npx,npy) real(R8), pointer :: ee1(:,:,:) real(R8), pointer :: ee2(:,:,:) real(R8), pointer :: ff1(:,:,:) real(R8), pointer :: ff2(:,:,:) end subroutine GetWeights end interface contains subroutine CubeLatLonDestroy( Trans, rc) type(T_CubeLatLonTransform), intent(inout) :: Trans integer, optional, intent( out) :: rc if(associated(Trans%weight)) deallocate(Trans%weight) if(associated(Trans%l2c )) deallocate(Trans%l2c ) if(associated(Trans%id1 )) deallocate(Trans%id1 ) if(associated(Trans%id2 )) deallocate(Trans%id2 ) if(associated(Trans%jdc )) deallocate(Trans%jdc ) if(associated(Trans%elon )) deallocate(Trans%elon ) if(associated(Trans%elat )) deallocate(Trans%elat ) if(associated(Trans%ee1 )) deallocate(Trans%ee1 ) if(associated(Trans%ee2 )) deallocate(Trans%ee2 ) if(associated(Trans%ff1 )) deallocate(Trans%ff1 ) if(associated(Trans%ff2 )) deallocate(Trans%ff2 ) Trans%Created = .false. RETURN_(SUCCESS) end subroutine CubeLatLonDestroy logical function CubeLatLonIsCreated(Trans) type(T_CubeLatLonTransform), intent(in ) :: Trans CubeLatLonIsCreated = Trans%Created end function CubeLatLonIsCreated function CubeLatLonCreate( npx, npy, nlon, nlat, lons, lats, rc ) result(Trans) integer, intent(in ) :: npx, npy integer, intent(in ) :: nlon, nlat real(R8), intent(in ) :: lons(:), lats(:) type(T_CubeLatLonTransform) :: Trans integer, optional, intent(out) :: rc ! Locals !------- integer :: npts, status, n, l integer :: i, j real(R8), allocatable :: slon(:), slat(:) real(R8), allocatable :: clon(:), clat(:) ! Real*8 are needed to make fv calls. !----------------------------------- ! Begin !------ ASSERT_(.not.Trans%Created) npts = npx + 1 write(Trans%name,'(i5.5,"x",i5.5,"_c2l_",i5.5,"x",i5.5)') npx,npy,nlon,nlat ! write(*,'(i5.5,"x",i5.5,"_c2l_",i5.5,"x",i5.5)') npx,npy,nlon,nlat Trans%npx = npx Trans%npy = npy Trans%nlon = nlon Trans%nlat = nlat ! allocate storage for weights and indeces for C2L !------------------------------------------------- if(associated(Trans%index )) deallocate(Trans%index ) if(associated(Trans%weight)) deallocate(Trans%weight) if(associated(Trans%l2c )) deallocate(Trans%l2c ) if(associated(Trans%id1 )) deallocate(Trans%id1 ) if(associated(Trans%id2 )) deallocate(Trans%id2 ) if(associated(Trans%jdc )) deallocate(Trans%jdc ) if(associated(Trans%elon )) deallocate(Trans%elon ) if(associated(Trans%elat )) deallocate(Trans%elat ) if(associated(Trans%ee1 )) deallocate(Trans%ee1 ) if(associated(Trans%ee2 )) deallocate(Trans%ee2 ) if(associated(Trans%ff1 )) deallocate(Trans%ff1 ) if(associated(Trans%ff2 )) deallocate(Trans%ff2 ) allocate(Trans%index(3,nlon,nlat),Trans%weight(4,nlon,nlat),stat=status) VERIFY_(STATUS) allocate( Trans%l2c(4,npx, npy), Trans%id1(npx, npy), & Trans%id2(npx, npy), Trans%jdc(npx, npy) ,stat=status) VERIFY_(STATUS) call GetWeights(npx, npy, nlat, nlon, Trans%index, Trans%weight, & Trans%id1, Trans%id2, Trans%jdc, Trans%l2c, & Trans%ee1, Trans%ee2, & Trans%ff1, Trans%ff2) allocate(Trans%elon(size(lons),size(lats),3)) allocate(Trans%elat(size(lons),size(lats),3)) allocate(slat(size(lats)),clat(size(lats))) allocate(slon(size(lons)),clon(size(lons))) do j=1,size(lats) SLAT(j) = SIN(lats(j)) CLAT(j) = COS(lats(j)) end do do I=1,size(lons) SLON(I) = sin(lons(i) - PI) CLON(I) = cos(lons(i) - PI) end DO do j=1,size(lats) do I=1,size(lons) Trans%elon(I,J,1) = -SLON(I) Trans%elon(I,J,2) = CLON(I) Trans%elon(I,J,3) = 0.0 Trans%elat(I,J,1) = -SLAT(J)*CLON(I) Trans%elat(I,J,2) = -SLAT(J)*SLON(I) Trans%elat(I,J,3) = CLAT(J) end do end do deallocate(slon,clon,slat,clat) Trans%Created=.true. RETURN_(SUCCESS) end function CubeLatLonCreate subroutine CubeToLatLonr8( Trans, data_cs, data_ll, transpose, misval, rc) type(T_CubeLatLonTransform), intent(in ) :: Trans real(R8), intent(inout) :: data_cs(:,:) real(R8), intent(inout) :: data_ll(:,:) logical, optional, intent(in ) :: transpose real, optional, intent(in ) :: misval integer, optional, intent(out) :: rc ! Locals !------- integer :: npx,npy,nlon,nlat integer :: nx,j1,j2,status,itile real(R8), allocatable :: var_cs(:,:,:) real(R8) :: misval_ ASSERT_(Trans%Created) if(present(misval)) then misval_ = misval else misval_ = 1.0 end if nx = Trans%npx !--------------------------------------------------------------------! ! perform interpolation ! !--------------------------------------------------------------------! allocate ( var_cs(0:nx+1,0:nx+1,ntiles),stat=status) VERIFY_(STATUS) ASSERT_(.not.transpose .or. misval_==1.0) var_cs=0.0 if(.not.transpose) then data_ll=0.0 do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx var_cs(1:nx,1:nx,itile)=data_cs(:,j1:j2) enddo end if call C2LInterp(var_cs, data_ll, Trans%index, Trans%weight,& misval_, transpose) if(transpose) then do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx data_cs(:,j1:j2) = var_cs(1:nx,1:nx,itile) enddo end if deallocate ( var_cs ,stat=status) VERIFY_(STATUS) RETURN_(SUCCESS) end subroutine CubeToLatLonr8 subroutine CubeToLatLonr4( Trans, data_cs, data_ll, transpose, misval, rc) type(T_CubeLatLonTransform), intent(in ) :: Trans real, intent(inout) :: data_cs(:,:) real, intent(inout) :: data_ll(:,:) logical, optional, intent(in ) :: transpose real, optional, intent(in ) :: misval integer, optional, intent(out) :: rc ! Locals !------- integer :: npx,npy,nlon,nlat integer :: nx,j1,j2,status,itile real(R8), allocatable :: var_cs(:,:,:), data_ll8(:,:) real(R8) :: misval_ ASSERT_(Trans%Created) if(present(misval)) then misval_ = misval else misval_ = 1.0 end if nx = Trans%npx !--------------------------------------------------------------------! ! perform interpolation ! !--------------------------------------------------------------------! allocate ( var_cs(0:nx+1,0:nx+1,ntiles),stat=status ) VERIFY_(STATUS) allocate ( data_ll8(size(data_ll,1),size(data_ll,2)),stat=status) VERIFY_(STATUS) ASSERT_(.not.transpose .or. misval_==1.0) var_cs=0.0 if(.not.transpose) then do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx var_cs(1:nx,1:nx,itile)=data_cs(:,j1:j2) enddo else data_ll8=data_ll end if call C2LInterp(var_cs, data_ll8, Trans%index, Trans%weight,& misval_, transpose) if(transpose) then do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx data_cs(:,j1:j2) = var_cs(1:nx,1:nx,itile) enddo else data_ll=data_ll8 end if deallocate ( var_cs, data_ll8 ) RETURN_(SUCCESS) end subroutine CubeToLatLonr4 subroutine LatLonToCuber8( Trans, data_ll, data_cs, transpose, misval, rc) type(T_CubeLatLonTransform), intent(in ) :: Trans real(R8), intent(inout) :: data_ll(:,:) real(R8), intent(inout) :: data_cs(:,:) logical, optional, intent(in ) :: transpose real, optional, intent(in ) :: misval integer, optional, intent(out) :: rc ! Locals !------- integer :: nx,j1,j2,status,itile real(R8), allocatable :: var_cs(:,:,:) real(R8) :: misval_ ASSERT_(Trans%Created) if(present(misval)) then misval_ = misval else misval_ = 1.0 end if nx = Trans%npx !--------------------------------------------------------------------! ! perform interpolation ! !--------------------------------------------------------------------! allocate ( var_cs(0:nx+1,0:nx+1,ntiles),stat=status) VERIFY_(STATUS) var_cs=0. if(transpose) then data_ll=0. do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx var_cs(1:nx,1:nx,itile) = data_cs(:,j1:j2) enddo end if call L2CInterp(data_ll, var_cs, Trans%id1, Trans%id2, Trans%jdc, & Trans%l2c, misval_, transpose) if(.not.transpose) then do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx data_cs(:,j1:j2) = var_cs(1:nx,1:nx,itile) enddo end if deallocate ( var_cs ,STAT=STATUS) VERIFY_(STATUS) RETURN_(SUCCESS) end subroutine LatLonToCuber8 subroutine LatLonToCuber4( Trans, data_ll, data_cs, transpose, misval, rc) type(T_CubeLatLonTransform), intent(in ) :: Trans real, intent(inout) :: data_ll(:,:) real, intent(inout) :: data_cs(:,:) logical, optional, intent(in ) :: transpose real, optional, intent(in ) :: misval integer, optional, intent(out) :: rc ! Locals !------- integer :: nx,j1,j2,status,itile real(R8) :: misval_ real(R8), allocatable :: data_cs8(:,:,:), data_ll8(:,:) ASSERT_(Trans%Created) if(present(misval)) then misval_ = misval else misval_ = 1.0 end if nx = Trans%npx allocate ( data_ll8(size(data_ll,1),size(data_ll,2)),stat=status) VERIFY_(STATUS) allocate ( data_cs8(0:nx+1,0:nx+1,ntiles),stat=status) VERIFY_(STATUS) !--------------------------------------------------------------------! ! perform interpolation ! !--------------------------------------------------------------------! data_cs8=0. if(.not.transpose) then data_ll8 = data_ll end if if(transpose) then data_ll8=0. do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx data_cs8(1:nx,1:nx,itile) = data_cs(:,j1:j2) enddo end if call L2CInterp(data_ll8, data_cs8, Trans%id1, Trans%id2, Trans%jdc, & Trans%l2c, misval_, transpose) if(.not.transpose) then do itile=1,ntiles j1 = nx*(itile-1) + 1 j2 = nx*(itile-1) + nx data_cs(:,j1:j2) = data_cs8(1:nx,1:nx,itile) enddo end if if(transpose) then data_ll = data_ll8 end if deallocate ( data_cs8 ) deallocate ( data_ll8 ) RETURN_(SUCCESS) end subroutine LatLonToCuber4 subroutine C2LInterp(cubsph, latlon, index, weight, misval, transpose) !------------------------------------------------------------------! ! do bilinear interpolation from cubed sphere to latlon grid ! ! using precalculated weights from get_weight ! !------------------------------------------------------------------! real(R8), dimension(0:,0:,:), intent(inout) :: cubsph real(R8), dimension(:,:), intent(inout) :: latlon real(R8), dimension(:,:,:), intent(in) :: weight integer, dimension(:,:,:), intent(in) :: index real(R8), intent(in) :: misval logical, intent(in) :: transpose !------------------------------------------------------------------! ! local variables ! !------------------------------------------------------------------! integer :: i, j, ic, jc, nx, ny, nlon, nlat, tile, ii real(R8) :: ww if(transpose .and. misval/=1) then print *, 'Trying to do C2L-transpose with missing value' return end if nx = size(cubsph,1)-2 ny = size(cubsph,2)-2 if(.not.transpose) then call GhostCube(cubsph) else cubsph = 0.0 endif ! not transpose nlon = size(latlon,1) nlat = size(latlon,2) if(mod(nlon,2)/=0) then print *, "NLON not even in cubetolatlon. Stopping." stop endif FACES: do tile = 1,ntiles ILOOP: do i=1,nlon ii = mod(i - 1 + nlon/2,nlon) + 1 JLOOP: do j=1,nlat HAVE_POINT: if (tile==index(3,i,j)) then ic=index(1,i,j) jc=index(2,i,j) ADJOINT: if(.not.transpose) then UNDEF: if(misval==1.D0) then latlon(ii,j) = weight(1,i,j)*cubsph(ic ,jc ,tile) & + weight(2,i,j)*cubsph(ic ,jc+1,tile) & + weight(3,i,j)*cubsph(ic+1,jc+1,tile) & + weight(4,i,j)*cubsph(ic+1,jc ,tile) else ww = 0.0 latlon(ii,j) = 0.0 if(cubsph(ic ,jc ,tile)/=misval) then latlon(ii,j) = latlon(ii,j) + weight(1,i,j)*cubsph(ic ,jc ,tile) ww = ww + weight(1,i,j) end if if(cubsph(ic ,jc+1,tile)/=misval) then latlon(ii,j) = latlon(ii,j) + weight(2,i,j)*cubsph(ic ,jc+1,tile) ww = ww + weight(2,i,j) end if if(cubsph(ic+1,jc+1,tile)/=misval) then latlon(ii,j) = latlon(ii,j) + weight(3,i,j)*cubsph(ic+1,jc+1,tile) ww = ww + weight(3,i,j) end if if(cubsph(ic+1,jc ,tile)/=misval) then latlon(ii,j) = latlon(ii,j) + weight(4,i,j)*cubsph(ic+1,jc ,tile) ww = ww + weight(4,i,j) end if if(ww==0.0) then latlon(ii,j) = misval else latlon(ii,j) = latlon(ii,j) / ww end if end if UNDEF else ! Transpose cubsph(ic ,jc ,tile)=cubsph(ic ,jc ,tile)+weight(1,i,j)*latlon(ii,j) cubsph(ic ,jc+1,tile)=cubsph(ic ,jc+1,tile)+weight(2,i,j)*latlon(ii,j) cubsph(ic+1,jc+1,tile)=cubsph(ic+1,jc+1,tile)+weight(3,i,j)*latlon(ii,j) cubsph(ic+1,jc ,tile)=cubsph(ic+1,jc ,tile)+weight(4,i,j)*latlon(ii,j) end if ADJOINT endif HAVE_POINT enddo JLOOP enddo ILOOP end do FACES if(transpose) then call GhostCubeT(cubsph) endif return end subroutine C2LInterp subroutine L2CInterp(latlon, cubsph, id1, id2, jdc, weight, misval, transpose) !------------------------------------------------------------------! ! do bilinear interpolation from cubed sphere to latlon grid ! ! using precalculated weights from get_weight ! !------------------------------------------------------------------! real(R8), dimension(0:,0:,:), intent(inout) :: cubsph real(R8), dimension(:,:), intent(inout) :: latlon real(R8), dimension(:,:,:), intent(in) :: weight integer, dimension(:,:), intent(in) :: id1, id2, jdc real(R8), intent(in) :: misval logical, intent(in) :: transpose !------------------------------------------------------------------! ! local variables ! !------------------------------------------------------------------! integer :: i, j, i1, i2, j1, nx, ny, nz, jx, k, nlon real(R8) :: ww nx = size(cubsph,1)-2 ny = size(cubsph,2)-2 nz = size(cubsph,3) nlon = size(latlon,1) if(transpose) latlon = 0.0 FACES: do k=1,nz FACE_Y: do jx=1,ny FACE_X: do i=1,nx j = (k-1)*ny + jx j1 = jdc(i,j) i1 = mod(id1(i,j)-1+nlon/2,nlon) + 1 i2 = mod(id2(i,j)-1+nlon/2,nlon) + 1 ADJOINT: if(.not.transpose) then UNDEF: if(misval==1.D0) then cubsph(i,jx,k) = weight(1,i,j)*latlon(i1,j1 ) & + weight(2,i,j)*latlon(i2,j1 ) & + weight(3,i,j)*latlon(i2,j1+1) & + weight(4,i,j)*latlon(i1,j1+1) else ww = 0.0 cubsph(i,jx,k) = 0.0 if(latlon(i1,j1 )/=misval) then cubsph(i,jx,k) = cubsph(i,jx,k) + weight(1,i,j)*latlon(i1,j1 ) ww = ww + weight(1,i,j) end if if(latlon(i2,j1 )/=misval) then cubsph(i,jx,k) = cubsph(i,jx,k) + weight(2,i,j)*latlon(i2,j1 ) ww = ww + weight(2,i,j) end if if(latlon(i1,j1+1)/=misval) then cubsph(i,jx,k) = cubsph(i,jx,k) + weight(3,i,j)*latlon(i1,j1+1) ww = ww + weight(3,i,j) end if if(latlon(i2,j1+1)/=misval) then cubsph(i,jx,k) = cubsph(i,jx,k) + weight(4,i,j)*latlon(i2,j1+1) ww = ww + weight(4,i,j) end if if(ww==0.0) then cubsph(i,jx,k) = misval else cubsph(i,jx,k) = cubsph(i,jx,k) / ww end if end if UNDEF else latlon(i1,j1 ) = latlon(i1,j1 ) + weight(1,i,j)*cubsph(i,jx,k) latlon(i2,j1 ) = latlon(i2,j1 ) + weight(2,i,j)*cubsph(i,jx,k) latlon(i2,j1+1) = latlon(i2,j1+1) + weight(3,i,j)*cubsph(i,jx,k) latlon(i1,j1+1) = latlon(i1,j1+1) + weight(4,i,j)*cubsph(i,jx,k) end if ADJOINT enddo FACE_X enddo FACE_Y end do FACES return end subroutine L2CInterp subroutine GhostCube(x) real(R8), intent(INOUT) :: x(0:,0:,:) integer :: nx, ny nx = size(x,1)-2 ny = nx x(1:nx,0, 1) = x(1:nx,ny, 6) x(1:nx,ny+1,1) = x(1,ny:1:-1, 3) x(0,1:ny, 1) = x(nx:1:-1,ny,5) x(nx+1,1:ny,1) = x(1,1:ny, 2) x(1:nx,0, 2) = x(nx,ny:1:-1,6) x(1:nx,ny+1,2) = x(1:nx,1, 3) x(0,1:ny, 2) = x(nx,1:ny, 1) x(nx+1,1:ny,2) = x(nx:1:-1,1, 4) x(1:nx,0, 3) = x(1:nx,ny, 2) x(1:nx,ny+1,3) = x(1,ny:1:-1, 5) x(0,1:ny, 3) = x(nx:1:-1,ny,1) x(nx+1,1:ny,3) = x(1,1:ny, 4) x(1:nx,0, 4) = x(nx,ny:1:-1,2) x(1:nx,ny+1,4) = x(1:nx,1, 5) x(0,1:ny, 4) = x(nx,1:ny, 3) x(nx+1,1:ny,4) = x(nx:1:-1,1, 6) x(1:nx,0, 5) = x(1:nx,ny, 4) x(1:nx,ny+1,5) = x(1,ny:1:-1, 1) x(0,1:ny, 5) = x(nx:1:-1,ny,3) x(nx+1,1:ny,5) = x(1,1:ny, 6) x(1:nx,0, 6) = x(nx,ny:1:-1,4) x(1:nx,ny+1,6) = x(1:nx,1, 1) x(0,1:ny, 6) = x(nx,1:ny, 5) x(nx+1,1:ny,6) = x(nx:1:-1,1, 2) ! Zero corners x(0, ny+1,:) = 0.0 x(0, 0,:) = 0.0 x(nx+1, 0,:) = 0.0 x(nx+1,ny+1,:) = 0.0 end subroutine GhostCube subroutine GhostCubeT(x) real(R8), intent(INOUT) :: x(0:,0:,:) integer :: nx, ny nx = size(x,1)-2 ny = nx x(1:nx,ny, 6) = x(1:nx,ny, 6) + x(1:nx,0, 1) x(1,ny:1:-1, 3) = x(1,ny:1:-1, 3) + x(1:nx,ny+1,1) x(nx:1:-1,ny,5) = x(nx:1:-1,ny,5) + x(0,1:ny, 1) x(1,1:ny, 2) = x(1,1:ny, 2) + x(nx+1,1:ny,1) x(nx,ny:1:-1,6) = x(nx,ny:1:-1,6) + x(1:nx,0, 2) x(1:nx,1, 3) = x(1:nx,1, 3) + x(1:nx,ny+1,2) x(nx,1:ny, 1) = x(nx,1:ny, 1) + x(0,1:ny, 2) x(nx:1:-1,1, 4) = x(nx:1:-1,1, 4) + x(nx+1,1:ny,2) x(1:nx,ny, 2) = x(1:nx,ny, 2) + x(1:nx,0, 3) x(1,ny:1:-1, 5) = x(1,ny:1:-1, 5) + x(1:nx,ny+1,3) x(nx:1:-1,ny,1) = x(nx:1:-1,ny,1) + x(0,1:ny, 3) x(1,1:ny, 4) = x(1,1:ny, 4) + x(nx+1,1:ny,3) x(nx,ny:1:-1,2) = x(nx,ny:1:-1,2) + x(1:nx,0, 4) x(1:nx,1, 5) = x(1:nx,1, 5) + x(1:nx,ny+1,4) x(nx,1:ny, 3) = x(nx,1:ny, 3) + x(0,1:ny, 4) x(nx:1:-1,1, 6) = x(nx:1:-1,1, 6) + x(nx+1,1:ny,4) x(1:nx,ny, 4) = x(1:nx,ny, 4) + x(1:nx,0, 5) x(1,ny:1:-1, 1) = x(1,ny:1:-1, 1) + x(1:nx,ny+1,5) x(nx:1:-1,ny,3) = x(nx:1:-1,ny,3) + x(0,1:ny, 5) x(1,1:ny, 6) = x(1,1:ny, 6) + x(nx+1,1:ny,5) x(nx,ny:1:-1,4) = x(nx,ny:1:-1,4) + x(1:nx,0, 6) x(1:nx,1, 1) = x(1:nx,1, 1) + x(1:nx,ny+1,6) x(nx,1:ny, 5) = x(nx,1:ny, 5) + x(0,1:ny, 6) x(nx:1:-1,1, 2) = x(nx:1:-1,1, 2) + x(nx+1,1:ny,6) ! Zero Halo x( :, 0,:) = 0.0 x( :,ny+1,:) = 0.0 x( 0, :,:) = 0.0 x(nx+1, :,:) = 0.0 end subroutine GhostCubeT ! Transforms a 2D vector on the surface of the sphere to ! a 3D Cartesian vector or 3D to 2D when (Inverse == True). ! It can also apply the transpose of either operation (Transpose==True). ! It can deal with 2D vectors defined either along Lat-Lon coordinates ! or along Cube-sphere grid lines (Cube==True) as defined by ! the FVCube dynamical core. The transforming coefficients for ! all of these operations are kept in Trans and are precomputed when ! it is initialized for a particular pair of grids when one grid is cube ! and the other is lat-lon. subroutine SphericalToCartesian(Tr, U, V, Uxyz, Transpose, SphIsLL) type(T_CubeLatLonTransform), intent(IN ) :: Tr real, intent(IN ) :: U(:,:,:), V(:,:,:) real, intent(OUT) :: Uxyz(:,:,:) logical, intent(IN ) :: Transpose logical, intent(IN ) :: SphIsLL integer :: K, LM real(R8), pointer :: e1(:,:,:), e2(:,:,:) if(SphIsLL) then e1=>Tr%elon e2=>Tr%elat else if(.not.Transpose) then e1=>Tr%ff1 e2=>Tr%ff2 else e1=>Tr%ee1 e2=>Tr%ee2 end if end if LM = size(U,3) do k=1,LM Uxyz(:,:,k ) = U(:,:,k)*e1(:,:,1) + V(:,:,k)*e2(:,:,1) Uxyz(:,:,k+ LM) = U(:,:,k)*e1(:,:,2) + V(:,:,k)*e2(:,:,2) Uxyz(:,:,k+2*LM) = U(:,:,k)*e1(:,:,3) + V(:,:,k)*e2(:,:,3) end do return end subroutine SphericalToCartesian subroutine CartesianToSpherical(Tr, Uxyz, U, V, Transpose, SphIsLL) type(T_CubeLatLonTransform), intent(IN ) :: Tr real, intent(OUT) :: U(:,:,:), V(:,:,:) real, intent(IN ) :: Uxyz(:,:,:) logical, intent(IN ) :: Transpose logical, intent(IN ) :: SphIsLL integer :: K, LM real(R8), pointer :: e1(:,:,:), e2(:,:,:) if(SphIsLL) then e1=>Tr%elon e2=>Tr%elat else if(Transpose) then e1=>Tr%ff1 e2=>Tr%ff2 else e1=>Tr%ee1 e2=>Tr%ee2 end if end if LM = size(U,3) do k=1,LM U(:,:,k) = Uxyz(:,:,k )*e1(:,:,1) + & Uxyz(:,:,k+ LM)*e1(:,:,2) + & Uxyz(:,:,k+2*LM)*e1(:,:,3) V(:,:,k) = Uxyz(:,:,k )*e2(:,:,1) + & Uxyz(:,:,k+ LM)*e2(:,:,2) + & Uxyz(:,:,k+2*LM)*e2(:,:,3) end do return end subroutine CartesianToSpherical end Module CubeLatLonTransformMod