LCOV - code coverage report
Current view: | top level - dynamics/se - interp_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 278 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module interp_mod
           2             :   use shr_kind_mod,        only: r8 => shr_kind_r8, r4 => shr_kind_r4
           3             :   use dimensions_mod,      only: nelemd, np, ne
           4             :   use interpolate_mod,     only: interpdata_t
           5             :   use interpolate_mod,     only: interp_lat => lat, interp_lon => lon
           6             :   use interpolate_mod,     only: interp_gweight => gweight
           7             :   use dyn_grid,            only: elem,fvm
           8             :   use spmd_utils,          only: iam
           9             :   use cam_history_support, only: fillvalue
          10             :   use hybrid_mod,          only: hybrid_t, config_thread_region
          11             :   use cam_abortutils,      only: endrun
          12             :
          13             :   implicit none
          14             :   private
          15             :   save
          16             :
          17             :   public :: setup_history_interpolation
          18             :   public :: set_interp_hfile
          19             :   public :: write_interpolated
          20             :
          21             :   interface write_interpolated
          22             :      module procedure write_interpolated_scalar
          23             :      module procedure write_interpolated_vector
          24             :   end interface
          25             :
          26             :   ! hybrid is created in setup_history_interpolation
          27             :   type(hybrid_t) :: hybrid
          28             :
          29             : ! A type to hold interpdata info for each interpolated history file
          30             : type cam_interpolate_t
          31             :   type(interpdata_t), pointer :: interpdata(:) => NULL()
          32             : end type cam_interpolate_t
          33             :
          34             : type(cam_interpolate_t), pointer :: interpdata_set(:) => NULL()  ! all files
          35             : type(interpdata_t),      pointer :: cam_interpolate(:) => NULL() ! curr. file
          36             :
          37             : CONTAINS
          38             :
          39           0 :   subroutine setup_history_interpolation(interp_ok, mtapes, interp_output,    &
          40           0 :        interp_info)
          41             :
          42             :     use cam_history_support, only: interp_info_t
          43             :     use cam_history_support, only: interp_type_native
          44             :     use cam_history_support, only: interp_type_bilinear
          45             :     use cam_history_support, only: interp_gridtype_equal_poles
          46             :     use cam_history_support, only: interp_gridtype_gauss
          47             :     use cam_history_support, only: interp_gridtype_equal_nopoles
          48             :     use cam_grid_support,    only: horiz_coord_t, horiz_coord_create, iMap
          49             :     use cam_grid_support,    only: cam_grid_register, cam_grid_attribute_register
          50             :     use cam_grid_support,    only: max_hcoordname_len
          51             :     use interpolate_mod,     only: get_interp_lat, get_interp_lon
          52             :     use interpolate_mod,     only: get_interp_parameter, set_interp_parameter
          53             :     use interpolate_mod,     only: get_interp_gweight, setup_latlon_interp
          54             :     use parallel_mod,        only: par
          55             :
          56             :     ! Dummy arguments
          57             :     logical,             intent(inout) :: interp_ok
          58             :     integer,             intent(in)    :: mtapes
          59             :     logical,             intent(in)    :: interp_output(:)
          60             :     type(interp_info_t), intent(inout) :: interp_info(:)
          61             :
          62             :     ! Local variables
          63             :     integer                            :: i
          64           0 :     real(r8),            pointer       :: w(:)
          65           0 :     integer(iMap),       pointer       :: grid_map(:,:)
          66             :     type(horiz_coord_t), pointer       :: lat_coord
          67             :     type(horiz_coord_t), pointer       :: lon_coord
          68             :     character(len=max_hcoordname_len)  :: gridname
          69             :     !---------------------------------------------------------------------------
          70             :
          71           0 :     nullify(grid_map)
          72             :
          73             :     ! For this dycore, interpolated output should be OK
          74           0 :     interp_ok = (iam < par%nprocs)
          75             :
          76           0 :     if (interp_ok) then
          77           0 :       hybrid = config_thread_region(par,'serial')
          78             :
          79           0 :       if(any(interp_output)) then
          80           0 :         allocate(interpdata_set(mtapes))
          81           0 :         do i = 1, mtapes
          82           0 :           if (interp_output(i)) then
          83           0 :             if ( (interp_info(i)%interp_nlon == 0) .or.                       &
          84             :                  (interp_info(i)%interp_nlat == 0)) then
          85             :               ! compute interpolation grid based on number of points around equator
          86           0 :               call set_interp_parameter('auto', (4 * ne * (np-1)))
          87           0 :               interp_info(i)%interp_nlat = get_interp_parameter('nlat')
          88           0 :               interp_info(i)%interp_nlon = get_interp_parameter('nlon')
          89             :             else
          90             :               call set_interp_parameter('nlat', interp_info(i)%interp_nlat)
          91           0 :               call set_interp_parameter('nlon', interp_info(i)%interp_nlon)
          92             :             end if
          93           0 :             call set_interp_parameter('itype', interp_info(i)%interp_type)
          94           0 :             call set_interp_parameter('gridtype', interp_info(i)%interp_gridtype)
          95             :
          96           0 :             allocate(interpdata_set(i)%interpdata(nelemd))
          97             :             ! Reset pointers in the interpolate module so they are not
          98             :             ! overwritten
          99           0 :             nullify(interp_lat)
         100           0 :             nullify(interp_lon)
         101           0 :             nullify(interp_gweight)
         102             :             call setup_latlon_interp(elem, interpdata_set(i)%interpdata, par)
         103             :             ! Create the grid coordinates
         104             :             lat_coord => horiz_coord_create('lat', '',                        &
         105           0 :                  interp_info(i)%interp_nlat, 'latitude', 'degrees_north',     &
         106           0 :                  1, interp_info(i)%interp_nlat, get_interp_lat())
         107             :             lon_coord => horiz_coord_create('lon', '',                        &
         108           0 :                  interp_info(i)%interp_nlon, 'longitude', 'degrees_east',     &
         109           0 :                  1, interp_info(i)%interp_nlon, get_interp_lon())
         110             :             ! Create a grid for this history file
         111           0 :             write(gridname, '(a,i0)') 'interp_out_', i
         112           0 :             interp_info(i)%grid_id = 200 + i
         113           0 :             call cam_grid_register(trim(gridname), interp_info(i)%grid_id,    &
         114           0 :                  lat_coord, lon_coord, grid_map, unstruct=.false.)
         115           0 :             interp_info(i)%gridname = trim(gridname)
         116             :             ! Add grid attributes
         117           0 :             allocate(w(get_interp_parameter('nlat')))
         118           0 :             w = get_interp_gweight()
         119           0 :             select case(interp_info(i)%interp_gridtype)
         120             :             case(interp_gridtype_equal_poles)
         121             :               call cam_grid_attribute_register(trim(gridname),                &
         122           0 :                    'interp_outputgridtype', 'equally spaced with poles')
         123             :               call cam_grid_attribute_register(trim(gridname), 'w',           &
         124           0 :                    'latitude weights', 'lat', w)
         125             :             case(interp_gridtype_equal_nopoles)
         126             :               call cam_grid_attribute_register(trim(gridname),                &
         127           0 :                    'interp_outputgridtype', 'equally spaced no poles')
         128             :               call cam_grid_attribute_register(trim(gridname), 'gw',          &
         129           0 :                    'latitude weights', 'lat', w)
         130             :             case(interp_gridtype_gauss)
         131             :               call cam_grid_attribute_register(trim(gridname),                &
         132           0 :                    'interp_outputgridtype', 'Gauss')
         133             :               call cam_grid_attribute_register(trim(gridname), 'gw',          &
         134           0 :                    'gauss weights', 'lat', w)
         135             :             case default
         136             :               call cam_grid_attribute_register(trim(gridname),                &
         137             :                    'interp_outputgridtype',                                   &
         138             :                    'Unknown interpolation output grid type',                  &
         139           0 :                    interp_info(i)%interp_gridtype)
         140             :             end select
         141           0 :             nullify(w) ! belongs to attribute
         142           0 :             if(interp_info(i)%interp_type == interp_type_native) then
         143             :               call cam_grid_attribute_register(trim(gridname),                &
         144           0 :                    'interp_type', 'se basis functions')
         145           0 :             else if(interp_info(i)%interp_type == interp_type_bilinear) then
         146             :               call cam_grid_attribute_register(trim(gridname),                &
         147           0 :                    'interp_type', 'bilinear')
         148             :             else
         149             :               call cam_grid_attribute_register(trim(gridname), 'interp_type', &
         150           0 :                    'Unknown interpolation type', interp_info(i)%interp_type)
         151             :             end if
         152             :             ! Store the data pointers for reuse later
         153           0 :             interp_info(i)%interp_lat => interp_lat
         154           0 :             interp_info(i)%interp_lon => interp_lon
         155           0 :             interp_info(i)%interp_gweight => interp_gweight
         156             :           end if
         157             :         end do
         158             :       end if
         159             :     end if
         160             :
         161           0 :   end subroutine setup_history_interpolation
         162             :
         163           0 :   subroutine set_interp_hfile(hfilenum, interp_info)
         164           0 :     use cam_history_support, only: interp_info_t
         165             :     use interpolate_mod,     only: set_interp_parameter
         166             :
         167             :     ! Dummy arguments
         168             :     integer,             intent(in)    :: hfilenum
         169             :     type(interp_info_t), intent(inout) :: interp_info(:)
         170             :
         171           0 :     if (.not. associated(interpdata_set)) then
         172           0 :       call endrun('SET_INTERP_HFILE: interpdata_set not allocated')
         173           0 :     else if ((hfilenum < 1) .or. (hfilenum > size(interpdata_set))) then
         174           0 :       call endrun('SET_INTERP_HFILE: hfilenum out of range')
         175           0 :     else if (hfilenum > size(interp_info)) then
         176           0 :       call endrun('SET_INTERP_HFILE: hfilenum out of range')
         177             :     else
         178           0 :       cam_interpolate => interpdata_set(hfilenum)%interpdata
         179           0 :       interp_lat => interp_info(hfilenum)%interp_lat
         180           0 :       interp_lon => interp_info(hfilenum)%interp_lon
         181           0 :       interp_gweight => interp_info(hfilenum)%interp_gweight
         182             :       call set_interp_parameter('nlat', interp_info(hfilenum)%interp_nlat)
         183             :       call set_interp_parameter('nlon', interp_info(hfilenum)%interp_nlon)
         184             :       call set_interp_parameter('itype', interp_info(hfilenum)%interp_type)
         185             :       call set_interp_parameter('gridtype', interp_info(hfilenum)%interp_gridtype)
         186             :     end if
         187           0 :   end subroutine set_interp_hfile
         188             :
         189           0 :   subroutine write_interpolated_scalar(File, varid, fld, numlev, data_type, decomp_type)
         190           0 :     use pio,              only: file_desc_t, var_desc_t
         191             :     use pio,              only: iosystem_desc_t
         192             :     use pio,              only: pio_initdecomp, pio_freedecomp
         193             :     use pio,              only: io_desc_t, pio_write_darray
         194             :     use pio,              only: pio_real
         195             :     use interpolate_mod,  only: interpolate_scalar
         196             :     use cam_instance,     only: atm_id
         197             :     use spmd_dyn,         only: local_dp_map
         198             :     use ppgrid,           only: begchunk
         199             :     use phys_grid,        only: get_dyn_col_p, columns_on_task, get_chunk_info_p
         200             :     use dimensions_mod,   only: fv_nphys, nc, nhc, nhc_phys
         201             :     use dof_mod,          only: PutUniquePoints
         202             :     use interpolate_mod,  only: get_interp_parameter
         203             :     use shr_pio_mod,      only: shr_pio_getiosys
         204             :     use edge_mod,         only: edgevpack, edgevunpack, initedgebuffer, freeedgebuffer
         205             :     use edgetype_mod,     only: EdgeBuffer_t
         206             :     use bndry_mod,        only: bndry_exchange
         207             :     use parallel_mod,     only: par
         208             :     use thread_mod,       only: horz_num_threads
         209             :     use cam_grid_support, only: cam_grid_id
         210             :     use hybrid_mod,       only: hybrid_t, config_thread_region, get_loop_ranges
         211             :     use fvm_mod,          only: fill_halo_and_extend_panel
         212             :
         213             :     type(file_desc_t), intent(inout) :: File
         214             :     type(var_desc_t) , intent(inout) :: varid
         215             :     real(r8),          intent(in)    :: fld(:,:,:)
         216             :     integer,           intent(in)    :: numlev, data_type, decomp_type
         217             :     !
         218             :     ! local variables
         219             :     !
         220             :     type(io_desc_t)                :: iodesc
         221             :     type(hybrid_t)                 :: hybrid
         222             :     type(iosystem_desc_t), pointer :: pio_subsystem
         223           0 :     type (EdgeBuffer_t)            :: edgebuf              ! edge buffer
         224             :
         225             :
         226             :     integer              :: lchnk, i, col_index, icol, ncols, ierr
         227             :     integer              :: nets, nete
         228             :     integer              :: phys_decomp, fvm_decomp, gll_decomp
         229             :
         230           0 :     real(r8), pointer     :: dest(:,:,:,:)
         231           0 :     real(r8), pointer     :: fldout(:,:)
         232           0 :     real(r8), allocatable :: fld_dyn(:,:,:), fld_tmp(:,:,:,:,:)
         233             :
         234             :     integer          :: st, en ! Start and end temporaries
         235             :     integer          :: ie, blk_ind(1), ncnt_out, k
         236           0 :     integer, pointer :: idof(:)
         237             :     integer          :: nlon, nlat, ncol, nsize, nhalo, nhcc
         238             :     logical          :: usefillvalues
         239             :     character(len=*), parameter :: subname = 'write_interpolated_scalar'
         240             :
         241           0 :     usefillvalues=.false.
         242             :
         243           0 :     phys_decomp = cam_grid_id('physgrid')
         244           0 :     fvm_decomp  = cam_grid_id('FVM')
         245           0 :     gll_decomp  = cam_grid_id('GLL')
         246             :     !
         247             :     ! There are 2 main scenarios regarding decomposition:
         248             :     !
         249             :     !   decomp_type==phys_decomp: we need to move data from physics decomposition to dynamics decomposition
         250             :     !   else                    : data is on dynamics decomposition
         251             :     !
         252           0 :     if (decomp_type==phys_decomp) then
         253           0 :        if(.not. local_dp_map) then
         254           0 :           call endrun(subname//': weak scaling does not support load balancing')
         255             :        end if
         256           0 :       if (fv_nphys > 0) then
         257             :         !
         258             :         ! note that even if fv_nphys<4 then SIZE(fld,DIM=1)=PCOLS
         259             :         !
         260           0 :         nsize = fv_nphys
         261           0 :         nhalo = 1!for bilinear only a halo of 1 is needed
         262           0 :         nhcc  = nhc_phys
         263             :       else
         264           0 :         nsize = np
         265           0 :         nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
         266           0 :         nhcc  = 0
         267             :       end if
         268           0 :     else if (decomp_type == fvm_decomp) then
         269             :       !
         270             :       ! CSLAM grid output
         271             :       !
         272           0 :       nsize = nc
         273           0 :       nhalo = 1!for bilinear only a halo of 1 is needed
         274           0 :       nhcc  = nhc
         275           0 :     else if (decomp_type == gll_decomp) then
         276           0 :       nsize = np
         277           0 :       nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
         278           0 :       nhcc  = 0
         279             :     else
         280           0 :       call endrun(subname//': unknown decomp_type')
         281             :     end if
         282           0 :     allocate(fld_tmp(1-nhcc:nsize+nhcc,1-nhcc:nsize+nhcc,numlev,1,nelemd))
         283           0 :     allocate(dest(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,numlev,nelemd))
         284             :
         285           0 :     nlon=get_interp_parameter('nlon')
         286           0 :     nlat=get_interp_parameter('nlat')
         287           0 :     pio_subsystem => shr_pio_getiosys(atm_id)
         288             :
         289           0 :     if(decomp_type == phys_decomp) then
         290             :
         291           0 :       allocate(fld_dyn(nsize*nsize,numlev,nelemd))
         292           0 :       fld_dyn = -999_R8
         293             :       !!$omp parallel do num_threads(horz_num_threads) private (col_index, lchnk, icol, ie, blk_ind, k)
         294           0 :       do col_index = 1, columns_on_task
         295           0 :          call get_dyn_col_p(col_index, ie, blk_ind)
         296           0 :          call get_chunk_info_p(col_index, lchnk, icol)
         297           0 :          do k = 1, numlev
         298           0 :             fld_dyn(blk_ind(1), k, ie) = fld(icol, k, lchnk-begchunk+1)
         299             :          end do
         300             :       end do
         301             :
         302           0 :       if (fv_nphys > 0) then
         303           0 :         do ie = 1, nelemd
         304           0 :           fld_tmp(1:nsize,1:nsize,:,1,ie) = RESHAPE(fld_dyn(:,:,ie),(/nsize,nsize,numlev/))
         305             :         end do
         306             :       else
         307           0 :         call initEdgeBuffer(par, edgebuf, elem, numlev,nthreads=1)
         308             :
         309           0 :         do ie = 1, nelemd
         310           0 :           ncols = elem(ie)%idxp%NumUniquePts
         311           0 :           call putUniquePoints(elem(ie)%idxP, numlev, fld_dyn(1:ncols,1:numlev,ie), fld_tmp(:,:,1:numlev,1,ie))
         312           0 :           call edgeVpack(edgebuf, fld_tmp(:,:,1:numlev,1,ie), numlev, 0, ie)
         313             :         end do
         314           0 :         if(iam < par%nprocs) then
         315           0 :           call bndry_exchange(par, edgebuf,location=subname)
         316             :         end if
         317           0 :         do ie = 1, nelemd
         318           0 :           call edgeVunpack(edgebuf, fld_tmp(:,:,1:numlev,1,ie), numlev, 0, ie)
         319             :         end do
         320           0 :         call freeEdgeBuffer(edgebuf)
         321             :         !check if fill values are present:
         322           0 :         usefillvalues = any(fld_tmp == fillvalue)
         323             :       end if
         324           0 :       deallocate(fld_dyn)
         325             :     else
         326             :       !
         327             :       ! not physics decomposition
         328             :       !
         329           0 :       do ie = 1, nelemd
         330           0 :         fld_tmp(1:nsize,1:nsize,1:numlev,1,ie) = RESHAPE(fld(1:nsize*nsize,1:numlev,ie),(/nsize,nsize,numlev/))
         331             :       end do
         332             :       !check if fillvalues are present:
         333           0 :       usefillvalues = any(fld_tmp == fillvalue)
         334             :     end if
         335             :     !
         336             :     ! code for non-GLL grids: need to fill halo and interpolate (if on panel edge/corner) for bilinear interpolation
         337             :     !
         338           0 :     if (decomp_type==fvm_decomp.or.(fv_nphys>0.and.decomp_type==phys_decomp)) then
         339             :       !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
         340             :       !JMD        hybrid = config_thread_region(par,'horizontal')
         341           0 :       hybrid = config_thread_region(par,'serial')
         342           0 :       call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
         343           0 :       call fill_halo_and_extend_panel(elem(nets:nete),fvm(nets:nete),&
         344           0 :            fld_tmp(:,:,:,:,nets:nete),hybrid,nets,nete,nsize,nhcc,nhalo,numlev,1,.true.,.true.)
         345             :
         346             :       !check if fill values are present:
         347           0 :       usefillvalues = any(fld_tmp(:,:,:,:,nets:nete) == fillvalue)
         348             :     end if
         349             :     !
         350             :     ! WARNING - 1:nelemd and nets:nete
         351             :     !
         352             :     !!$OMP MASTER   !JMD
         353           0 :     dest(:,:,:,1:nelemd) = fld_tmp(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,:,1,1:nelemd)
         354             :     !!$OMP END MASTER
         355           0 :     deallocate(fld_tmp)
         356             :     !
         357             :     !***************************************************************************
         358             :     !
         359             :     ! now data is on dynamics decomposition
         360             :     !
         361             :     !***************************************************************************
         362             :     !
         363           0 :     ncnt_out = sum(cam_interpolate(1:nelemd)%n_interp)
         364           0 :     allocate(fldout(ncnt_out,numlev))
         365           0 :     allocate(idof(ncnt_out*numlev))
         366           0 :     fldout = -999_r8
         367           0 :     idof = 0
         368           0 :     st = 1
         369             :
         370           0 :     do ie=1,nelemd
         371           0 :       ncol = cam_interpolate(ie)%n_interp
         372           0 :       do k=0,numlev-1
         373           0 :         do i=1,ncol
         374           0 :           idof(st+i-1+k*ncnt_out)=cam_interpolate(ie)%ilon(i)+nlon*(cam_interpolate(ie)%ilat(i)-1)+nlon*nlat*k
         375             :         enddo
         376             :       enddo
         377             :       ! Now that we have the field on the dyn grid we need to interpolate
         378           0 :       en = st+cam_interpolate(ie)%n_interp-1
         379           0 :       if(usefillvalues) then
         380           0 :         call interpolate_scalar(cam_interpolate(ie),dest(:,:,:,ie), nsize, nhalo, numlev, fldout(st:en,:), fillvalue)
         381             :       else
         382           0 :         call interpolate_scalar(cam_interpolate(ie),dest(:,:,:,ie), nsize, nhalo, numlev, fldout(st:en,:))
         383             :       end if
         384           0 :       st = en+1
         385             :     end do
         386             :
         387           0 :     if(numlev==1) then
         388           0 :        call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat/), idof, iodesc)
         389             :     else
         390           0 :        call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat,numlev/), idof, iodesc)
         391             :     end if
         392             :
         393           0 :     if(data_type == pio_real) then
         394           0 :       call pio_write_darray(File, varid, iodesc, real(fldout, r4), ierr)
         395             :     else
         396           0 :       call pio_write_darray(File, varid, iodesc, fldout, ierr)
         397             :     end if
         398             :
         399           0 :     deallocate(dest)
         400             :
         401           0 :     deallocate(fldout)
         402           0 :     deallocate(idof)
         403           0 :     call pio_freedecomp(file,iodesc)
         404             :
         405           0 :   end subroutine write_interpolated_scalar
         406             :
         407           0 :   subroutine write_interpolated_vector(File, varidu, varidv, fldu, fldv, numlev, data_type, decomp_type)
         408           0 :     use pio,              only: file_desc_t, var_desc_t
         409             :     use pio,              only: iosystem_desc_t
         410             :     use pio,              only: pio_initdecomp, pio_freedecomp
         411             :     use pio,              only: io_desc_t, pio_write_darray
         412             :     use pio,              only: pio_real
         413             :     use cam_instance,     only: atm_id
         414             :     use interpolate_mod,  only: interpolate_scalar, vec_latlon_to_contra,get_interp_parameter
         415             :     use spmd_dyn,         only: local_dp_map
         416             :     use ppgrid,           only: begchunk
         417             :     use phys_grid,        only: get_dyn_col_p, columns_on_task, get_chunk_info_p
         418             :     use dimensions_mod,   only: fv_nphys,nc,nhc,nhc_phys
         419             :     use dof_mod,          only: PutUniquePoints
         420             :     use shr_pio_mod,      only: shr_pio_getiosys
         421             :     use edge_mod,         only: edgevpack, edgevunpack, initedgebuffer, freeedgebuffer
         422             :     use edgetype_mod,     only: EdgeBuffer_t
         423             :     use bndry_mod,        only: bndry_exchange
         424             :     use parallel_mod,     only: par
         425             :     use thread_mod,       only: horz_num_threads
         426             :     use cam_grid_support, only: cam_grid_id
         427             :     use hybrid_mod,       only: hybrid_t,config_thread_region, get_loop_ranges
         428             :     use fvm_mod,          only: fill_halo_and_extend_panel
         429             :     use control_mod,      only: cubed_sphere_map
         430             :     use cube_mod,         only: dmap
         431             :
         432             :     type(file_desc_t), intent(inout) :: File
         433             :     type(var_desc_t),  intent(inout) :: varidu, varidv
         434             :     real(r8),          intent(in)    :: fldu(:,:,:), fldv(:,:,:)
         435             :     integer,           intent(in)    :: numlev, data_type, decomp_type
         436             :
         437             :     type(hybrid_t)                 :: hybrid
         438             :     type(io_desc_t)                :: iodesc
         439             :     type(iosystem_desc_t), pointer :: pio_subsystem
         440           0 :     type (EdgeBuffer_t)            :: edgebuf              ! edge buffer
         441             :
         442             :     integer              :: lchnk, i, col_index, icol, ncols, ierr
         443             :     integer              :: nets, nete
         444             :
         445           0 :     real(r8), allocatable :: dest(:,:,:,:,:)
         446           0 :     real(r8), pointer     :: fldout(:,:,:)
         447           0 :     real(r8), allocatable :: fld_dyn(:,:,:,:), fld_tmp(:,:,:,:,:)
         448             :
         449             :     integer          :: st, en ! Start and end temporaries
         450             :     integer          :: ie, blk_ind(1), ncnt_out, k
         451           0 :     integer, pointer :: idof(:)
         452             :     integer          :: nlon, nlat, ncol,nsize,nhalo,nhcc
         453             :     logical          :: usefillvalues
         454             :     integer          :: phys_decomp, fvm_decomp,gll_decomp
         455             :     real (r8)        :: D(2,2)   ! derivative of gnomonic mapping
         456             :     real (r8)        :: v1,v2
         457             :     character(len=*), parameter :: subname = 'write_interpolated_vector'
         458             :
         459           0 :     usefillvalues=.false.
         460             :
         461           0 :     phys_decomp = cam_grid_id('physgrid')
         462           0 :     fvm_decomp  = cam_grid_id('FVM')
         463           0 :     gll_decomp  = cam_grid_id('GLL')
         464             :     !
         465             :     ! There are 2 main scenarios regarding decomposition:
         466             :     !
         467             :     !   decomp_type==phys_decomp: we need to move data from physics decomposition to dynamics decomposition
         468             :     !   else                    : data is on dynamics decomposition
         469             :     !
         470           0 :     if (decomp_type==phys_decomp) then
         471           0 :        if(.not. local_dp_map) then
         472           0 :           call endrun(subname//': weak scaling does not support load balancing')
         473             :        end if
         474           0 :       if (fv_nphys > 0) then
         475             :         !
         476             :         ! note that even if fv_nphys<4 then SIZE(fld,DIM=1)=npsq
         477             :         !
         478           0 :         nsize = fv_nphys
         479           0 :         nhalo = 1!for bilinear only a halo of 1 is needed
         480           0 :         nhcc  = nhc_phys
         481             :       else
         482           0 :         nsize = np
         483           0 :         nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
         484           0 :         nhcc  = 0
         485             :       end if
         486           0 :     else if (decomp_type == fvm_decomp) then
         487             :       !
         488             :       ! CSLAM grid output
         489             :       !
         490           0 :       nsize = nc
         491           0 :       nhalo = 1!for bilinear only a halo of 1 is needed
         492           0 :       nhcc  = nhc
         493           0 :     else if (decomp_type == gll_decomp) then
         494           0 :       nsize = np
         495           0 :       nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
         496           0 :       nhcc  = 0
         497             :     else
         498           0 :       call endrun(subname//': unknown decomp_type')
         499             :     end if
         500           0 :     allocate(fld_tmp(1-nhcc:nsize+nhcc,1-nhcc:nsize+nhcc,2,numlev,nelemd))
         501           0 :     allocate(dest(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,2,numlev,nelemd))
         502             :
         503           0 :     nlon=get_interp_parameter('nlon')
         504           0 :     nlat=get_interp_parameter('nlat')
         505           0 :     pio_subsystem => shr_pio_getiosys(atm_id)
         506           0 :     if(decomp_type == phys_decomp) then
         507           0 :       allocate(fld_dyn(nsize*nsize,2,numlev,nelemd))
         508           0 :       fld_dyn = -999_R8
         509             :        !!$omp parallel do num_threads(horz_num_threads) private (col_index, lchnk, icol, ie, blk_ind, k)
         510           0 :       do col_index = 1, columns_on_task
         511           0 :          call get_dyn_col_p(col_index, ie, blk_ind)
         512           0 :          call get_chunk_info_p(col_index, lchnk, icol)
         513           0 :          do k = 1, numlev
         514           0 :             fld_dyn(blk_ind(1), 1, k, ie) = fldu(icol, k, lchnk-begchunk+1)
         515           0 :             fld_dyn(blk_ind(1), 2, k, ie) = fldv(icol, k, lchnk-begchunk+1)
         516             :          end do
         517             :       end do
         518           0 :       if (fv_nphys > 0) then
         519           0 :         do ie = 1, nelemd
         520           0 :           fld_tmp(1:nsize,1:nsize,:,:,ie) = RESHAPE(fld_dyn(:,:,:,ie),(/nsize,nsize,2,numlev/))
         521             :         end do
         522             :       else
         523           0 :         call initEdgeBuffer(par, edgebuf, elem, 2*numlev,nthreads=1)
         524             :
         525           0 :         do ie = 1, nelemd
         526           0 :           ncols = elem(ie)%idxp%NumUniquePts
         527           0 :           call putUniquePoints(elem(ie)%idxP, 2, numlev, fld_dyn(1:ncols,:,1:numlev,ie), fld_tmp(:,:,:,1:numlev,ie))
         528           0 :           call edgeVpack(edgebuf, fld_tmp(:,:,:,:,ie), 2*numlev, 0, ie)
         529             :         end do
         530           0 :         if(iam < par%nprocs) then
         531           0 :           call bndry_exchange(par, edgebuf,location=subname)
         532             :         end if
         533             :
         534           0 :         do ie = 1, nelemd
         535           0 :           call edgeVunpack(edgebuf, fld_tmp(:,:,:,:,ie), 2*numlev, 0, ie)
         536             :         end do
         537           0 :         call freeEdgeBuffer(edgebuf)
         538             :         !check if fill values are present:
         539           0 :         usefillvalues = any(fld_tmp==fillvalue)
         540             :       end if
         541           0 :       deallocate(fld_dyn)
         542             :     else
         543             :       !
         544             :       ! not physics decomposition
         545             :       !
         546             :       !check if fill values are present:
         547           0 :       usefillvalues = (any(fldu(1:nsize:1,nsize,:)==fillvalue) .or. any(fldv(1:nsize:1,nsize,:)==fillvalue))
         548           0 :       do ie = 1, nelemd
         549           0 :         fld_tmp(1:nsize,1:nsize,1,1:numlev,ie) = RESHAPE(fldu(1:nsize*nsize,1:numlev,ie),(/nsize,nsize,numlev/))
         550           0 :         fld_tmp(1:nsize,1:nsize,2,1:numlev,ie) = RESHAPE(fldv(1:nsize*nsize,1:numlev,ie),(/nsize,nsize,numlev/))
         551             :       end do
         552             :     endif
         553             :     !
         554             :     !***************************************************************************
         555             :     !
         556             :     ! now data is on dynamics decomposition
         557             :     !
         558             :     !***************************************************************************
         559             :     !
         560           0 :     if (decomp_type==fvm_decomp.or.(fv_nphys>0.and.decomp_type==phys_decomp)) then
         561             :       !
         562             :       !***************************************************************************
         563             :       !
         564             :       ! code for non-GLL grids: need to fill halo and interpolate
         565             :       ! (if on panel edge/corner) for bilinear interpolation
         566             :       !
         567             :       !***************************************************************************
         568             :       !
         569             :
         570             :       !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
         571             :       !JMD        hybrid = config_thread_region(par,'horizontal')
         572           0 :       hybrid = config_thread_region(par,'serial')
         573           0 :       call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
         574           0 :       call fill_halo_and_extend_panel(elem(nets:nete),fvm(nets:nete),&
         575           0 :            fld_tmp(:,:,:,:,nets:nete),hybrid,nets,nete,nsize,nhcc,nhalo,2,numlev,.true.,.false.)
         576           0 :       do ie=nets,nete
         577           0 :         call vec_latlon_to_contra(elem(ie),nsize,nhcc,numlev,fld_tmp(:,:,:,:,ie),fvm(ie))
         578             :       end do
         579           0 :       call fill_halo_and_extend_panel(elem(nets:nete),fvm(nets:nete),&
         580           0 :            fld_tmp(:,:,:,:,nets:nete),hybrid,nets,nete,nsize,nhcc,nhalo,2,numlev,.false.,.true.)
         581             :
         582             :       !check if fill values are present:
         583           0 :       usefillvalues = any(fld_tmp(:,:,:,:,nets:nete) == fillvalue)
         584             :     else
         585           0 :       do ie=1,nelemd
         586           0 :         call vec_latlon_to_contra(elem(ie),nsize,nhcc,numlev,fld_tmp(:,:,:,:,ie))
         587             :       end do
         588             :     end if
         589             :     !
         590             :     ! WARNING - 1:nelemd and nets:nete
         591             :     !
         592             :     !!$OMP MASTER   !JMD
         593           0 :     dest(:,:,:,:,1:nelemd) = fld_tmp(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,:,:,1:nelemd)
         594             :     !!$OMP END MASTER
         595           0 :     deallocate(fld_tmp)
         596             :     !
         597             :     !***************************************************************************
         598             :     !
         599             :     ! do mapping from source grid to latlon grid
         600             :     !
         601             :     !***************************************************************************
         602             :     !
         603           0 :     ncnt_out = sum(cam_interpolate(1:nelemd)%n_interp)
         604           0 :     allocate(fldout(ncnt_out,numlev,2))
         605           0 :     allocate(idof(ncnt_out*numlev))
         606             :
         607           0 :     fldout = -999_r8
         608           0 :     idof = 0
         609           0 :     st = 1
         610           0 :     do ie=1,nelemd
         611           0 :       ncol = cam_interpolate(ie)%n_interp
         612           0 :       do k=0,numlev-1
         613           0 :         do i=1,ncol
         614           0 :           idof(st+i-1+k*ncnt_out)=cam_interpolate(ie)%ilon(i)+nlon*(cam_interpolate(ie)%ilat(i)-1)+nlon*nlat*k
         615             :         enddo
         616             :       enddo
         617             :       ! Now that we have the field on the dyn grid we need to interpolate
         618           0 :       en = st+cam_interpolate(ie)%n_interp-1
         619           0 :       if(usefillvalues) then
         620           0 :         call interpolate_scalar(cam_interpolate(ie),dest(:,:,1,:,ie), nsize, nhalo, numlev, fldout(st:en,:,1), fillvalue)
         621           0 :         call interpolate_scalar(cam_interpolate(ie),dest(:,:,2,:,ie), nsize, nhalo, numlev, fldout(st:en,:,2), fillvalue)
         622             :       else
         623           0 :         call interpolate_scalar(cam_interpolate(ie),dest(:,:,1,:,ie), nsize, nhalo, numlev, fldout(st:en,:,1))
         624           0 :         call interpolate_scalar(cam_interpolate(ie),dest(:,:,2,:,ie), nsize, nhalo, numlev, fldout(st:en,:,2))
         625             :       end if
         626             :       !
         627             :       ! convert from contravariant components to lat-lon
         628             :       !
         629           0 :       do i=1,cam_interpolate(ie)%n_interp
         630             :         ! convert fld from contra->latlon
         631           0 :         call dmap(D,cam_interpolate(ie)%interp_xy(i)%x,cam_interpolate(ie)%interp_xy(i)%y,&
         632           0 :              elem(ie)%corners3D,cubed_sphere_map,elem(ie)%corners,elem(ie)%u2qmap,elem(ie)%facenum)
         633             :         ! convert fld from contra->latlon
         634           0 :         do k=1,numlev
         635           0 :           v1 = fldout(st+i-1,k,1)
         636           0 :           v2 = fldout(st+i-1,k,2)
         637           0 :           fldout(st+i-1,k,1)=D(1,1)*v1 + D(1,2)*v2
         638           0 :           fldout(st+i-1,k,2)=D(2,1)*v1 + D(2,2)*v2
         639             :         end do
         640             :       end do
         641           0 :       st = en+1
         642             :     end do
         643             :
         644           0 :     if(numlev==1) then
         645           0 :        call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat/), idof, iodesc)
         646             :     else
         647           0 :        call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat,numlev/), idof, iodesc)
         648             :     end if
         649             :
         650           0 :     if(data_type == pio_real) then
         651           0 :       call pio_write_darray(File, varidu, iodesc, real(fldout(:,:,1), r4), ierr)
         652           0 :       call pio_write_darray(File, varidv, iodesc, real(fldout(:,:,2), r4), ierr)
         653             :     else
         654           0 :       call pio_write_darray(File, varidu, iodesc, fldout(:,:,1), ierr)
         655           0 :       call pio_write_darray(File, varidv, iodesc, fldout(:,:,2), ierr)
         656             :     end if
         657             :
         658             :
         659           0 :     deallocate(fldout)
         660           0 :     deallocate(idof)
         661           0 :     deallocate(dest)
         662           0 :     call pio_freedecomp(file,iodesc)
         663             :
         664           0 :   end subroutine write_interpolated_vector
         665             :
         666           0 : end module interp_mod

Generated by: LCOV version 1.14