LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - interpolate_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 72 | 665 | 10.8 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 2 | 28 | 7.1 %

              Line data    Source code
           1             : module interpolate_mod
           2             :   use shr_kind_mod,           only: r8=>shr_kind_r8
           3             :   use element_mod,            only: element_t
           4             :   use dimensions_mod,         only: np, ne, nelemd, nc, nhe, nhc
           5             :   use quadrature_mod,         only: quadrature_t, legendre, quad_norm
           6             :   use coordinate_systems_mod, only: spherical_polar_t, cartesian2d_t,         &
           7             :        cartesian3D_t, sphere2cubedsphere, spherical_to_cart,                  &
           8             :        cubedsphere2cart, distance, change_coordinates, projectpoint
           9             :   use physconst,              only: PI
          10             :   use quadrature_mod,         only: quadrature_t, gauss, gausslobatto
          11             :   use parallel_mod,           only: syncmp, parallel_t
          12             :   use cam_abortutils,         only: endrun
          13             :   use spmd_utils,             only: MPI_MAX, MPI_SUM, MPI_MIN, mpi_real8, MPI_integer
          14             :   use cube_mod,               only: convert_gbl_index, dmap, ref2sphere
          15             :   use mesh_mod,               only: MeshUseMeshFile
          16             :   use control_mod,            only: cubed_sphere_map
          17             :   use cam_logfile,            only: iulog
          18             :   use string_utils,           only: int2str
          19             :
          20             :   implicit none
          21             :   private
          22             :   save
          23             :
          24             :   logical   :: debug=.false.
          25             :
          26             :   type, public :: interpolate_t
          27             :      real (kind=r8), dimension(:,:), pointer :: Imat  ! P_k(xj)*wj/gamma(k)
          28             :      real (kind=r8), dimension(:)  , pointer :: rk    ! 1/k
          29             :      real (kind=r8), dimension(:)  , pointer :: vtemp ! temp results
          30             :      real (kind=r8), dimension(:)  , pointer :: glp   ! GLL pts (nair)
          31             :   end type interpolate_t
          32             :
          33             :   type, public :: interpdata_t
          34             :      ! Output Interpolation points.  Used to output data on lat-lon (or other grid)
          35             :      ! with native element interpolation.  Each element keeps a list of points from the
          36             :      ! interpolation grid that are in this element
          37             :      type (cartesian2D_t),pointer,dimension(:):: interp_xy      ! element coordinate
          38             :      integer, pointer,dimension(:)            :: ilat,ilon   ! position of interpolation point in lat-lon grid
          39             :      integer                                  :: n_interp
          40             :      integer                                  :: nlat
          41             :      integer                                  :: nlon
          42             :      logical                                  :: first_entry = .TRUE.
          43             :   end type interpdata_t
          44             :
          45             :   real (kind=r8), private :: delta  = 1.0e-9_r8  ! move tiny bit off center to
          46             :   ! avoid landing on element edges
          47             :
          48             :
          49             :   ! static data for interp_tracers
          50             :   logical                           :: interp_tracers_init=.false.
          51             :   real (kind=r8      )       :: interp_c(np,np)
          52             :   real (kind=r8      )       :: interp_gll(np)
          53             :
          54             :   public :: interp_init
          55             :   public :: setup_latlon_interp
          56             :   public :: interpolate_scalar
          57             :   public :: interpolate_ce
          58             :
          59             :   public :: interpol_phys_latlon
          60             :   public :: interpolate_vector
          61             :   public :: set_interp_parameter
          62             :   public :: get_interp_parameter
          63             :   public :: get_interp_gweight
          64             :   public :: get_interp_lat
          65             :   public :: get_interp_lon
          66             :   public :: cube_facepoint_ne
          67             :   public :: cube_facepoint_unstructured
          68             :   public :: parametric_coordinates
          69             :
          70             :   public :: interpolate_tracers
          71             :   public :: interpolate_tracers_init
          72             :   public :: minmax_tracers
          73             :   public :: interpolate_2d
          74             :   public :: interpolate_create
          75             :   public :: point_inside_quad
          76             :   public :: vec_latlon_to_contra
          77             :
          78             :
          79             :   interface interpolate_scalar
          80             :      module procedure interpolate_scalar2d
          81             :      module procedure interpolate_scalar3d
          82             :   end interface
          83             :   interface interpolate_vector
          84             :      module procedure interpolate_vector2d
          85             :      module procedure interpolate_vector3d
          86             :   end interface
          87             :
          88             :   type (interpolate_t), target ::  interp_p
          89             :
          90             :   ! store the  lat-lon grid
          91             :   ! gridtype = 1       equally spaced, including poles (FV scalars output grid)
          92             :   ! gridtype = 2       Gauss grid (CAM Eulerian)
          93             :   ! gridtype = 3       equally spaced, no poles (FV staggered velocity)
          94             :   ! Seven possible history files, last one is inithist and should be native grid
          95             :   integer :: nlat,nlon
          96             :   real (kind=r8), pointer, public   :: lat(:)     => NULL()
          97             :   real (kind=r8), pointer, public   :: lon(:)     => NULL()
          98             :   real (kind=r8), pointer, public   :: gweight(:) => NULL()
          99             :   integer :: gridtype = 1        !
         100             :   integer :: itype = 1           ! 0 = native high order
         101             :                                  ! 1 = bilinear
         102             :
         103             :   integer :: auto_grid = 0        ! 0 = interpolation grid set by namelist
         104             :                                   ! 1 = grid set via mesh resolution
         105             :
         106             :
         107             :   ! static data, used by bilin_phys2gll()
         108             :   ! shared by all threads.  only allocate if subroutine will be used
         109             : !JMD  integer  :: nphys_init=0
         110             : !JMD  integer  :: index_l(np),index_r(np)
         111             : !JMD  real(kind=r8),allocatable :: weights(:,:,:,:,:) ! np,np,2,2,nelemd
         112             :
         113             : !JMD  public :: bilin_phys2gll
         114             : !JMD  public :: bilin_phys2gll_init
         115             : contains
         116             :
         117             :
         118           0 :   subroutine set_interp_parameter(parm_name, value)
         119             :     character*(*), intent(in) :: parm_name
         120             :     character(len=80) :: msg
         121             :     integer :: value,power
         122             :     real (kind=r8) :: value_target
         123             :
         124           0 :     if(parm_name .eq. 'itype') then
         125           0 :        itype=value
         126           0 :     else if(parm_name .eq. 'nlon') then
         127           0 :        nlon=value
         128           0 :     else if(parm_name .eq. 'nlat') then
         129           0 :        nlat=value
         130           0 :     else if(parm_name.eq. 'gridtype') then
         131           0 :        gridtype=value
         132           0 :     else if(parm_name.eq. 'auto') then
         133           0 :        auto_grid=1
         134             :        ! compute recommended nlat,nlon which has slightly higher
         135             :        ! resolution than the specifed number of points around equator given in "value"
         136             :        ! computed recommended lat-lon grid.
         137             :        ! nlon > peq   peq = points around equator cubed sphere grid
         138             :        ! take nlon power of 2, and at most 1 power of 3
         139           0 :        if (value.eq.0) then
         140             :            ! If reading in unstructured mesh, ne = 0
         141             :            ! This makes it hard to guess how many interpolation points to use
         142             :            ! So We'll set the default as 720 x 360
         143             :            ! BUT if you're running with an unstructured mesh, set interp_nlon and interp_nlat
         144           0 :            nlon = 1536
         145           0 :            nlat = 768
         146             :        else
         147           0 :            value_target=value*1.25_r8
         148           0 :            power = nint(0.5_r8 +  log( value_target)/log(2.0_r8) )
         149           0 :            power = max(power,7) ! min grid: 64x128
         150           0 :            if ( 3*2**(power-2) > value_target) then
         151           0 :                nlon=3*2**(power-2)   ! use 1 power of 3
         152             :            else
         153           0 :                nlon=2**power
         154             :            endif
         155             :        endif
         156           0 :        nlat=nlon/2
         157           0 :        if (gridtype==1) nlat=nlat+1
         158             :     else
         159           0 :        write(msg,*) 'Did not recognize parameter named ',parm_name,' in interpolate_mod:set_interp_parameter'
         160           0 :        call endrun(msg)
         161             :     end if
         162           0 :   end subroutine set_interp_parameter
         163           0 :   function get_interp_parameter(parm_name) result(value)
         164             :     character*(*), intent(in) :: parm_name
         165             :     integer :: value
         166             :     character(len=80) :: msg
         167           0 :     if(parm_name .eq. 'itype') then
         168           0 :        value=itype
         169           0 :     else if(parm_name .eq. 'nlon') then
         170           0 :        value=nlon
         171           0 :     else if(parm_name .eq. 'nlat') then
         172           0 :        value=nlat
         173           0 :     else if(parm_name.eq. 'gridtype') then
         174           0 :        value=gridtype
         175           0 :     else if(parm_name.eq. 'auto_grid') then
         176           0 :        value=auto_grid
         177             :     else
         178           0 :        write(msg,*) 'Did not recognize parameter named ',parm_name,' in interpolate_mod:get_interp_parameter'
         179           0 :        value=-1
         180           0 :        call endrun(msg)
         181             :     end if
         182             :     return
         183             :   end function get_interp_parameter
         184           0 :   function get_interp_gweight() result(gw)
         185             :     real(kind=r8) :: gw(nlat)
         186           0 :     gw=gweight
         187           0 :     return
         188           0 :   end function get_interp_gweight
         189           0 :   function get_interp_lat() result(thislat)
         190             :     real(kind=r8) :: thislat(nlat)
         191           0 :     thislat=lat*180.0_r8/PI
         192           0 :     return
         193           0 :   end function get_interp_lat
         194           0 :   function get_interp_lon() result(thislon)
         195             :     real(kind=r8) :: thislon(nlon)
         196           0 :     thislon=lon*180.0_r8/PI
         197           0 :     return
         198           0 :   end function get_interp_lon
         199             :
         200     5216400 :   subroutine interpolate_create(gquad,interp)
         201             :     type (quadrature_t) , intent(in)   :: gquad
         202             :     type (interpolate_t), intent(out)  :: interp
         203             :
         204             :
         205             :     ! Local variables
         206             :
         207             :     integer k,j
         208             :     integer npts
         209     5216400 :     real (kind=r8), dimension(:), allocatable :: gamma
         210     5216400 :     real (kind=r8), dimension(:), allocatable :: leg
         211             :
         212     5216400 :     npts = size(gquad%points)
         213             :
         214    20865600 :     allocate(interp%Imat(npts,npts))
         215    15649200 :     allocate(interp%rk(npts))
         216    10432800 :     allocate(interp%vtemp(npts))
         217    10432800 :     allocate(interp%glp(npts))
         218    10432800 :     allocate(gamma(npts))
         219    10432800 :     allocate(leg(npts))
         220             :
         221     5216400 :     gamma = quad_norm(gquad,npts)
         222             :
         223    26082000 :     do k=1,npts
         224    20865600 :        interp%rk(k) = 1.0_r8/k
         225    26082000 :        interp%glp(k) = gquad%points(k)    !nair
         226             :     end do
         227             :
         228    26082000 :     do j=1,npts
         229    20865600 :        leg=legendre(gquad%points(j),npts-1)
         230   109544400 :        do k=1,npts
         231   104328000 :           interp%Imat(j,k)=leg(k)*gquad%weights(j)/gamma(k)
         232             :        end do
         233             :     end do
         234             :
         235     5216400 :     deallocate(gamma)
         236     5216400 :     deallocate(leg)
         237             :
         238     5216400 :   end subroutine interpolate_create
         239             :
         240             :
         241           0 :   subroutine interpolate_tracers_init()
         242             :     use dimensions_mod, only : np, qsize
         243             :     use quadrature_mod, only : quadrature_t, gausslobatto
         244             :
         245             :
         246             :     implicit none
         247             :
         248             :     type (quadrature_t        )       :: gll
         249             :     real (kind=r8      )       :: dp    (np)
         250             :     integer                           :: i,j
         251             :
         252           0 :       gll=gausslobatto(np)
         253           0 :       dp = 1
         254           0 :       do i=1,np
         255           0 :         do j=1,np
         256           0 :           if (i /= j) then
         257           0 :             dp(i) = dp(i) * (gll%points(i) - gll%points(j))
         258             :           end if
         259             :         end do
         260             :       end do
         261           0 :       do i=1,np
         262           0 :         do j=1,np
         263           0 :           interp_c(i,j) = 1/(dp(i)*dp(j))
         264             :         end do
         265             :       end do
         266           0 :       interp_gll(:) = gll%points(:)
         267           0 :       interp_tracers_init = .true.
         268             :
         269           0 :       deallocate(gll%points)
         270           0 :       deallocate(gll%weights)
         271             :
         272             :
         273           0 :   end subroutine interpolate_tracers_init
         274             :
         275             :
         276             :
         277             :
         278           0 :   subroutine interpolate_tracers(r, tracers, f)
         279             :     use dimensions_mod, only : np, qsize
         280             :
         281             :
         282             :     implicit none
         283             :     type (cartesian2D_t), intent(in)  :: r
         284             :     real (kind=r8),intent(in)  :: tracers(np*np,qsize)
         285             :     real (kind=r8),intent(out) :: f(qsize)
         286             :
         287             :     real (kind=r8      )       :: x     (np)
         288             :     real (kind=r8      )       :: y     (np)
         289             :     real (kind=r8      )       :: xy    (np*np)
         290             :
         291             :     integer                           :: i,j
         292             :
         293             :
         294           0 :     if (.not. interp_tracers_init   ) then
         295           0 :        call endrun('ERROR: interpolate_tracers() was not initialized')
         296             :     endif
         297             :
         298           0 :     x = 1
         299           0 :     y = 1
         300           0 :     do i=1,np
         301           0 :       do j=1,np
         302           0 :         if (i /= j) then
         303           0 :           x(i) = x(i) * (r%x - interp_gll(j))
         304           0 :           y(i) = y(i) * (r%y - interp_gll(j))
         305             :         end if
         306             :       end do
         307             :     end do
         308             :
         309           0 :     do j=1,np
         310           0 :       do i=1,np
         311           0 :         xy(i + (j-1)*np) = x(i)*y(j)*interp_c(i,j)
         312             :       end do
         313             :     end do
         314           0 :     f = MATMUL(xy,tracers)
         315           0 :   end subroutine interpolate_tracers
         316             :
         317             :
         318             :
         319             :   subroutine linear_interpolate_2d(x,y,s,v)
         320             :     use dimensions_mod, only : np, qsize
         321             :
         322             :     real(kind=r8) ,       intent(in) :: x(np)
         323             :     real(kind=r8),        intent(in) :: y(np,np,qsize)
         324             :     type (cartesian2D_t), intent(in) :: s
         325             :     real(kind=r8),        intent(inout) :: v(qsize)
         326             :
         327             :     integer                           :: i,j,q
         328             :     real (kind=r8)  dx, dy(qsize), dydx(qsize)
         329             :     real (kind=r8)  y0(qsize), y1(qsize)
         330             :     type (cartesian2D_t)              :: r
         331             :
         332             :     r = s
         333             :     if (r%x < -1) r%x = -1
         334             :     if (r%y < -1) r%y = -1
         335             :     if ( 1 < r%x) r%x =  1
         336             :     if ( 1 < r%y) r%y =  1
         337             :     do i=1,np
         338             :       if (r%x < x(i)) exit
         339             :     end do
         340             :     do j=1,np
         341             :       if (r%y < x(j)) exit
         342             :     end do
         343             :     if (1 < i) i = i-1
         344             :     if (1 < j) j = j-1
         345             :     if (np==i) i = i-1
         346             :     if (np==j) j = j-1
         347             :
         348             :     dx = x(i+1)     - x(i)
         349             :     dy = y(i+1,j,:) - y(i,j,:)
         350             :     dydx = dy/dx
         351             :     y0 = y(i,j,:) + (r%x-x(i))*dydx
         352             :
         353             :     dy = y(i+1,j+1,:) - y(i,j+1,:)
         354             :     dydx = dy/dx
         355             :     y1 = y(i,j+1,:) + (r%x-x(i))*dydx
         356             :
         357             :     dx = x(j+1)     - x(j)
         358             :     dy = y1         - y0
         359             :     dydx = dy/dx
         360             :     v  = y0         + (r%y-x(j))*dydx
         361             :
         362             :   end subroutine linear_interpolate_2d
         363             :
         364           0 :   subroutine minmax_tracers(r, tracers, mint, maxt)
         365             :     use dimensions_mod, only : np, qsize
         366             :     use quadrature_mod, only : quadrature_t, gausslobatto
         367             :
         368             :
         369             :     implicit none
         370             :
         371             :     type (cartesian2D_t), intent(in)  :: r
         372             :     real (kind=r8),intent(in)  :: tracers(np,np,qsize)
         373             :     real (kind=r8),intent(out) :: mint         (qsize)
         374             :     real (kind=r8),intent(out) :: maxt         (qsize)
         375             :
         376             :     type (quadrature_t), save         :: gll
         377             :     integer                           :: i,j
         378             :     logical            , save         :: first_time=.true.
         379             :     real (kind=r8)             :: y1           (qsize)
         380             :     real (kind=r8)             :: y2           (qsize)
         381           0 :     real (kind=r8)             :: q_interp     (4,qsize)
         382             :     type (cartesian2D_t)              :: s
         383             :     real (kind=r8)             :: delta
         384             :     integer :: q
         385             :
         386           0 :     do q=1,qsize
         387           0 :        mint(q) = minval(tracers(:,:,q))
         388           0 :        maxt(q) = maxval(tracers(:,:,q))
         389             :     enddo
         390             :     return
         391             :
         392             :     delta = 1._r8/8._r8
         393             :
         394             :     if (first_time) then
         395             :       first_time = .false.
         396             :       gll=gausslobatto(np)
         397             :     end if
         398             :
         399             :     do i=1,np
         400             :       if (r%x < gll%points(i)) exit
         401             :     end do
         402             :     do j=1,np
         403             :       if (r%y < gll%points(j)) exit
         404             :     end do
         405             :     if (1 < i) i = i-1
         406             :     if (1 < j) j = j-1
         407             :     if (np==i) i = i-1
         408             :     if (np==j) j = j-1
         409             :
         410             : !   mint(:) = minval(minval(tracers(i:i+1,j:j+1,:),1),1)
         411             : !   maxt(:) = maxval(maxval(tracers(i:i+1,j:j+1,:),1),1)
         412             :
         413             : ! Or check this out:
         414             :     s   = r
         415             :     s%x = s%x - delta
         416             :     s%y = s%y - delta
         417             :     call linear_interpolate_2d(gll%points,tracers,s,q_interp(1,:))
         418             :     s   = r
         419             :     s%x = s%x + delta
         420             :     s%y = s%y - delta
         421             :     call linear_interpolate_2d(gll%points,tracers,s,q_interp(2,:))
         422             :     s   = r
         423             :     s%x = s%x - delta
         424             :     s%y = s%y + delta
         425             :     call linear_interpolate_2d(gll%points,tracers,s,q_interp(3,:))
         426             :     s   = r
         427             :     s%x = s%x + delta
         428             :     s%y = s%y + delta
         429             :     call linear_interpolate_2d(gll%points,tracers,s,q_interp(4,:))
         430             :
         431             :     mint(:) = minval(q_interp(:,:),1)
         432             :     maxt(:) = maxval(q_interp(:,:),1)
         433             :   end subroutine minmax_tracers
         434             :
         435  8732253600 :   function interpolate_2d(cart, f, interp, npts, fillvalue) result(fxy)
         436             :     integer, intent(in)               :: npts
         437             :     type (cartesian2D_t), intent(in)  :: cart
         438             :     real (kind=r8), intent(in) :: f(npts,npts)
         439             :     type (interpolate_t)              :: interp
         440             :     real (kind=r8)             :: fxy     ! value of f interpolated to (x,y)
         441             :     real (kind=r8), intent(in), optional :: fillvalue
         442             :     ! local variables
         443             :
         444             :     real (kind=r8)             :: tmp_1,tmp_2
         445             :     real (kind=r8)             :: fk0,fk1
         446             :     real (kind=r8)             :: pk
         447             :
         448             :     integer                           :: l,j,k
         449             :
         450  8732253600 :     if(present(fillvalue)) then
         451           0 :        if (any(f==fillvalue)) then
         452  8732253600 :           fxy = fillvalue
         453             :           return
         454             :        endif
         455             :     endif
         456             :
         457             :
         458 26196760800 :     do l=1,npts,2
         459             :
         460             :        ! Compute Pk(cart%x) for Legendre order 0
         461             :
         462 17464507200 :        pk = 1.0_r8
         463             :
         464 17464507200 :        fk0=0.0_r8
         465 17464507200 :        fk1=0.0_r8
         466 87322536000 :        do j=1,npts
         467 69858028800 :           fk0 = fk0 + interp%Imat(j,1)*f(j,l  )
         468 87322536000 :           fk1 = fk1 + interp%Imat(j,1)*f(j,l+1)
         469             :        end do
         470 17464507200 :        interp%vtemp(l  ) = pk*fk0
         471 17464507200 :        interp%vtemp(l+1) = pk*fk1
         472             :
         473             :        ! Compute Pk(cart%x) for Legendre order 1
         474             :
         475 17464507200 :        tmp_2 = pk
         476 17464507200 :        pk    = cart%x
         477             :
         478 17464507200 :        fk0=0.0_r8
         479 17464507200 :        fk1=0.0_r8
         480 87322536000 :        do j=1,npts
         481 69858028800 :           fk0 = fk0 + interp%Imat(j,2)*f(j,l  )
         482 87322536000 :           fk1 = fk1 + interp%Imat(j,2)*f(j,l+1)
         483             :        end do
         484 17464507200 :        interp%vtemp(l  ) = interp%vtemp(l  ) + pk*fk0
         485 17464507200 :        interp%vtemp(l+1) = interp%vtemp(l+1) + pk*fk1
         486             :
         487             :        ! Compute Pk(cart%x) for Legendre order 2 to npts-1
         488             :
         489 61125775200 :        do k = 2,npts-1
         490             :
         491 34929014400 :           tmp_1  = tmp_2
         492 34929014400 :           tmp_2  = pk
         493 34929014400 :           pk = ( (2*k-1)*cart%x*tmp_2 - (k-1)*tmp_1 )*interp%rk(k)
         494             :
         495 34929014400 :           fk0=0.0_r8
         496 34929014400 :           fk1=0.0_r8
         497 >17464*10^7 :           do j=1,npts
         498 >13971*10^7 :              fk0 = fk0 + interp%Imat(j,k+1)*f(j,l  )
         499 >17464*10^7 :              fk1 = fk1 + interp%Imat(j,k+1)*f(j,l+1)
         500             :           end do
         501 34929014400 :           interp%vtemp(l  ) = interp%vtemp(l  ) + pk*fk0
         502 52393521600 :           interp%vtemp(l+1) = interp%vtemp(l+1) + pk*fk1
         503             :
         504             :        end do
         505             :
         506             :     end do
         507             :
         508             :     ! Compute Pk(cart%y) for Legendre order 0
         509             :
         510  8732253600 :     pk = 1.0_r8
         511             :
         512  8732253600 :     fk0 = 0.0_r8
         513 43661268000 :     do j=1,npts
         514 43661268000 :        fk0 = fk0 + interp%Imat(j,1)*interp%vtemp(j)
         515             :     end do
         516  8732253600 :     fxy = pk*fk0
         517             :
         518             :     ! Compute Pk(cart%y) for Legendre order 1
         519             :
         520  8732253600 :     tmp_2 = pk
         521  8732253600 :     pk    = cart%y
         522             :
         523  8732253600 :     fk0=0.0_r8
         524 43661268000 :     do j=1,npts
         525 43661268000 :        fk0 = fk0 + interp%Imat(j,2)*interp%vtemp(j)
         526             :     end do
         527  8732253600 :     fxy = fxy + pk*fk0
         528             :
         529             :     ! Compute Pk(cart%y) for Legendre order 2, npts-1
         530             :
         531 26196760800 :     do k = 2,npts-1
         532 17464507200 :        tmp_1  = tmp_2
         533 17464507200 :        tmp_2  = pk
         534 17464507200 :        pk = ( (2*k-1)*cart%y*tmp_2 - (k-1)*tmp_1 )*interp%rk(k)
         535             :
         536 17464507200 :        fk0 = 0.0_r8
         537 87322536000 :        do j=1,npts
         538 87322536000 :           fk0 = fk0 + interp%Imat(j,k+1)*interp%vtemp(j)
         539             :        end do
         540             :
         541 26196760800 :        fxy = fxy + pk*fk0
         542             :
         543             :     end do
         544             :
         545             :   end function interpolate_2d
         546             :
         547             :   !===============================
         548             :   !(Nair) Bilinear interpolation for every GLL grid cell
         549             :   !===============================
         550             :
         551           0 :   function interpol_bilinear(cart, f, xoy, imin, imax, fillvalue) result(fxy)
         552             :     integer, intent(in)               :: imin,imax
         553             :     type (cartesian2D_t), intent(in)  :: cart
         554             :     real (kind=r8), intent(in) :: f(imin:imax,imin:imax)
         555             :     real (kind=r8)             :: xoy(imin:imax)
         556             :     real (kind=r8)             :: fxy     ! value of f interpolated to (x,y)
         557             :     real (kind=r8), intent(in), optional :: fillvalue
         558             :     ! local variables
         559             :
         560             :     real (kind=r8) :: p,q,xp,yp ,y4(4)
         561             :     integer        :: l,j,k, ii, jj, na,nb,nm
         562             :
         563           0 :     xp = cart%x
         564           0 :     yp = cart%y
         565             :
         566             :     ! Search index along "x"  (bisection method)
         567             :
         568           0 :     na = imin
         569           0 :     nb = imax
         570             :     do
         571           0 :        if  ((nb-na) <=  1)  exit
         572           0 :        nm = (nb + na)/2
         573           0 :        if (xp  >  xoy(nm)) then
         574             :           na = nm
         575             :        else
         576           0 :           nb = nm
         577             :        endif
         578             :     enddo
         579             :     ii = na
         580             :
         581             :     ! Search index along "y"
         582             :
         583             :     na = imin
         584             :     nb = imax
         585             :     do
         586           0 :        if  ((nb-na) <=  1)  exit
         587           0 :        nm = (nb + na)/2
         588           0 :        if (yp  >  xoy(nm)) then
         589             :           na = nm
         590             :        else
         591           0 :           nb = nm
         592             :        endif
         593             :     enddo
         594           0 :     jj = na
         595             :
         596             :     ! GLL cell containing (xp,yp)
         597             :
         598           0 :     y4(1) = f(ii,jj)
         599           0 :     y4(2) = f(ii+1,jj)
         600           0 :     y4(3) = f(ii+1,jj+1)
         601           0 :     y4(4) = f(ii,jj+1)
         602             :
         603           0 :     if(present(fillvalue)) then
         604           0 :        if (any(y4==fillvalue)) then
         605           0 :           fxy = fillvalue
         606             :           return
         607             :        endif
         608             :     endif
         609             :
         610           0 :     p = (xp - xoy(ii))/(xoy(ii+1) - xoy(ii))
         611           0 :     q = (yp - xoy(jj))/(xoy(jj+1) - xoy(jj))
         612             :
         613             :     fxy = (1.0_r8 - p)*(1.0_r8 - q)* y4(1) + p*(1.0_r8 - q) * y4(2)   &
         614           0 :          + p*q* y4(3) + (1.0_r8 - p)*q * y4(4)
         615           0 :   end function interpol_bilinear
         616             :
         617             :   ! ----------------------------------------------------------------------------------!
         618             :   !FUNCTION   interpol_phys_latlon----------------------------------------CE-for fvm!
         619             :   ! AUTHOR: CHRISTOPH ERATH, 23. May 2012                                             !
         620             :   ! DESCRIPTION: evaluation of the reconstruction for every physics grid cell         !
         621             :   !                                                                                   !
         622             :   ! CALLS:
         623             :   ! INPUT:
         624             :   !
         625             :   ! OUTPUT:
         626             :   !-----------------------------------------------------------------------------------!
         627           0 :   subroutine interpol_phys_latlon(interpdata,f, fvm, corners, desc, flatlon,lmono)
         628             :     use fvm_control_volume_mod, only : fvm_struct
         629             :     !  use fvm_reconstruction_mod, only: reconstruction_gradient, recons_val_cart
         630             :     use edgetype_mod, only : edgedescriptor_t
         631             :
         632             :     type (interpdata_t), intent(in)     :: interpdata
         633             :     real (kind=r8), intent(inout)   :: f(1-nhc:nc+nhc,1-nhc:nc+nhc)
         634             :     type (fvm_struct), intent(in)       :: fvm
         635             :     type (cartesian2d_t), intent(in)    :: corners(:)
         636             :     type (edgedescriptor_t),intent(in)  :: desc
         637             :     logical, intent(in) :: lmono
         638             :
         639             :     real (kind=r8)             :: flatlon(:)
         640             :     ! local variables
         641             :     real (kind=r8)             :: xp,yp, tmpval
         642             :     real (kind=r8)             :: tmpaxp,tmpaxm, tmpayp, tmpaym
         643             :     integer                           :: i, ix, jy, starti,endi,tmpi
         644             :     real (kind=r8), dimension(1-nhe:nc+nhe,1-nhe:nc+nhe,6)      :: recons
         645             :
         646             :     real (kind=r8), dimension(nc+1) :: x, y
         647             :
         648             :     !  call reconstruction_gradient(f, fvm,recons,6,lmono)
         649             :     !  recons=0.0 ! PCoM
         650             :
         651           0 :     x(1:nc) = fvm%vtx_cart(1,1,1:nc,1   )
         652           0 :     y(1:nc) = fvm%vtx_cart(1,2,1   ,1:nc)
         653           0 :     x(nc+1) = fvm%vtx_cart(2,1,nc,1     )
         654           0 :     y(nc+1) = fvm%vtx_cart(3,2,1   ,nc  )
         655             :
         656           0 :     tmpaxp=(corners(1)%x+corners(2)%x)/2
         657           0 :     tmpaxm=(corners(2)%x-corners(1)%x)/2
         658           0 :     tmpayp=(corners(1)%y+corners(4)%y)/2
         659           0 :     tmpaym=(corners(4)%y-corners(1)%y)/2
         660           0 :     do i=1,interpdata%n_interp
         661             :       ! caculation phys grid coordinate of xp point, note the interp_xy are on the reference [-1,1]x[-1,1]
         662           0 :       xp=tan(tmpaxp+interpdata%interp_xy(i)%x*tmpaxm)
         663           0 :       yp=tan(tmpayp+interpdata%interp_xy(i)%y*tmpaym)
         664             :
         665             :       ! Search index along "x"  (bisection method)
         666           0 :       starti = 1
         667           0 :       endi = nc+1
         668             :       do
         669           0 :         if  ((endi-starti) <=  1)  exit
         670           0 :         tmpi = (endi + starti)/2
         671           0 :         if (xp  >  x(tmpi)) then
         672             :           starti = tmpi
         673             :         else
         674           0 :           endi = tmpi
         675             :         endif
         676             :       enddo
         677             :       ix = starti
         678             :
         679             :       ! Search index along "y"
         680             :       starti = 1
         681             :       endi = nc+1
         682             :       do
         683           0 :         if  ((endi-starti) <=  1)  exit
         684           0 :         tmpi = (endi + starti)/2
         685           0 :         if (yp  >  y(tmpi)) then
         686             :           starti = tmpi
         687             :         else
         688           0 :           endi = tmpi
         689             :         endif
         690             :       enddo
         691           0 :       jy = starti
         692             :
         693             :       !    call recons_val_cart(f(ix,jy), xp,yp, fvm%spherecentroid(ix,jy,:), fvm%recons_metrics(ix,jy,:), &
         694             :       !         recons(ix,jy,:), tmpval)
         695           0 :       tmpval=f(ix,jy)
         696           0 :       flatlon(i)=tmpval
         697             :       !phl PCoM
         698             :       !    flatlon(i)=f(ix,jy)
         699             :     end do
         700           0 :   end subroutine interpol_phys_latlon
         701             :
         702           0 :   function parametric_coordinates(sphere, corners3D,ref_map_in, corners,u2qmap,facenum) result (ref)
         703             :     implicit none
         704             :     type (spherical_polar_t), intent(in) :: sphere
         705             :     type (cartesian2D_t) :: ref
         706             :
         707             :     type (cartesian3D_t)   :: corners3D(4)  !x,y,z coords of element corners
         708             :     integer,optional  :: ref_map_in    ! default is global variable 'cubed_sphere_map'
         709             :     ! optional arguments, only needed for ref_map=1 (equi-angle gnomonic projection):
         710             :     type (cartesian2D_t),optional   :: corners(4)    ! gnomonic coords of element corners
         711             :     real (kind=r8),optional  :: u2qmap(4,2)
         712             :     integer,optional  :: facenum
         713             :
         714             :
         715             :     ! local
         716             :     integer               :: i, MAX_NR_ITER=10
         717             :     real(kind=r8)  :: D(2,2),Dinv(2,2),detD,a,b,resa,resb,dela,delb,costh
         718             :     real(kind=r8)  :: tol_sq = 1.0e-26_r8
         719             :     type (spherical_polar_t) :: sphere1, sphere_tmp
         720             :     integer  :: ref_map
         721             :
         722             :     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         723             :     ! newton iteration on: ref=ref - df^-1 (ref2sphere(ref) - sphere)
         724             :     !
         725             :     ! Generic version written in terms of HOMME's 'ref2sphere' and 'Dmap' operaters,
         726             :     ! with no assumption as to the type of map (gnomonic, equi-angular, parametric)
         727             :     !
         728             :     ! Note that the coordinate increment from newton iterations is not a direction and thus
         729             :     ! should not be converted into motion along a great circle arc - this routine
         730             :     ! correclty applies the increment by just adding it to the coordintes
         731             :     !
         732             :     ! f = ref2sphere(xvec) - sphere
         733             :     ! df = d(ref2sphere)
         734             :     !
         735             :     ! D = diag(cos(theta),1) * d(ref2sphere)       d(ref2sphere) = diag(1/cos(theta),1)*D
         736             :     !
         737             :     ! df = diag(1/cos(theta),1)*D
         738             :     ! df^-1 =  D^-1 *  diag(cos(theta),1)
         739             :     !
         740             :     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         741           0 :     if (present(ref_map_in)) then
         742           0 :        ref_map=ref_map_in
         743             :     else
         744           0 :        ref_map=cubed_sphere_map
         745             :     endif
         746           0 :     costh=cos(sphere%lat)
         747           0 :     a=0
         748           0 :     b=0
         749           0 :     i=0
         750             :     do
         751           0 :        sphere1 = ref2sphere(a,b,corners3D,ref_map,corners,facenum)
         752           0 :        resa = sphere1%lon - sphere%lon
         753           0 :        if (resa >  pi) resa= resa - 2*pi
         754           0 :        if (resa < -pi) resa= resa + 2*pi
         755             :
         756           0 :        resb = sphere1%lat - sphere%lat
         757             :
         758           0 :        call Dmap(D,a,b,corners3D,ref_map,corners,u2qmap,facenum)
         759           0 :        detD = D(1,1)*D(2,2) - D(1,2)*D(2,1)
         760           0 :        Dinv(1,1) =  D(2,2)/detD
         761           0 :        Dinv(1,2) = -D(1,2)/detD
         762           0 :        Dinv(2,1) = -D(2,1)/detD
         763           0 :        Dinv(2,2) =  D(1,1)/detD
         764             :
         765           0 :        dela =  Dinv(1,1)*costh*resa + Dinv(1,2)*resb
         766           0 :        delb =  Dinv(2,1)*costh*resa + Dinv(2,2)*resb
         767           0 :        a = a - dela
         768           0 :        b = b - delb
         769           0 :        i=i+1
         770           0 :        if ( (costh*resa)**2 + resb**2 < tol_sq .or. MAX_NR_ITER < i) exit
         771             :     end do
         772           0 :     ref%x=a
         773           0 :     ref%y=b
         774             :
         775           0 :   end function parametric_coordinates
         776             :
         777             :
         778             :
         779             :
         780             : !
         781             : ! find element containing given point, useing HOMME's standard
         782             : ! equi-angular gnomonic map.
         783             : ! note that with this map, only coordinate lines are great circle arcs
         784             : !
         785           0 :   function point_inside_equiangular(elem, sphere, sphere_xyz) result(inside)
         786             :     implicit none
         787             :     type (spherical_polar_t), intent(in)     :: sphere
         788             :     type (cartesian3D_t),     intent(in)    :: sphere_xyz
         789             :     type (element_t)        , intent(in)     :: elem
         790             :     logical                              :: inside, inside2
         791             :     integer               :: i,j
         792             :     type (cartesian2D_t) :: corners(4),sphere_xy,cart
         793             :     type (cartesian3D_t) :: corners_xyz(4),center,a,b,cross(4)
         794             :     real (kind=r8) :: yp(4), y, elem_diam,dotprod
         795             :     real (kind=r8) :: xp(4), x, xc,yc
         796             :     real (kind=r8) :: tol_inside
         797             :     real (kind=r8) :: d1,d2
         798             :
         799             :     type (spherical_polar_t)    :: sphere_tmp
         800             :
         801           0 :     inside = .false.
         802             :
         803             :
         804             :     ! first check if point is near the element:
         805           0 :     corners_xyz(:) = elem%corners3D(:)
         806             :     elem_diam = max( distance(corners_xyz(1),corners_xyz(3)), &
         807           0 :          distance(corners_xyz(2),corners_xyz(4)) )
         808             :
         809           0 :     center%x = sum(corners_xyz(1:4)%x)/4
         810           0 :     center%y = sum(corners_xyz(1:4)%y)/4
         811           0 :     center%z = sum(corners_xyz(1:4)%z)/4
         812           0 :     if ( distance(center,sphere_xyz) > 1.0_r8*elem_diam ) return
         813             :
         814           0 :     tol_inside = 1.0e-10_r8*elem_diam**2
         815             :     ! the point is close to the element, so project both to cubed sphere
         816             :     ! and perform contour integral
         817           0 :     sphere_xy=sphere2cubedsphere(sphere,elem%FaceNum)
         818           0 :     x = sphere_xy%x
         819           0 :     y = sphere_xy%y
         820           0 :     do i=1,4
         821           0 :       xp(i) = elem%corners(i)%x
         822           0 :       yp(i) = elem%corners(i)%y
         823             :     end do
         824             :
         825             :
         826           0 :     if (debug) then
         827           0 :        print *,'point: ',x,y,elem%FaceNum
         828           0 :        print *,'element:'
         829           0 :        write(*,'(a,4e16.8,a)') 'x=[',xp(1:4),']'
         830           0 :        write(*,'(a,4e16.8,a)') 'y=[',yp(1:4),']'
         831             :
         832             :        ! first check if centroid is in this element (sanity check)
         833           0 :        sphere_tmp=change_coordinates(center)
         834           0 :        sphere_xy=sphere2cubedsphere(sphere_tmp,elem%FaceNum)
         835           0 :        xc=sphere_xy%x
         836           0 :        yc=sphere_xy%y
         837           0 :        print *,'cross product with centroid: all numbers should be negative'
         838           0 :        j = 4
         839           0 :        do i=1,4
         840           0 :           print *,i,(xc-xp(j))*(yp(i)-yp(j))  - (yc-yp(j))*(xp(i)-xp(j))
         841           0 :           j = i  ! within this loopk j = i-1
         842             :        end do
         843             :
         844           0 :        print *,'cross product with search point'
         845           0 :        j = 4
         846           0 :        do i=1,4
         847           0 :           print *,i,(x-xp(j))*(yp(i)-yp(j))  - (y-yp(j))*(xp(i)-xp(j))
         848           0 :           j = i  ! within this loopk j = i-1
         849             :        end do
         850             :     endif
         851             :
         852             :
         853           0 :     j = 4
         854           0 :     do i=1,4
         855             :       ! a = x-xp(j), y-yp(j)
         856             :       ! b = xp(i)-xp(j), yp(i)-yp(j)
         857             :       ! compute a cross b:
         858           0 :       if ( -( (x-xp(j))*(yp(i)-yp(j))  - (y-yp(j))*(xp(i)-xp(j))) > tol_inside ) then
         859             :          return
         860             :       endif
         861           0 :       j = i  ! within this loopk j = i-1
         862             :     end do
         863             :     ! all cross products were negative, must be inside:
         864           0 :     inside=.true.
         865             :   end function point_inside_equiangular
         866             :
         867             :
         868             : !
         869             : ! find if quad contains given point, with quad edges assumed to be great circle arcs
         870             : ! this will work with any map where straight lines are mapped to great circle arcs.
         871             : ! (thus it will fail on unstructured grids using the equi-angular gnomonic map)
         872             : !
         873           0 :   function point_inside_quad(corners_xyz, sphere_xyz) result(inside)
         874             :     implicit none
         875             :     type (cartesian3D_t),     intent(in)    :: sphere_xyz
         876             :     type (cartesian3D_t)    , intent(in)    :: corners_xyz(4)
         877             :     logical                              :: inside, inside2
         878             :     integer               :: i,j,ii
         879             :     type (cartesian2D_t) :: corners(4),sphere_xy,cart
         880             :     type (cartesian3D_t) :: center,a,b,cross(4)
         881             :     real (kind=r8) :: yp(4), y, elem_diam,dotprod
         882             :     real (kind=r8) :: xp(4), x
         883             :     real (kind=r8) :: d1,d2, tol_inside = 1.0e-12_r8
         884             :
         885             :     type (spherical_polar_t)   :: sphere  ! debug
         886             :
         887           0 :     inside = .false.
         888             :
         889             :     ! first check if point is near the corners:
         890             :     elem_diam = max( distance(corners_xyz(1),corners_xyz(3)), &
         891           0 :          distance(corners_xyz(2),corners_xyz(4)) )
         892             :
         893           0 :     center%x = sum(corners_xyz(1:4)%x)/4
         894           0 :     center%y = sum(corners_xyz(1:4)%y)/4
         895           0 :     center%z = sum(corners_xyz(1:4)%z)/4
         896           0 :     if ( distance(center,sphere_xyz) > 1.0_r8*elem_diam ) return
         897             :
         898             :     j = 4
         899           0 :     do i=1,4
         900             :       ! outward normal to plane containing j->i edge:  corner(i) x corner(j)
         901             :       ! sphere dot (corner(i) x corner(j) ) = negative if inside
         902           0 :        cross(i)%x =  corners_xyz(i)%y*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%y
         903           0 :        cross(i)%y =-(corners_xyz(i)%x*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%x)
         904           0 :        cross(i)%z =  corners_xyz(i)%x*corners_xyz(j)%y - corners_xyz(i)%y*corners_xyz(j)%x
         905             :        dotprod = cross(i)%x*sphere_xyz%x + cross(i)%y*sphere_xyz%y +&
         906           0 :                cross(i)%z*sphere_xyz%z
         907           0 :        j = i  ! within this loopk j = i-1
         908             :
         909             :        ! dot product is proportional to elem_diam. positive means outside,
         910             :        ! but allow machine precision tolorence:
         911           0 :        if (dotprod > tol_inside*elem_diam) return
         912             :        !if (dotprod > 0) return
         913             :     end do
         914           0 :     inside=.true.
         915             :     return
         916             :   end function point_inside_quad
         917             :
         918             : !
         919             : ! find element containing given point, with element edges assumed to be great circle arcs
         920             : ! this will work with any map where straight lines are mapped to great circle arcs.
         921             : ! (thus it will fail on unstructured grids using the equi-angular gnomonic map)
         922             : !
         923           0 :   function point_inside_gc(elem, sphere_xyz) result(inside)
         924             :     implicit none
         925             :     type (cartesian3D_t),     intent(in)    :: sphere_xyz
         926             :     type (element_t)        , intent(in)     :: elem
         927             :     logical                              :: inside, inside2
         928             :     integer               :: i,j,ii
         929             :     type (cartesian2D_t) :: corners(4),sphere_xy,cart
         930             :     type (cartesian3D_t) :: corners_xyz(4),center,a,b,cross(4)
         931             :     real (kind=r8) :: yp(4), y, elem_diam,dotprod
         932             :     real (kind=r8) :: xp(4), x
         933             :     real (kind=r8) :: d1,d2, tol_inside = 1.0e-12_r8
         934             :
         935             :     type (spherical_polar_t)   :: sphere  ! debug
         936             :
         937           0 :     inside = .false.
         938             :
         939             :     ! first check if point is near the element:
         940           0 :     corners_xyz(:) = elem%corners3D(:)
         941             :     elem_diam = max( distance(corners_xyz(1),corners_xyz(3)), &
         942           0 :          distance(corners_xyz(2),corners_xyz(4)) )
         943             :
         944           0 :     center%x = sum(corners_xyz(1:4)%x)/4
         945           0 :     center%y = sum(corners_xyz(1:4)%y)/4
         946           0 :     center%z = sum(corners_xyz(1:4)%z)/4
         947           0 :     if ( distance(center,sphere_xyz) > 1.0_r8*elem_diam ) return
         948             :
         949             :     j = 4
         950           0 :     do i=1,4
         951             :       ! outward normal to plane containing j->i edge:  corner(i) x corner(j)
         952             :       ! sphere dot (corner(i) x corner(j) ) = negative if inside
         953           0 :        cross(i)%x =  corners_xyz(i)%y*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%y
         954           0 :        cross(i)%y =-(corners_xyz(i)%x*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%x)
         955           0 :        cross(i)%z =  corners_xyz(i)%x*corners_xyz(j)%y - corners_xyz(i)%y*corners_xyz(j)%x
         956             :        dotprod = cross(i)%x*sphere_xyz%x + cross(i)%y*sphere_xyz%y +&
         957           0 :                cross(i)%z*sphere_xyz%z
         958           0 :        j = i  ! within this loopk j = i-1
         959             :
         960             :        !if (dotprod>0 .and. dotprod/elem_diam < 1e-5) print *,dotprod/elem_diam
         961             :
         962             :        ! dot product is proportional to elem_diam. positive means outside,
         963             :        ! but allow machine precision tolorence:
         964           0 :        if (dotprod > tol_inside*elem_diam) return
         965             :        !if (dotprod > 0) return
         966             :     end do
         967           0 :     inside=.true.
         968             :     return
         969             :   end function point_inside_gc
         970             :
         971             :
         972             :   !================================================
         973             :   !  (Nair) Cube face index and local coordinates
         974             :   !================================================
         975             :
         976           0 :   subroutine cube_facepoint_ne(sphere, ne, cart, number)
         977             :     use coordinate_systems_mod, only : cube_face_number_from_sphere, sphere2cubedsphere
         978             :
         979             :     type(spherical_polar_t), intent(in)  :: sphere
         980             :     integer,                 intent(in)  :: ne
         981             :     type(cartesian2D_t),     intent(out) :: cart
         982             :     integer,                 intent(out) :: number
         983             :
         984             :     real(kind=r8)       :: xp, yp
         985             :     type(cartesian2D_t) :: cube
         986             :     integer             :: ie, je, face_no
         987             :     real(kind=r8)       :: x1, x2
         988             :     real(kind=r8)       :: dx
         989             :
         990           0 :     face_no = cube_face_number_from_sphere(sphere)
         991           0 :     cube    = sphere2cubedsphere(sphere, face_no)
         992           0 :     xp      = cube%x
         993           0 :     yp      = cube%y
         994             :
         995             :     ! MNL: for uniform grids (on cube face), analytic solution is fine
         996           0 :     x1 = xp + 0.25_r8*PI
         997           0 :     x2 = yp + 0.25_r8*PI
         998             :
         999           0 :     dx = (0.5_r8*PI)/ne
        1000           0 :     ie = INT(ABS(x1)/dx)
        1001           0 :     je = INT(ABS(x2)/dx)
        1002             :     ! if we are exactly on an element edge, we can put the point in
        1003             :     ! either the ie or ie+1 element, EXCEPT if ie==ne.
        1004           0 :     if ( ABS(x1) < ne*dx  ) then
        1005           0 :       ie = ie + 1
        1006             :     end if
        1007           0 :     if ( ABS(x2) < ne*dx  ) then
        1008           0 :       je = je + 1
        1009             :     end if
        1010           0 :     if ((ie > ne) .or. (je > ne)) then
        1011           0 :       write(iulog, *) 'ERROR: ',ie,je,ne
        1012           0 :       write(iulog, *) 'lat,lon=',sphere%lat,sphere%lon
        1013           0 :       write(iulog, *) 'face no=',face_no
        1014           0 :       write(iulog, *) x1,x2,x1/dx,x2/dx
        1015           0 :       call endrun('interpolate_mod: bad argument')
        1016             :     endif
        1017             :
        1018             :     ! bug fix MT 1/2009.  This was creating a plotting error at
        1019             :     ! the row of elements in iface=2 at 50 degrees (NE=16 128x256 lat/lon grid)
        1020             :     ! For point on element edge, we can have ie=2, but x1=dx
        1021             :     ! but if ie>1, we must execute this statement.
        1022             :     ! The only time we can skip this statement is if ie=1, but then
        1023             :     ! the statement has no effect, so lets never skip it:
        1024             :     !    if (x1 > dx ) then
        1025           0 :     x1 = x1 - dble(ie-1)*dx
        1026             :     !    endif
        1027             :
        1028           0 :     x1 = 2.0_r8*(x1/dx)-1.0_r8
        1029             :
        1030             :     !    if (x2 > dx ) then    ! removed MT 1/2009, see above
        1031           0 :     x2 = x2 - dble(je-1)*dx
        1032             :     !    endif
        1033             :
        1034           0 :     x2 = 2.0_r8*(x2/dx)-1.0_r8
        1035             :
        1036             :     ! coordinates within an element [-1,1]
        1037           0 :     cart%x = x1
        1038           0 :     cart%y = x2
        1039           0 :     number = ie + (je-1)*ne + (face_no-1)*ne*ne
        1040           0 :   end subroutine cube_facepoint_ne
        1041             :   !================================================
        1042             :   !  (Nair) Cube face index and local coordinates
        1043             :   !================================================
        1044             :
        1045             :
        1046           0 :   subroutine cube_facepoint_unstructured(sphere,cart, number, elem)
        1047             :     use coordinate_systems_mod, only : cube_face_number_from_sphere, &
        1048             :                                        sphere2cubedsphere,change_coordinates,cube_face_number_from_cart
        1049             :     implicit none
        1050             :
        1051             :     type (element_t)     , intent(in), target :: elem(:)
        1052             :     type (spherical_polar_t), intent (in) :: sphere
        1053             :     type (cartesian2D_t), intent(out)     :: cart
        1054             :     integer             , intent(out)     :: number
        1055             :
        1056             :     integer               :: ii
        1057             :     Logical               :: found
        1058             :     type (cartesian3D_t)       :: sphere_xyz
        1059             :     type (cartesian2D_t)  :: cube
        1060           0 :     sphere_xyz=spherical_to_cart(sphere)
        1061             :
        1062           0 :     number=-1
        1063             : !    print *,'WARNING: using GC map'
        1064           0 :     do ii = 1,nelemd
        1065             :        ! for equiangular gnomonic map:
        1066             :        ! unstructed grid element edges are NOT great circles
        1067           0 :        if (cubed_sphere_map==0) then
        1068           0 :           found = point_inside_equiangular(elem(ii), sphere, sphere_xyz)
        1069             :        else
        1070             :           ! assume element edges are great circle arcs:
        1071           0 :           found = point_inside_gc(elem(ii), sphere_xyz)
        1072             :        endif
        1073             :
        1074           0 :        if (found) then
        1075           0 :           number = ii
        1076           0 :           cart = parametric_coordinates(sphere, elem(ii)%corners3D,&
        1077           0 :                cubed_sphere_map,elem(ii)%corners,elem(ii)%u2qmap,elem(ii)%facenum)
        1078           0 :           exit
        1079             :        end if
        1080             :     end do
        1081           0 :   end subroutine cube_facepoint_unstructured
        1082             :
        1083             :
        1084           0 :   subroutine interp_init()
        1085             :     type (quadrature_t)   :: gp
        1086             :
        1087           0 :     gp = gausslobatto(np)
        1088           0 :     call interpolate_create(gp,interp_p)
        1089           0 :   end subroutine interp_init
        1090             :
        1091             :
        1092           0 :   subroutine setup_latlon_interp(elem,interpdata,par)
        1093             :     !
        1094             :     ! initialize interpolation data structures to interpolate to a lat-lon grid
        1095             :     !
        1096             :     !
        1097             :
        1098             :     implicit none
        1099             :     type (element_t)     , intent(in), target :: elem(:)
        1100             :     type (parallel_t)      , intent(in)       :: par
        1101             :     type (interpdata_t)  , intent(out)        :: interpdata(:)
        1102             :
        1103             :     ! local
        1104             :     integer i,j,ii,count_total,n_interp,count_max
        1105             :     integer ngrid, number, elem_num, plat
        1106             :     integer countx, missing_pts,ierr
        1107             :     integer :: npts_mult_claims,max_claims
        1108             :
        1109           0 :     real (kind=r8)    ::  dp,latdeg(nlat+1),clat(nlat+1),w(nlat+1),w_staggered(nlat)
        1110           0 :     real (kind=r8)    ::  clat_staggered(nlat),latdeg_st(nlat),err,err2
        1111             :
        1112             :     type (spherical_polar_t) :: sphere
        1113             :     type (cartesian2D_t)     :: cart
        1114             :     type (cartesian3D_t)     :: sphere_xyz,sphere2_xyz
        1115             :
        1116             :     type (quadrature_t)       :: gp
        1117             :
        1118             :
        1119             :     ! Array to make sure each interp point is on exactly one process
        1120           0 :     type (cartesian2D_t),allocatable    :: cart_vec(:,:)
        1121             :     integer :: k
        1122           0 :     integer, allocatable :: global_elem_gid(:,:),local_elem_gid(:,:), local_elem_num(:,:)
        1123             :
        1124             :     ! these arrays often are too large for stack, so lets make sure
        1125             :     ! they go on the heap:
        1126           0 :     allocate(local_elem_num(nlat,nlon))
        1127           0 :     allocate(local_elem_gid(nlat,nlon))
        1128           0 :     allocate(global_elem_gid(nlat,nlon))
        1129           0 :     allocate(cart_vec(nlat,nlon))
        1130             :
        1131           0 :     if (par%masterproc) then
        1132           0 :        write(iulog,'(a,i4,a,i4,a)') 'Initializing ',nlat,' x ',nlon,' lat-lon interpolation grid: '
        1133             :     endif
        1134             :
        1135           0 :     do ii=1,nelemd
        1136           0 :        interpdata(ii)%n_interp=0  ! reset counter
        1137             :     enddo
        1138             :
        1139           0 :     if (associated(lat))then
        1140           0 :        deallocate(lat)
        1141             :        nullify(lat)
        1142             :     endif
        1143           0 :     if (associated(gweight))then
        1144           0 :        deallocate(gweight)
        1145             :        nullify(gweight)
        1146             :     endif
        1147             :
        1148           0 :     if (associated(lon))then
        1149           0 :        deallocate(lon)
        1150             :        nullify(lon)
        1151             :     endif
        1152             :
        1153           0 :     allocate(lat(nlat))
        1154           0 :     allocate(gweight(nlat))
        1155           0 :     allocate(lon(nlon))
        1156           0 :     call interp_init()
        1157           0 :     gweight=0
        1158           0 :     do i=1,nlon
        1159           0 :        lon(i)=2*pi*(i-1)/nlon
        1160             :     enddo
        1161           0 :     if (gridtype==1) then
        1162           0 :        do j=1,nlat
        1163           0 :           lat(j) = -pi/2 + pi*(j-1)/(nlat-1)
        1164             :        end do
        1165             :        plat=nlat
        1166             :     endif
        1167           0 :     if (gridtype==2) then
        1168           0 :        gp=gauss(nlat)
        1169           0 :        do j=1,nlat
        1170           0 :           lat(j) = asin(gp%points(j))
        1171           0 :           gweight(j) = gp%weights(j)
        1172             :        end do
        1173             :     endif
        1174           0 :     if (gridtype==3) then
        1175           0 :        do j=1,nlat
        1176           0 :           lat(j) = -pi/2 + pi*(j-.5_r8)/nlat
        1177             :        end do
        1178           0 :        plat=nlat+1
        1179             :     endif
        1180             :
        1181           0 :     if (gridtype==1 .or. gridtype==3) then
        1182             :        ! gridtype=1    plat=nlat    gweight(1:nlat)=w(1:plat)
        1183             :        ! gridtype=3    plat=nlat+1  gweight(1:nlat)=w_staggered(1:plat-1)
        1184             :
        1185             :        ! L-R dynamics uses a regular latitude distribution (not gausian).
        1186             :        ! The algorithm below is a bastardized version of LSM: map.F.
        1187           0 :        dp = 180.0_r8/(plat-1)
        1188           0 :        do j = 1, plat
        1189           0 :           latdeg(j) = -90.0_r8 + (j-1)*dp
        1190           0 :           clat(j) = latdeg(j)*pi/180.0_r8
        1191             :        end do
        1192             :
        1193             :        ! Calculate latitudes for the staggered grid
        1194             :
        1195           0 :        do j = 1, plat-1
        1196           0 :           clat_staggered(j) = (clat(j) + clat(j+1)) / 2
        1197           0 :           latdeg_st     (j) = clat_staggered(j)*180.0_r8/pi
        1198             :        end do
        1199             :
        1200             :        ! Weights are defined as cos(phi)*(delta-phi)
        1201             :        ! For a sanity check, the sum of w across all lats should be 2, or 1 across
        1202             :        ! half of the latitudes.
        1203             :
        1204           0 :        do j = 2, plat-1
        1205           0 :           w(j) = sin(clat_staggered(j)) - sin(clat_staggered(j-1))
        1206             :        end do
        1207           0 :        w(1) = sin(clat_staggered(1)) + 1
        1208           0 :        w(plat) = w(1)
        1209             :
        1210             :        ! with nlat=2048, this error was 4e-16
        1211           0 :        if (abs(sum(w(1:plat)) - 2) > 1.0e-8_r8) then
        1212           0 :           write(iulog,*) 'interpolate_mod: w weights do not sum to 2. sum=',sum(w(1:plat))
        1213           0 :           call endrun('interpolate_mod: weights do not sum to 2.')
        1214             :        end if
        1215             :
        1216           0 :        dp = pi / (plat-1)
        1217           0 :        do j = 1, plat-1
        1218           0 :           w_staggered(j) = sin(clat(j+1)) - sin(clat(j))
        1219             :        end do
        1220             :
        1221             :
        1222           0 :        if (abs(sum(w_staggered(1:plat-1)) - 2) > 1.0e-8_r8) then
        1223           0 :           write(iulog,*) 'interpolate_mod: staggered weights do not sum to 2. sum=',sum(w_staggered(1:plat-1))
        1224           0 :           call endrun('interpolate_mod: weights do not sum to 2.')
        1225             :        end if
        1226             :
        1227           0 :        if (gridtype==1) then
        1228           0 :           gweight(1:nlat)=w(1:plat)
        1229             :        endif
        1230           0 :        if (gridtype==3) then
        1231           0 :           gweight(1:nlat)=w_staggered(1:plat-1)
        1232             :        endif
        1233             :     endif
        1234             :
        1235             :
        1236             :     ! go through once, counting the number of points on each element
        1237           0 :     sphere%r=1
        1238           0 :     local_elem_num  = -1
        1239           0 :     local_elem_gid  = -1
        1240           0 :     global_elem_gid = -1
        1241           0 :     err=0
        1242           0 :     do j=1,nlat
        1243           0 :        do i=1,nlon
        1244           0 :           sphere%lat=lat(j)
        1245           0 :           sphere%lon=lon(i)
        1246             :
        1247           0 :           number = -1
        1248           0 :           if ( (cubed_sphere_map /= 0) .or. MeshUseMeshFile) then
        1249           0 :              call cube_facepoint_unstructured(sphere, cart, number, elem)
        1250           0 :              if (number /= -1) then
        1251             :                 ! If points are outside element but within tolerance, move to boundary
        1252           0 :                 if (cart%x + 1.0_r8.le.0.0_r8) cart%x = -1.0_r8
        1253           0 :                 if (cart%x - 1.0_r8.ge.0.0_r8) cart%x = 1.0_r8
        1254           0 :                 if (cart%y + 1.0_r8.le.0.0_r8) cart%y = -1.0_r8
        1255           0 :                 if (cart%y - 1.0_r8.ge.0.0_r8) cart%y = 1.0_r8
        1256             :
        1257           0 :                 local_elem_num(j,i) = number
        1258           0 :                 local_elem_gid(j,i) = elem(number)%vertex%number
        1259           0 :                 cart_vec(j,i)    = cart  ! local element coordiante of interpolation point
        1260             :              endif
        1261             :           else
        1262           0 :              call cube_facepoint_ne(sphere, ne, cart, number)
        1263             :              ! the sphere point belongs to the element number on face = face_no.
        1264             :              ! do I own this element?
        1265           0 :              if (number /= -1) then
        1266           0 :                 do ii=1,nelemd
        1267           0 :                    if (number == elem(ii)%vertex%number) then
        1268           0 :                       local_elem_gid(j,i) = number
        1269           0 :                       local_elem_num(j,i) = ii
        1270           0 :                       cart_vec(j,i)        = cart   ! local element coordinate found above
        1271           0 :                       exit
        1272             :                    endif
        1273             :                 enddo
        1274             :              endif
        1275             :           endif
        1276           0 :           ii=local_elem_num(j,i)
        1277           0 :           if (ii /= -1) then
        1278             :              ! compute error: map 'cart' back to sphere and compare with original
        1279             :              ! interpolation point:
        1280             :              sphere2_xyz = spherical_to_cart( ref2sphere(cart%x,cart%y,     &
        1281           0 :                   elem(ii)%corners3D,cubed_sphere_map,elem(ii)%corners,elem(ii)%facenum ))
        1282           0 :              sphere_xyz = spherical_to_cart(sphere)
        1283           0 :              err=max(err,distance(sphere2_xyz,sphere_xyz))
        1284             :           endif
        1285             :        enddo
        1286           0 :        if (par%masterproc) then
        1287           0 :           if ((MOD(j,64).eq.1).or.(j.eq.nlat)) then
        1288           0 :              print *,'finished latitude ',j,' of ',nlat
        1289             :           endif
        1290             :        endif
        1291             :     enddo
        1292           0 :     err2=err
        1293           0 :     call MPI_Allreduce(err, err2, 1, MPI_real8, MPI_MAX, par%comm, ierr)
        1294           0 :     if (par%masterproc) then
        1295           0 :        write(iulog,'(a,e12.4)') 'Max interpolation point search error: ',err2
        1296             :     endif
        1297             :
        1298             :     ! if multile elements claim a interpolation point, take the one with largest gid:
        1299           0 :     global_elem_gid = local_elem_gid
        1300           0 :     call MPI_Allreduce(local_elem_gid, global_elem_gid, nlat*nlon, MPI_integer, MPI_MAX, par%comm,ierr)
        1301             :
        1302           0 :     missing_pts=0
        1303           0 :     do j=1,nlat
        1304           0 :        do i=1,nlon
        1305           0 :           if (global_elem_gid(j,i) == -1 ) then
        1306           0 :              missing_pts = missing_pts + 1
        1307           0 :              if (par%masterproc) &
        1308           0 :                   print *,'Error: point not claimed by any element j,i,lat(j),lon(i)=',j,i,lat(j),lon(i)
        1309           0 :           else if (local_elem_gid(j,i) == global_elem_gid(j,i)  ) then
        1310           0 :              ii = local_elem_num(j,i)
        1311           0 :              interpdata(ii)%n_interp = interpdata(ii)%n_interp + 1
        1312             :           endif
        1313             :        end do
        1314             :     end do
        1315             :
        1316           0 :     countx=maxval(interpdata(1:nelemd)%n_interp)
        1317           0 :     count_max = countx
        1318           0 :     call MPI_Allreduce(countx,count_max,1,MPI_integer,MPI_MAX,par%comm,ierr)
        1319             :
        1320           0 :     if (par%masterproc) then
        1321           0 :        write(iulog,'(a,i6)') 'Maximum number of interpolation points claimed by an element: ',count_max
        1322             :     endif
        1323             :
        1324             :     ! allocate storage
        1325           0 :     do ii=1,nelemd
        1326           0 :        ngrid = interpdata(ii)%n_interp
        1327           0 :        if(interpdata(ii)%first_entry)then
        1328           0 :           NULLIFY(interpdata(ii)%interp_xy)
        1329           0 :           NULLIFY(interpdata(ii)%ilat)
        1330           0 :           NULLIFY(interpdata(ii)%ilon)
        1331             :
        1332           0 :           interpdata(ii)%first_entry=.FALSE.
        1333             :        endif
        1334           0 :        if(associated(interpdata(ii)%interp_xy))then
        1335           0 :           if(size(interpdata(ii)%interp_xy)>0)deallocate(interpdata(ii)%interp_xy)
        1336             :        endif
        1337           0 :        if(associated(interpdata(ii)%ilat))then
        1338           0 :           if(size(interpdata(ii)%ilat)>0)deallocate(interpdata(ii)%ilat)
        1339             :        endif
        1340             :
        1341           0 :        if (associated(interpdata(ii)%ilon))then
        1342           0 :           if(size(interpdata(ii)%ilon)>0)deallocate(interpdata(ii)%ilon)
        1343             :        endif
        1344           0 :        allocate(interpdata(ii)%interp_xy( ngrid ) )
        1345           0 :        allocate(interpdata(ii)%ilat( ngrid ) )
        1346           0 :        allocate(interpdata(ii)%ilon( ngrid ) )
        1347           0 :        interpdata(ii)%n_interp=0  ! reset counter
        1348             :     enddo
        1349           0 :     do j=1,nlat
        1350           0 :        do i=1,nlon
        1351           0 :           if (local_elem_gid(j,i) == global_elem_gid(j,i) .and. &
        1352           0 :                local_elem_gid(j,i) /= -1 ) then
        1353           0 :              ii = local_elem_num(j,i)
        1354           0 :              ngrid = interpdata(ii)%n_interp + 1
        1355           0 :              interpdata(ii)%n_interp = ngrid
        1356           0 :              interpdata(ii)%interp_xy( ngrid )   = cart_vec(j,i)
        1357           0 :              interpdata(ii)%ilon( ngrid ) = i
        1358           0 :              interpdata(ii)%ilat( ngrid ) = j
        1359             :           endif
        1360             :        enddo
        1361             :     enddo
        1362             :
        1363             :     ! now lets compute the number of points that were claimed by
        1364             :     ! more than one element:
        1365           0 :     do j=1,nlat
        1366           0 :        do i=1,nlon
        1367           0 :           if (local_elem_gid(j,i) == -1) then
        1368           0 :              local_elem_gid(j,i)=0
        1369             :           else
        1370           0 :              local_elem_gid(j,i)=1
        1371             :           endif
        1372             :        enddo
        1373             :     enddo
        1374           0 :     global_elem_gid = local_elem_gid
        1375           0 :     call MPI_Allreduce(local_elem_gid, global_elem_gid, nlat*nlon, MPI_integer, MPI_SUM, par%comm,ierr)
        1376           0 :     if (par%masterproc) then
        1377           0 :        countx=0
        1378           0 :        do j=1,nlat
        1379           0 :           do i=1,nlon
        1380           0 :              if (global_elem_gid(j,i)>1) countx=countx+1
        1381             :           enddo
        1382             :        enddo
        1383           0 :        npts_mult_claims=countx
        1384           0 :        max_claims=maxval(global_elem_gid)
        1385             :     endif
        1386             :
        1387           0 :     if (par%masterproc) then
        1388           0 :        print *,'Number of interpolation points claimed by more than one element: ',npts_mult_claims
        1389           0 :        print *,'max number of elements which claimed the same interpolation point:',max_claims
        1390             :     endif
        1391             :
        1392           0 :     deallocate(global_elem_gid)
        1393           0 :     deallocate(local_elem_num)
        1394           0 :     deallocate(local_elem_gid)
        1395           0 :     deallocate(cart_vec)
        1396             :
        1397             :     ! check if every point in interpolation grid was claimed by an element:
        1398           0 :     if (missing_pts>0) then
        1399           0 :        count_total = nlat*nlon
        1400           0 :        if(par%masterproc) then
        1401           0 :           write(iulog,"(3A,I4,A,I7,a,i5)")"Error:",__FILE__," ",__LINE__," count_total:",count_total," missing:",missing_pts
        1402             :        end if
        1403           0 :        call syncmp(par)
        1404           0 :        call endrun('Error: interpolation points not claimed by any element')
        1405             :     endif
        1406             :
        1407             :
        1408           0 :   end subroutine setup_latlon_interp
        1409             :
        1410             :
        1411             :
        1412             : ! interpolate_scalar
        1413             : !
        1414             : ! Interpolate a scalar field given in an element (fld_cube) to the points in
        1415             : ! interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
        1416             : !
        1417             : ! Note that it is possible the given element contains none of the interpolation points
        1418             : ! =======================================
        1419           0 : subroutine interpolate_ce(cart,fld_cube,npts,fld, fillvalue)
        1420             :   type (cartesian2D_t) :: cart
        1421             :   integer                  ::  npts
        1422             :   real (kind=r8)    ::  fld_cube(npts,npts) ! cube field
        1423             :   real (kind=r8)    ::  fld          ! field at new grid lat,lon coordinates
        1424             :   real (kind=r8), intent(in), optional :: fillvalue
        1425             :   ! Local variables
        1426             :   type (interpolate_t), pointer  ::  interp          ! interpolation structure
        1427             :
        1428             :   integer :: ne
        1429             :   integer :: i
        1430             :
        1431           0 :   if (npts==np) then
        1432             :      interp => interp_p
        1433             :   else
        1434           0 :      call endrun('Error in interpolate_scalar(): must be called with p or v grid data')
        1435             :   endif
        1436             :
        1437           0 :   fld=interpolate_2d(cart,fld_cube,interp,npts,fillvalue)
        1438             :
        1439           0 : end subroutine interpolate_ce
        1440             :
        1441             :
        1442             :
        1443             :   ! =======================================
        1444             :   ! interpolate_scalar
        1445             :   !
        1446             :   ! Interpolate a scalar field given in an element (fld_cube) to the points in
        1447             :   ! interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
        1448             :   !
        1449             :   ! Note that it is possible the given element contains none of the interpolation points
        1450             :   ! =======================================
        1451           0 :   subroutine interpolate_scalar2d(interpdata,fld_cube,nsize,nhalo,fld, fillvalue)
        1452             :     use dimensions_mod, only: npsq, fv_nphys,nc
        1453             :     integer,             intent(in) ::  nsize,nhalo
        1454             :     real (kind=r8),      intent(in) ::  fld_cube(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo) ! cube field
        1455             :     real (kind=r8),      intent(out)::  fld(:)          ! field at new grid lat,lon coordinates
        1456             :     type (interpdata_t), intent(in) ::  interpdata
        1457             :     real (kind=r8),      intent(in), optional :: fillvalue
        1458             :     ! Local variables
        1459             :     type (interpolate_t), pointer  ::  interp          ! interpolation structure
        1460             :
        1461             :     integer :: i,imin,imax,ne
        1462           0 :     real (kind=r8):: xoy(1-nhalo:nsize+nhalo),dx
        1463             :     type (cartesian2D_t) :: cart
        1464             :
        1465           0 :     if (nsize==np.and.nhalo==0) then
        1466             :       !
        1467             :       ! GLL grid
        1468             :       !
        1469           0 :       interp => interp_p
        1470           0 :       xoy = interp%glp(:)
        1471           0 :       imin = 1
        1472           0 :       imax = np
        1473           0 :     else if (nhalo>0.and.(nsize==fv_nphys.or.nsize==nc)) then
        1474             :       !
        1475             :       ! finite-volume grid
        1476             :       !
        1477           0 :       if (itype.ne.1) then
        1478           0 :         call endrun('itype must be 1 for latlon output from finite-volume (non-GLL) grids')
        1479             :       end if
        1480           0 :       imin = 1-nhalo
        1481           0 :       imax = nsize+nhalo
        1482             :       !
        1483             :       ! create normalized coordinates
        1484             :       !
        1485           0 :       dx      = 2.0_r8/REAL(nsize,KIND=r8)
        1486           0 :       do i=imin,imax
        1487           0 :         xoy(i) = -1.0_r8+(i-0.5_r8)*dx
        1488             :       end do
        1489             :     else
        1490           0 :       call endrun('interpolate_scalar2d: resolution not supported')
        1491             :     endif
        1492             :
        1493             :     ! Choice for Native (high-order) or Bilinear interpolations
        1494           0 :     if (itype == 0) then
        1495           0 :       do i=1,interpdata%n_interp
        1496           0 :         fld(i)=interpolate_2d(interpdata%interp_xy(i),fld_cube,interp,nsize,fillvalue)
        1497             :       end do
        1498           0 :     else if (itype == 1) then
        1499           0 :       do i=1,interpdata%n_interp
        1500           0 :         fld(i)=interpol_bilinear(interpdata%interp_xy(i),fld_cube,xoy,imin,imax,fillvalue)
        1501             :       end do
        1502             :     else
        1503           0 :       call endrun("interpolate_scalar2d: wrong interpolation type: "//int2str(itype))
        1504             :     end if
        1505             :
        1506           0 :   end subroutine interpolate_scalar2d
        1507           0 :   subroutine interpolate_scalar3d(interpdata,fld_cube,nsize,nhalo,nlev,fld, fillvalue)
        1508             :     use dimensions_mod, only: npsq, fv_nphys,nc
        1509             :     integer ,            intent(in)  ::  nsize, nhalo, nlev
        1510             :     real (kind=r8),      intent(in)  ::  fld_cube(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,nlev) ! cube field
        1511             :     real (kind=r8),      intent(out) ::  fld(:,:)          ! field at new grid lat,lon coordinates
        1512             :     type (interpdata_t), intent(in)  ::  interpdata
        1513             :     real (kind=r8), intent(in), optional :: fillvalue
        1514             :     ! Local variables
        1515             :     type (interpolate_t), pointer  ::  interp          ! interpolation structure
        1516             :
        1517             :     integer :: ne
        1518             :
        1519             :     integer :: i, k, imin, imax
        1520           0 :     real (kind=r8) :: xoy(1-nhalo:nsize+nhalo),dx
        1521             :
        1522             :     type (cartesian2D_t) :: cart
        1523             :
        1524           0 :     if (nsize==np.and.nhalo==0) then
        1525             :       !
        1526             :       ! GLL grid
        1527             :       !
        1528           0 :       interp => interp_p
        1529           0 :       xoy = interp%glp(:)
        1530           0 :       imin = 1
        1531           0 :       imax = np
        1532           0 :     else if (nhalo>0.and.(nsize==fv_nphys.or.nsize==nc)) then
        1533             :       !
        1534             :       ! finite-volume grid
        1535             :       !
        1536           0 :       if (itype.ne.1) then
        1537           0 :         call endrun('itype must be 1 for latlon output from finite-volume (non-GLL) grids')
        1538             :       end if
        1539           0 :       imin = 1-nhalo
        1540           0 :       imax = nsize+nhalo
        1541             :       !
        1542             :       ! create normalized coordinates
        1543             :       !
        1544           0 :       dx      = 2.0_r8/REAL(nsize,KIND=r8)
        1545           0 :       do i=imin,imax
        1546           0 :         xoy(i) = -1.0_r8+(i-0.5_r8)*dx
        1547             :       end do
        1548             :     else
        1549           0 :       call endrun('interpolate_scalar3d: resolution not supported')
        1550             :     endif
        1551             :
        1552             :     ! Choice for Native (high-order) or Bilinear interpolations
        1553           0 :     if (itype == 0) then
        1554           0 :       do k=1,nlev
        1555           0 :         do i=1,interpdata%n_interp
        1556           0 :           fld(i,k)=interpolate_2d(interpdata%interp_xy(i),fld_cube(:,:,k),interp,nsize,fillvalue)
        1557             :         end do
        1558             :       end do
        1559           0 :     elseif (itype == 1) then
        1560           0 :       do k=1,nlev
        1561           0 :         do i=1,interpdata%n_interp
        1562           0 :           fld(i,k)=interpol_bilinear(interpdata%interp_xy(i),fld_cube(:,:,k),xoy,imin,imax,fillvalue)
        1563             :         end do
        1564             :       end do
        1565             :     else
        1566           0 :       call endrun("interpolate_scalar3d: wrong interpolation type: "//int2str(itype))
        1567             :     endif
        1568           0 :   end subroutine interpolate_scalar3d
        1569             :
        1570             :
        1571             :   ! =======================================
        1572             :   ! interpolate_vector
        1573             :   !
        1574             :   ! Interpolate a vector field given in an element (fld_cube)
        1575             :   ! to the points in interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
        1576             :   !
        1577             :   ! input_coords = 0    fld_cube given in lat-lon
        1578             :   ! input_coords = 1    fld_cube given in contravariant
        1579             :   !
        1580             :   ! Note that it is possible the given element contains none of the interpolation points
        1581             :   ! =======================================
        1582           0 :   subroutine interpolate_vector2d(interpdata,elem,fld_cube,npts,fld,input_coords, fillvalue)
        1583             :     implicit none
        1584             :     integer                  ::  npts
        1585             :     real (kind=r8)    ::  fld_cube(npts,npts,2) ! vector field
        1586             :     real (kind=r8)    ::  fld(:,:)          ! field at new grid lat,lon coordinates
        1587             :     type (interpdata_t)      ::  interpdata
        1588             :     type (element_t), intent(in)         :: elem
        1589             :     real (kind=r8), intent(in), optional :: fillvalue
        1590             :     integer                  ::  input_coords
        1591             :
        1592             :
        1593             :     ! Local variables
        1594           0 :     real (kind=r8)    ::  fld_contra(npts,npts,2) ! vector field
        1595             :     type (interpolate_t), pointer  ::  interp          ! interpolation structure
        1596             :
        1597             :     real (kind=r8)    ::  v1,v2
        1598             :     real (kind=r8)    ::  D(2,2)   ! derivative of gnomonic mapping
        1599             :     real (kind=r8)    ::  JJ(2,2), tmpD(2,2)   ! derivative of gnomonic mapping
        1600             :
        1601             :     integer :: i,j
        1602             :
        1603             :     type (cartesian2D_t) :: cart
        1604             :
        1605           0 :     if(present(fillvalue)) then
        1606           0 :        if (any(fld_cube==fillvalue)) then
        1607           0 :           fld = fillvalue
        1608             :           return
        1609             :        end if
        1610             :     end if
        1611             :
        1612           0 :     if (input_coords==0 ) then
        1613             :        ! convert to contra
        1614           0 :        do j=1,npts
        1615           0 :           do i=1,npts
        1616             :              ! latlon->contra
        1617           0 :              fld_contra(i,j,1) = elem%Dinv(i,j,1,1)*fld_cube(i,j,1) + elem%Dinv(i,j,1,2)*fld_cube(i,j,2)
        1618           0 :              fld_contra(i,j,2) = elem%Dinv(i,j,2,1)*fld_cube(i,j,1) + elem%Dinv(i,j,2,2)*fld_cube(i,j,2)
        1619             :           enddo
        1620             :        enddo
        1621             :     else
        1622           0 :        fld_contra=fld_cube
        1623             :     endif
        1624             :
        1625             :
        1626           0 :     if (npts==np) then
        1627             :        interp => interp_p
        1628             :     else
        1629           0 :        call endrun('interpolate_vector2d: Error in interpolate_vector(): input must be on GLL grid')
        1630             :     endif
        1631             :
        1632             :
        1633             :        ! Choice for Native (high-order) or Bilinear interpolations
        1634             :
        1635           0 :     if (itype == 0) then
        1636           0 :        do i=1,interpdata%n_interp
        1637           0 :           fld(i,1)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,1),interp,npts)
        1638           0 :           fld(i,2)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,2),interp,npts)
        1639             :        end do
        1640           0 :     elseif (itype == 1) then
        1641           0 :        do i=1,interpdata%n_interp
        1642           0 :           fld(i,1)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,1),interp%glp(:),1,np)
        1643           0 :           fld(i,2)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,2),interp%glp(:),1,np)
        1644             :        end do
        1645             :     else
        1646           0 :        call endrun("interpolate_vector2d: wrong interpolation type: "//int2str(itype))
        1647             :     endif
        1648           0 :     do i=1,interpdata%n_interp
        1649             :        ! convert fld from contra->latlon
        1650           0 :        call dmap(D,interpdata%interp_xy(i)%x,interpdata%interp_xy(i)%y,&
        1651           0 :             elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
        1652             :        ! convert fld from contra->latlon
        1653           0 :        v1 = fld(i,1)
        1654           0 :        v2 = fld(i,2)
        1655             :
        1656           0 :        fld(i,1)=D(1,1)*v1 + D(1,2)*v2
        1657           0 :        fld(i,2)=D(2,1)*v1 + D(2,2)*v2
        1658             :     end do
        1659             :
        1660             :   end subroutine interpolate_vector2d
        1661             :
        1662             :   ! =======================================
        1663             :   ! interpolate_vector
        1664             :   !
        1665             :   ! Interpolate a vector field given in an element (fld_cube)
        1666             :   ! to the points in interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
        1667             :   !
        1668             :   ! input_coords = 0    fld_cube given in lat-lon
        1669             :   ! input_coords = 1    fld_cube given in contravariant
        1670             :   !
        1671             :   ! Note that it is possible the given element contains none of the interpolation points
        1672             :   ! =======================================
        1673           0 :   subroutine interpolate_vector3d(interpdata,elem,fld_cube,npts,nlev,fld,input_coords, fillvalue)
        1674             :     implicit none
        1675             :     type (interpdata_t),intent(in)       ::  interpdata
        1676             :     type (element_t), intent(in)         :: elem
        1677             :     integer, intent(in)                  ::  npts, nlev
        1678             :     real (kind=r8), intent(in)    ::  fld_cube(npts,npts,2,nlev) ! vector field
        1679             :     real (kind=r8), intent(out)   ::  fld(:,:,:)          ! field at new grid lat,lon coordinates
        1680             :     real (kind=r8), intent(in),optional :: fillvalue
        1681             :     integer, intent(in)                  ::  input_coords
        1682             :
        1683             :     ! Local variables
        1684           0 :     real (kind=r8)    ::  fld_contra(npts,npts,2,nlev) ! vector field
        1685             :     type (interpolate_t), pointer  ::  interp          ! interpolation structure
        1686             :
        1687             :     real (kind=r8)    ::  v1,v2
        1688             :     real (kind=r8)    ::  D(2,2)   ! derivative of gnomonic mapping
        1689             :     real (kind=r8)    ::  JJ(2,2), tmpD(2,2)   ! derivative of gnomonic mapping
        1690             :
        1691             :
        1692             :     integer :: i,j,k
        1693             :
        1694             :     type (cartesian2D_t) :: cart
        1695           0 :     if(present(fillvalue)) then
        1696           0 :        if (any(fld_cube==fillvalue)) then
        1697           0 :           fld = fillvalue
        1698             :           return
        1699             :        end if
        1700             :     end if
        1701           0 :     if (input_coords==0 ) then
        1702             :        ! convert to contra
        1703           0 :        do k=1,nlev
        1704           0 :           do j=1,npts
        1705           0 :              do i=1,npts
        1706             :                 ! latlon->contra
        1707           0 :                 fld_contra(i,j,1,k) = elem%Dinv(i,j,1,1)*fld_cube(i,j,1,k) + elem%Dinv(i,j,1,2)*fld_cube(i,j,2,k)
        1708           0 :                 fld_contra(i,j,2,k) = elem%Dinv(i,j,2,1)*fld_cube(i,j,1,k) + elem%Dinv(i,j,2,2)*fld_cube(i,j,2,k)
        1709             :              enddo
        1710             :           enddo
        1711             :        end do
        1712             :     else
        1713           0 :        fld_contra=fld_cube
        1714             :     endif
        1715             :
        1716           0 :     if (npts==np) then
        1717             :        interp => interp_p
        1718             :     else
        1719           0 :        call endrun('interpolate_vector3d: Error in interpolate_vector(): input must be on GLL grid')
        1720             :     endif
        1721             :
        1722             :
        1723             :        ! Choice for Native (high-order) or Bilinear interpolations
        1724             :
        1725           0 :     if (itype == 0) then
        1726           0 :        do k=1,nlev
        1727           0 :           do i=1,interpdata%n_interp
        1728           0 :              fld(i,k,1)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,1,k),interp,npts)
        1729           0 :              fld(i,k,2)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,2,k),interp,npts)
        1730             :           end do
        1731             :        end do
        1732           0 :     elseif (itype == 1) then
        1733           0 :        do k=1,nlev
        1734           0 :           do i=1,interpdata%n_interp
        1735           0 :              fld(i,k,1)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,1,k),interp%glp(:),1,np)
        1736           0 :              fld(i,k,2)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,2,k),interp%glp(:),1,np)
        1737             :           end do
        1738             :        end do
        1739             :     else
        1740           0 :        call endrun("interpolate_vector3d: wrong interpolation type: "//int2str(itype))
        1741             :     endif
        1742             :
        1743             :
        1744           0 :     do i=1,interpdata%n_interp
        1745             :        ! compute D(:,:) at the point elem%interp_cube(i)
        1746           0 :        call dmap(D,interpdata%interp_xy(i)%x,interpdata%interp_xy(i)%y,&
        1747           0 :             elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
        1748           0 :        do k=1,nlev
        1749             :           ! convert fld from contra->latlon
        1750           0 :           v1 = fld(i,k,1)
        1751           0 :           v2 = fld(i,k,2)
        1752             :
        1753           0 :           fld(i,k,1)=D(1,1)*v1 + D(1,2)*v2
        1754           0 :           fld(i,k,2)=D(2,1)*v1 + D(2,2)*v2
        1755             :        end do
        1756             :     end do
        1757             :   end subroutine interpolate_vector3d
        1758             :
        1759           0 :   subroutine vec_latlon_to_contra(elem,nphys,nhalo,nlev,fld,fvm)
        1760             :     use fvm_control_volume_mod, only: fvm_struct
        1761             :     use dimensions_mod,         only: fv_nphys
        1762             :     integer      , intent(in)   :: nphys,nhalo,nlev
        1763             :     real(kind=r8), intent(inout):: fld(1-nhalo:nphys+nhalo,1-nhalo:nphys+nhalo,2,nlev)
        1764             :     type (element_t), intent(in)           :: elem
        1765             :     type(fvm_struct), intent(in), optional :: fvm
        1766             :     !
        1767             :     ! local variables
        1768             :     !
        1769             :     integer :: i,j,k
        1770             :     real(r8):: v1,v2
        1771             :
        1772           0 :     if (nhalo==0.and.nphys==np) then
        1773           0 :       do k=1,nlev
        1774           0 :         do j=1,nphys
        1775           0 :           do i=1,nphys
        1776             :             ! latlon->contra
        1777           0 :             v1 = fld(i,j,1,k)
        1778           0 :             v2 = fld(i,j,2,k)
        1779           0 :             fld(i,j,1,k) = elem%Dinv(i,j,1,1)*v1 + elem%Dinv(i,j,1,2)*v2
        1780           0 :             fld(i,j,2,k) = elem%Dinv(i,j,2,1)*v1 + elem%Dinv(i,j,2,2)*v2
        1781             :           enddo
        1782             :         enddo
        1783             :       end do
        1784           0 :     else if (nphys==fv_nphys.and.nhalo.le.fv_nphys) then
        1785           0 :       do k=1,nlev
        1786           0 :         do j=1-nhalo,nphys+nhalo
        1787           0 :           do i=1-nhalo,nphys+nhalo
        1788             :             ! latlon->contra
        1789           0 :             v1 = fld(i,j,1,k)
        1790           0 :             v2 = fld(i,j,2,k)
        1791           0 :             fld(i,j,1,k) = fvm%Dinv_physgrid(i,j,1,1)*v1 + fvm%Dinv_physgrid(i,j,1,2)*v2
        1792           0 :             fld(i,j,2,k) = fvm%Dinv_physgrid(i,j,2,1)*v1 + fvm%Dinv_physgrid(i,j,2,2)*v2
        1793             :           enddo
        1794             :         enddo
        1795             :       end do
        1796             :     else
        1797           0 :       call endrun('ERROR: vec_latlon_to_contra - grid not supported or halo too large')
        1798             :     end if
        1799           0 :   end subroutine vec_latlon_to_contra
        1800           0 : end module interpolate_mod

Generated by: LCOV version 1.14