metgrid/src/remapper.F

   1 module remapper
   2
   3     use mpas_mesh
   4     use scan_input, only : input_field_type, FIELD_TYPE_REAL, FIELD_TYPE_DOUBLE, FIELD_TYPE_INTEGER
   5     use target_mesh
   6
   7     integer, parameter :: max_queue_length    = 2700    ! should be at least (earth circumference / minimum grid distance)
   8     integer, parameter :: max_dictionary_size = 82000   ! should be at least (nCells/32)
   9
  10     integer :: queue_head = 0
  11     integer :: queue_tail = 0
  12     integer :: queue_size = 0
  13     integer, dimension(0:max_queue_length-1) :: queue_array
  14
  15     integer :: int_size = 32
  16     integer, dimension(max_dictionary_size) :: dictionary_array
  17
  18
  19     type remap_info_type
  20         integer :: method = -1
  21         type (mpas_mesh_type), pointer :: src_mesh
  22         type (target_mesh_type), pointer :: dst_mesh
  23
  24         ! For nearest-neighbor
  25         integer, dimension(:,:), pointer :: nearestCell => null()
  26         integer, dimension(:,:), pointer :: nearestVertex => null()
  27         integer, dimension(:,:), pointer :: nearestEdge => null()
  28
  29         ! For Wachspress interpolation
  30         real, dimension(:,:,:), pointer :: cellWeights => null()
  31         real, dimension(:,:,:), pointer :: vertexWeights => null()
  32         real, dimension(:,:,:), pointer :: edgeWeights => null()
  33         integer, dimension(:,:,:), pointer :: sourceCells => null()
  34         integer, dimension(:,:,:), pointer :: sourceVertices => null()
  35         integer, dimension(:,:,:), pointer :: sourceEdges => null()
  36
  37         ! For masked interpolation
  38         real, dimension(:,:,:), pointer :: cellMaskedWeights => null()
  39         integer, dimension(:,:,:), pointer :: sourceMaskedCells => null()
  40     end type remap_info_type
  41
  42
  43     type target_field_type
  44         character (len=64) :: name
  45         integer :: ndims = -1
  46         integer :: xtype = -1
  47         logical :: isTimeDependent = .false.
  48         integer, dimension(:), pointer :: dimlens => null()
  49         character (len=64), dimension(:), pointer :: dimnames => null()
  50
  51         !  Members to store field data
  52         real :: array0r
  53         real, dimension(:), pointer :: array1r => null()
  54         real, dimension(:,:), pointer :: array2r => null()
  55         real, dimension(:,:,:), pointer :: array3r => null()
  56         real, dimension(:,:,:,:), pointer :: array4r => null()
  57         double precision :: array0d
  58         double precision, dimension(:), pointer :: array1d => null()
  59         double precision, dimension(:,:), pointer :: array2d => null()
  60         double precision, dimension(:,:,:), pointer :: array3d => null()
  61         double precision, dimension(:,:,:,:), pointer :: array4d => null()
  62         integer :: array0i
  63         integer, dimension(:), pointer :: array1i => null()
  64         integer, dimension(:,:), pointer :: array2i => null()
  65         integer, dimension(:,:,:), pointer :: array3i => null()
  66         integer, dimension(:,:,:,:), pointer :: array4i => null()
  67     end type target_field_type
  68
  69
  70     private :: nearest_cell, &
  71                nearest_vertex, &
  72                sphere_distance, &
  73                mpas_arc_length, &
  74                mpas_triangle_signed_area_sphere, &
  75                mpas_wachspress_coordinates, &
  76                convert_lx, &
  77                index2d
  78
  79
  80     contains
  81
  82
  83     integer function remap_info_setup(src_mesh, dst_mesh, remap_info) result(stat)
  84
  85         implicit none
  86
  87         type (mpas_mesh_type), intent(in), target :: src_mesh
  88         type (target_mesh_type), intent(in), target :: dst_mesh
  89         type (remap_info_type), intent(out) :: remap_info
  90
  91         integer :: idx
  92         integer :: j
  93         integer :: nn
  94         integer :: ix, iy
  95         integer :: irank
  96         real :: sumWeights
  97         real, dimension(:,:), allocatable :: vertCoords
  98         real, dimension(3) :: pointInterp
  99
 100         stat = 0
 101
 102         remap_info % method = 1
 103         remap_info % src_mesh => src_mesh
 104         remap_info % dst_mesh => dst_mesh
 105
 106         !
 107         ! For nearest neighbor
 108         !
 109         allocate(remap_info % nearestCell(dst_mesh % nlon, dst_mesh % nlat))
 110         allocate(remap_info % nearestVertex(dst_mesh % nlon, dst_mesh % nlat))
 111         allocate(remap_info % nearestEdge(dst_mesh % nlon, dst_mesh % nlat))
 112
 113         irank = dst_mesh % irank
 114
 115         idx = 1
 116         do iy=1,dst_mesh % nlat
 117         do ix=1,dst_mesh % nlon
 118             idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 119                                src_mesh % nCells, src_mesh % maxEdges, &
 120                                src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell)
 121             remap_info % nearestCell(ix, iy) = idx
 122         end do
 123         end do
 124
 125         idx = 1
 126         do iy=1,dst_mesh % nlat
 127         do ix=1,dst_mesh % nlon
 128             idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 129                                  src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, &
 130                                  src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, &
 131                                  src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, &
 132                                  src_mesh % latVertex, src_mesh % lonVertex )
 133             remap_info % nearestVertex(ix, iy) = idx
 134         end do
 135         end do
 136
 137         idx = 1
 138         do iy=1,dst_mesh % nlat
 139         do ix=1,dst_mesh % nlon
 140             idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 141                                  src_mesh % nCells, src_mesh % nEdges, src_mesh % maxEdges, 2, &
 142                                  src_mesh % nEdgesOnCell, src_mesh % edgesOnCell, &
 143                                  src_mesh % cellsOnEdge, src_mesh % latCell, src_mesh % lonCell, &
 144                                  src_mesh % latEdge, src_mesh % lonEdge )
 145             remap_info % nearestEdge(ix, iy) = idx
 146         end do
 147         end do
 148
 149
 150         !
 151         ! For Wachspress interpolation
 152         !
 153         allocate(vertCoords(3,3))
 154
 155         allocate(remap_info % cellWeights(3, dst_mesh % nlon, dst_mesh % nlat))
 156         allocate(remap_info % sourceCells(3, dst_mesh % nlon, dst_mesh % nlat))
 157
 158         idx = 1
 159         do iy=1,dst_mesh % nlat
 160         do ix=1,dst_mesh % nlon
 161             idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 162                                  src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, &
 163                                  src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, &
 164                                  src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, &
 165                                  src_mesh % latVertex, src_mesh % lonVertex )
 166             remap_info % sourceCells(:,ix,iy) = src_mesh % cellsOnVertex(:,idx)
 167             pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0)
 168             do j=1,3
 169                 vertCoords(:,j) = convert_lx(src_mesh % latCell(src_mesh % cellsOnVertex(j,idx)), &
 170                                              src_mesh % lonCell(src_mesh % cellsOnVertex(j,idx)), &
 171                                              6371229.0)
 172             end do
 173             remap_info % cellWeights(:,ix,iy) = mpas_wachspress_coordinates(3, vertCoords, pointInterp)
 174         end do
 175         end do
 176
 177         deallocate(vertCoords)
 178
 179
 180         allocate(vertCoords(3,src_mesh % maxEdges))
 181
 182         allocate(remap_info % vertexWeights(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat))
 183         allocate(remap_info % sourceVertices(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat))
 184
 185         idx = 1
 186         do iy=1,dst_mesh % nlat
 187         do ix=1,dst_mesh % nlon
 188             idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 189                                src_mesh % nCells, src_mesh % maxEdges, &
 190                                src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell)
 191             nn = src_mesh % nEdgesOnCell(idx)
 192             remap_info % sourceVertices(:,ix,iy) = 1
 193             remap_info % sourceVertices(1:nn,ix,iy) = src_mesh % verticesOnCell(1:nn,idx)
 194             pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0)
 195             do j=1,nn
 196                 vertCoords(:,j) = convert_lx(src_mesh % latVertex(src_mesh % verticesOnCell(j,idx)), &
 197                                              src_mesh % lonVertex(src_mesh % verticesOnCell(j,idx)), &
 198                                              6371229.0)
 199             end do
 200             remap_info % vertexWeights(:,ix,iy) = 0.0
 201             remap_info % vertexWeights(1:nn,ix,iy) = mpas_wachspress_coordinates(3, vertCoords(:,1:nn), pointInterp)
 202         end do
 203         end do
 204
 205         deallocate(vertCoords)
 206
 207
 208         allocate(vertCoords(3,src_mesh % maxEdges))
 209
 210         allocate(remap_info % edgeWeights(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat))
 211         allocate(remap_info % sourceEdges(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat))
 212
 213         idx = 1
 214         do iy=1,dst_mesh % nlat
 215         do ix=1,dst_mesh % nlon
 216             idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 217                                src_mesh % nCells, src_mesh % maxEdges, &
 218                                src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell)
 219             nn = src_mesh % nEdgesOnCell(idx)
 220             remap_info % sourceEdges(:,ix,iy) = 1
 221             remap_info % sourceEdges(1:nn,ix,iy) = src_mesh % edgesOnCell(1:nn,idx)
 222             pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0)
 223             do j=1,nn
 224                 vertCoords(:,j) = convert_lx(src_mesh % latEdge(src_mesh % edgesOnCell(j,idx)), &
 225                                              src_mesh % lonEdge(src_mesh % edgesOnCell(j,idx)), &
 226                                              6371229.0)
 227             end do
 228             remap_info % edgeWeights(:,ix,iy) = 0.0
 229             remap_info % edgeWeights(1:nn,ix,iy) = mpas_wachspress_coordinates(3, vertCoords(:,1:nn), pointInterp)
 230         end do
 231         end do
 232
 233         deallocate(vertCoords)
 234
 235
 236         !
 237         ! For masked interpolation
 238         !
 239         allocate(vertCoords(3,3))
 240
 241         allocate(remap_info % cellMaskedWeights(3, dst_mesh % nlon, dst_mesh % nlat))
 242         allocate(remap_info % sourceMaskedCells(3, dst_mesh % nlon, dst_mesh % nlat))
 243
 244         idx = 1
 245         do iy=1,dst_mesh % nlat
 246         do ix=1,dst_mesh % nlon
 247             idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, &
 248                                  src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, &
 249                                  src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, &
 250                                  src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, &
 251                                  src_mesh % latVertex, src_mesh % lonVertex )
 252             remap_info % sourceMaskedCells(:,ix,iy) = src_mesh % cellsOnVertex(:,idx)
 253             pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0)
 254             do j=1,3
 255                 vertCoords(:,j) = convert_lx(src_mesh % latCell(src_mesh % cellsOnVertex(j,idx)), &
 256                                              src_mesh % lonCell(src_mesh % cellsOnVertex(j,idx)), &
 257                                              6371229.0)
 258             end do
 259             remap_info % cellMaskedWeights(:,ix,iy) = mpas_wachspress_coordinates(3, vertCoords, pointInterp)
 260             sumWeights = 0.0
 261             do j=1,3
 262                 if (src_mesh % landmask(remap_info % sourceMaskedCells(j,ix,iy)) == 1) then
 263                    sumWeights = sumWeights + remap_info % cellMaskedWeights(j,ix,iy)
 264                 else
 265                    remap_info % cellMaskedWeights(j,ix,iy) = 0.0
 266                 end if
 267             end do
 268             if (sumWeights > 0.0) then
 269                 remap_info % cellMaskedWeights(:,ix,iy) = remap_info % cellMaskedWeights(:,ix,iy) / sumWeights
 270             else
 271                 idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), &
 272                                    remap_info % sourceMaskedCells(1,ix,iy), &
 273                                    src_mesh % nCells, src_mesh % maxEdges, &
 274                                    src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell)
 275                 call search_for_cells(src_mesh % nCells, src_mesh % maxEdges, src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, &
 276                                       src_mesh % landmask, dst_mesh % lats(index2d(irank,ix),iy), &
 277                                       dst_mesh % lons(ix,index2d(irank,iy)), &
 278                                       src_mesh % latCell, src_mesh % lonCell, idx, remap_info % sourceMaskedCells(:,ix,iy), &
 279                                       remap_info % cellMaskedWeights(:,ix,iy))
 280             end if
 281         end do
 282         end do
 283
 284         deallocate(vertCoords)
 285
 286     end function remap_info_setup
 287
 288
 289     integer function remap_info_free(remap_info) result(stat)
 290
 291         implicit none
 292
 293         type (remap_info_type), intent(inout) :: remap_info
 294
 295
 296         stat = 0
 297
 298         remap_info % method = -1
 299         nullify(remap_info % src_mesh)
 300         nullify(remap_info % dst_mesh)
 301
 302         if (associated(remap_info % nearestCell)) then
 303             deallocate(remap_info % nearestCell)
 304         end if
 305         if (associated(remap_info % nearestVertex)) then
 306             deallocate(remap_info % nearestVertex)
 307         end if
 308         if (associated(remap_info % nearestEdge)) then
 309             deallocate(remap_info % nearestEdge)
 310         end if
 311         if (associated(remap_info % cellWeights)) then
 312             deallocate(remap_info % cellWeights)
 313         end if
 314         if (associated(remap_info % vertexWeights)) then
 315             deallocate(remap_info % vertexWeights)
 316         end if
 317         if (associated(remap_info % edgeWeights)) then
 318             deallocate(remap_info % edgeWeights)
 319         end if
 320         if (associated(remap_info % sourceCells)) then
 321             deallocate(remap_info % sourceCells)
 322         end if
 323         if (associated(remap_info % sourceVertices)) then
 324             deallocate(remap_info % sourceVertices)
 325         end if
 326         if (associated(remap_info % sourceEdges)) then
 327             deallocate(remap_info % sourceEdges)
 328         end if
 329         if (associated(remap_info % cellMaskedWeights)) then
 330             deallocate(remap_info % cellMaskedWeights)
 331         end if
 332         if (associated(remap_info % sourceMaskedCells)) then
 333             deallocate(remap_info % sourceMaskedCells)
 334         end if
 335
 336     end function remap_info_free
 337
 338
 339     logical function can_remap_field(field)
 340
 341         implicit none
 342
 343         type (input_field_type), intent(in) :: field
 344
 345         integer :: decompDim
 346
 347
 348         can_remap_field = .true.
 349
 350         if (field % xtype /= FIELD_TYPE_INTEGER .and. &
 351             field % xtype /= FIELD_TYPE_REAL .and. &
 352             field % xtype /= FIELD_TYPE_DOUBLE) then
 353             can_remap_field = .false.
 354             return
 355         end if
 356
 357         if (field % ndims == 0 .or. &
 358             (field % ndims == 1 .and. field % isTimeDependent)) then
 359             can_remap_field = .false.
 360             return
 361         end if
 362
 363         decompDim = field % ndims
 364         if (field % isTimeDependent) then
 365             decompDim = decompDim - 1
 366         end if
 367
 368         if (trim(field % dimnames(decompDim)) /= 'nCells' .and. &
 369             trim(field % dimnames(decompDim)) /= 'nVertices' .and. &
 370             trim(field % dimnames(decompDim)) /= 'nEdges') then
 371
 372             can_remap_field = .false.
 373             return
 374         end if
 375
 376     end function can_remap_field
 377
 378
 379     integer function remap_field_dryrun(remap_info, src_field, dst_field) result(stat)
 380
 381         implicit none
 382
 383         type (remap_info_type), intent(in) :: remap_info
 384         type (input_field_type), intent(in) :: src_field
 385         type (target_field_type), intent(out) :: dst_field
 386
 387         integer :: idim
 388
 389         stat = 0
 390
 391         dst_field % name = src_field % name
 392         dst_field % xtype = src_field % xtype
 393         if (src_field % isTimeDependent) then
 394             ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon,
 395             ! but the time dimension is not counted in the target field
 396             dst_field % ndims = src_field % ndims
 397         else
 398             ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon
 399             dst_field % ndims = src_field % ndims + 1
 400         end if
 401         allocate(dst_field % dimnames(dst_field % ndims))
 402         allocate(dst_field % dimlens(dst_field % ndims))
 403         dst_field % isTimeDependent = src_field % isTimeDependent
 404         do idim=1,dst_field % ndims-2
 405             dst_field % dimlens(idim) = src_field % dimlens(idim)
 406             dst_field % dimnames(idim) = src_field % dimnames(idim)
 407         end do
 408
 409         dst_field % dimlens(dst_field % ndims-1) = remap_info % dst_mesh % nlon
 410         dst_field % dimnames(dst_field % ndims-1) = 'nLons'
 411         dst_field % dimlens(dst_field % ndims) = remap_info % dst_mesh % nlat
 412         dst_field % dimnames(dst_field % ndims) = 'nLats'
 413
 414     end function remap_field_dryrun
 415
 416
 417     integer function remap_field(remap_info, src_field, dst_field, masked) result(stat)
 418
 419         implicit none
 420
 421         type (remap_info_type), intent(in) :: remap_info
 422         type (input_field_type), intent(in) :: src_field
 423         type (target_field_type), intent(out) :: dst_field
 424         logical, intent(in), optional :: masked
 425
 426         integer :: decompDim
 427         integer :: idim
 428         integer :: j
 429         integer :: iy, ix
 430         logical :: local_masked
 431         integer, dimension(:,:), pointer :: nearestIndex
 432         integer, dimension(:,:,:), pointer :: sourceNodes
 433         real, dimension(:,:,:), pointer :: nodeWeights
 434
 435         stat = 0
 436
 437         decompDim = src_field % ndims
 438         if (src_field % isTimeDependent) then
 439             decompDim = decompDim - 1
 440         end if
 441
 442         local_masked = .false.
 443         if (present(masked)) then
 444             local_masked = masked
 445         end if
 446
 447         if (trim(src_field % dimnames(decompDim)) == 'nCells') then
 448             nearestIndex => remap_info % nearestCell
 449             if (local_masked) then
 450                 sourceNodes => remap_info % sourceMaskedCells
 451                 nodeWeights => remap_info % cellMaskedWeights
 452             else
 453                 sourceNodes => remap_info % sourceCells
 454                 nodeWeights => remap_info % cellWeights
 455             end if
 456         else if (trim(src_field % dimnames(decompDim)) == 'nVertices') then
 457             nearestIndex => remap_info % nearestVertex
 458             sourceNodes => remap_info % sourceVertices
 459             nodeWeights => remap_info % vertexWeights
 460         else if (trim(src_field % dimnames(decompDim)) == 'nEdges') then
 461             nearestIndex => remap_info % nearestEdge
 462             sourceNodes => remap_info % sourceEdges
 463             nodeWeights => remap_info % edgeWeights
 464         else
 465             write(0,*) 'Remap exception: unhandled decomposed dim'
 466             stat = 1
 467             return
 468         end if
 469
 470         dst_field % name = src_field % name
 471         dst_field % xtype = src_field % xtype
 472         if (src_field % isTimeDependent) then
 473             ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon,
 474             ! but the time dimension is not counted in the target field
 475             dst_field % ndims = src_field % ndims
 476         else
 477             ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon
 478             dst_field % ndims = src_field % ndims + 1
 479         end if
 480         allocate(dst_field % dimnames(dst_field % ndims))
 481         allocate(dst_field % dimlens(dst_field % ndims))
 482         dst_field % isTimeDependent = src_field % isTimeDependent
 483         do idim=1,dst_field % ndims-2
 484             dst_field % dimlens(idim) = src_field % dimlens(idim)
 485             dst_field % dimnames(idim) = src_field % dimnames(idim)
 486         end do
 487         dst_field % dimlens(dst_field % ndims-1) = remap_info % dst_mesh % nlon
 488         dst_field % dimnames(dst_field % ndims-1) = 'nLons'
 489         dst_field % dimlens(dst_field % ndims) = remap_info % dst_mesh % nlat
 490         dst_field % dimnames(dst_field % ndims) = 'nLats'
 491
 492         if (src_field % xtype == FIELD_TYPE_REAL) then
 493             if (dst_field % ndims == 2) then
 494                 allocate(dst_field % array2r(dst_field % dimlens(1), &
 495                                              dst_field % dimlens(2)))
 496 #ifdef NEAREST_NEIGHBOR
 497                 do iy=1,size(nearestIndex, dim=2)
 498                 do ix=1,size(nearestIndex, dim=1)
 499                     dst_field % array2r(ix,iy) = src_field % array1r(nearestIndex(ix,iy))
 500                 end do
 501                 end do
 502 #else
 503                 do iy=1,size(nearestIndex, dim=2)
 504                 do ix=1,size(nearestIndex, dim=1)
 505                     dst_field % array2r(ix,iy) = 0.0
 506                     do j=1,size(sourceNodes, dim=1)
 507                         dst_field % array2r(ix,iy) = dst_field % array2r(ix,iy) + &
 508                                                       nodeWeights(j,ix,iy) &
 509                                                       * src_field % array1r(sourceNodes(j,ix,iy))
 510                     end do
 511                 end do
 512                 end do
 513 #endif
 514             else if (dst_field % ndims == 3) then
 515                 allocate(dst_field % array3r(dst_field % dimlens(1), &
 516                                              dst_field % dimlens(2), &
 517                                              dst_field % dimlens(3)))
 518 #ifdef NEAREST_NEIGHBOR
 519                 do iy=1,size(nearestIndex, dim=2)
 520                 do ix=1,size(nearestIndex, dim=1)
 521                     dst_field % array3r(:,ix,iy) = src_field % array2r(:,nearestIndex(ix,iy))
 522                 end do
 523                 end do
 524 #else
 525                 do iy=1,size(nearestIndex, dim=2)
 526                 do ix=1,size(nearestIndex, dim=1)
 527                     dst_field % array3r(:,ix,iy) = 0.0
 528                     do j=1,3
 529                         dst_field % array3r(:,ix,iy) = dst_field % array3r(:,ix,iy) + &
 530                                                         nodeWeights(j,ix,iy) &
 531                                                         * src_field % array2r(:,sourceNodes(j,ix,iy))
 532                     end do
 533                 end do
 534                 end do
 535 #endif
 536             else if (dst_field % ndims == 4) then
 537                 allocate(dst_field % array4r(dst_field % dimlens(1), &
 538                                              dst_field % dimlens(2), &
 539                                              dst_field % dimlens(3), &
 540                                              dst_field % dimlens(4)))
 541 #ifdef NEAREST_NEIGHBOR
 542                 do iy=1,size(nearestIndex, dim=2)
 543                 do ix=1,size(nearestIndex, dim=1)
 544                     dst_field % array4r(:,:,ix,iy) = src_field % array3r(:,:,nearestIndex(ix,iy))
 545                 end do
 546                 end do
 547 #endif
 548                 do iy=1,size(nearestIndex, dim=2)
 549                 do ix=1,size(nearestIndex, dim=1)
 550                     dst_field % array4r(:,:,ix,iy) = 0.0
 551                     do j=1,3
 552                         dst_field % array4r(:,:,ix,iy) = dst_field % array4r(:,:,ix,iy) + &
 553                                                           remap_info % cellWeights(j,ix,iy) &
 554                                                           * src_field % array3r(:,:,remap_info % sourceCells(j,ix,iy))
 555                     end do
 556                 end do
 557                 end do
 558             else
 559                 write(0,*) 'Remap exception: unhandled dimension for real ', dst_field % ndims
 560             end if
 561         else if (src_field % xtype == FIELD_TYPE_DOUBLE) then
 562             if (dst_field % ndims == 2) then
 563                 allocate(dst_field % array2d(dst_field % dimlens(1), &
 564                                              dst_field % dimlens(2)))
 565 #ifdef NEAREST_NEIGHBOR
 566                 do iy=1,size(nearestIndex, dim=2)
 567                 do ix=1,size(nearestIndex, dim=1)
 568                     dst_field % array2d(ix,iy) = src_field % array1d(nearestIndex(ix,iy))
 569                 end do
 570                 end do
 571 #endif
 572                 do iy=1,size(nearestIndex, dim=2)
 573                 do ix=1,size(nearestIndex, dim=1)
 574                     dst_field % array2d(ix,iy) = 0.0
 575                     do j=1,size(sourceNodes, dim=1)
 576                         dst_field % array2d(ix,iy) = dst_field % array2d(ix,iy) + &
 577                                                       nodeWeights(j,ix,iy) &
 578                                                       * src_field % array1d(sourceNodes(j,ix,iy))
 579                     end do
 580                 end do
 581                 end do
 582             else if (dst_field % ndims == 3) then
 583                 allocate(dst_field % array3d(dst_field % dimlens(1), &
 584                                              dst_field % dimlens(2), &
 585                                              dst_field % dimlens(3)))
 586 #ifdef NEAREST_NEIGHBOR
 587                 do iy=1,size(nearestIndex, dim=2)
 588                 do ix=1,size(nearestIndex, dim=1)
 589                     dst_field % array3d(:,ix,iy) = src_field % array2d(:,nearestIndex(ix,iy))
 590                 end do
 591                 end do
 592 #endif
 593                 do iy=1,size(nearestIndex, dim=2)
 594                 do ix=1,size(nearestIndex, dim=1)
 595                     dst_field % array3d(:,ix,iy) = 0.0
 596                     do j=1,3
 597                         dst_field % array3d(:,ix,iy) = dst_field % array3d(:,ix,iy) + &
 598                                                         remap_info % cellWeights(j,ix,iy) &
 599                                                         * src_field % array2d(:,remap_info % sourceCells(j,ix,iy))
 600                     end do
 601                 end do
 602                 end do
 603             else if (dst_field % ndims == 4) then
 604                 allocate(dst_field % array4d(dst_field % dimlens(1), &
 605                                              dst_field % dimlens(2), &
 606                                              dst_field % dimlens(3), &
 607                                              dst_field % dimlens(4)))
 608 #ifdef NEAREST_NEIGHBOR
 609                 do iy=1,size(nearestIndex, dim=2)
 610                 do ix=1,size(nearestIndex, dim=1)
 611                     dst_field % array4d(:,:,ix,iy) = src_field % array3d(:,:,nearestIndex(ix,iy))
 612                 end do
 613                 end do
 614 #endif
 615                 do iy=1,size(nearestIndex, dim=2)
 616                 do ix=1,size(nearestIndex, dim=1)
 617                     dst_field % array4d(:,:,ix,iy) = 0.0
 618                     do j=1,3
 619                         dst_field % array4d(:,:,ix,iy) = dst_field % array4d(:,:,ix,iy) + &
 620                                                           remap_info % cellWeights(j,ix,iy) &
 621                                                           * src_field % array3d(:,:,remap_info % sourceCells(j,ix,iy))
 622                     end do
 623                 end do
 624                 end do
 625             else
 626                 write(0,*) 'Remap exception: unhandled dimension for dbl ', dst_field % ndims
 627             end if
 628         else if (src_field % xtype == FIELD_TYPE_INTEGER) then
 629             if (dst_field % ndims == 2) then
 630                 allocate(dst_field % array2i(dst_field % dimlens(1), &
 631                                              dst_field % dimlens(2)))
 632                 do iy=1,size(nearestIndex, dim=2)
 633                 do ix=1,size(nearestIndex, dim=1)
 634                     dst_field % array2i(ix,iy) = src_field % array1i(nearestIndex(ix,iy))
 635                 end do
 636                 end do
 637             else if (dst_field % ndims == 3) then
 638                 allocate(dst_field % array3i(dst_field % dimlens(1), &
 639                                              dst_field % dimlens(2), &
 640                                              dst_field % dimlens(3)))
 641                 do iy=1,size(nearestIndex, dim=2)
 642                 do ix=1,size(nearestIndex, dim=1)
 643                     dst_field % array3i(:,ix,iy) = src_field % array2i(:,nearestIndex(ix,iy))
 644                 end do
 645                 end do
 646             else if (dst_field % ndims == 4) then
 647                 allocate(dst_field % array4i(dst_field % dimlens(1), &
 648                                              dst_field % dimlens(2), &
 649                                              dst_field % dimlens(3), &
 650                                              dst_field % dimlens(4)))
 651                 do iy=1,size(nearestIndex, dim=2)
 652                 do ix=1,size(nearestIndex, dim=1)
 653                     dst_field % array4i(:,:,ix,iy) = src_field % array3i(:,:,nearestIndex(ix,iy))
 654                 end do
 655                 end do
 656             else
 657                 write(0,*) 'Remap exception: unhandled dimension for int ', dst_field % ndims
 658             end if
 659         else
 660             write(0,*) 'Remap exception: unhandled type'
 661         end if
 662
 663     end function remap_field
 664
 665
 666     integer function remap_get_target_latitudes(remap_info, lat_field) result(stat)
 667
 668         implicit none
 669
 670         type (remap_info_type), intent(in) :: remap_info
 671         type (target_field_type), intent(out) :: lat_field
 672
 673         real, parameter :: rad2deg = 90.0 / asin(1.0)
 674
 675         stat = 0
 676
 677
 678         lat_field % name = 'lat'
 679         lat_field % xtype = FIELD_TYPE_REAL
 680         lat_field % ndims = 1
 681         lat_field % isTimeDependent = .false.
 682
 683         allocate(lat_field % dimnames(lat_field % ndims))
 684         allocate(lat_field % dimlens(lat_field % ndims))
 685
 686         lat_field % dimnames(1) = 'nLats'
 687         lat_field % dimlens(1) = remap_info % dst_mesh % nlat
 688
 689         allocate(lat_field % array1r(lat_field % dimlens(1)))
 690         lat_field % array1r(:) = remap_info % dst_mesh % lats(1,:) * rad2deg
 691
 692     end function remap_get_target_latitudes
 693
 694
 695     integer function remap_get_target_longitudes(remap_info, lon_field) result(stat)
 696
 697         implicit none
 698
 699         type (remap_info_type), intent(in) :: remap_info
 700         type (target_field_type), intent(out) :: lon_field
 701
 702         real, parameter :: rad2deg = 90.0 / asin(1.0)
 703
 704         stat = 0
 705
 706
 707         lon_field % name = 'lon'
 708         lon_field % xtype = FIELD_TYPE_REAL
 709         lon_field % ndims = 1
 710         lon_field % isTimeDependent = .false.
 711
 712         allocate(lon_field % dimnames(lon_field % ndims))
 713         allocate(lon_field % dimlens(lon_field % ndims))
 714
 715         lon_field % dimnames(1) = 'nLons'
 716         lon_field % dimlens(1) = remap_info % dst_mesh % nlon
 717
 718         allocate(lon_field % array1r(lon_field % dimlens(1)))
 719         lon_field % array1r(:) = remap_info % dst_mesh % lons(:,1) * rad2deg
 720
 721     end function remap_get_target_longitudes
 722
 723
 724     integer function free_target_field(field) result(stat)
 725
 726         implicit none
 727
 728         type (target_field_type), intent(inout) :: field
 729
 730         stat = 0
 731
 732         if (associated(field % dimlens)) then
 733             deallocate(field % dimlens)
 734         end if
 735         if (associated(field % dimnames)) then
 736             deallocate(field % dimnames)
 737         end if
 738
 739         if (associated(field % array1r)) then
 740             deallocate(field % array1r)
 741         end if
 742         if (associated(field % array2r)) then
 743             deallocate(field % array2r)
 744         end if
 745         if (associated(field % array3r)) then
 746             deallocate(field % array3r)
 747         end if
 748         if (associated(field % array4r)) then
 749             deallocate(field % array4r)
 750         end if
 751
 752         if (associated(field % array1d)) then
 753             deallocate(field % array1d)
 754         end if
 755         if (associated(field % array2d)) then
 756             deallocate(field % array2d)
 757         end if
 758         if (associated(field % array3d)) then
 759             deallocate(field % array3d)
 760         end if
 761         if (associated(field % array4d)) then
 762             deallocate(field % array4d)
 763         end if
 764
 765         if (associated(field % array1i)) then
 766             deallocate(field % array1i)
 767         end if
 768         if (associated(field % array2i)) then
 769             deallocate(field % array2i)
 770         end if
 771         if (associated(field % array3i)) then
 772             deallocate(field % array3i)
 773         end if
 774         if (associated(field % array4i)) then
 775             deallocate(field % array4i)
 776         end if
 777
 778     end function free_target_field
 779
 780
 781     integer function nearest_cell(target_lat, target_lon, start_cell, nCells, maxEdges, &
 782                                   nEdgesOnCell, cellsOnCell, latCell, lonCell)
 783
 784         implicit none
 785
 786         real, intent(in) :: target_lat, target_lon
 787         integer, intent(in) :: start_cell
 788         integer, intent(in) :: nCells, maxEdges
 789         integer, dimension(nCells), intent(in) :: nEdgesOnCell
 790         integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell
 791         real, dimension(nCells), intent(in) :: latCell, lonCell
 792
 793         integer :: i
 794         integer :: iCell
 795         integer :: current_cell
 796         real :: current_distance, d
 797         real :: nearest_distance
 798
 799         nearest_cell = start_cell
 800         current_cell = -1
 801
 802         do while (nearest_cell /= current_cell)
 803             current_cell = nearest_cell
 804             current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, &
 805                                                target_lon, 1.0)
 806             nearest_cell = current_cell
 807             nearest_distance = current_distance
 808             do i = 1, nEdgesOnCell(current_cell)
 809                 iCell = cellsOnCell(i,current_cell)
 810                 if (iCell <= nCells) then
 811                     d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0)
 812                     if (d < nearest_distance) then
 813                         nearest_cell = iCell
 814                         nearest_distance = d
 815                     end if
 816                 end if
 817             end do
 818         end do
 819
 820     end function nearest_cell
 821
 822
 823     integer function nearest_vertex( target_lat, target_lon, &
 824                                      start_vertex, &
 825                                      nCells, nVertices, maxEdges, vertexDegree, &
 826                                      nEdgesOnCell, verticesOnCell, &
 827                                      cellsOnVertex, latCell, lonCell, &
 828                                      latVertex, lonVertex )
 829
 830         implicit none
 831
 832         real, intent(in) :: target_lat, target_lon
 833         integer, intent(in) :: start_vertex
 834         integer, intent(in) :: nCells, nVertices, maxEdges, vertexDegree
 835         integer, dimension(nCells), intent(in) :: nEdgesOnCell
 836         integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell
 837         integer, dimension(vertexDegree,nVertices), intent(in) :: cellsOnVertex
 838         real, dimension(nCells), intent(in) :: latCell, lonCell
 839         real, dimension(nVertices), intent(in) :: latVertex, lonVertex
 840
 841
 842         integer :: i, cell1, cell2, cell3, iCell
 843         integer :: iVtx
 844         integer :: current_vertex
 845         real :: cell1_dist, cell2_dist, cell3_dist
 846         real :: current_distance, d
 847         real :: nearest_distance
 848
 849         nearest_vertex = start_vertex
 850         current_vertex = -1
 851
 852         do while (nearest_vertex /= current_vertex)
 853             current_vertex = nearest_vertex
 854             current_distance = sphere_distance(latVertex(current_vertex), lonVertex(current_vertex), &
 855                                                target_lat,                target_lon,                1.0)
 856             nearest_vertex = current_vertex
 857             nearest_distance = current_distance
 858             cell1 = cellsOnVertex(1,current_vertex)
 859             cell1_dist = sphere_distance(latCell(cell1), lonCell(cell1), target_lat, target_lon, 1.0)
 860             cell2 = cellsOnVertex(2,current_vertex)
 861             cell2_dist = sphere_distance(latCell(cell2), lonCell(cell2), target_lat, target_lon, 1.0)
 862             if (vertexDegree == 3) then
 863                cell3 = cellsOnVertex(3,current_vertex)
 864                cell3_dist = sphere_distance(latCell(cell3), lonCell(cell3), target_lat, target_lon, 1.0)
 865             end if
 866             if (vertexDegree == 3) then
 867                 if (cell1_dist < cell2_dist) then
 868                     if (cell1_dist < cell3_dist) then
 869                         iCell = cell1
 870                     else
 871                         iCell = cell3
 872                     end if
 873                 else
 874                     if (cell2_dist < cell3_dist) then
 875                         iCell = cell2
 876                     else
 877                         iCell = cell3
 878                     end if
 879                 end if
 880             else
 881                 if (cell1_dist < cell2_dist) then
 882                     iCell = cell1
 883                 else
 884                     iCell = cell2
 885                 end if
 886             end if
 887             do i = 1, nEdgesOnCell(iCell)
 888                 iVtx = verticesOnCell(i,iCell)
 889                 d = sphere_distance(latVertex(iVtx), lonVertex(iVtx), target_lat, target_lon, 1.0)
 890                 if (d < nearest_distance) then
 891                     nearest_vertex = iVtx
 892                     nearest_distance = d
 893                 end if
 894             end do
 895         end do
 896
 897     end function nearest_vertex
 898
 899
 900     real function sphere_distance(lat1, lon1, lat2, lon2, radius)
 901     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 902     ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
 903     !   sphere with given radius.
 904     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 905
 906         implicit none
 907
 908         real, intent(in) :: lat1, lon1, lat2, lon2, radius
 909
 910         real :: arg1
 911
 912         arg1 = sqrt( sin(0.5*(lat2-lat1))**2 +  &
 913                      cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
 914         sphere_distance = 2.*radius*asin(arg1)
 915
 916     end function sphere_distance
 917
 918
 919     !***********************************************************************
 920     !
 921     !  function mpas_wachspress_coordinates
 922     !
 923     !> \brief Compute the barycentric Wachspress coordinates for a polygon
 924     !> \author  Phillip Wolfram
 925     !> \date    01/26/2015
 926     !> \details
 927     !>  Computes the barycentric Wachspress coordinates for a polygon with nVertices
 928     !>  points in R3, vertCoords for a particular pointInterp with normalized radius.
 929     !>  Follows Gillette, A., Rand, A., Bajaj, C., 2011.
 930     !>  Error estimates for generalized barycentric interpolation.
 931     !>  Advances in computational mathematics 37 (3), 417–439.
 932     !>  Optimized version of mpas_wachspress_coordinates uses optional cached B_i areas
 933     !------------------------------------------------------------------------
 934     function mpas_wachspress_coordinates(nVertices, vertCoords, pointInterp, areaBin)
 935
 936         implicit none
 937
 938         ! input points
 939         integer, intent(in) :: nVertices
 940         real, dimension(3, nVertices), intent(in) :: vertCoords
 941         real, dimension(3), intent(in) :: pointInterp
 942         real, dimension(nVertices), optional, intent(in) :: areaBin
 943
 944         ! output
 945         real, dimension(nVertices) :: mpas_wachspress_coordinates
 946
 947         ! computational intermediates
 948         real, dimension(nVertices) :: wach       ! The wachpress area-product
 949         real :: wach_total                       ! The wachpress total weight
 950         integer :: i, j                                       ! Loop indices
 951         integer :: im1, i0, ip1                               ! im1 = (i-1), i0 = i, ip1 = (i+1)
 952
 953         ! triangle areas to compute wachspress coordinate
 954         real, dimension(nVertices) :: areaA
 955         real, dimension(nVertices) :: areaB
 956
 957         real :: radiusLocal
 958
 959         radiusLocal = sqrt(sum(vertCoords(:,1)**2))
 960
 961         if (.not. present(areaBin)) then
 962             ! compute areas
 963             do i = 1, nVertices
 964                 ! compute first area B_i
 965                 ! get vertex indices
 966                 im1 = mod(nVertices + i - 2, nVertices) + 1
 967                 i0  = mod(nVertices + i - 1, nVertices) + 1
 968                 ip1 = mod(nVertices + i    , nVertices) + 1
 969
 970                 ! precompute B_i areas
 971                 ! always the same because B_i independent of xp,yp,zp
 972                 ! (COULD CACHE AND USE RESULT FROM ARRAY FOR FURTHER OPTIMIZATION)
 973                 areaB(i) = mpas_triangle_signed_area_sphere(vertCoords(:, im1), vertCoords(:, i0), vertCoords(:, ip1), radiusLocal)
 974             end do
 975         else
 976             ! assign areas
 977             do i = 1, nVertices
 978                 areaB(i) = areaBin(i)
 979             end do
 980         end if
 981
 982         ! compute areas
 983         do i = 1, nVertices
 984             ! compute first area B_i
 985             ! get vertex indices
 986             im1 = mod(nVertices + i - 2, nVertices) + 1
 987             i0  = mod(nVertices + i - 1, nVertices) + 1
 988             ip1 = mod(nVertices + i    , nVertices) + 1
 989
 990             ! compute A_ij areas
 991             ! must be computed each time
 992             areaA(i0) = mpas_triangle_signed_area_sphere(pointInterp, vertCoords(:, i0), vertCoords(:, ip1), radiusLocal)
 993
 994             ! precomputed B_i areas, cached
 995         end do
 996
 997
 998         ! for each vertex compute wachpress coordinate
 999         do i = 1, nVertices
1000             wach(i) = areaB(i)
1001             do j = (i + 1), (i + nVertices - 2)
1002                 i0  = mod(nVertices + j - 1, nVertices) + 1
1003                 ! accumulate products for A_ij subareas
1004                 wach(i) = wach(i) * areaA(i0)
1005             end do
1006         end do
1007
1008         ! get summed weights for normalization
1009         wach_total = 0
1010         do i = 1, nVertices
1011             wach_total = wach_total + wach(i)
1012         end do
1013
1014         ! compute lambda
1015         mpas_wachspress_coordinates= 0.0
1016         do i = 1, nVertices
1017             mpas_wachspress_coordinates(i) = wach(i)/wach_total
1018         end do
1019
1020     end function mpas_wachspress_coordinates
1021
1022
1023     !***********************************************************************
1024     !
1025     !  routine mpas_triangle_signed_area_sphere
1026     !
1027     !> \brief   Calculates area of a triangle on a sphere
1028     !> \author  Matthew Hoffman
1029     !> \date    13 January 2015
1030     !> \details
1031     !>  This routine calculates the area of a triangle on the surface of a sphere.
1032     !>  Uses the spherical analog of Heron's formula.
1033     !>  Copied from mesh generator.  A CCW winding angle is positive.
1034     !-----------------------------------------------------------------------
1035     real function mpas_triangle_signed_area_sphere(a, b, c, radius)
1036
1037         implicit none
1038
1039         !-----------------------------------------------------------------
1040         ! input variables
1041         !-----------------------------------------------------------------
1042         real, dimension(3), intent(in) :: a, b, c  !< Input: 3d (x,y,z) points forming the triangle in which to calculate the bary weights
1043         real, intent(in) :: radius  !< sphere radius
1044
1045         !-----------------------------------------------------------------
1046         ! local variables
1047         !-----------------------------------------------------------------
1048         real :: ab, bc, ca, semiperim, tanqe
1049         real, dimension(3) :: ablen, aclen, Dlen
1050
1051         ab = mpas_arc_length(a(1), a(2), a(3), b(1), b(2), b(3))/radius
1052         bc = mpas_arc_length(b(1), b(2), b(3), c(1), c(2), c(3))/radius
1053         ca = mpas_arc_length(c(1), c(2), c(3), a(1), a(2), a(3))/radius
1054         semiperim = 0.5 * (ab + bc + ca)
1055
1056         tanqe = sqrt(max(0.0,tan(0.5 * semiperim) * tan(0.5 * (semiperim - ab)) &
1057                      * tan(0.5 * (semiperim - bc)) * tan(0.5 * (semiperim - ca))))
1058
1059         mpas_triangle_signed_area_sphere = 4.0 * radius * radius * atan(tanqe)
1060
1061         ! computing correct signs (in similar fashion to mpas_sphere_angle)
1062         ablen(1) = b(1) - a(1)
1063         ablen(2) = b(2) - a(2)
1064         ablen(3) = b(3) - a(3)
1065
1066         aclen(1) = c(1) - a(1)
1067         aclen(2) = c(2) - a(2)
1068         aclen(3) = c(3) - a(3)
1069
1070         dlen(1) =   (ablen(2) * aclen(3)) - (ablen(3) * aclen(2))
1071         dlen(2) = -((ablen(1) * aclen(3)) - (ablen(3) * aclen(1)))
1072         dlen(3) =   (ablen(1) * aclen(2)) - (ablen(2) * aclen(1))
1073
1074         if ((Dlen(1)*a(1) + Dlen(2)*a(2) + Dlen(3)*a(3)) < 0.0) then
1075             mpas_triangle_signed_area_sphere = -mpas_triangle_signed_area_sphere
1076         end if
1077
1078     end function mpas_triangle_signed_area_sphere
1079
1080
1081     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1082     ! FUNCTION MPAS_ARC_LENGTH
1083     !
1084     ! Returns the length of the great circle arc from A=(ax, ay, az) to
1085     !    B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
1086     !    same sphere centered at the origin.
1087     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1088     real function mpas_arc_length(ax, ay, az, bx, by, bz)
1089
1090         implicit none
1091
1092         real, intent(in) :: ax, ay, az, bx, by, bz
1093
1094         real :: r, c
1095         real :: cx, cy, cz
1096
1097         cx = bx - ax
1098         cy = by - ay
1099         cz = bz - az
1100
1101         r = sqrt(ax*ax + ay*ay + az*az)
1102         c = sqrt(cx*cx + cy*cy + cz*cz)
1103
1104         mpas_arc_length = r * 2.0 * asin(c/(2.0*r))
1105
1106     end function mpas_arc_length
1107
1108
1109     function convert_lx(lat, lon, radius) result(vec)
1110
1111         implicit none
1112
1113         real, intent(in) :: lat, lon, radius
1114
1115         real, dimension(3) :: vec
1116
1117         vec(1) = radius * cos(lon) * cos(lat)
1118         vec(2) = radius * sin(lon) * cos(lat)
1119         vec(3) = radius * sin(lat)
1120
1121     end function convert_lx
1122
1123
1124     subroutine search_for_cells(nCells, maxEdges, nEdgesOnCell, cellsOnCell, cellMask, &
1125                                 targetLat, targetLon, latCell, lonCell, &
1126                                 startIdx, sourceMaskedCells, cellMaskedWeights)
1127
1128         implicit none
1129
1130         integer, intent(in) :: nCells
1131         integer, intent(in) :: maxEdges
1132         integer, dimension(nCells), intent(in) :: nEdgesOnCell
1133         integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell
1134         integer, dimension(nCells), intent(in) :: cellMask
1135         real, intent(in) :: targetLat
1136         real, intent(in) :: targetLon
1137         real, dimension(nCells), intent(in) :: latCell
1138         real, dimension(nCells), intent(in) :: lonCell
1139         integer, intent(in) :: startIdx
1140         integer, dimension(3), intent(inout) :: sourceMaskedCells
1141         real, dimension(3), intent(inout) :: cellMaskedWeights
1142
1143         integer :: i
1144         integer :: scan_cell
1145         integer :: neighbor
1146         integer :: unscanned
1147         logical :: no_more_queueing
1148         real :: d, d_min
1149
1150         !
1151         ! Reset data structures
1152         !
1153         call queue_reset()
1154         call dictionary_reset()
1155
1156         no_more_queueing = .false.
1157         unscanned = 0
1158
1159         d_min = 1000.0
1160
1161         sourceMaskedCells(:) = 1
1162         cellMaskedWeights(:) = 0.0
1163
1164         !
1165         ! Insert the origin cell into the queue for processing
1166         !
1167         call queue_insert(startIdx)
1168         call dictionary_insert(startIdx)
1169
1170         do while (queue_size > 0)
1171             scan_cell = queue_remove()
1172
1173             !
1174             ! Each cell index removed from the queue represents a unique grid cell
1175             !    that falls within the specified radius of the origin point
1176             ! Here, we can do any processing we like for these cells
1177             !
1178             if (cellMask(scan_cell) == 1) then
1179                d = sphere_distance(targetLat, targetLon, latCell(scan_cell), lonCell(scan_cell), 1.0)
1180                if (d < d_min) then
1181                    d_min = d
1182                    sourceMaskedCells(1) = scan_cell
1183                    cellMaskedWeights(1) = 1.0
1184                    if (.not. no_more_queueing) then
1185                       unscanned = queue_size
1186                    end if
1187                    no_more_queueing = .true.
1188                end if
1189             end if
1190
1191             !
1192             ! Add any neighbors of scan_cell within specified radius to the queue for processing
1193             !
1194             if (.not. no_more_queueing .or. unscanned > 0) then
1195                 do i=1,nEdgesOnCell(scan_cell)
1196                     neighbor = cellsOnCell(i,scan_cell)
1197                     if (.not. dictionary_search(neighbor)) then
1198                         call queue_insert(neighbor)
1199                         call dictionary_insert(neighbor)
1200                     end if
1201                 end do
1202             end if
1203
1204             if (no_more_queueing) then
1205                unscanned = unscanned - 1
1206             end if
1207
1208         end do
1209
1210     end subroutine search_for_cells
1211
1212
1213     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1214     ! Insert a new integer into the queue
1215     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1216     subroutine queue_insert(i)
1217
1218         implicit none
1219
1220         integer, intent(in) :: i
1221
1222         if (queue_size == max_queue_length) then
1223             write(0,*) 'Error: queue overrun'
1224             return
1225         end if
1226         queue_size = queue_size + 1
1227         queue_array(queue_head) = i
1228         queue_head = mod(queue_head + 1, max_queue_length)
1229
1230     end subroutine queue_insert
1231
1232
1233     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1234     ! Remove the oldest integer from the queue
1235     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1236     integer function queue_remove()
1237
1238         implicit none
1239
1240         if (queue_size <= 0) then
1241             write(0,*) 'Error: queue underrun'
1242             queue_remove = -1
1243             return
1244         end if
1245         queue_size = queue_size - 1
1246         queue_remove = queue_array(queue_tail)
1247         queue_tail = mod(queue_tail + 1, max_queue_length)
1248
1249     end function queue_remove
1250
1251
1252     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1253     ! Reset the queue to an empty state
1254     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1255     subroutine queue_reset()
1256
1257         implicit none
1258
1259         queue_head = 0
1260         queue_tail = 0
1261         queue_size = 0
1262
1263     end subroutine queue_reset
1264
1265
1266     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1267     ! Insert an integer into the dictionary
1268     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1269     subroutine dictionary_insert(i)
1270
1271         implicit none
1272
1273         integer, intent(in) :: i
1274
1275         integer :: n_integer
1276         integer :: n_bit
1277
1278         n_integer = ((i-1) / int_size) + 1
1279         n_bit = mod((i-1), int_size)
1280
1281         if (n_integer > max_dictionary_size) then
1282             write(0,*) 'Error: dictionary insert out of bounds'
1283             return
1284         end if
1285
1286         dictionary_array(n_integer) = ibset(dictionary_array(n_integer), n_bit)
1287
1288     end subroutine dictionary_insert
1289
1290
1291     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1292     ! Search for an integer in the dictionary
1293     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1294     logical function dictionary_search(i)
1295
1296         implicit none
1297
1298         integer, intent(in) :: i
1299
1300         integer :: n_integer
1301         integer :: n_bit
1302
1303         n_integer = ((i-1) / int_size) + 1
1304         n_bit = mod((i-1), int_size)
1305
1306         if (n_integer > max_dictionary_size) then
1307             write(0,*) 'Error: dictionary search out of bounds'
1308             dictionary_search = .false.
1309             return
1310         end if
1311
1312         dictionary_search = btest(dictionary_array(n_integer), n_bit)
1313
1314     end function dictionary_search
1315
1316
1317     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1318     ! Reset the dictionary to an empty state
1319     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1320     subroutine dictionary_reset()
1321
1322         implicit none
1323
1324         dictionary_array(:) = 0
1325
1326     end subroutine dictionary_reset
1327
1328
1329     function index2d(irank, idx) result(i)
1330
1331         implicit none
1332
1333         integer, intent(in) :: irank, idx
1334
1335         integer :: i
1336
1337         i = irank * (idx - 1) + 1
1338
1339     end function index2d
1340
1341 end module remapper