LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - fvm_mapping.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 297 | 584 | 50.9 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 13 | 21 | 61.9 %

              Line data    Source code
           1             : !
           2             : ! pg3->GLL and GLL->pg3 mapping algorithm described in:
           3             : !
           4             : ! Adam R. Herrington, Peter H. Lauritzen, Mark A. Taylor, Steve Goldhaber, Brian Eaton, Kevin A  Reed and Paul A. Ullrich, 2018:
           5             : ! Physics-dynamics coupling with element-based high-order Galerkin methods: quasi equal-area physics grid:
           6             : ! Mon. Wea. Rev., DOI:MWR-D-18-0136.1
           7             : !
           8             : ! pg2->pg3 mapping algorithm described in:
           9             : !
          10             : ! Adam R. Herrington, Peter H. Lauritzen, Kevin A  Reed, Steve Goldhaber, and Brian Eaton, 2019:
          11             : ! Exploring a lower resolution physics grid in CAM-SE-CSLAM. J. Adv. Model. Earth Syst.
          12             : !
          13             : !#define PCoM !replace PPM with PCoM for mass variables for fvm2phys and phys2fvm
          14             : !#define skip_high_order_fq_map !do mass and correlation preserving phys2fvm mapping but no high-order pre-mapping of fq
          15             : #define mass_fix
          16             : module fvm_mapping
          17             :   use shr_kind_mod,           only: r8=>shr_kind_r8
          18             :   use dimensions_mod,         only: irecons_tracer
          19             :   use element_mod,            only: element_t
          20             :   use fvm_control_volume_mod, only: fvm_struct
          21             :   use perf_mod,               only: t_startf, t_stopf
          22             :   use cam_abortutils,         only: endrun
          23             :   use cam_logfile,            only: iulog
          24             :   implicit none
          25             :   private
          26             :
          27             :   public :: phys2dyn_forcings_fvm, dyn2phys, dyn2phys_vector, dyn2phys_all_vars,dyn2fvm_mass_vars
          28             :   public :: phys2dyn,fvm2dyn,dyn2fvm,cslam2gll
          29             :   save
          30             :   integer                                            :: save_max_overlap
          31             :   real(kind=r8), allocatable, dimension(:,:,:,:,:)   :: save_air_mass_overlap
          32             :   real(kind=r8), allocatable, dimension(:,:,:,:,:,:) :: save_q_overlap
          33             :   real(kind=r8), allocatable, dimension(:,:,:,:,:)   :: save_q_phys
          34             :   real(kind=r8), allocatable, dimension(:,:,:,:)     :: save_dp_phys
          35             :   real(kind=r8), allocatable, dimension(:,:,:,:)     :: save_overlap_area
          36             :   integer      , allocatable, dimension(:,:,:,:,:)   :: save_overlap_idx
          37             :   integer      , allocatable, dimension(:,:,:,:)     :: save_num_overlap
          38             :
          39             :   interface fvm2dyn
          40             :      module procedure fvm2dynt1
          41             :      module procedure fvm2dyntn
          42             :   end interface fvm2dyn
          43             :
          44             : contains
          45             :   !
          46             :   ! map all mass variables from gll to fvm
          47             :   !
          48      369408 :   subroutine phys2dyn_forcings_fvm(elem, fvm, hybrid,nets,nete,no_cslam, tl_f, tl_qdp)
          49             :     use dimensions_mod,         only: np, nc,nlev
          50             :     use dimensions_mod,         only: fv_nphys, nhc_phys,ntrac,nhc,ksponge_end, nu_scale_top
          51             :     use hybrid_mod,             only: hybrid_t
          52             :     use air_composition,        only: thermodynamic_active_species_num, thermodynamic_active_species_idx
          53             :     type (element_t), intent(inout):: elem(:)
          54             :     type(fvm_struct), intent(inout):: fvm(:)
          55             :
          56             :     type (hybrid_t), intent(in)    :: hybrid  ! distributed parallel structure (shared)
          57             :     logical, intent(in)            :: no_cslam
          58             :     integer, intent(in)            :: nets, nete, tl_f, tl_qdp
          59             :
          60             :     integer                                             :: ie,i,j,k,m_cnst,nq
          61      369408 :     real (kind=r8), dimension(:,:,:,:,:)  , allocatable :: fld_phys, fld_gll
          62             :     real (kind=r8)  :: element_ave
          63             :     !
          64             :     ! for tensor product Lagrange interpolation
          65             :     !
          66             :     integer              :: nflds
          67      369408 :     logical, allocatable :: llimiter(:)
          68             :
          69             :     integer :: ierr
          70             :
          71      369408 :     if (no_cslam) then
          72           0 :       call endrun("phys2dyn_forcings_fvm: no cslam case: NOT SUPPORTED")
          73      369408 :     else if (nc.ne.fv_nphys) then
          74             :       !
          75             :       !***********************************************************
          76             :       !
          77             :       ! using cslam and different resolution physics grid
          78             :       !
          79             :       !***********************************************************
          80             :       !
          81           0 :       call t_startf('p2d-pg2:copying')
          82           0 :       nflds = 4+ntrac
          83           0 :       allocate(fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev,nflds,nets:nete), stat=ierr)
          84           0 :       if( ierr /= 0 ) then
          85           0 :          write(iulog,*) 'phys2dyn_forcings_fvm:  fld_phys allocation error = ',ierr
          86           0 :          call endrun('phys2dyn_forcings_fvm: failed to allocate fld_phys array')
          87             :       end if
          88           0 :       allocate(fld_gll(np,np,nlev,3,nets:nete), stat=ierr)
          89           0 :       if( ierr /= 0 ) then
          90           0 :          write(iulog,*) 'phys2dyn_forcings_fvm:  fld_gll allocation error = ',ierr
          91           0 :          call endrun('phys2dyn_forcings_fvm: failed to allocate fld_gll array')
          92             :       end if
          93           0 :       allocate(llimiter(3), stat=ierr)
          94           0 :       if( ierr /= 0 ) then
          95           0 :          write(iulog,*) 'phys2dyn_forcings_fvm:  llimiter allocation error = ',ierr
          96           0 :          call endrun('phys2dyn_forcings_fvm: failed to allocate llimiter array')
          97             :       end if
          98           0 :       fld_phys = -9.99E99_r8!xxx necessary?
          99             :
         100           0 :       llimiter          = .false.
         101             :
         102           0 :       do ie=nets,nete
         103             :         !
         104             :         ! pack fields that need to be interpolated
         105             :         !
         106           0 :         fld_phys(1:fv_nphys,1:fv_nphys,:,1,ie)       = fvm(ie)%ft(1:fv_nphys,1:fv_nphys,:)
         107           0 :         fld_phys(1:fv_nphys,1:fv_nphys,:,2,ie)       = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,1,:)
         108           0 :         fld_phys(1:fv_nphys,1:fv_nphys,:,3,ie)       = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,2,:)
         109           0 :         fld_phys(1:fv_nphys,1:fv_nphys,:,4,ie)       = fvm(ie)%dp_phys(1:fv_nphys,1:fv_nphys,:)
         110           0 :         do m_cnst=1,ntrac
         111           0 :           fld_phys(1:fv_nphys,1:fv_nphys,:,4+m_cnst,ie) = &
         112           0 :                fvm(ie)%fc_phys(1:fv_nphys,1:fv_nphys,:,m_cnst)
         113             :         end do
         114             :       end do
         115           0 :       call t_stopf('p2d-pg2:copying')
         116           0 :       call t_startf('p2d-pg2:fill_halo_phys')
         117           0 :       call fill_halo_phys(fld_phys,hybrid,nets,nete,nlev,nflds)
         118             :       !
         119             :       ! do mapping of fu,fv,ft
         120             :       !
         121           0 :       call phys2dyn(hybrid,elem,fld_phys(:,:,:,1:3,:),fld_gll,nets,nete,nlev,3,fvm,llimiter, &
         122           0 :                              istart_vector=2,halo_filled=.true.)
         123           0 :       do ie=nets,nete
         124           0 :         elem(ie)%derived%fT(:,:,:)   = fld_gll(:,:,:,1,ie)
         125           0 :         elem(ie)%derived%fM(:,:,1,:) = fld_gll(:,:,:,2,ie)
         126           0 :         elem(ie)%derived%fM(:,:,2,:) = fld_gll(:,:,:,3,ie)
         127             :       end do
         128           0 :       call t_stopf('p2d-pg2:fill_halo_phys')
         129             :
         130           0 :       deallocate(fld_gll)
         131             :       !
         132             :       ! map fq from phys to fvm
         133             :       !
         134           0 :       call t_startf('p2d-pg2:phys2fvm')
         135             :
         136           0 :       do ie=nets,nete
         137           0 :          do k=1,nlev
         138           0 :            call phys2fvm(ie,k,fvm(ie),&
         139           0 :                 fld_phys(:,:,k,5:4+ntrac,ie),fvm(ie)%fc(:,:,k,1:ntrac),ntrac)
         140             :          end do
         141             :        end do
         142           0 :        call t_stopf('p2d-pg2:phys2fvm')
         143             :        !deallocate arrays allocated in dyn2phys_all_vars
         144           0 :        deallocate(save_air_mass_overlap,save_q_phys,save_q_overlap,&
         145           0 :             save_overlap_area,save_num_overlap,save_overlap_idx,save_dp_phys)
         146             :      else
         147             :        !
         148             :        !
         149             :        !*****************************************************************************************
         150             :        !
         151             :        ! using cslam with same physics grid resolution as cslam resolution
         152             :        !
         153             :        !*****************************************************************************************
         154             :        !
         155             :        ! nflds is ft, fu, fv, + thermo species
         156      369408 :        nflds = 3
         157     2585856 :        allocate(fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev,nflds,nets:nete))
         158     1108224 :        allocate(fld_gll(np,np,nlev,nflds,nets:nete))
         159      369408 :        allocate(llimiter(nflds))
         160     1477632 :        llimiter(1:nflds) = .false.
         161     2966808 :        do ie=nets,nete
         162             :          !
         163             :          ! pack fields that need to be interpolated
         164             :          !
         165  3142854000 :          fld_phys(1:fv_nphys,1:fv_nphys,:,1,ie)       = fvm(ie)%ft(1:fv_nphys,1:fv_nphys,:)
         166  3142854000 :          fld_phys(1:fv_nphys,1:fv_nphys,:,2,ie)       = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,1,:)
         167  3143223408 :          fld_phys(1:fv_nphys,1:fv_nphys,:,3,ie)       = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,2,:)
         168             :        end do
         169             :        !
         170             :        ! do mapping
         171             :        !
         172      369408 :        call phys2dyn(hybrid,elem,fld_phys,fld_gll,nets,nete,nlev,nflds,fvm,llimiter,2)
         173     2966808 :        do ie=nets,nete
         174  5075319600 :          elem(ie)%derived%fT(:,:,:)   = fld_gll(:,:,:,1,ie)
         175  5075319600 :          elem(ie)%derived%fM(:,:,1,:) = fld_gll(:,:,:,2,ie)
         176  5075689008 :          elem(ie)%derived%fM(:,:,2,:) = fld_gll(:,:,:,3,ie)
         177             :        end do
         178      369408 :        deallocate(fld_gll)
         179     2966808 :        do ie=nets,nete
         180    31538208 :          do m_cnst = 1,ntrac
         181 34573991400 :            fvm(ie)%fc(1:nc,1:nc,:,m_cnst) = fvm(ie)%fc_phys(1:nc,1:nc,:,m_cnst)*fvm(ie)%dp_fvm(1:nc,1:nc,:)
         182             :          end do
         183             :        end do
         184             :      end if
         185      369408 :      deallocate(fld_phys,llimiter)
         186      369408 :   end subroutine phys2dyn_forcings_fvm
         187             :
         188             :   ! for multiple fields
         189     1477632 :   subroutine fvm2dyntn(fld_fvm,fld_gll,hybrid,nets,nete,numlev,num_flds,fvm,llimiter,halo_filled)
         190             :     use dimensions_mod, only: np, nhc, nc
         191             :     use hybrid_mod    , only: hybrid_t
         192             :     use bndry_mod     , only: ghost_exchange
         193             :     use edge_mod      , only: ghostpack,ghostunpack
         194             :     use fvm_mod       , only: ghostBufQnhc_s
         195             :     !
         196             :     integer              , intent(in)    :: nets,nete,num_flds,numlev
         197             :     real (kind=r8), intent(inout) :: fld_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,numlev,num_flds,nets:nete)
         198             :     real (kind=r8), intent(out)   :: fld_gll(np,np,numlev,num_flds,nets:nete)
         199             :     type (hybrid_t)      , intent(in)    :: hybrid
         200             :     type(fvm_struct)     , intent(in)    :: fvm(nets:nete)
         201             :     logical              , intent(in)    :: llimiter(num_flds)
         202             :     logical, optional    , intent(in)    :: halo_filled !optional if boundary exchange for fld_fvm has already been called
         203             :
         204             :     integer                              :: ie, iwidth
         205             :     logical                              :: fill_halo
         206             :     !
         207             :     !*********************************************
         208             :     !
         209             :     ! halo exchange
         210             :     !
         211             :     !*********************************************
         212             :     !
         213      738816 :     fill_halo = .true.
         214      738816 :     if (present(halo_filled)) then
         215      738816 :        fill_halo = .not. halo_filled
         216             :     end if
         217             :
         218      738816 :     if (fill_halo) then
         219           0 :       do ie=nets,nete
         220           0 :         call ghostpack(ghostBufQnhc_s, fld_fvm(:,:,:,:,ie),numlev*num_flds,0,ie)
         221             :       end do
         222           0 :       call ghost_exchange(hybrid,ghostbufQnhc_s,location='fvm2dyntn')
         223           0 :       do ie=nets,nete
         224           0 :         call ghostunpack(ghostbufQnhc_s, fld_fvm(:,:,:,:,ie),numlev*num_flds,0,ie)
         225             :       end do
         226             :     end if
         227             :     !
         228             :     ! mapping
         229             :     !
         230      738816 :     iwidth=2
         231             : !    iwidth=1 !low-order mapping
         232     5933616 :     do ie=nets,nete
         233     5194800 :       call tensor_lagrange_interp(fvm(ie)%cubeboundary,np,nc,nhc,numlev,num_flds,fld_fvm(:,:,:,:,ie),&
         234    11128416 :            fld_gll(:,:,:,:,ie),llimiter,iwidth,fvm(ie)%norm_elem_coord)
         235             :     end do
         236      738816 :   end subroutine fvm2dyntn
         237             :
         238             :   ! for single field
         239           0 :   subroutine fvm2dynt1(fld_fvm,fld_gll,hybrid,nets,nete,numlev,fvm,llimiter,halo_filled)
         240      738816 :     use dimensions_mod, only: np, nhc, nc
         241             :     use hybrid_mod    , only: hybrid_t
         242             :     use bndry_mod     , only: ghost_exchange
         243             :     use edge_mod      , only: ghostpack,ghostunpack
         244             :     use fvm_mod       , only: ghostBufQnhc_t1
         245             :     !
         246             :     integer              , intent(in)    :: nets,nete,numlev
         247             :     real (kind=r8), intent(inout) :: fld_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,numlev,1,nets:nete)
         248             :     real (kind=r8), intent(out)   :: fld_gll(np,np,numlev,1,nets:nete)
         249             :     type (hybrid_t)      , intent(in)    :: hybrid
         250             :     type(fvm_struct)     , intent(in)    :: fvm(nets:nete)
         251             :     logical              , intent(in)    :: llimiter(1)
         252             :     logical, optional    , intent(in)    :: halo_filled!optional if boundary exchange for fld_fvm has already been called
         253             :
         254             :     integer                              :: ie, iwidth
         255             :     logical                              :: fill_halo
         256             :     !
         257             :     !*********************************************
         258             :     !
         259             :     ! halo exchange
         260             :     !
         261             :     !*********************************************
         262             :     !
         263           0 :     fill_halo = .true.
         264           0 :     if (present(halo_filled)) then
         265           0 :        fill_halo = .not. halo_filled
         266             :     end if
         267             :
         268           0 :     if (fill_halo) then
         269           0 :       do ie=nets,nete
         270           0 :         call ghostpack(ghostBufQnhc_t1, fld_fvm(:,:,:,1,ie),numlev,0,ie)
         271             :       end do
         272           0 :       call ghost_exchange(hybrid,ghostbufQnhc_t1,location='fvm2dynt1')
         273           0 :       do ie=nets,nete
         274           0 :         call ghostunpack(ghostbufQnhc_t1, fld_fvm(:,:,:,1,ie),numlev,0,ie)
         275             :       end do
         276             :     end if
         277             :     !
         278             :     ! mapping
         279             :     !
         280           0 :     iwidth=2
         281           0 :     do ie=nets,nete
         282           0 :       call tensor_lagrange_interp(fvm(ie)%cubeboundary,np,nc,nhc,numlev,1,fld_fvm(:,:,:,:,ie),&
         283           0 :            fld_gll(:,:,:,:,ie),llimiter,iwidth,fvm(ie)%norm_elem_coord)
         284             :     end do
         285           0 :   end subroutine fvm2dynt1
         286             :
         287      738816 :   subroutine fill_halo_phys(fld_phys,hybrid,nets,nete,num_lev,num_flds)
         288           0 :     use dimensions_mod, only: nhc_phys, fv_nphys
         289             :     use hybrid_mod    , only: hybrid_t
         290             :     use bndry_mod     , only: ghost_exchange
         291             :     use edge_mod      , only: ghostpack, ghostunpack
         292             :     use fvm_mod       , only: ghostBufPG_s
         293             :
         294             :     integer              , intent(in)    :: nets,nete,num_lev,num_flds
         295             :     real (kind=r8), intent(inout) :: fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_lev,num_flds, &
         296             :          nets:nete)
         297             :     type (hybrid_t)      , intent(in)    :: hybrid  ! distributed parallel structure (shared)
         298             :
         299             :     integer                 :: ie
         300             :     !
         301             :     !*********************************************
         302             :     !
         303             :     ! halo exchange
         304             :     !
         305             :     !*********************************************
         306             :     !
         307      369408 :     call t_startf('fvm:fill_halo_phys')
         308     2966808 :     do ie=nets,nete
         309     2966808 :        call ghostpack(ghostBufPG_s, fld_phys(:,:,:,:,ie),num_lev*num_flds,0,ie)
         310             :     end do
         311             :
         312      369408 :     call ghost_exchange(hybrid,ghostBufPG_s,location='fill_halo_phys')
         313             :
         314     2966808 :     do ie=nets,nete
         315     2966808 :        call ghostunpack(ghostBufPG_s, fld_phys(:,:,:,:,ie),num_lev*num_flds,0,ie)
         316             :     end do
         317             :     !
         318      369408 :     call t_stopf('fvm:fill_halo_phys')
         319      369408 :   end subroutine fill_halo_phys
         320             :   !
         321             :   ! must call fill_halo_phys before calling this subroutine
         322             :   !
         323      369408 :   subroutine phys2dyn(hybrid,elem,fld_phys,fld_gll,nets,nete,num_lev,num_flds,fvm,llimiter,istart_vector,halo_filled)
         324      369408 :     use dimensions_mod, only: np, nhc_phys, fv_nphys
         325             :     use hybrid_mod, only : hybrid_t
         326             :     type (hybrid_t), intent(in)   :: hybrid  ! distributed parallel structure (shared)
         327             :     integer       , intent(in)    :: nets,nete,num_flds,num_lev
         328             :     real (kind=r8), intent(inout) :: fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_lev,num_flds, &
         329             :          nets:nete)
         330             :     real (kind=r8), intent(out)   :: fld_gll(np,np,num_lev,num_flds,nets:nete)
         331             :     type (element_t)     , intent(inout) :: elem(:)
         332             :     type(fvm_struct)     , intent(in)    :: fvm(:)
         333             :     integer, optional    , intent(in)    :: istart_vector
         334             :     logical              , intent(in)    :: llimiter(num_flds)
         335             :     logical, optional    , intent(in)    :: halo_filled!optional if boundary exchange for fld_fvm has already been called
         336             :
         337             :     integer                 :: i, j, ie, k, iwidth
         338             :     real (kind=r8)   :: v1,v2
         339             :
         340      369408 :     if (present(halo_filled)) then
         341           0 :       if (.not.halo_filled) call fill_halo_phys(fld_phys,hybrid,nets,nete,num_lev,num_flds)
         342             :     else
         343      369408 :       call fill_halo_phys(fld_phys,hybrid,nets,nete,num_lev,num_flds)
         344             :     end if
         345      369408 :     if (present(istart_vector)) then
         346     2966808 :       do ie=nets,nete
         347   244525008 :         do k=1,num_lev
         348  2418179400 :           do j=1-nhc_phys,fv_nphys+nhc_phys
         349 21981796200 :             do i=1-nhc_phys,fv_nphys+nhc_phys
         350             :               !
         351             :               ! convert lat-lon vectors to contra-variant gnomonic
         352             :               !
         353 19566214200 :               v1 = fld_phys(i,j,k,istart_vector  ,ie)
         354 19566214200 :               v2 = fld_phys(i,j,k,istart_vector+1,ie)
         355 19566214200 :               fld_phys(i,j,k,istart_vector  ,ie)=fvm(ie)%Dinv_physgrid(i,j,1,1)*v1 + fvm(ie)%Dinv_physgrid(i,j,1,2)*v2
         356 21740238000 :               fld_phys(i,j,k,istart_vector+1,ie)=fvm(ie)%Dinv_physgrid(i,j,2,1)*v1 + fvm(ie)%Dinv_physgrid(i,j,2,2)*v2
         357             :             end do
         358             :           end do
         359             :         end do
         360             :       end do
         361             :     end if
         362             :     !
         363             :     ! mapping
         364             :     !
         365      369408 :     iwidth=2
         366             : !    iwidth=1
         367      369408 :     if (fv_nphys==1) iwidth=1
         368     2966808 :     do ie=nets,nete
         369     2597400 :       call tensor_lagrange_interp(fvm(ie)%cubeboundary,np,fv_nphys,nhc_phys,num_lev,num_flds,fld_phys(:,:,:,:,ie),&
         370     5564208 :            fld_gll(:,:,:,:,ie),llimiter,iwidth,fvm(ie)%norm_elem_coord_physgrid)
         371             :     end do
         372             :
         373      369408 :     if (present(istart_vector)) then
         374             :       !
         375             :       ! convert contra-variant to lat-lon
         376             :       !
         377     2966808 :       do ie=nets,nete
         378   244525008 :         do k=1,num_lev
         379  1210388400 :           do j=1,np
         380  5072722200 :             do i=1,np
         381  3864931200 :               v1 = fld_gll(i,j,k,istart_vector  ,ie)
         382  3864931200 :               v2 = fld_gll(i,j,k,istart_vector+1,ie)
         383  3864931200 :               fld_gll(i,j,k,istart_vector  ,ie) = elem(ie)%D(i,j,1,1)*v1 + elem(ie)%D(i,j,1,2)*v2
         384  4831164000 :               fld_gll(i,j,k,istart_vector+1,ie) = elem(ie)%D(i,j,2,1)*v1 + elem(ie)%D(i,j,2,2)*v2
         385             :             end do
         386             :           end do
         387             :         end do
         388             :       end do
         389             :     end if
         390      369408 :   end subroutine phys2dyn
         391             :   !
         392             :   ! map all mass variables from gll to fvm
         393             :   !
         394        5400 :   subroutine dyn2fvm_mass_vars(dp_gll,ps_gll,q_gll,&
         395        5400 :        dp_fvm,ps_fvm,q_fvm,num_trac,metdet,inv_area)
         396             :     use dimensions_mod, only: np, nc,nlev
         397             :     integer, intent(in) :: num_trac
         398             :     real (kind=r8), dimension(np,np,nlev)         , intent(in) :: dp_gll
         399             :     real (kind=r8), dimension(np,np,nlev,num_trac), intent(in) :: q_gll
         400             :     real (kind=r8), dimension(np,np)              , intent(in) :: ps_gll
         401             :
         402             :
         403             :     real (kind=r8), dimension(nc,nc,nlev)         , intent(inout) :: dp_fvm
         404             :     real (kind=r8), dimension(nc,nc,nlev,num_trac), intent(inout) :: q_fvm
         405             :     real (kind=r8), dimension(nc,nc)     , intent(inout)        :: ps_fvm
         406             :     real (kind=r8), dimension(nc,nc)     , intent(out)          :: inv_area
         407             :
         408             :     real (kind=r8), intent(in)           :: metdet(np,np)
         409             :
         410             :     real (kind=r8) :: se_area_sphere(nc,nc), tmp(np,np)
         411             :     real (kind=r8) :: inv_darea_dp_fvm(nc,nc)
         412             :     integer        :: k,m_cnst
         413             :
         414      113400 :     tmp = 1.0_r8
         415        5400 :     se_area_sphere = dyn2fvm(tmp,metdet)
         416       70200 :     inv_area = 1.0_r8/se_area_sphere
         417             :
         418       70200 :     ps_fvm(:,:) = dyn2fvm(ps_gll,metdet,inv_area)
         419      507600 :     do k=1,nlev
         420     6528600 :       dp_fvm(:,:,k) = dyn2fvm(dp_gll(:,:,k),metdet,inv_area)
         421     6528600 :       inv_darea_dp_fvm = inv_area/dp_fvm(:,:,k)
         422     6031800 :       do m_cnst=1,num_trac
         423             :         q_fvm(:,:,k,m_cnst) = &
         424     5524200 :              dyn2fvm(q_gll(:,:,k,m_cnst)*dp_gll(:,:,k),metdet,&
         425   188325000 :              inv_darea_dp_fvm,q_gll(:,:,k,m_cnst))
         426             :       end do
         427             :     end do
         428        5400 :   end subroutine dyn2fvm_mass_vars
         429             :
         430             :   !
         431             :   ! this subroutine assumes that the fvm halo has already been filled
         432             :   ! (if nc/=fv_nphys)
         433             :   !
         434             :
         435      370944 :   subroutine dyn2phys_all_vars(nets,nete,elem,fvm,&
         436             :        num_trac,ptop,tl,&
         437      370944 :        dp3d_phys,ps_phys,q_phys,T_phys,omega_phys,phis_phys)
         438             :     use dimensions_mod, only: np, nc,nlev,fv_nphys
         439             :     use dp_mapping,     only: nphys_pts
         440             :     use element_mod,    only: element_t
         441             :     integer, intent(in) :: nets,nete,num_trac,tl
         442             :
         443             :     type(fvm_struct), dimension(nets:nete), intent(inout):: fvm
         444             :     type(element_t),  dimension(nets:nete), intent(in)   :: elem
         445             :
         446             :     real (kind=r8), intent(in)           :: ptop
         447             :
         448             :     real (kind=r8), dimension(nphys_pts,nets:nete)               , intent(out) :: ps_phys,phis_phys
         449             :     real (kind=r8), dimension(nphys_pts,nlev,nets:nete)          , intent(out) :: dp3d_phys,T_phys,omega_phys
         450             :     real (kind=r8), dimension(nphys_pts,nlev,num_trac,nets:nete) , intent(out) :: q_phys
         451             :
         452             :
         453             :     real (kind=r8) :: tmp(np,np)
         454      741888 :     real (kind=r8), dimension(fv_nphys,fv_nphys)          :: inv_area,inv_darea_dp_phys,dp3d_tmp
         455      741888 :     real (kind=r8), dimension(fv_nphys,fv_nphys,num_trac) :: q_phys_tmp
         456             :     real (kind=r8), dimension(nc,nc)                      :: inv_darea_dp_fvm
         457             :     integer :: k,m_cnst,ie
         458             :
         459             :
         460             :
         461             :     !OMP BARRIER OMP MASTER needed
         462      370944 :     if (nc.ne.fv_nphys) then
         463           0 :       save_max_overlap = 4 !max number of mass overlap areas between phys and fvm grids
         464           0 :       allocate(save_air_mass_overlap(save_max_overlap,fv_nphys,fv_nphys,nlev,nets:nete))
         465           0 :       allocate(save_q_overlap(save_max_overlap,fv_nphys,fv_nphys,nlev,num_trac,nets:nete))
         466           0 :       allocate(save_q_phys(fv_nphys,fv_nphys,nlev,num_trac,nets:nete))
         467           0 :       allocate(save_dp_phys(fv_nphys,fv_nphys,nlev,nets:nete))
         468           0 :       allocate(save_overlap_area(save_max_overlap,fv_nphys,fv_nphys,nets:nete))
         469           0 :       allocate(save_num_overlap(fv_nphys,fv_nphys,nlev,nets:nete))
         470           0 :       save_num_overlap = 0
         471           0 :       allocate(save_overlap_idx(2,save_max_overlap,fv_nphys,fv_nphys,nets:nete))
         472             :     end if
         473             :
         474     2979144 :     do ie=nets,nete
         475    54772200 :       tmp = 1.0_r8
         476    33906600 :       inv_area  = 1.0_r8/dyn2phys(tmp,elem(ie)%metdet(:,:))
         477     5216400 :       phis_phys(:,ie) = RESHAPE(dyn2phys(elem(ie)%state%phis(:,:),elem(ie)%metdet(:,:),inv_area),SHAPE(phis_phys(:,ie)))
         478    26082000 :       ps_phys(:,ie) = ptop
         479     2979144 :       if (nc.ne.fv_nphys) then
         480           0 :         tmp = 1.0_r8
         481           0 :         do k=1,nlev
         482           0 :           inv_darea_dp_fvm = dyn2fvm(elem(ie)%state%dp3d(:,:,k,tl),elem(ie)%metdet(:,:))
         483           0 :           inv_darea_dp_fvm = 1.0_r8/inv_darea_dp_fvm
         484           0 :           T_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%state%T(:,:,k,tl),elem(ie)%metdet(:,:),inv_area),SHAPE(T_phys(:,k,ie)))
         485           0 :           Omega_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%derived%omega(:,:,k),elem(ie)%metdet(:,:),inv_area), &
         486           0 :                                        SHAPE(Omega_phys(:,k,ie)))
         487           0 :           call fvm2phys(ie,k,fvm(ie),fvm(ie)%c(:,:,k,:),q_phys_tmp,num_trac)
         488           0 :           dp3d_phys(:,k,ie) = RESHAPE(save_dp_phys(:,:,k,ie),SHAPE(dp3d_phys(:,k,ie)))
         489           0 :           ps_phys(:,ie) = ps_phys(:,ie)+RESHAPE(save_dp_phys(:,:,k,ie),SHAPE(ps_phys(:,ie)))
         490           0 :           do m_cnst=1,num_trac
         491           0 :             q_phys(:,k,m_cnst,ie) = RESHAPE(q_phys_tmp(:,:,m_cnst),SHAPE(q_phys(:,k,m_cnst,ie)))
         492             :           end do
         493             :         end do
         494             :       else
         495   245170800 :         do k=1,nlev
         496   242562600 :           dp3d_tmp       = dyn2phys(elem(ie)%state%dp3d(:,:,k,tl),elem(ie)%metdet(:,:),inv_area)
         497  3153313800 :           inv_darea_dp_phys = inv_area/dp3d_tmp
         498   727687800 :           T_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%state%T(:,:,k,tl)*elem(ie)%state%dp3d(:,:,k,tl),elem(ie)%metdet(:,:),&
         499  6064065000 :                inv_darea_dp_phys),SHAPE(T_phys(:,k,ie)))
         500   242562600 :           Omega_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%derived%OMEGA(:,:,k),elem(ie)%metdet(:,:),inv_area), &
         501   727687800 :                                        SHAPE(Omega_phys(:,k,ie)))
         502             :           !
         503             :           ! no mapping needed - just copy fields into physics structure
         504             :           !
         505   485125200 :           dp3d_phys(:,k,ie) = RESHAPE(fvm(ie)%dp_fvm(1:nc,1:nc,k),SHAPE(dp3d_phys(:,k,ie)))
         506  2668188600 :           ps_phys(:,ie) = ps_phys(:,ie)+RESHAPE(fvm(ie)%dp_fvm(1:nc,1:nc,k),SHAPE(ps_phys(:,ie)))
         507  2913359400 :           do m_cnst=1,num_trac
         508  5578939800 :             q_phys(:,k,m_cnst,ie) = RESHAPE(fvm(ie)%c(1:nc,1:nc,k,m_cnst),SHAPE(q_phys(:,k,m_cnst,ie)))
         509             :           end do
         510             :         end do
         511             :       end if
         512             :     end do
         513      370944 :   end subroutine dyn2phys_all_vars
         514             :
         515             :
         516   978075000 :   function dyn2phys(qdp_gll,metdet,inv_dp_darea_phys) result(qdp_phys)
         517             :     use dimensions_mod, only: np, nc, fv_nphys
         518             :     use derivative_mod, only: subcell_integration
         519             :     real (kind=r8), intent(in)           :: qdp_gll(np,np)
         520             :     real (kind=r8)                       :: qdp_phys(fv_nphys,fv_nphys)
         521             :     real (kind=r8), intent(in)           :: metdet(np,np)
         522             :     real (kind=r8), intent(in), optional :: inv_dp_darea_phys(fv_nphys,fv_nphys)
         523             :
         524   978075000 :     call subcell_integration(qdp_gll(:,:), np, fv_nphys, metdet,qdp_phys,nc.ne.fv_nphys)
         525   978075000 :     if (present(inv_dp_darea_phys)) &
         526 13620020400 :       qdp_phys = qdp_phys*inv_dp_darea_phys ! convert qdp to q
         527  2934225000 :   end function dyn2phys
         528             :
         529             :
         530     6037200 :   function dyn2fvm(qdp_gll,metdet,inv_dp_darea_phys,q_gll) result(qdp_phys)
         531   978075000 :     use dimensions_mod, only: np, nc
         532             :     use derivative_mod, only: subcell_integration
         533             :     real (kind=r8), intent(in)           :: qdp_gll(np,np)
         534             :     real (kind=r8), intent(in)           :: metdet(np,np)
         535             :     real (kind=r8), intent(in), optional :: inv_dp_darea_phys(nc,nc)
         536             :     real (kind=r8), intent(in), optional :: q_gll(np,np)
         537             :     real (kind=r8)                       :: qdp_phys(nc,nc), min_val, max_val
         538             :     integer                              :: i,j
         539             :
         540     6037200 :     call subcell_integration(qdp_gll(:,:), np, nc, metdet,qdp_phys)
         541     6037200 :     if (present(inv_dp_darea_phys)) then
         542             :       !
         543             :       ! convert qdp to q
         544             :       !
         545    78413400 :       qdp_phys = qdp_phys*inv_dp_darea_phys
         546             :       !
         547             :       ! simple limiter
         548             :       !
         549     6031800 :       if (present(q_gll)) then
         550   116008200 :         min_val = minval(q_gll)
         551   116008200 :         max_val = maxval(q_gll)
         552    22096800 :         do j = 1, nc
         553    71814600 :           do i = 1, nc
         554             :             !
         555             :             ! simple limiter: only coded for nc=3 and np4
         556             :             !
         557    66290400 :             qdp_phys(i,j) = max(min_val,min(max_val,qdp_phys(i,j)))
         558             :           end do
         559             :         end do
         560             :       end if
         561             :     end if
         562    12074400 :   end function dyn2fvm
         563             :
         564    26082000 :   function dyn2phys_vector(v_gll,elem) result(v_phys)
         565     6037200 :     use dimensions_mod, only: np, nlev,  fv_nphys
         566             :     use interpolate_mod,only: interpdata_t,interpolate_2d,interpolate_t
         567             :     use cube_mod       ,only: dmap
         568             :     use control_mod    ,only: cubed_sphere_map
         569             :
         570             :     type (interpdata_t):: interpdata
         571             :     type (element_t), intent(in)   :: elem
         572             :     type (interpolate_t) , target :: interp_p
         573             :     real (kind=r8), intent(in)           :: v_gll(np,np,2,nlev)
         574             :     real (kind=r8)                       :: v_phys(fv_nphys*fv_nphys,2,nlev)
         575             :
         576             :     integer :: i,j,k
         577             :
         578             :     ! Local variables
         579             :     real (kind=r8)    ::  fld_contra(np,np,2,nlev) ! vector field
         580             :
         581             :     real (kind=r8)    ::  v1,v2
         582    10432800 :     real (kind=r8)    ::  D(2,2,fv_nphys*fv_nphys)   ! derivative of gnomonic mapping
         583             :     !
         584             :     ! this could be done at initialization and does not need to be repeated
         585             :     !
         586     5216400 :     call setup_interpdata_for_gll_to_phys_vec_mapping(interpdata, interp_p)
         587             :     ! convert to contra
         588   490341600 :     do k=1,nlev
         589  2430842400 :       do j=1,np
         590 10187629200 :         do i=1,np
         591             :           ! latlon->contra
         592  7762003200 :           fld_contra(i,j,1,k) = elem%Dinv(i,j,1,1)*v_gll(i,j,1,k) + elem%Dinv(i,j,1,2)*v_gll(i,j,2,k)
         593  9702504000 :           fld_contra(i,j,2,k) = elem%Dinv(i,j,2,1)*v_gll(i,j,1,k) + elem%Dinv(i,j,2,2)*v_gll(i,j,2,k)
         594             :         enddo
         595             :       enddo
         596             :     end do
         597             :
         598   490341600 :     do k=1,nlev
         599  4856468400 :       do i=1,interpdata%n_interp
         600  4366126800 :         v_phys(i,1,k)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,1,k),interp_p,np)
         601  4851252000 :         v_phys(i,2,k)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,2,k),interp_p,np)
         602             :       end do
         603             :     end do
         604    52164000 :     do i=1,interpdata%n_interp
         605             :       ! convert fld from contra->latlon
         606           0 :       call dmap(D(:,:,i),interpdata%interp_xy(i)%x,interpdata%interp_xy(i)%y,&
         607    52164000 :            elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
         608             :     end do
         609   490341600 :     do k=1,nlev
         610  4856468400 :       do i=1,interpdata%n_interp
         611             :         ! convert fld from contra->latlon
         612  4366126800 :         v1 = v_phys(i,1,k)
         613  4366126800 :         v2 = v_phys(i,2,k)
         614             :
         615  4366126800 :         v_phys(i,1,k)=D(1,1,i)*v1 + D(1,2,i)*v2
         616  4851252000 :         v_phys(i,2,k)=D(2,1,i)*v1 + D(2,2,i)*v2
         617             :       end do
         618             :     end do
         619     5216400 :   end function dyn2phys_vector
         620             :
         621     5216400 :   subroutine setup_interpdata_for_gll_to_phys_vec_mapping(interpdata,interp_p)
         622             :     !
         623             :     ! initialize interpolation data structures to interpolate to phys grid
         624             :     ! using interpolate_mod subroutines
         625             :     !
         626             :     use interpolate_mod, only: interpolate_t, interpdata_t, interpolate_create
         627             :     use dimensions_mod, only : np
         628             :     use quadrature_mod, only : quadrature_t, gausslobatto
         629             :     use dimensions_mod, only : fv_nphys
         630             :     type (interpdata_t)  , intent(out)         :: interpdata
         631             :     type (interpolate_t) , intent(out), target :: interp_p
         632             :
         633             :     ! local
         634             :     type (quadrature_t)   :: gp_quadrature
         635             :     integer i,j,ioff,ngrid
         636             :     real (kind=r8) ::  dx
         637             :
         638     5216400 :     ngrid = fv_nphys*fv_nphys
         639     5216400 :     interpdata%n_interp=ngrid
         640             :     !
         641             :     ! initialize interpolation stuff related to basis functions
         642             :     !
         643     5216400 :     gp_quadrature = gausslobatto(np)
         644     5216400 :     call interpolate_create(gp_quadrature,interp_p)
         645    15649200 :     allocate(interpdata%interp_xy(ngrid))
         646    15649200 :     allocate(interpdata%ilat(ngrid) )
         647    10432800 :     allocate(interpdata%ilon(ngrid) )
         648             :     !
         649             :     !WARNING: THIS CODE INTERFERES WITH LAT-LON OUTPUT
         650             :     !         OF REGULAR SE IF nc>0
         651             :     !
         652     5216400 :     ioff=1
         653     5216400 :     dx = 2.0_r8/dble(fv_nphys)
         654    20865600 :     do j=1,fv_nphys
         655    67813200 :       do i=1,fv_nphys
         656    46947600 :         interpdata%interp_xy(ioff)%x = -1_r8+(i-0.5_r8)*dx
         657    46947600 :         interpdata%interp_xy(ioff)%y = -1_r8+(j-0.5_r8)*dx
         658    46947600 :         interpdata%ilon(ioff) = i
         659    46947600 :         interpdata%ilat(ioff) = j
         660    62596800 :         ioff=ioff+1
         661             :       enddo
         662             :     enddo
         663     5216400 :   end subroutine setup_interpdata_for_gll_to_phys_vec_mapping
         664             :
         665             :
         666 >15958*10^7 :   function lagrange_1d(src_grid,src_val,ngrid,dst_point,iwidth) result(val)
         667             :     integer              , intent(in)  :: ngrid,iwidth
         668             :     real (kind=r8), intent(in)  :: src_grid(ngrid), src_val(ngrid)
         669             :     real (kind=r8)              :: val
         670             :
         671             :     real (kind=r8), intent(in)  :: dst_point
         672             :
         673             :     integer :: iref, j,k
         674 >31917*10^7 :     real (kind=r8)              :: w(ngrid)
         675             :
         676 >15958*10^7 :     if (dst_point.LE.src_grid(1)) then
         677             :       iref=1
         678             :     else
         679             :       iref=1
         680 >81931*10^7 :       do while (dst_point>src_grid(iref))
         681 >65974*10^7 :         iref = iref + 1
         682 >81931*10^7 :         if (iref>ngrid) then
         683             :           exit
         684             :         end if
         685             :       end do
         686 >15957*10^7 :       iref=iref-1
         687             :     end if
         688             :
         689 >15958*10^7 :     iref=MIN(MAX(iref,iwidth),ngrid-iwidth)
         690             :
         691 >14790*10^8 :     w = 1.0_r8
         692 >79794*10^7 :     do j=iref-(iwidth-1),iref+iwidth
         693 >33513*10^8 :       do k=iref-(iwidth-1),iref+iwidth
         694 >31917*10^8 :         if (k.ne.j) then
         695 >19150*10^8 :           w(j)=w(j)*(dst_point-src_grid(k))/(src_grid(j)-src_grid(k))
         696             :         end if
         697             :       end do
         698             :     end do
         699             :
         700             :     val=0.0_r8
         701 >79794*10^7 :     do j=iref-(iwidth-1),iref+iwidth
         702 >79794*10^7 :       val=val+w(j)*src_val(j)
         703             :     end do
         704 >15958*10^7 :   end function lagrange_1d
         705             :
         706     7792200 :   subroutine tensor_lagrange_interp(cubeboundary,np,nc,nhc,num_lev,nflds,psi,interp_value,llimiter,iwidth,norm_elem_coord)
         707             :     use control_mod, only : north, south, east, west, neast, nwest, seast, swest
         708             :     implicit none
         709             :
         710             :     integer              , intent(in)    :: cubeboundary,nc, np, iwidth,nhc,num_lev,nflds
         711             :     logical              , intent(in)    :: llimiter(nflds)                   !apply limiter
         712             :     real (kind=r8), intent(inout) :: psi(1-nhc:nc+nhc,1-nhc:nc+nhc,num_lev,nflds) !fvm grid values with filled halo
         713             :     real (kind=r8), intent(out)   :: interp_value(np,np,num_lev,nflds)            !interpolated field
         714             :     real (kind=r8), intent(in)    :: norm_elem_coord(2,1-nhc:nc+nhc,1-nhc:nc+nhc)
         715    15584400 :     integer :: which_nc_cell(np)
         716             :
         717    15584400 :     real (kind=r8):: dx,gll_points(np)
         718    15584400 :     real (kind=r8):: nc_points(1-nc:nc+nc)
         719             :
         720    15584400 :     real (kind=r8):: value(1-iwidth:nc+iwidth)
         721    15584400 :     real (kind=r8):: val_tmp(1-nhc:nc+nhc,1-nhc:nc+nhc)
         722             :
         723    15584400 :     real (kind=r8):: min_value(np,np,num_lev,nflds), max_value(np,np,num_lev,nflds)
         724             :
         725    15584400 :     integer :: imin(1-nhc:nc+nhc), imax(1-nhc:nc+nhc)
         726             :     integer :: k,i,j,isearch,igll,jgll,jrow,h,irow,itr
         727             :
         728     7792200 :     gll_points(1) = -1.0_r8
         729     7792200 :     gll_points(2) = -sqrt(1.0_r8/5.0_r8)
         730     7792200 :     gll_points(3) =  sqrt(1.0_r8/5.0_r8)
         731     7792200 :     gll_points(4) =  1.0_r8
         732             :
         733     7792200 :     dx = 2_r8/dble(nc)
         734    77922000 :     do k=1-nc,2*nc
         735    77922000 :       nc_points(k) = -1.0_r8+dx*0.5_r8+dble(k-1)*dx
         736             :     end do
         737             :     !
         738             :     ! find fvm point surrounding gll points for simple limiter
         739             :     !
         740    38961000 :     do k=1,np
         741    77922000 :       do isearch=0,nc+1
         742    77922000 :         if (nc_points(isearch)<gll_points(k).and.nc_points(isearch+1).ge.gll_points(k)) exit
         743             :       end do
         744    38961000 :       which_nc_cell(k)=isearch
         745             :     end do
         746    46753200 :     do itr=1,nflds
         747    46753200 :       if (llimiter(itr)) then
         748             :         !
         749             :         ! fill non-existent halo cells for limiter
         750             :         !
         751           0 :         if (cubeboundary>4) then
         752           0 :           h=1
         753           0 :           select case(cubeboundary)
         754             :           case (nwest)
         755           0 :             psi(0,nc+h  ,:,itr) = psi(1-h,nc  ,:,itr)
         756           0 :             psi(1-h,nc+1,:,itr) = psi(1  ,nc+h,:,itr)
         757             :           case (swest)
         758           0 :             psi(1-h,0,:,itr) = psi(1,1-h,:,itr)
         759           0 :             psi(0,1-h,:,itr) = psi(1-h,1,:,itr)
         760             :           case (seast)
         761           0 :             psi(nc+h,0,:,itr) = psi(nc,1-h,:,itr)
         762           0 :             psi(nc+1,1-h,:,itr) = psi(nc+h,1,:,itr)
         763             :           case (neast)
         764           0 :             psi(nc+h,nc+1,:,itr) = psi(nc,nc+h,:,itr)
         765           0 :             psi(nc+1,nc+h,:,itr) = psi(nc+h,nc,:,itr)
         766             :           end select
         767             :         end if
         768           0 :         do k=1,num_lev
         769           0 :           do j=1,np
         770           0 :             do i=1,np
         771           0 :               max_value(i,j,k,itr) = max(&
         772           0 :                    psi(which_nc_cell(i)  ,which_nc_cell(j)  ,k,itr),&
         773           0 :                    psi(which_nc_cell(i)+1,which_nc_cell(j)  ,k,itr),&
         774           0 :                    psi(which_nc_cell(i)  ,which_nc_cell(j)+1,k,itr),&
         775             :                    psi(which_nc_cell(i)+1,which_nc_cell(j)+1,k,itr) &
         776           0 :                    )
         777             :               min_value(i,j,k,itr) = min(&
         778           0 :                    psi(which_nc_cell(i)  ,which_nc_cell(j)  ,k,itr),&
         779           0 :                    psi(which_nc_cell(i)+1,which_nc_cell(j)  ,k,itr),&
         780           0 :                    psi(which_nc_cell(i)  ,which_nc_cell(j)+1,k,itr),&
         781             :                    psi(which_nc_cell(i)+1,which_nc_cell(j)+1,k,itr) &
         782           0 :                    )
         783             :             end do
         784             :           end do
         785             :         end do
         786             :       end if
         787             :     end do
         788             :
         789    77922000 :     imin=1-nhc
         790    77922000 :     imax=nc+nhc
         791             :     !
         792             :     ! special corner treatment
         793             :     !
         794     7792200 :     if (cubeboundary==swest) then
         795       51948 :       do itr=1,nflds
         796     4077918 :         do k=1,num_lev
         797    24155820 :           do jrow=1,nc+iwidth
         798             :             !
         799             :             ! cubic along constant x (i=irow) in west halo to fvm points in halo
         800             :             !
         801    64415520 :             do irow=1-iwidth,0
         802   161038800 :               val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1:nc+nhc),psi(irow,1:nc+nhc,k,itr),nc+nhc,&
         803   704544750 :                    norm_elem_coord(2,1,jrow),iwidth)
         804             :             end do
         805             :           end do
         806    64458810 :           psi(1-iwidth:0,1:nc+iwidth,k,itr) = val_tmp(1-iwidth:0,1:nc+iwidth)
         807             :         enddo
         808             :       end do
         809       34632 :       imin(1-nhc:0) = 1
         810             :     end if
         811     7792200 :     if (cubeboundary==nwest) then
         812       51948 :       do itr=1,nflds
         813     4077918 :         do k=1,num_lev
         814    24155820 :           do jrow=1-iwidth,nc
         815             :             !
         816             :             ! cubic along constant x (i=irow) in west halo to fvm points in halo
         817             :             !
         818    64415520 :             do irow=1-iwidth,0
         819   161038800 :               val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1-nhc:nc),psi(irow,1-nhc:nc,k,itr),nc+nhc,&
         820   704544750 :                    norm_elem_coord(2,1,jrow),iwidth)
         821             :             end do
         822             :           end do
         823    64458810 :           psi(1-iwidth:0,1-iwidth:nc,k,itr) = val_tmp(1-iwidth:0,1-iwidth:nc)
         824             :         end do
         825             :       end do
         826       34632 :       imin(nc+1:nc+nhc) = 1
         827             :     end if
         828             :
         829     7792200 :     if (cubeboundary==seast) then
         830       51948 :       do itr=1,nflds
         831     4077918 :         do k=1,num_lev
         832    24155820 :           do jrow=1,nc+iwidth
         833   161038800 :             value=0.0_r8
         834             :             !
         835             :             ! cubic along constant y in ease halo to fvm points in halo
         836             :             !
         837    64415520 :             do irow=nc+1,nc+iwidth
         838   161038800 :               val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1:nc+nhc),psi(irow,1:nc+nhc,k,itr),nc+nhc,&
         839   704544750 :                    norm_elem_coord(2,1,jrow),iwidth)
         840             :             end do
         841             :           end do
         842    64458810 :           psi(nc+1:nc+iwidth,1:nc+iwidth,k,itr) = val_tmp(nc+1:nc+iwidth,1:nc+iwidth)
         843             :         end do
         844             :       end do
         845       34632 :       imax(1-nhc:0) = nc
         846             :     end if
         847             :
         848     7792200 :     if (cubeboundary==neast) then
         849       51948 :       do itr=1,nflds
         850     4077918 :         do k=1,num_lev
         851    24155820 :           do jrow=1-iwidth,nc
         852             :             !
         853             :             ! cubic along constant y in ease halo to fvm points in halo
         854             :             !
         855    64415520 :             do irow=nc+1,nc+iwidth
         856   161038800 :               val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1-nhc:nc),psi(irow,1-nhc:nc,k,itr),nc+nhc,&
         857   704544750 :                    norm_elem_coord(2,1,jrow),iwidth)
         858             :             end do
         859             :           end do
         860    64458810 :           psi(nc+1:nc+iwidth,1-iwidth:nc,k,itr) = val_tmp(nc+1:nc+iwidth,1-iwidth:nc)
         861             :         end do
         862             :       end do
         863       34632 :       imax(nc+1:nc+nhc) = nc
         864             :     end if
         865             :     !
         866             :     ! mapping
         867             :     !
         868             :     !
         869             :     if (cubeboundary==0.or.cubeboundary==north.or.cubeboundary==south.or.&
         870             :          cubeboundary==swest.or.cubeboundary==nwest.or.&
         871     7792200 :          cubeboundary==seast.or.cubeboundary==neast) then
         872    43844112 :       do itr=1,nflds
         873  3441762792 :         do k=1,num_lev
         874 17026130160 :           do igll=1,np
         875             :             !
         876             :             ! cubic along constant y (j=jrow)
         877             :             !
         878 >10873*10^7 :             do jrow=1-iwidth,nc+iwidth
         879 95141723040 :               value(jrow) = lagrange_1d(norm_elem_coord(1,imin(jrow):imax(jrow),jrow),psi(imin(jrow):imax(jrow),jrow,k,itr),&
         880 >10597*10^8 :                    imax(jrow)-imin(jrow)+1,gll_points(igll),iwidth)
         881             :             end do
         882 71356292280 :             do jgll=1,np
         883 >27183*10^7 :               interp_value(igll,jgll,k,itr) = lagrange_1d(norm_elem_coord(2,1,1-iwidth:nc+iwidth),value,nc+2*iwidth,&
         884 >72035*10^7 :                    gll_points(jgll),iwidth)
         885             :             end do
         886             :           end do
         887             :         end do
         888             :       end do
         889      484848 :     else if (cubeboundary==east.or.cubeboundary==west) then
         890     2909088 :       do itr=1,nflds
         891   228363408 :         do k=1,num_lev
         892  1129695840 :           do jgll=1,np
         893             :             !
         894             :             ! cubic along constant x (i=irow)
         895             :             !
         896  7214538240 :             do irow=1-iwidth,nc+iwidth
         897 31563604800 :               value(irow) = lagrange_1d(norm_elem_coord(2,irow,1-nhc:nc+nhc),psi(irow,1-nhc:nc+nhc,k,itr),nc+2*nhc,&
         898 >15240*10^7 :                    gll_points(jgll),iwidth)
         899             :             end do
         900  4734540720 :             do igll=1,np
         901 18036345600 :               interp_value(igll,jgll,k,itr) = lagrange_1d(norm_elem_coord(1,1-iwidth:nc+iwidth,1),value,nc+2*iwidth,&
         902 47796315840 :                    gll_points(igll),iwidth)
         903             :             end do
         904             :           end do
         905             :         end do
         906             :       end do
         907             :     end if
         908    46753200 :     do itr=1,nflds
         909    46753200 :       if (llimiter(itr)) then
         910           0 :         do k=1,num_lev
         911           0 :           do j=1,np
         912           0 :             do i=1,np
         913           0 :               interp_value(i,j,k,itr)=max(min_value(i,j,k,itr),min(max_value(i,j,k,itr),interp_value(i,j,k,itr)))
         914             :             end do
         915             :           enddo
         916             :         end do
         917             :       end if
         918             :     end do
         919     7792200 :   end subroutine tensor_lagrange_interp
         920             :
         921           0 :   subroutine fvm2phys(ie,k,fvm,q_fvm,q_phys,num_trac)
         922             :     use dimensions_mod, only: nc,nhc,fv_nphys
         923             :     !
         924             :     ! weights must be initialized in fvm2phys_init before using these functions
         925             :     !
         926             :     type(fvm_struct)     , intent(inout)        :: fvm
         927             :     integer              , intent(in)           :: ie,k
         928             :     integer              , intent(in)           :: num_trac
         929             :
         930             :     real (kind=r8), intent(inout)        :: q_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,num_trac)
         931             :     real (kind=r8), intent(out)          :: q_phys(fv_nphys,fv_nphys,num_trac)
         932             :
         933             :     real (kind=r8)                       :: recons    (irecons_tracer,1:nc,1:nc,1)
         934             :
         935             :     integer                              :: jx,jy
         936             :
         937           0 :     call get_dp_overlap_save(ie,k,fvm,recons)
         938           0 :     do jy=1,fv_nphys
         939           0 :       do jx=1,fv_nphys
         940           0 :         save_dp_phys(jx,jy,k,ie) = SUM(save_air_mass_overlap(1:save_num_overlap(jx,jy,k,ie),jx,jy,k,ie))
         941             :       end do
         942             :     end do
         943           0 :     call get_q_overlap_save(ie,k,fvm,q_fvm,num_trac,q_phys)
         944           0 :     save_dp_phys(:,:,k,ie) = save_dp_phys(:,:,k,ie)/fvm%area_sphere_physgrid
         945           0 :   end subroutine fvm2phys
         946             :
         947             :
         948           0 :   subroutine phys2fvm(ie,k,fvm,fq_phys,fqdp_fvm,num_trac)
         949             :     use dimensions_mod, only: nhc_phys,fv_nphys,nc
         950             :     integer              , intent(in)           :: ie,k
         951             :     type(fvm_struct)     , intent(inout)        :: fvm
         952             :     integer              , intent(in)           :: num_trac
         953             :     real (kind=r8), intent(inout)        :: fq_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_trac)
         954             :     real (kind=r8), intent(out)          :: fqdp_fvm (nc,nc,num_trac)
         955             :
         956             :
         957             :     integer                              :: h,jx,jy,jdx,jdy,m_cnst
         958             :
         959           0 :     real(kind=r8), dimension(fv_nphys,fv_nphys) :: phys_cdp_max, phys_cdp_min
         960             :
         961             :     integer       :: num
         962             :     real(kind=r8) :: tmp,sum_dq_min,sum_dq_max,fq
         963             :
         964           0 :     real(kind=r8) :: mass_phys(fv_nphys,fv_nphys)
         965             :     real(kind=r8) :: min_patch,max_patch,gamma
         966             :     real (kind=r8):: q_prev,mass_forcing,mass_forcing_phys
         967             :
         968           0 :     real(kind=r8), allocatable, dimension(:,:,:) :: dq_min_overlap,dq_max_overlap
         969           0 :     real(kind=r8), allocatable, dimension(:,:,:) :: dq_overlap
         970           0 :     real(kind=r8), allocatable, dimension(:,:,:) :: fq_phys_overlap
         971             :
         972           0 :     allocate(dq_min_overlap       (save_max_overlap,fv_nphys,fv_nphys))
         973           0 :     allocate(dq_max_overlap       (save_max_overlap,fv_nphys,fv_nphys))
         974           0 :     allocate(dq_overlap           (save_max_overlap,fv_nphys,fv_nphys))
         975           0 :     allocate(fq_phys_overlap      (save_max_overlap,fv_nphys,fv_nphys))
         976             :
         977           0 :     do m_cnst=1,num_trac
         978           0 :       fqdp_fvm(:,:,m_cnst) = 0.0_r8
         979             :       call get_fq_overlap(ie,k,fvm,&
         980           0 :            fq_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,m_cnst),save_max_overlap,&
         981           0 :            fq_phys_overlap,1)
         982           0 :       mass_phys(1:fv_nphys,1:fv_nphys) = fq_phys(1:fv_nphys,1:fv_nphys,m_cnst)*&
         983           0 :            (save_dp_phys(1:fv_nphys,1:fv_nphys,k,ie)*fvm%area_sphere_physgrid)
         984             :
         985           0 :       min_patch = MINVAL(fvm%c(0:nc+1,0:nc+1,k,m_cnst))
         986           0 :       max_patch = MAXVAL(fvm%c(0:nc+1,0:nc+1,k,m_cnst))
         987           0 :       do jy=1,fv_nphys
         988           0 :         do jx=1,fv_nphys
         989           0 :           num = save_num_overlap(jx,jy,k,ie)
         990             : #ifdef debug_coupling
         991             :           save_q_overlap(:,jx,jy,k,m_cnst,ie) = 0.0_r8
         992             :           save_q_phys(jx,jy,k,m_cnst,ie)      = 0.0_r8
         993             :           tmp = save_q_phys(jx,jy,k,m_cnst,ie)+fq_phys(jx,jy,m_cnst) !updated physics grid mixing ratio
         994             :           phys_cdp_max(jx,jy)= MAX(max_patch,tmp)
         995             :           phys_cdp_min(jx,jy)= MIN(min_patch,tmp)
         996             : #else
         997           0 :           tmp = save_q_phys(jx,jy,k,m_cnst,ie)+fq_phys(jx,jy,m_cnst) !updated physics grid mixing ratio
         998           0 :           phys_cdp_max(jx,jy)= MAX(MAX(MAXVAL(save_q_overlap(1:num,jx,jy,k,m_cnst,ie)),max_patch),tmp)
         999           0 :           phys_cdp_min(jx,jy)= MIN(MIN(MINVAL(save_q_overlap(1:num,jx,jy,k,m_cnst,ie)),min_patch),tmp)
        1000             : #endif
        1001             :           !
        1002             :           ! add high-order fq change when it does not violate monotonicity
        1003             :           !
        1004           0 :           mass_forcing_phys = 0.0_r8
        1005           0 :           do h=1,num
        1006           0 :             jdx = save_overlap_idx(1,h,jx,jy,ie); jdy = save_overlap_idx(2,h,jx,jy,ie)
        1007           0 :             q_prev = save_q_overlap(h,jx,jy,k,m_cnst,ie)
        1008             : #ifndef skip_high_order_fq_map
        1009           0 :             save_q_overlap(h,jx,jy,k,m_cnst,ie) = save_q_overlap(h,jx,jy,k,m_cnst,ie)+fq_phys_overlap(h,jx,jy)
        1010           0 :             save_q_overlap(h,jx,jy,k,m_cnst,ie) = MIN(save_q_overlap(h,jx,jy,k,m_cnst,ie),phys_cdp_max(jx,jy))
        1011           0 :             save_q_overlap(h,jx,jy,k,m_cnst,ie) = MAX(save_q_overlap(h,jx,jy,k,m_cnst,ie),phys_cdp_min(jx,jy))
        1012           0 :             mass_forcing = (save_q_overlap(h,jx,jy,k,m_cnst,ie)-q_prev)*save_air_mass_overlap(h,jx,jy,k,ie)
        1013           0 :             mass_forcing_phys = mass_forcing_phys + mass_forcing
        1014           0 :             fqdp_fvm(jdx,jdy,m_cnst) = fqdp_fvm(jdx,jdy,m_cnst)+mass_forcing
        1015             : #endif
        1016             :             !
        1017             :             ! prepare for mass fixing algorithm
        1018             :             !
        1019           0 :             dq_min_overlap(h,jx,jy)   = save_q_overlap(h,jx,jy,k,m_cnst,ie)-phys_cdp_min(jx,jy)
        1020           0 :             dq_max_overlap  (h,jx,jy) = save_q_overlap(h,jx,jy,k,m_cnst,ie)-phys_cdp_max(jx,jy)
        1021             :           end do
        1022           0 :           mass_phys(jx,jy) = mass_phys(jx,jy) -mass_forcing_phys
        1023             :         end do
        1024             :       end do
        1025             :       !
        1026             :       ! let physics mass tendency remove excess mass (as defined above) first proportional to how much is availabe
        1027             :       !
        1028             : #ifdef mass_fix
        1029           0 :       do jy=1,fv_nphys
        1030           0 :         do jx=1,fv_nphys
        1031             :           !
        1032             :           ! total mass change from physics on physics grid
        1033             :           !
        1034           0 :           num = save_num_overlap(jx,jy,k,ie)
        1035           0 :           fq = mass_phys(jx,jy)/(fvm%area_sphere_physgrid(jx,jy)*save_dp_phys(jx,jy,k,ie))
        1036           0 :           if (fq<0.0_r8) then
        1037           0 :             sum_dq_min = SUM(dq_min_overlap(1:num,jx,jy)*save_air_mass_overlap(1:num,jx,jy,k,ie))
        1038           0 :             if (sum_dq_min>1.0E-14_r8) then
        1039           0 :               gamma=mass_phys(jx,jy)/sum_dq_min
        1040           0 :               do h=1,num
        1041           0 :                 jdx = save_overlap_idx(1,h,jx,jy,ie); jdy = save_overlap_idx(2,h,jx,jy,ie)
        1042           0 :                 fqdp_fvm(jdx,jdy,m_cnst) = fqdp_fvm(jdx,jdy,m_cnst)&
        1043           0 :                      +gamma*dq_min_overlap(h,jx,jy)*save_air_mass_overlap(h,jx,jy,k,ie)
        1044             :               end do
        1045             :             end if
        1046             :           end if
        1047             :
        1048           0 :           if (fq>0.0_r8) then
        1049           0 :             sum_dq_max = SUM(dq_max_overlap(1:num,jx,jy)*save_air_mass_overlap(1:num,jx,jy,k,ie))
        1050           0 :             if (sum_dq_max<-1.0E-14_r8) then
        1051           0 :               gamma=mass_phys(jx,jy)/sum_dq_max
        1052           0 :               do h=1,num
        1053           0 :                 jdx = save_overlap_idx(1,h,jx,jy,ie); jdy = save_overlap_idx(2,h,jx,jy,ie)
        1054           0 :                 fqdp_fvm(jdx,jdy,m_cnst) = fqdp_fvm(jdx,jdy,m_cnst)&
        1055           0 :                      +gamma*dq_max_overlap(h,jx,jy)*save_air_mass_overlap(h,jx,jy,k,ie)
        1056             :               end do
        1057             :             end if
        1058             :           end if
        1059             :         end do
        1060             :       end do
        1061             : #endif
        1062             :       !
        1063             :       ! convert to mass per unit area
        1064             :       !
        1065           0 :       fqdp_fvm(:,:,m_cnst) = fqdp_fvm(:,:,m_cnst)*fvm%inv_area_sphere(:,:)
        1066             :     end do
        1067           0 :     deallocate(dq_min_overlap)
        1068           0 :     deallocate(dq_max_overlap)
        1069           0 :     deallocate(fq_phys_overlap)
        1070           0 :   end subroutine phys2fvm
        1071             :
        1072             :
        1073           0 :   subroutine get_dp_overlap_save(ie,k,fvm,recons)
        1074             :     use dimensions_mod, only: nc,nhr,nhc
        1075             :     !
        1076             :     ! weights must be initialized in fvm2phys_init before using these functions
        1077             :     !
        1078             :     use dp_mapping, only: weights_all_fvm2phys, weights_eul_index_all_fvm2phys
        1079             :     use dp_mapping, only: weights_lgr_index_all_fvm2phys, jall_fvm2phys
        1080             :     !
        1081             :     ! setting nhe=0 because we do not need reconstruction outside of element
        1082             :     !
        1083             :     integer, parameter :: nh = nhr!+(nhe-1) ! = 2 (nhr=2; nhe_local=1),! = 3 (nhr=2; nhe_local=2)
        1084             :
        1085             :     type(fvm_struct)                                         , intent(inout):: fvm
        1086             :     integer                                                  , intent(in)   :: ie, k
        1087             :
        1088             :     real (kind=r8)                                           , intent(out)  :: recons    (irecons_tracer,nc,nc)
        1089             :     logical                              :: llimiter(1)
        1090             :     integer                              :: h,jx,jy,jdx,jdy,idx
        1091           0 :     llimiter=.false.
        1092           0 :     call get_fvm_recons(fvm,fvm%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,k),recons,1,llimiter)
        1093             :
        1094           0 :      do h=1,jall_fvm2phys(ie)
        1095           0 :        jx  = weights_lgr_index_all_fvm2phys(h,1,ie); jy  = weights_lgr_index_all_fvm2phys(h,2,ie)
        1096           0 :        jdx = weights_eul_index_all_fvm2phys(h,1,ie); jdy = weights_eul_index_all_fvm2phys(h,2,ie)
        1097           0 :        save_num_overlap(jx,jy,k,ie) = save_num_overlap(jx,jy,k,ie)+1!could be pre-computed
        1098           0 :        idx = save_num_overlap(jx,jy,k,ie)
        1099           0 :        save_overlap_idx(1,idx,jx,jy,ie) = jdx; save_overlap_idx(2,idx,jx,jy,ie) = jdy;
        1100           0 :        save_overlap_area(idx,jx,jy,ie)  = weights_all_fvm2phys(h,1,ie)
        1101           0 :        save_air_mass_overlap(idx,jx,jy,k,ie) = SUM(weights_all_fvm2phys(h,:,ie)*recons(:,jdx,jdy))
        1102             : #ifdef PCoM
        1103             :        save_air_mass_overlap(idx,jx,jy,k,ie) = fvm%dp_fvm(jdx,jdy,k)*weights_all_fvm2phys(h,1,ie)!PCoM
        1104             : #endif
        1105             :      end do
        1106             :
        1107           0 :    end subroutine get_dp_overlap_save
        1108             :
        1109             :
        1110           0 :   subroutine get_fq_overlap(ie,k,fvm,fq_phys,max_overlap,fq_phys_overlap,num_trac)
        1111             :     use dimensions_mod, only: fv_nphys, nhc_phys, nc
        1112             :     use dp_mapping, only: weights_lgr_index_all_fvm2phys, jall_fvm2phys
        1113             :     use dp_mapping, only: weights_eul_index_all_fvm2phys
        1114             :     use dp_mapping, only: weights_lgr_index_all_phys2fvm, weights_eul_index_all_phys2fvm,jall_phys2fvm
        1115             :     use dp_mapping, only: weights_all_phys2fvm
        1116             :
        1117             :     integer              , intent(in)           :: ie,k
        1118             :     type(fvm_struct)     , intent(in)           :: fvm
        1119             :     integer              , intent(in)           :: num_trac, max_overlap
        1120             :     real(kind=r8), dimension(max_overlap,fv_nphys,fv_nphys,num_trac),intent(out) :: fq_phys_overlap
        1121             :
        1122             :     real (kind=r8), intent(inout) :: fq_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_trac)
        1123           0 :     real (kind=r8)                :: recons_q  (irecons_tracer,fv_nphys,fv_nphys,num_trac)
        1124           0 :     integer                       :: num_overlap(fv_nphys,fv_nphys)
        1125           0 :     logical                       :: llimiter_q(num_trac)
        1126             :     integer                       :: h,jx,jy,m_cnst,jdx,jdy
        1127             :     real (kind=r8)                :: dp_tmp
        1128             :     integer                       :: idx,jxx,jyy,jdxx,jdyy,hh
        1129           0 :     real(kind=r8)                 :: weights_all_phys2fvm_local((nc+fv_nphys)**2,irecons_tracer)
        1130             :     !
        1131             :     ! could be pre-computed
        1132             :     !
        1133           0 :     do h=1,jall_fvm2phys(ie)
        1134           0 :       jx   = weights_lgr_index_all_fvm2phys(h,1,ie)
        1135           0 :       jy   = weights_lgr_index_all_fvm2phys(h,2,ie)
        1136           0 :       jdx  = weights_eul_index_all_fvm2phys(h,1,ie)
        1137           0 :       jdy  = weights_eul_index_all_fvm2phys(h,2,ie)
        1138           0 :       do hh=1,jall_phys2fvm(ie)
        1139           0 :         jxx  = weights_lgr_index_all_phys2fvm(hh,1,ie)
        1140           0 :         jyy  = weights_lgr_index_all_phys2fvm(hh,2,ie)
        1141           0 :         jdxx  = weights_eul_index_all_phys2fvm(hh,1,ie)
        1142           0 :         jdyy  = weights_eul_index_all_phys2fvm(hh,2,ie)
        1143           0 :         if (jx==jdxx.and.jy==jdyy.and.jdx==jxx.and.jdy==jyy) then
        1144           0 :           weights_all_phys2fvm_local(h,:) = weights_all_phys2fvm(hh,:,ie)
        1145             :           exit
        1146             :         end if
        1147             :       end do
        1148             :     end do
        1149             :
        1150           0 :     llimiter_q=.false.
        1151           0 :     call get_physgrid_recons(fvm,fq_phys,recons_q,num_trac,llimiter_q)
        1152             :     !
        1153             :     ! q-dp coupling as described in equation (55) in Appendinx B of
        1154             :     ! Nair and Lauritzen, 2010: A Class of Deformational Flow Test Cases for Linear Transport Problems on the Sphere.
        1155             :     ! J. Comput. Phys.: Vol. 229, Issue 23, pp. 8868-8887, DOI:10.1016/j.jcp.2010.08.014.
        1156             :     !
        1157           0 :     num_overlap = 0
        1158           0 :     do h=1,jall_fvm2phys(ie)
        1159           0 :        jx  = weights_lgr_index_all_fvm2phys(h,1,ie)
        1160           0 :        jy  = weights_lgr_index_all_fvm2phys(h,2,ie)
        1161           0 :        jdx  = weights_eul_index_all_fvm2phys(h,1,ie)
        1162           0 :        jdy  = weights_eul_index_all_fvm2phys(h,2,ie)
        1163           0 :        num_overlap(jx,jy) = num_overlap(jx,jy)+1
        1164           0 :        idx = num_overlap(jx,jy)
        1165           0 :        dp_tmp = save_air_mass_overlap(idx,jx,jy,k,ie)-fvm%dp_fvm(jdx,jdy,k)*weights_all_phys2fvm_local(h,1)
        1166           0 :        do m_cnst=1,num_trac
        1167           0 :          fq_phys_overlap(idx,jx,jy,m_cnst) = &
        1168             :               (fvm%dp_fvm(jdx,jdy,k)*SUM(weights_all_phys2fvm_local(h,:)*recons_q(:,jx,jy,m_cnst))+&
        1169           0 :               fq_phys(jx,jy,m_cnst)*dp_tmp)/save_air_mass_overlap(idx,jx,jy,k,ie)
        1170             :        end do
        1171             :      end do
        1172           0 :   end subroutine get_fq_overlap
        1173             :
        1174           0 :   subroutine get_physgrid_recons(fvm,field_phys,recons_phys,num_trac,llimiter)
        1175             :     use dimensions_mod, only: fv_nphys,nhr_phys,nhc_phys,ns_phys
        1176             :     use fvm_reconstruction_mod, only: reconstruction
        1177             :     type(fvm_struct), intent(in)           :: fvm
        1178             :     integer,          intent(in)           :: num_trac
        1179             :     real (kind=r8),   intent(inout)        :: field_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,1,num_trac)
        1180             :     real (kind=r8),   intent(out)          :: recons_phys(irecons_tracer,1:fv_nphys,1:fv_nphys,num_trac)
        1181             :     logical,          intent(in)           :: llimiter(num_trac)
        1182             :
        1183             :     integer, dimension(3)                  :: jx_min_local, jx_max_local, jy_min_local, jy_max_local
        1184             :
        1185           0 :     jx_min_local(1) = 1            ; jx_max_local(1) = fv_nphys+1
        1186           0 :     jy_min_local(1) = 1            ; jy_max_local(1) = fv_nphys+1
        1187           0 :     jx_min_local(2) = 0            ; jx_max_local(2) = -1
        1188           0 :     jy_min_local(2) = 0            ; jy_max_local(2) = -1
        1189           0 :     jx_min_local(3) = 0            ; jx_max_local(3) = -1
        1190           0 :     jy_min_local(3) = 0            ; jy_max_local(3) = -1
        1191             :
        1192             :     call reconstruction(field_phys,1,1,recons_phys,irecons_tracer,llimiter,num_trac,&
        1193             :        fv_nphys,0,nhr_phys,nhc_phys,nhr_phys,ns_phys,nhr_phys,&
        1194             :        jx_min_local,jx_max_local,jy_min_local,jy_max_local,&
        1195           0 :        fvm%cubeboundary,fvm%halo_interp_weight_physgrid(1:ns_phys,1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,:),&
        1196           0 :        fvm%ibase_physgrid(1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,:),&
        1197           0 :        fvm%spherecentroid_physgrid(:,1:fv_nphys,1:fv_nphys),&
        1198           0 :        fvm%recons_metrics_physgrid(:,1:fv_nphys,1:fv_nphys),&
        1199           0 :        fvm%recons_metrics_integral_physgrid(:,1:fv_nphys,1:fv_nphys)    ,&
        1200             :        fvm%rot_matrix_physgrid,&
        1201           0 :        fvm%centroid_stretch_physgrid(1:7,1:fv_nphys,1:fv_nphys),&
        1202           0 :        fvm%vertex_recons_weights_physgrid(:,1:irecons_tracer-1,1:fv_nphys,1:fv_nphys),&
        1203           0 :        fvm%vtx_cart_physgrid(:,:,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
        1204           0 :   end subroutine get_physgrid_recons
        1205             :
        1206           0 :   subroutine get_fvm_recons(fvm,field_fvm,recons_fvm,num_trac,llimiter)
        1207           0 :     use dimensions_mod, only: nc,nhr,nhc,ns
        1208             :     use fvm_reconstruction_mod, only: reconstruction
        1209             :
        1210             :     type(fvm_struct), intent(in)   :: fvm
        1211             :     integer,          intent(in)   :: num_trac
        1212             :     logical,          intent(in)   :: llimiter(num_trac)
        1213             :     real (kind=r8),   intent(inout):: field_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,1,num_trac)
        1214             :     real (kind=r8),   intent(out)  :: recons_fvm(irecons_tracer,nc,nc,num_trac)
        1215             :     integer                        :: jx_min_local(3), jx_max_local(3), jy_min_local(3), jy_max_local(3)
        1216             :
        1217           0 :     jx_min_local(1) = 1            ; jx_max_local(1) = nc+1
        1218           0 :     jy_min_local(1) = 1            ; jy_max_local(1) = nc+1
        1219           0 :     jx_min_local(2) = 0            ; jx_max_local(2) = -1
        1220           0 :     jy_min_local(2) = 0            ; jy_max_local(2) = -1
        1221           0 :     jx_min_local(3) = 0            ; jx_max_local(3) = -1
        1222           0 :     jy_min_local(3) = 0            ; jy_max_local(3) = -1
        1223             :
        1224             :     call reconstruction(field_fvm,1,1,recons_fvm,irecons_tracer,&
        1225             :          llimiter,num_trac,nc,0,nhr,nhc,nhr,ns,nhr,&
        1226             :          jx_min_local,jx_max_local,jy_min_local,jy_max_local,&
        1227             :          fvm%cubeboundary,fvm%halo_interp_weight(1:ns,1-nhr:nc+nhr,1:nhr,:),fvm%ibase(1-nhr:nc+nhr,1:nhr,:),&
        1228             :          fvm%spherecentroid(:,1:nc,1:nc),&
        1229             :          fvm%recons_metrics(:,1:nc,1:nc),&
        1230             :          fvm%recons_metrics_integral(:,1:nc,1:nc)    ,&
        1231             :          fvm%rot_matrix,fvm%centroid_stretch(1:7,1:nc,1:nc),&
        1232             :          fvm%vertex_recons_weights(:,1:irecons_tracer-1,1:nc,1:nc),&
        1233           0 :          fvm%vtx_cart(:,:,1-nhc:nc+nhc,1-nhc:nc+nhc))
        1234           0 :   end subroutine get_fvm_recons
        1235             :
        1236           0 :   subroutine get_q_overlap_save(ie,k,fvm,q_fvm,num_trac,q_phys)
        1237           0 :     use dimensions_mod, only: nc,nhr,nhc,fv_nphys
        1238             :     !
        1239             :     ! weights must be initialized in fvm2phys_init before using these functions
        1240             :     !
        1241             :     use dp_mapping, only: weights_all_fvm2phys, weights_eul_index_all_fvm2phys
        1242             :     use dp_mapping, only: weights_lgr_index_all_fvm2phys, jall_fvm2phys
        1243             :     !
        1244             :     ! setting nhe=0 because we do not need reconstruction outside of element
        1245             :     !
        1246             :     integer, parameter :: nhe_local=0
        1247             :     integer, parameter :: nh = nhr!+(nhe-1) ! = 2 (nhr=2; nhe_local=1),! = 3 (nhr=2; nhe_local=2)
        1248             :
        1249             :     type(fvm_struct)     , intent(inout)        :: fvm
        1250             :     integer              , intent(in)           :: ie, k
        1251             :     integer              , intent(in)           :: num_trac
        1252             :
        1253             :     real(kind=r8), dimension(1-nhc:nc+nhc,1-nhc:nc+nhc,1:num_trac)      :: q_fvm
        1254             :     real(kind=r8), dimension(fv_nphys,fv_nphys, num_trac),   intent(out):: q_phys
        1255             :
        1256           0 :     real (kind=r8)                       :: recons_q  (irecons_tracer,1:nc,1:nc,num_trac)
        1257           0 :     logical                              :: llimiter_q(num_trac)
        1258             :     integer                              :: h,jx,jy,jdx,jdy,m_cnst,idx
        1259             :     real (kind=r8)                       :: dp_tmp, dp_fvm_tmp, tmp
        1260           0 :     real (kind=r8)                       :: dp_phys_inv(fv_nphys,fv_nphys)
        1261           0 :     integer, dimension(fv_nphys,fv_nphys):: num_overlap
        1262             :
        1263           0 :     llimiter_q=.true.
        1264           0 :     call get_fvm_recons(fvm,q_fvm,recons_q,num_trac,llimiter_q)
        1265           0 :     num_overlap(:,:) = 0
        1266           0 :     q_phys = 0.0_r8
        1267           0 :     do h=1,jall_fvm2phys(ie)
        1268           0 :        jx  = weights_lgr_index_all_fvm2phys(h,1,ie); jy  = weights_lgr_index_all_fvm2phys(h,2,ie)
        1269           0 :        jdx = weights_eul_index_all_fvm2phys(h,1,ie); jdy = weights_eul_index_all_fvm2phys(h,2,ie)
        1270             :
        1271           0 :        num_overlap(jx,jy) = num_overlap(jx,jy)+1
        1272           0 :        idx = num_overlap(jx,jy)
        1273             :
        1274           0 :        dp_fvm_tmp = fvm%dp_fvm(jdx,jdy,k)
        1275           0 :        dp_tmp = save_air_mass_overlap(idx,jx,jy,k,ie)-dp_fvm_tmp*weights_all_fvm2phys(h,1,ie)
        1276             : #ifdef PCoM
        1277             :        dp_tmp = save_air_mass_overlap(idx,jx,jy,k,ie)
        1278             : #endif
        1279           0 :        do m_cnst=1,num_trac
        1280           0 :          tmp = dp_fvm_tmp*SUM(weights_all_fvm2phys(h,:,ie)*recons_q(:,jdx,jdy,m_cnst))+q_fvm(jdx,jdy,m_cnst)*dp_tmp
        1281             : #ifdef PCoM
        1282             :          tmp = dp_fvm_tmp*weights_all_fvm2phys(h,1,ie)*q_fvm(jdx,jdy,m_cnst)
        1283             : #endif
        1284           0 :          save_q_overlap(idx,jx,jy,k,m_cnst,ie) = tmp/save_air_mass_overlap(idx,jx,jy,k,ie)
        1285           0 :          q_phys(jx,jy,m_cnst) = q_phys(jx,jy,m_cnst)+tmp
        1286             :        end do
        1287             :      end do
        1288             :      !
        1289             :      ! q_phys holds mass - convert to mixing ratio
        1290             :      !
        1291           0 :      dp_phys_inv = 1.0_r8/save_dp_phys(:,:,k,ie)!*fvm%area_sphere_physgrid)
        1292           0 :      do m_cnst=1,num_trac
        1293           0 :        q_phys(:,:,m_cnst) = q_phys(:,:,m_cnst)*dp_phys_inv
        1294           0 :        save_q_phys(:,:,k,m_cnst,ie) = q_phys(:,:,m_cnst)
        1295             :      end do
        1296           0 :    end subroutine get_q_overlap_save
        1297             :    !
        1298             :    ! Routine to overwrite thermodynamic active tracers on the GLL grid with CSLAM values
        1299             :    ! by Lagrange interpolation from 3x3 CSLAM grid to GLL grid.
        1300             :    !
        1301      738816 :    subroutine cslam2gll(elem, fvm, hybrid,nets,nete, tl_f, tl_qdp)
        1302             :      use dimensions_mod,  only: nc,nlev,np,nhc
        1303             :      use hybrid_mod,      only: hybrid_t
        1304             :      use air_composition, only: thermodynamic_active_species_num, thermodynamic_active_species_idx
        1305             :      use fvm_mod,         only: ghostBuf_cslam2gll
        1306             :      use bndry_mod,       only: ghost_exchange
        1307             :      use edge_mod,        only: ghostpack,ghostunpack
        1308             :
        1309             :      type (element_t), intent(inout):: elem(:)
        1310             :      type(fvm_struct), intent(inout):: fvm(:)
        1311             :
        1312             :      type (hybrid_t), intent(in)    :: hybrid  ! distributed parallel structure (shared)
        1313             :      integer, intent(in)            :: nets, nete, tl_f, tl_qdp
        1314             :
        1315             :      integer                                             :: ie,i,j,k,m_cnst,nq,ierr
        1316      738816 :      real (kind=r8), dimension(:,:,:,:,:)  , allocatable :: fld_fvm, fld_gll
        1317             :      !
        1318             :      ! for tensor product Lagrange interpolation
        1319             :      !
        1320             :      integer              :: nflds
        1321      738816 :      logical, allocatable :: llimiter(:)
        1322      738816 :      call t_startf('cslam2gll')
        1323      738816 :      nflds = thermodynamic_active_species_num
        1324             :
        1325             :      !Allocate variables
        1326             :      !------------------
        1327     2955264 :      allocate(fld_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev,nflds,nets:nete), stat=ierr)
        1328      738816 :      if( ierr /= 0 ) then
        1329           0 :          write(iulog,*) 'cslam2gll: fld_fvm allocation error = ', ierr
        1330           0 :          call endrun('cslam2gll: failed to allocate fld_fvm array')
        1331             :      end if
        1332             :
        1333     2955264 :      allocate(fld_gll(np,np,nlev,thermodynamic_active_species_num,nets:nete),stat=ierr)
        1334      738816 :      if( ierr /= 0 ) then
        1335           0 :          write(iulog,*) 'cslam2gll: fld_gll allocation error = ', ierr
        1336           0 :          call endrun('cslam2gll: failed to allocate fld_gll array')
        1337             :      end if
        1338             :
        1339     2216448 :      allocate(llimiter(nflds), stat=ierr)
        1340      738816 :      if( ierr /= 0 ) then
        1341           0 :          write(iulog,*) 'cslam2gll: llimiter allocation error = ', ierr
        1342           0 :          call endrun('cslam2gll: failed to allocate llimiter array')
        1343             :      end if
        1344             :      !------------------
        1345             :
        1346     5910528 :      llimiter(1:nflds) = .false.
        1347     5933616 :      do ie=nets,nete
        1348    37102416 :        do m_cnst=1,thermodynamic_active_species_num
        1349  2935062000 :          do k=1,nlev
        1350  5797396800 :            fld_fvm(1:nc,1:nc,k,m_cnst,ie) = &
        1351 43511644800 :                   fvm(ie)%c(1:nc,1:nc,k,thermodynamic_active_species_idx(m_cnst))
        1352             :          end do
        1353             :        end do
        1354             :      end do
        1355      738816 :      call t_startf('fvm:fill_halo_cslam2gll')
        1356     5933616 :      do ie=nets,nete
        1357     5933616 :        call ghostpack(ghostBuf_cslam2gll, fld_fvm(:,:,:,:,ie),nlev*nflds,0,ie)
        1358             :      end do
        1359             :
        1360      738816 :      call ghost_exchange(hybrid,ghostBuf_cslam2gll,location='cslam2gll')
        1361             :
        1362     5933616 :      do ie=nets,nete
        1363     5933616 :        call ghostunpack(ghostBuf_cslam2gll, fld_fvm(:,:,:,:,ie),nlev*nflds,0,ie)
        1364             :      end do
        1365      738816 :      call t_stopf('fvm:fill_halo_cslam2gll')
        1366             :      !
        1367             :      ! do mapping
        1368             :      !
        1369      738816 :      call fvm2dyn(fld_fvm,fld_gll,hybrid,nets,nete,nlev,nflds,fvm,llimiter,halo_filled=.true.)
        1370             :
        1371     5933616 :      do ie=nets,nete
        1372    37102416 :        do m_cnst=1,thermodynamic_active_species_num
        1373    62337600 :          elem(ie)%state%qdp(:,:,:,m_cnst,tl_qdp) = fld_gll(:,:,:,m_cnst,ie)*&
        1374 60971367600 :                elem(ie)%state%dp3d(:,:,:,tl_f)
        1375             :        end do
        1376             :      end do
        1377      738816 :      deallocate(fld_fvm, fld_gll, llimiter)
        1378      738816 :      call t_stopf('cslam2gll')
        1379      738816 :    end subroutine cslam2gll
        1380             : end module fvm_mapping

Generated by: LCOV version 1.14