LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - vertremap_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 31 | 246 | 12.6 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 1 | 8 | 12.5 %

              Line data    Source code
           1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           2             : !! Begin GPU remap module  !!
           3             : !! by Rick Archibald, 2010  !!
           4             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           5             : module vertremap_mod
           6             :
           7             :   !**************************************************************************************
           8             :   !
           9             :   !  Purpose:
          10             :   !        Construct sub-grid-scale polynomials using piecewise spline method with
          11             :   !        monotone filters.
          12             :   !
          13             :   !  References: PCM - Zerroukat et al., Q.J.R. Meteorol. Soc., 2005. (ZWS2005QJR)
          14             :   !              PSM - Zerroukat et al., Int. J. Numer. Meth. Fluids, 2005. (ZWS2005IJMF)
          15             :   !
          16             :   !**************************************************************************************
          17             :
          18             :   use shr_kind_mod,           only: r8=>shr_kind_r8
          19             :   use dimensions_mod,         only: np,nlev,qsize,nlevp,npsq,nc
          20             :   use element_mod,            only: element_t
          21             :   use fvm_control_volume_mod, only: fvm_struct
          22             :   use perf_mod,               only: t_startf, t_stopf ! _EXTERNAL
          23             :   use parallel_mod,           only: parallel_t
          24             :   use cam_abortutils,         only: endrun
          25             :
          26             :   implicit none
          27             :
          28             :   public remap1                  ! remap any field, splines, monotone
          29             :   public remap1_nofilter         ! remap any field, splines, no filter
          30             : ! todo: tweak interface to match remap1 above, rename remap1_ppm:
          31             :   public remap_q_ppm             ! remap state%Q, PPM, monotone
          32             :
          33             :   contains
          34             :
          35             : !=======================================================================================================!
          36             :
          37    25974000 :     subroutine remap1(Qdp,nx,qstart,qstop,qsize,dp1,dp2,ptop,identifier,Qdp_mass,kord)
          38             :       use fv_mapz,      only: map1_ppm
          39             :       ! remap 1 field
          40             :       ! input:  Qdp   field to be remapped (NOTE: MASS, not MIXING RATIO)
          41             :       !         dp1   layer thickness (source)
          42             :       !         dp2   layer thickness (target)
          43             :       !
          44             :       ! output: remaped Qdp, conserving mass, monotone on Q=Qdp/dp
          45             :       !
          46             :       integer, intent(in) :: nx,qstart,qstop,qsize
          47             :       real (kind=r8), intent(inout) :: Qdp(nx,nx,nlev,qsize)
          48             :       real (kind=r8), intent(in) :: dp1(nx,nx,nlev),dp2(nx,nx,nlev)
          49             :       integer, intent(in) :: identifier  !0: tracers, 1: T, -1: u,v
          50             :       real (kind=r8), intent(in) :: ptop
          51             :       logical, intent(in) :: Qdp_mass
          52             :       integer, intent(in) :: kord(qsize)
          53             :       ! ========================
          54             :       ! Local Variables
          55             :       ! ========================
          56    51948000 :       real (kind=r8) :: pe1(nx,nlev+1),pe2(nx,nlev+1),inv_dp(nx,nx,nlev),dp2_local(nx,nlev)
          57    51948000 :       real (kind=r8) :: tmp(nx,nlev), gz(nx)
          58             :       integer        :: i,j,k,itrac
          59             :       logical        :: logp
          60    25974000 :       integer        :: kord_local(qsize)
          61             :
          62    72727200 :       kord_local = kord
          63             :
          64    72727200 :       if (any(kord(:) >= 0)) then
          65           0 :         if (.not.qdp_mass) then
          66           0 :           do itrac=1,qsize
          67           0 :             if (kord(itrac) >= 0) then
          68           0 :               Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)*dp1(:,:,:)
          69             :             end if
          70             :           end do
          71             :         end if
          72           0 :         call remap_Q_ppm(qdp,nx,qstart,qstop,qsize,dp1,dp2,kord)
          73           0 :         if (.not.qdp_mass) then
          74           0 :           do itrac=1,qsize
          75           0 :             if (kord(itrac) >= 0) then
          76           0 :               Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)/dp2(:,:,:)
          77             :             end if
          78             :           end do
          79             :         end if
          80             :       endif
          81    25974000 :       if (any(kord(:)<0)) then
          82             :         !
          83             :         ! check if remapping over p or log(p)
          84             :         !
          85             :         ! can not mix and match here (all kord's must >-20 or <=-20)
          86             :         !
          87    25974000 :         if (any(kord(:)>-20)) then
          88    72727200 :           kord_local = abs(kord)
          89             :           logp    = .false.
          90             :         else
          91           0 :           kord_local = abs(kord/10)
          92           0 :           if (identifier==1) then
          93             :             logp    = .true.
          94             :           else
          95           0 :             logp    = .false.
          96             :           end if
          97             :         end if
          98             :         !
          99             :         ! modified FV3 vertical remapping
         100             :         !
         101    25974000 :         if (qdp_mass) then
         102  5683111200 :           inv_dp = 1.0_r8/dp1
         103    31168800 :           do itrac=1,qsize
         104    31168800 :             if (kord(itrac)<0) then
         105 11366222400 :               Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)*inv_dp(:,:,:)
         106             :             end if
         107             :           end do
         108             :         end if
         109    25974000 :         if (logp) then
         110           0 :           do j=1,nx
         111           0 :             pe1(:,1) = ptop
         112           0 :             pe2(:,1) = ptop
         113           0 :             do k=1,nlev
         114           0 :               do i=1,nx
         115           0 :                 pe1(i,k+1) = pe1(i,k)+dp1(i,j,k)
         116           0 :                 pe2(i,k+1) = pe2(i,k)+dp2(i,j,k)
         117             :               end do
         118             :             end do
         119           0 :             pe1(:,nlev+1) = pe2(:,nlev+1)
         120           0 :             do k=1,nlev+1
         121           0 :               do i=1,nx
         122           0 :                 pe1(i,k) = log(pe1(i,k))
         123           0 :                 pe2(i,k) = log(pe2(i,k))
         124             :               end do
         125             :             end do
         126             :
         127           0 :             do itrac=1,qsize
         128           0 :               if (kord(itrac)<0) then
         129             :                 call map1_ppm( nlev, pe1(:,:),   Qdp(:,:,:,itrac),   gz,   &
         130             :                      nlev, pe2(:,:),    Qdp(:,:,:,itrac),               &
         131           0 :                      1, nx, j, 1, nx, 1, nx, identifier, kord_local(itrac))
         132             :               end if
         133             :             end do
         134             :             !      call mapn_tracer(qsize, nlev, pe1, pe2, Qdp, dp2_local, kord, j,     &
         135             :             !           1, nx, 1, nx, 1, nx, 0.0_r8, fill)
         136             :           end do
         137             :         else
         138   124675200 :           do j=1,nx
         139   477921600 :             pe1(:,1) = ptop
         140   477921600 :             pe2(:,1) = ptop
         141  2664932400 :             do k=1,nlev
         142 12524662800 :               do i=1,nx
         143  9859730400 :                 pe1(i,k+1) = pe1(i,k)+dp1(i,j,k)
         144 12425961600 :                 pe2(i,k+1) = pe2(i,k)+dp2(i,j,k)
         145             :               end do
         146             :             end do
         147   477921600 :             pe1(:,nlev+1) = pe2(:,nlev+1)
         148   296103600 :             do itrac=1,qsize
         149   270129600 :               if (kord(itrac)<0) then
         150             :                 call map1_ppm( nlev, pe1(:,:),   Qdp(:,:,:,itrac),   gz,   &!phl
         151             :                      nlev, pe2(:,:),    Qdp(:,:,:,itrac),               &
         152   171428400 :                      1, nx, j, 1, nx, 1, nx, identifier, kord_local(itrac))
         153             :               end if
         154             :             end do
         155             :             !      call mapn_tracer(qsize, nlev, pe1, pe2, Qdp, dp2_local, kord, j,     &
         156             :             !           1, nx, 1, nx, 1, nx, 0.0_r8, fill)
         157             :           end do
         158             :         end if
         159    25974000 :         if (qdp_mass) then
         160    31168800 :           do itrac=1,qsize
         161    31168800 :             if (kord(itrac)<0) then
         162 11366222400 :               Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)*dp2(:,:,:)
         163             :             end if
         164             :           end do
         165             :         end if
         166             :       end if
         167    25974000 :     end subroutine remap1
         168             :
         169           0 : subroutine remap1_nofilter(Qdp,nx,qsize,dp1,dp2)
         170             :   ! remap 1 field
         171             :   ! input:  Qdp   field to be remapped (NOTE: MASS, not MIXING RATIO)
         172             :   !         dp1   layer thickness (source)
         173             :   !         dp2   layer thickness (target)
         174             :   !
         175             :   ! output: remaped Qdp, conserving mass
         176             :   !
         177             :   implicit none
         178             :   integer, intent(in) :: nx,qsize
         179             :   real (kind=r8), intent(inout) :: Qdp(nx,nx,nlev,qsize)
         180             :   real (kind=r8), intent(in) :: dp1(nx,nx,nlev),dp2(nx,nx,nlev)
         181             :   ! ========================
         182             :   ! Local Variables
         183             :   ! ========================
         184             :
         185             :   real (kind=r8), dimension(nlev+1)    :: rhs,lower_diag,diag,upper_diag,q_diag,zgam,z1c,z2c,zv
         186             :   real (kind=r8), dimension(nlev)      :: h,Qcol,za0,za1,za2,zarg,zhdp
         187             :   real (kind=r8)  :: tmp_cal,zv1,zv2
         188             :   integer :: zkr(nlev+1),i,ilev,j,jk,k,q
         189             :   logical :: abort=.false.
         190             :   !   call t_startf('remap1_nofilter')
         191             :
         192             : #if (defined COLUMN_OPENMP)
         193             :   !$omp parallel do num_threads(tracer_num_threads) &
         194             :   !$omp    private(q,i,j,z1c,z2c,zv,k,Qcol,zkr,ilev) &
         195             :   !$omp    private(jk,zgam,zhdp,h,zarg,rhs,lower_diag,diag,upper_diag,q_diag,tmp_cal) &
         196             :   !$omp    private(za0,za1,za2) &
         197             :   !$omp    private(ip2,zv1,zv2)
         198             : #endif
         199           0 :   do q=1,qsize
         200           0 :     do i=1,nx
         201           0 :       do j=1,nx
         202             :
         203           0 :         z1c(1)=0 ! source grid
         204           0 :         z2c(1)=0 ! target grid
         205           0 :         do k=1,nlev
         206           0 :           z1c(k+1)=z1c(k)+dp1(i,j,k)
         207           0 :           z2c(k+1)=z2c(k)+dp2(i,j,k)
         208             :         enddo
         209             :
         210           0 :         zv(1)=0
         211           0 :         do k=1,nlev
         212           0 :           Qcol(k)=Qdp(i,j,k,q)!  *(z1c(k+1)-z1c(k)) input is mass
         213           0 :           zv(k+1) = zv(k)+Qcol(k)
         214             :         enddo
         215             :
         216           0 :         if (ABS(z2c(nlev+1)-z1c(nlev+1)) >= 0.000001_r8) then
         217           0 :           write(6,*) 'SURFACE PRESSURE IMPLIED BY ADVECTION SCHEME'
         218           0 :           write(6,*) 'NOT CORRESPONDING TO SURFACE PRESSURE IN    '
         219           0 :           write(6,*) 'DATA FOR MODEL LEVELS'
         220           0 :           write(6,*) 'PLEVMODEL=',z2c(nlev+1)
         221           0 :           write(6,*) 'PLEV     =',z1c(nlev+1)
         222           0 :           write(6,*) 'DIFF     =',z2c(nlev+1)-z1c(nlev+1)
         223           0 :           abort=.true.
         224             :         endif
         225             :
         226             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         227             :         !! quadratic splies with UK met office monotonicity constraints  !!
         228             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         229             :
         230           0 :         zkr  = 99
         231           0 :         ilev = 2
         232           0 :         zkr(1) = 1
         233           0 :         zkr(nlev+1) = nlev
         234           0 :         kloop: do k = 2,nlev
         235           0 :           do jk = ilev,nlev+1
         236           0 :             if (z1c(jk) >= z2c(k)) then
         237           0 :               ilev      = jk
         238           0 :               zkr(k)   = jk-1
         239           0 :               cycle kloop
         240             :             endif
         241             :           enddo
         242             :         enddo kloop
         243             :
         244           0 :         zgam  = (z2c(1:nlev+1)-z1c(zkr)) / (z1c(zkr+1)-z1c(zkr))
         245           0 :         zgam(1)      = 0.0_r8
         246           0 :         zgam(nlev+1) = 1.0_r8
         247           0 :         zhdp = z1c(2:nlev+1)-z1c(1:nlev)
         248             :
         249             :
         250           0 :         h = 1/zhdp
         251           0 :         zarg = Qcol * h
         252           0 :         rhs = 0
         253           0 :         lower_diag = 0
         254           0 :         diag = 0
         255           0 :         upper_diag = 0
         256             :
         257           0 :         rhs(1)=3*zarg(1)
         258           0 :         rhs(2:nlev) = 3*(zarg(2:nlev)*h(2:nlev) + zarg(1:nlev-1)*h(1:nlev-1))
         259           0 :         rhs(nlev+1)=3*zarg(nlev)
         260             :
         261           0 :         lower_diag(1)=1
         262           0 :         lower_diag(2:nlev) = h(1:nlev-1)
         263           0 :         lower_diag(nlev+1)=1
         264             :
         265           0 :         diag(1)=2
         266           0 :         diag(2:nlev) = 2*(h(2:nlev) + h(1:nlev-1))
         267           0 :         diag(nlev+1)=2
         268             :
         269           0 :         upper_diag(1)=1
         270           0 :         upper_diag(2:nlev) = h(2:nlev)
         271           0 :         upper_diag(nlev+1)=0
         272             :
         273           0 :         q_diag(1)=-upper_diag(1)/diag(1)
         274           0 :         rhs(1)= rhs(1)/diag(1)
         275             :
         276           0 :         do k=2,nlev+1
         277           0 :           tmp_cal    =  1/(diag(k)+lower_diag(k)*q_diag(k-1))
         278           0 :           q_diag(k) = -upper_diag(k)*tmp_cal
         279           0 :           rhs(k) =  (rhs(k)-lower_diag(k)*rhs(k-1))*tmp_cal
         280             :         enddo
         281           0 :         do k=nlev,1,-1
         282           0 :           rhs(k)=rhs(k)+q_diag(k)*rhs(k+1)
         283             :         enddo
         284             :
         285           0 :         za0 = rhs(1:nlev)
         286           0 :         za1 = -4*rhs(1:nlev) - 2*rhs(2:nlev+1) + 6*zarg
         287           0 :         za2 =  3*rhs(1:nlev) + 3*rhs(2:nlev+1) - 6*zarg
         288             :
         289             :
         290             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         291             :         !! start iteration from top to bottom of atmosphere !!
         292             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         293             :
         294             :         zv1 = 0
         295           0 :         do k=1,nlev
         296           0 :           if (zgam(k+1)>1_r8) then
         297           0 :             WRITE(*,*) 'r not in [0:1]', zgam(k+1)
         298           0 :             abort=.true.
         299             :           endif
         300           0 :           zv2 = zv(zkr(k+1))+(za0(zkr(k+1))*zgam(k+1)+(za1(zkr(k+1))/2)*(zgam(k+1)**2)+ &
         301           0 :                (za2(zkr(k+1))/3)*(zgam(k+1)**3))*zhdp(zkr(k+1))
         302           0 :           Qdp(i,j,k,q) = (zv2 - zv1) ! / (z2c(k+1)-z2c(k) ) dont convert back to mixing ratio
         303           0 :           zv1 = zv2
         304             :         enddo
         305             :       enddo
         306             :     enddo
         307             :   enddo ! q loop
         308           0 :   if (abort) then
         309           0 :     call endrun('Bad levels in remap1_nofilter.  usually CFL violatioin')
         310             :   end if
         311           0 : end subroutine remap1_nofilter
         312             :
         313             : !=============================================================================!
         314             :
         315             : !This uses the exact same model and reference grids and data as remap_Q, but it interpolates
         316             : !using PPM instead of splines.
         317           0 : subroutine remap_Q_ppm(Qdp,nx,qstart,qstop,qsize,dp1,dp2,kord)
         318             :   ! remap 1 field
         319             :   ! input:  Qdp   field to be remapped (NOTE: MASS, not MIXING RATIO)
         320             :   !         dp1   layer thickness (source)
         321             :   !         dp2   layer thickness (target)
         322             :   !
         323             :   ! output: remaped Qdp, conserving mass
         324             :   !
         325             :   implicit none
         326             :   integer,intent(in) :: nx,qstart,qstop,qsize
         327             :   real (kind=r8), intent(inout) :: Qdp(nx,nx,nlev,qsize)
         328             :   real (kind=r8), intent(in) :: dp1(nx,nx,nlev),dp2(nx,nx,nlev)
         329             :   integer       , intent(in) :: kord(qsize)
         330             :   ! Local Variables
         331             :   integer, parameter :: gs = 2                              !Number of cells to place in the ghost region
         332             :   real(kind=r8), dimension(       nlev+2 ) :: pio    !Pressure at interfaces for old grid
         333             :   real(kind=r8), dimension(       nlev+1 ) :: pin    !Pressure at interfaces for new grid
         334             :   real(kind=r8), dimension(       nlev+1 ) :: masso  !Accumulate mass up to each interface
         335             :   real(kind=r8), dimension(  1-gs:nlev+gs) :: ao     !Tracer value on old grid
         336             :   real(kind=r8), dimension(  1-gs:nlev+gs) :: dpo    !change in pressure over a cell for old grid
         337             :   real(kind=r8), dimension(  1-gs:nlev+gs) :: dpn    !change in pressure over a cell for old grid
         338             :   real(kind=r8), dimension(3,     nlev   ) :: coefs  !PPM coefficients within each cell
         339             :   real(kind=r8), dimension(       nlev   ) :: z1, z2
         340             :   real(kind=r8) :: ppmdx(10,0:nlev+1)  !grid spacings
         341             :   real(kind=r8) :: massn1, massn2, ext(2)
         342             :   integer :: i, j, k, q, kk, kid(nlev)
         343             :
         344           0 :   do j = 1 , nx
         345           0 :     do i = 1 , nx
         346             :
         347           0 :       pin(1)=0
         348           0 :       pio(1)=0
         349           0 :       do k=1,nlev
         350           0 :          dpn(k)=dp2(i,j,k)
         351           0 :          dpo(k)=dp1(i,j,k)
         352           0 :          pin(k+1)=pin(k)+dpn(k)
         353           0 :          pio(k+1)=pio(k)+dpo(k)
         354             :       enddo
         355             :
         356             :
         357             :
         358           0 :       pio(nlev+2) = pio(nlev+1) + 1._r8  !This is here to allow an entire block of k threads to run in the remapping phase.
         359             :                                       !It makes sure there's an old interface value below the domain that is larger.
         360           0 :       pin(nlev+1) = pio(nlev+1)       !The total mass in a column does not change.
         361             :                                       !Therefore, the pressure of that mass cannot either.
         362             :       !Fill in the ghost regions with mirrored values. if vert_remap_q_alg is defined, this is of no consequence.
         363           0 :       do k = 1 , gs
         364           0 :         dpo(1   -k) = dpo(       k)
         365           0 :         dpo(nlev+k) = dpo(nlev+1-k)
         366             :       enddo
         367             :
         368             :       !Compute remapping intervals once for all tracers. Find the old grid cell index in which the
         369             :       !k-th new cell interface resides. Then integrate from the bottom of that old cell to the new
         370             :       !interface location. In practice, the grid never deforms past one cell, so the search can be
         371             :       !simplified by this. Also, the interval of integration is usually of magnitude close to zero
         372             :       !or close to dpo because of minimial deformation.
         373             :       !Numerous tests confirmed that the bottom and top of the grids match to machine precision, so
         374             :       !I set them equal to each other.
         375           0 :       do k = 1 , nlev
         376           0 :         kk = k  !Keep from an order n^2 search operation by assuming the old cell index is close.
         377             :         !Find the index of the old grid cell in which this new cell's bottom interface resides.
         378             : !        do while ( pio(kk) <= pin(k+1) )
         379             : !          kk = kk + 1
         380             : !          if(kk==nlev+2) exit
         381             : !        enddo
         382             :         !        kk = kk - 1                   !kk is now the cell index we're integrating over.
         383             :
         384           0 :         if (pio(kk) <= pin(k+1)) then
         385           0 :           do while ( pio(kk) <= pin(k+1) )
         386           0 :             kk = kk + 1
         387             :           enddo
         388           0 :           kk = kk - 1                   !kk is now the cell index we're integrating over.
         389             :         else
         390             :           call binary_search(pio, pin(k+1), kk)
         391             :         end if
         392           0 :         if (kk == nlev+1) kk = nlev   !This is to keep the indices in bounds.
         393             :                                       !Top bounds match anyway, so doesn't matter what coefficients are used
         394           0 :         kid(k) = kk                   !Save for reuse
         395           0 :         z1(k) = -0.5_R8                !This remapping assumes we're starting from the left interface of an old grid cell
         396             :                                       !In fact, we're usually integrating very little or almost all of the cell in question
         397           0 :         z2(k) = ( pin(k+1) - ( pio(kk) + pio(kk+1) ) * 0.5_r8 ) / dpo(kk)  !PPM interpolants are normalized to an independent
         398             :                                                                         !coordinate domain [-0.5,0.5].
         399             :       enddo
         400             :
         401             :       !This turned out a big optimization, remembering that only parts of the PPM algorithm depends on the data, namely the
         402             :       !limiting. So anything that depends only on the grid is pre-computed outside the tracer loop.
         403           0 :       ppmdx(:,:) = compute_ppm_grids( dpo)
         404             :
         405             :       !From here, we loop over tracers for only those portions which depend on tracer data, which includes PPM limiting and
         406             :       !mass accumulation
         407           0 :       do q = qstart, qstop
         408           0 :         if (kord(q) >= 0) then
         409             :         !Accumulate the old mass up to old grid cell interface locations to simplify integration
         410             :         !during remapping. Also, divide out the grid spacing so we're working with actual tracer
         411             :         !values and can conserve mass. The option for ifndef ZEROHORZ I believe is there to ensure
         412             :         !tracer consistency for an initially uniform field. I copied it from the old remap routine.
         413           0 :         masso(1) = 0._r8
         414             :
         415           0 :         do k = 1 , nlev
         416           0 :           ao(k) = Qdp(i,j,k,q)
         417           0 :           masso(k+1) = masso(k) + ao(k) !Accumulate the old mass. This will simplify the remapping
         418           0 :           ao(k) = ao(k) / dpo(k)        !Divide out the old grid spacing because we want the tracer mixing ratio, not mass.
         419             :         enddo
         420             :         !Fill in ghost values. Ignored if kord == 2
         421           0 :         if (kord(q) == 10) then
         422           0 :           ext(1) = minval(ao(1:nlev))
         423           0 :           ext(2) = maxval(ao(1:nlev))
         424           0 :           call linextrap(dpo(2), dpo(1), dpo(0), dpo(-1), ao(2), ao(1), ao(0), ao(-1), ext(1), ext(2))
         425             :           call linextrap(dpo(nlev-1), dpo(nlev), dpo(nlev+1), dpo(nlev+2), &
         426           0 :                ao(nlev-1), ao(nlev), ao(nlev+1), ao(nlev+2), ext(1), ext(2))
         427             :         else
         428           0 :           do k = 1 , gs
         429           0 :             ao(1   -k) = ao(       k)
         430           0 :             ao(nlev+k) = ao(nlev+1-k)
         431             :           enddo
         432             :         end if
         433             :         !Compute monotonic and conservative PPM reconstruction over every cell
         434           0 :         coefs(:,:) = compute_ppm( ao , ppmdx, kord(q) )
         435             :         !Compute tracer values on the new grid by integrating from the old cell bottom to the new
         436             :         !cell interface to form a new grid mass accumulation. Taking the difference between
         437             :         !accumulation at successive interfaces gives the mass inside each cell. Since Qdp is
         438             :         !supposed to hold the full mass this needs no normalization.
         439           0 :         massn1 = 0._r8
         440           0 :         do k = 1 , nlev
         441           0 :           kk = kid(k)
         442           0 :           massn2 = masso(kk) + integrate_parabola( coefs(:,kk) , z1(k) , z2(k) ) * dpo(kk)
         443           0 :           Qdp(i,j,k,q) = massn2 - massn1
         444           0 :           massn1 = massn2
         445             :         enddo
         446             :         end if
         447             :       enddo
         448             :     enddo
         449             :   enddo
         450             : ! call t_stopf('remap_Q_ppm')
         451           0 : end subroutine remap_Q_ppm
         452             :
         453             : ! Find k such that pio(k) <= pivot < pio(k+1). Provide a reasonable input
         454             : ! value for k.
         455           0 : subroutine binary_search(pio, pivot, k)
         456             :   real(kind=r8), intent(in)    :: pio(nlev+2), pivot
         457             :   integer,       intent(inout) :: k
         458             :   integer :: lo, hi, mid
         459             :
         460           0 :   if (pio(k) > pivot) then
         461             :     lo = 1
         462             :     hi = k
         463             :   else
         464           0 :     lo = k
         465           0 :     hi = nlev+2
         466             :   end if
         467           0 :   do while (hi > lo + 1)
         468           0 :     k = (lo + hi)/2
         469           0 :     if (pio(k) > pivot) then
         470             :       hi = k
         471             :     else
         472           0 :       lo = k
         473             :     end if
         474             :   end do
         475           0 :   k = lo
         476           0 : end subroutine binary_search
         477             : !=======================================================================================================!
         478             :
         479             : !This compute grid-based coefficients from Collela & Woodward 1984.
         480           0 : function compute_ppm_grids( dx )   result(rslt)
         481             :   implicit none
         482             :   real(kind=r8), intent(in) :: dx(-1:nlev+2)  !grid spacings
         483             :   real(kind=r8)             :: rslt(10,0:nlev+1)  !grid spacings
         484             :   integer :: j
         485             :
         486             :   !Calculate grid-based coefficients for stage 1 of compute_ppm
         487           0 :   do j = 0 , nlev+1
         488           0 :     rslt( 1,j) = dx(j) / ( dx(j-1) + dx(j) + dx(j+1) )
         489           0 :     rslt( 2,j) = ( 2._r8*dx(j-1) + dx(j) ) / ( dx(j+1) + dx(j) )
         490           0 :     rslt( 3,j) = ( dx(j) + 2._r8*dx(j+1) ) / ( dx(j-1) + dx(j) )
         491             :   enddo
         492             :
         493             :   !Caculate grid-based coefficients for stage 2 of compute_ppm
         494           0 :   do j = 0 , nlev
         495           0 :     rslt( 4,j) = dx(j) / ( dx(j) + dx(j+1) )
         496           0 :     rslt( 5,j) = 1._r8 / sum( dx(j-1:j+2) )
         497           0 :     rslt( 6,j) = ( 2._r8 * dx(j+1) * dx(j) ) / ( dx(j) + dx(j+1 ) )
         498           0 :     rslt( 7,j) = ( dx(j-1) + dx(j  ) ) / ( 2._r8 * dx(j  ) + dx(j+1) )
         499           0 :     rslt( 8,j) = ( dx(j+2) + dx(j+1) ) / ( 2._r8 * dx(j+1) + dx(j  ) )
         500           0 :     rslt( 9,j) = dx(j  ) * ( dx(j-1) + dx(j  ) ) / ( 2._r8*dx(j  ) +    dx(j+1) )
         501           0 :     rslt(10,j) = dx(j+1) * ( dx(j+1) + dx(j+2) ) / (    dx(j  ) + 2._r8*dx(j+1) )
         502             :   enddo
         503           0 : end function compute_ppm_grids
         504             :
         505             :
         506             : !=======================================================================================================!
         507             :
         508             :
         509             :
         510             : !This computes a limited parabolic interpolant using a net 5-cell stencil, but the stages of computation are broken up into 3 stages
         511           0 : function compute_ppm( a , dx , kord)    result(coefs)
         512             :   implicit none
         513             :   real(kind=r8), intent(in) :: a    (    -1:nlev+2)  !Cell-mean values
         514             :   real(kind=r8), intent(in) :: dx   (10,  0:nlev+1)  !grid spacings
         515             :   integer,       intent(in) :: kord
         516             :   real(kind=r8) ::             coefs(0:2,   nlev  )  !PPM coefficients (for parabola)
         517             :   real(kind=r8) :: ai (0:nlev  )                     !fourth-order accurate, then limited interface values
         518             :   real(kind=r8) :: dma(0:nlev+1)                     !An expression from Collela's '84 publication
         519             :   real(kind=r8) :: da                                !Ditto
         520             :   ! Hold expressions based on the grid (which are cumbersome).
         521             :   real(kind=r8) :: al, ar                            !Left and right interface values for cell-local limiting
         522             :   integer :: j
         523             :   integer :: indB, indE
         524             :
         525             :   ! Stage 1: Compute dma for each cell, allowing a 1-cell ghost stencil below and above the domain
         526           0 :   do j = 0 , nlev+1
         527           0 :     da = dx(1,j) * ( dx(2,j) * ( a(j+1) - a(j) ) + dx(3,j) * ( a(j) - a(j-1) ) )
         528           0 :     dma(j) = minval( (/ abs(da) , 2._r8 * abs( a(j) - a(j-1) ) , 2._r8 * abs( a(j+1) - a(j) ) /) ) * sign(1._r8,da)
         529           0 :     if ( ( a(j+1) - a(j) ) * ( a(j) - a(j-1) ) <= 0._r8 ) dma(j) = 0._r8
         530             :   enddo
         531             :
         532             :   ! Stage 2: Compute ai for each cell interface in the physical domain (dimension nlev+1)
         533           0 :   do j = 0 , nlev
         534           0 :     ai(j) = a(j) + dx(4,j) * ( a(j+1) - a(j) ) + dx(5,j) * ( dx(6,j) * ( dx(7,j) - dx(8,j) ) &
         535           0 :          * ( a(j+1) - a(j) ) - dx(9,j) * dma(j+1) + dx(10,j) * dma(j) )
         536             :   enddo
         537             :
         538             :   ! Stage 3: Compute limited PPM interpolant over each cell in the physical domain
         539             :   ! (dimension nlev) using ai on either side and ao within the cell.
         540           0 :   do j = 1 , nlev
         541           0 :     al = ai(j-1)
         542           0 :     ar = ai(j  )
         543           0 :     if ( (ar - a(j)) * (a(j) - al) <= 0._r8 ) then
         544           0 :       al = a(j)
         545           0 :       ar = a(j)
         546             :     endif
         547           0 :     if ( (ar - al) * (a(j) - (al + ar)/2._r8) >  (ar - al)**2/6._r8 ) al = 3._r8*a(j) - 2._r8 * ar
         548           0 :     if ( (ar - al) * (a(j) - (al + ar)/2._r8) < -(ar - al)**2/6._r8 ) ar = 3._r8*a(j) - 2._r8 * al
         549             :     !Computed these coefficients from the edge values and cell mean in Maple. Assumes normalized coordinates: xi=(x-x0)/dx
         550           0 :     coefs(0,j) = 1.5_r8 * a(j) - ( al + ar ) / 4._r8
         551           0 :     coefs(1,j) = ar - al
         552           0 :     coefs(2,j) = 3._r8 * (-2._r8 * a(j) + ( al + ar ))
         553             :   enddo
         554             :
         555             :   !If kord == 2, use piecewise constant in the boundaries, and don't use ghost cells.
         556           0 :   if (kord == 2) then
         557           0 :     coefs(0,1:2) = a(1:2)
         558           0 :     coefs(1:2,1:2) = 0._r8
         559           0 :     coefs(0,nlev-1:nlev) = a(nlev-1:nlev)
         560           0 :     coefs(1:2,nlev-1:nlev) = 0._r8
         561             :   endif
         562           0 : end function compute_ppm
         563             :
         564             :
         565             : !=======================================================================================================!
         566             :
         567             :
         568             : !Simple function computes the definite integral of a parabola in normalized coordinates, xi=(x-x0)/dx,
         569             : !given two bounds. Make sure this gets inlined during compilation.
         570           0 : function integrate_parabola( a , x1 , x2 )    result(mass)
         571             :   implicit none
         572             :   real(kind=r8), intent(in) :: a(0:2)  !Coefficients of the parabola
         573             :   real(kind=r8), intent(in) :: x1      !lower domain bound for integration
         574             :   real(kind=r8), intent(in) :: x2      !upper domain bound for integration
         575             :   real(kind=r8)             :: mass
         576           0 :   mass = a(0) * (x2 - x1) + a(1) * (x2 ** 2 - x1 ** 2) / 0.2D1 + a(2) * (x2 ** 3 - x1 ** 3) / 0.3D1
         577           0 : end function integrate_parabola
         578             :
         579             :
         580             : !=============================================================================================!
         581           0 :   subroutine linextrap(dx1,dx2,dx3,dx4,y1,y2,y3,y4,lo,hi)
         582             :     real(kind=r8), intent(in) :: dx1,dx2,dx3,dx4,y1,y2,lo,hi
         583             :     real(kind=r8), intent(out) :: y3,y4
         584             :
         585             :     real(kind=r8), parameter :: half = 0.5_r8
         586             :
         587             :     real(kind=r8) :: x1,x2,x3,x4,a
         588             :
         589           0 :     x1 = half*dx1
         590           0 :     x2 = x1 + half*(dx1 + dx2)
         591           0 :     x3 = x2 + half*(dx2 + dx3)
         592           0 :     x4 = x3 + half*(dx3 + dx4)
         593           0 :     a  = (x3-x1)/(x2-x1)
         594           0 :     y3 = (1.0_r8-a)*y1 + a*y2
         595           0 :     a  = (x4-x1)/(x2-x1)
         596           0 :     y4 = (1.0_r8-a)*y1 + a*y2
         597           0 :     y3 = max(lo, min(hi, y3))
         598           0 :     y4 = max(lo, min(hi, y4))
         599           0 :   end subroutine linextrap
         600             : end module vertremap_mod
         601             :
         602             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         603             : !! End GPU remap module    !!
         604             : !! by Rick Archibald, 2010  !!
         605             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Generated by: LCOV version 1.14