LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - fvm_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 289 | 391 | 73.9 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 5 | 7 | 71.4 %

              Line data    Source code
           1             : #define FVM_TIMERS .FALSE.
           2             : !-----------------------------------------------------------------------------!
           3             : !MODULE FVM_MOD-----------------------------------------------------CE-for FVM!
           4             : ! FVM_MOD File for the fvm project in HOMME                                   !
           5             : ! Author: Christoph Erath                                                     !
           6             : ! Date: 25.January 2011                                                       !
           7             : ! MAIN module to run fvm on HOMME                                             !
           8             : ! 14.November 2011: reorganisation done                                       !
           9             : ! 7.Februar 2012: cslam_run and cslam_runair                                  !
          10             : !-----------------------------------------------------------------------------!
          11             :
          12             : module fvm_mod
          13             :   use shr_kind_mod,           only: r8=>shr_kind_r8
          14             :   use edge_mod,               only: initghostbuffer, freeghostbuffer, ghostpack, ghostunpack
          15             :   use edgetype_mod,           only: edgebuffer_t
          16             :   use bndry_mod,              only: ghost_exchange
          17             :   use thread_mod,             only: horz_num_threads, vert_num_threads
          18             :
          19             :   use element_mod,            only: element_t
          20             :   use fvm_control_volume_mod, only: fvm_struct
          21             :   use hybrid_mod,             only: hybrid_t
          22             :
          23             :   implicit none
          24             :   private
          25             :   save
          26             :
          27             :   type (EdgeBuffer_t) :: edgeveloc
          28             :   type (EdgeBuffer_t), public  :: ghostBufQnhc_s
          29             :   type (EdgeBuffer_t), public  :: ghostBufQnhc_t1
          30             :   type (EdgeBuffer_t), public  :: ghostBufQnhc_vh
          31             :   type (EdgeBuffer_t), public  :: ghostBufQnhc_h
          32             :   type (EdgeBuffer_t), public  :: ghostBufQ1_h
          33             :   type (EdgeBuffer_t), public  :: ghostBufQ1_vh
          34             : !  type (EdgeBuffer_t), private  :: ghostBufFlux_h
          35             :   type (EdgeBuffer_t), public  :: ghostBufFlux_vh
          36             :   type (EdgeBuffer_t), public  :: ghostBufQnhcJet_h
          37             :   type (EdgeBuffer_t), public  :: ghostBufFluxJet_h
          38             :   type (EdgeBuffer_t), public  :: ghostBufPG_s
          39             :   type (EdgeBuffer_t), public  :: ghostBuf_cslam2gll
          40             :
          41             :   interface fill_halo_fvm
          42             :      module procedure fill_halo_fvm_noprealloc
          43             :      module procedure fill_halo_fvm_prealloc
          44             :   end interface
          45             :
          46             :
          47             :   public :: edgeveloc, fvm_init1,fvm_init2, fill_halo_fvm, fvm_pg_init,fvm_init3,fill_halo_and_extend_panel
          48             :
          49             : contains
          50             :
          51           0 :   subroutine fill_halo_fvm_noprealloc(elem,fvm,hybrid,nets,nete,ndepth,kmin,kmax,ksize)
          52             :     use perf_mod, only : t_startf, t_stopf ! _EXTERNAL
          53             :     use dimensions_mod, only: nc, ntrac, nlev
          54             :     implicit none
          55             :     type (element_t),intent(inout)            :: elem(:)
          56             :     type (fvm_struct),intent(inout)           :: fvm(:)
          57             :     type (hybrid_t),intent(in)                :: hybrid
          58             :
          59           0 :     type (edgeBuffer_t)                      :: cellghostbuf
          60             :
          61             :     integer,intent(in)                        :: nets,nete
          62             :     integer,intent(in)                        :: ndepth     ! depth of halo
          63             :     integer,intent(in)                        :: kmin,kmax  ! min and max vertical level
          64             :     integer,intent(in)                        :: ksize      ! the total number of vertical
          65             :
          66             :     integer                                   :: ie,i1,i2,kblk,kptr,q
          67             :     !
          68             :     !
          69             :
          70           0 :     if(kmin .ne. 1 .or. kmax .ne. nlev) then
          71           0 :        print *,'WARNING: fill_halo_fvm_noprealloc does not support the passing of non-contigous arrays'
          72           0 :        print *,'WARNING:   incorrect answers are likely'
          73             :     endif
          74             :     if(FVM_TIMERS) call t_startf('FVM:initbuf')
          75           0 :     i1=1-ndepth
          76           0 :     i2=nc+ndepth
          77           0 :     kblk = kmax-kmin+1
          78           0 :     call initghostbuffer(hybrid%par,cellghostbuf,elem,kblk*(ntrac+1),ndepth,nc)
          79             :     if(FVM_TIMERS) call t_stopf('FVM:initbuf')
          80             :     if(FVM_TIMERS) call t_startf('FVM:pack')
          81           0 :     do ie=nets,nete
          82           0 :        kptr = kmin-1
          83           0 :        call ghostpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk,   kptr,ie)
          84           0 :        do q=1,ntrac
          85           0 :           kptr = kptr + ksize
          86           0 :           call ghostpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,q)   ,kblk,kptr,ie)
          87             :        enddo
          88             :     end do
          89             :     if(FVM_TIMERS) call t_stopf('FVM:pack')
          90             :     if(FVM_TIMERS) call t_startf('FVM:Communication')
          91           0 :     call ghost_exchange(hybrid,cellghostbuf,location='fill_halo_fvm_noprealloc')
          92             :     if(FVM_TIMERS) call t_stopf('FVM:Communication')
          93             :     !-----------------------------------------------------------------------------------!
          94             :     if(FVM_TIMERS) call t_startf('FVM:Unpack')
          95           0 :     do ie=nets,nete
          96           0 :        kptr = kmin-1
          97           0 :        call ghostunpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk   ,kptr,ie)
          98           0 :        do q=1,ntrac
          99           0 :           kptr = kptr + ksize
         100           0 :           call ghostunpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,:),   kblk,kptr,ie)
         101             :        enddo
         102             :     enddo
         103             :     if(FVM_TIMERS) call t_stopf('FVM:Unpack')
         104             :     if(FVM_TIMERS) call t_startf('FVM:freebuf')
         105           0 :     call freeghostbuffer(cellghostbuf)
         106             :     if(FVM_TIMERS) call t_stopf('FVM:freebuf')
         107           0 :   end subroutine fill_halo_fvm_noprealloc
         108             :
         109      738816 : subroutine fill_halo_fvm_prealloc(cellghostbuf,elem,fvm,hybrid,nets,nete,ndepth,kmin,kmax,ksize,active)
         110           0 :     use perf_mod, only : t_startf, t_stopf ! _EXTERNAL
         111             :     use dimensions_mod, only: nc, ntrac, nlev
         112             :     implicit none
         113             :     type (EdgeBuffer_t), intent(inout)       :: cellghostbuf
         114             :     type (element_t),intent(inout)            :: elem(:)
         115             :     type (fvm_struct),intent(inout)           :: fvm(:)
         116             :     type (hybrid_t),intent(in)                :: hybrid
         117             :
         118             :
         119             :     integer,intent(in)                        :: nets,nete
         120             :     integer,intent(in)                        :: ndepth     ! depth of halo
         121             :     integer,intent(in)                        :: kmin,kmax  ! min and max vertical level
         122             :     integer,intent(in)                        :: ksize      ! the total number of vertical
         123             :     logical, optional                         :: active     ! indicates if te current thread is active
         124             :     integer                                   :: ie,i1,i2,kblk,q,kptr
         125             :     !
         126             :     !
         127             :     logical :: lactive
         128             :
         129      738816 :     if(present(active)) then
         130      738816 :        lactive = active
         131             :     else
         132             :        lactive = .true.
         133             :     endif
         134             : !    call t_startf('FVM:initbuf')
         135      738816 :     i1=1-ndepth
         136      738816 :     i2=nc+ndepth
         137      738816 :     kblk = kmax-kmin+1
         138             :     if(FVM_TIMERS) call t_startf('FVM:pack')
         139      738816 :     if(lactive) then
         140     5933616 :     do ie=nets,nete
         141     5194800 :        kptr = kmin-1
         142  4192203600 :        call ghostpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk, kptr,ie)
         143    21518016 :        do q=1, ntrac
         144    15584400 :           kptr = kptr + ksize
         145 12581805600 :           call ghostpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,q) ,kblk,kptr,ie)
         146             :        enddo
         147             :     end do
         148             :     endif
         149             :     if(FVM_TIMERS) call t_stopf('FVM:pack')
         150             :     if(FVM_TIMERS) call t_startf('FVM:Communication')
         151      738816 :     call ghost_exchange(hybrid,cellghostbuf,location='fill_halo_fvm_prealloc')
         152             :     if(FVM_TIMERS) call t_stopf('FVM:Communication')
         153             :     !-----------------------------------------------------------------------------------!
         154             :     if(FVM_TIMERS) call t_startf('FVM:Unpack')
         155      738816 :     if(lactive) then
         156     5933616 :     do ie=nets,nete
         157     5194800 :        kptr = kmin-1
         158  8379212400 :        call ghostunpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk, kptr,ie)
         159    21518016 :        do q=1, ntrac
         160    15584400 :           kptr = kptr + ksize
         161 25142832000 :           call ghostunpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,q), kblk,kptr,ie)
         162             :        enddo
         163             :     enddo
         164             :     endif
         165             :     if(FVM_TIMERS) call t_stopf('FVM:Unpack')
         166             :
         167      738816 :   end subroutine fill_halo_fvm_prealloc
         168             :
         169             :  subroutine PrintArray(i1,i2,array)
         170             :    ! debug routine potentially called from any MPI rank
         171             :    integer :: i1,i2
         172             :    real(kind=r8) :: array(i1:i2,i1:i2)
         173             :    integer :: sz,i,ub
         174             :
         175             :    sz = size(array,dim=1)
         176             :
         177             :    if (sz == 9) then
         178             :      do i=i2,i1,-1
         179             :         write(6,9) array(-2,i),array(-1,i), array(0,i), &
         180             :                    array( 1,i), array(2,i), array(3,i), &
         181             :                    array( 4,i), array(5,i), array(6,i)
         182             :      enddo
         183             :    endif
         184             :
         185             :  9      format('|',9(f10.1,'|'))
         186             :
         187             :
         188      738816 :  end subroutine
         189             :
         190             :
         191           0 :   subroutine fill_halo_and_extend_panel(elem,fvm,fld,hybrid,nets,nete,nphys,nhcc, ndepth,numlev,num_flds,lfill_halo,lextend_panel)
         192             :     use hybrid_mod,             only: hybrid_t
         193             :     use edge_mod,               only: initghostbuffer, freeghostbuffer, ghostpack, ghostunpack
         194             :
         195             :     use fvm_reconstruction_mod, only: extend_panel_interpolate
         196             :     use cam_abortutils,         only: endrun
         197             :     use dimensions_mod,         only: fv_nphys,nhr,nhr_phys,nhc,nhc_phys,ns,ns_phys,nhe_phys,nc
         198             :     use perf_mod, only : t_startf, t_stopf ! _EXTERNAL
         199             :
         200             :     integer              , intent(in)    :: nets,nete,nphys,ndepth,numlev,num_flds,nhcc
         201             :     real (kind=r8)       , intent(inout) :: fld(1-nhcc:nphys+nhcc,1-nhcc:nphys+nhcc,numlev,num_flds,nets:nete)
         202             :     type (hybrid_t)      , intent(in)    :: hybrid  ! distributed parallel structure (shared)
         203             :     type (element_t)     , intent(inout) :: elem(:)
         204             :     type(fvm_struct)     , intent(in)    :: fvm(:)
         205             :     logical              , intent(in)    :: lfill_halo,lextend_panel
         206             : !    real (kind=r8)       , allocatable   :: ftmp(:,:)
         207             : !    real (kind=r8)               :: ftmp(1-nhcc:nphys+nhcc,1-nhcc:nphys+nhcc,numlev,num_flds,nets:nete)
         208           0 :     real (kind=r8), allocatable  :: fld_tmp(:,:)
         209             :
         210             :     integer                 :: ie,k,itr,nht_phys,nh_phys
         211           0 :     type (edgeBuffer_t)   :: cellghostbuf
         212             :
         213           0 :     if (lfill_halo) then
         214             :       !
         215             :       !*********************************************
         216             :       !
         217             :       ! halo exchange
         218             :       !
         219             :       !*********************************************
         220             :       !
         221           0 :       call t_startf('fill_halo_and_extend_panel initbuffer')
         222           0 :       call initghostbuffer(hybrid%par,cellghostbuf,elem,numlev*num_flds,nhcc,nphys)
         223           0 :       call t_stopf('fill_halo_and_extend_panel initbuffer')
         224           0 :       do ie=nets,nete
         225           0 :         call ghostpack(cellghostbuf, fld(:,:,:,:,ie),numlev*num_flds,0,ie)
         226             :       end do
         227           0 :       call ghost_exchange(hybrid,cellghostbuf,location='fill_halo_and_extend_panel')
         228           0 :       do ie=nets,nete
         229           0 :         call ghostunpack(cellghostbuf, fld(:,:,:,:,ie),numlev*num_flds,0,ie)
         230             :       end do
         231           0 :       call freeghostbuffer(cellghostbuf)
         232             :     end if
         233           0 :     if (lextend_panel) then
         234             :       !
         235             :       !*********************************************
         236             :       !
         237             :       ! extend panel
         238             :       !
         239             :       !*********************************************
         240             :       !
         241           0 :       if (nphys==fv_nphys) then
         242           0 :         if (ndepth>nhr_phys) &
         243           0 :              call endrun("fill_halo_and_extend_panel: ndepth>nhr_phys")
         244           0 :         nht_phys = nhe_phys+nhr_phys
         245           0 :         nh_phys  = nhr_phys
         246           0 :         allocate(fld_tmp(1-nht_phys:nphys+nht_phys,1-nht_phys:nphys+nht_phys))
         247           0 :         do ie=nets,nete
         248           0 :           do itr=1,num_flds
         249           0 :             do k=1,numlev
         250             :               call extend_panel_interpolate(fv_nphys,nhc_phys,nhr_phys,nht_phys,ns_phys,nh_phys,&
         251           0 :                    fld(:,:,k,itr,ie),fvm(ie)%cubeboundary,&
         252           0 :                    fvm(ie)%halo_interp_weight_physgrid(1:ns_phys,1-nh_phys:fv_nphys+nh_phys,1:nhr_phys,:),&
         253           0 :                    fvm(ie)%ibase_physgrid(1-nh_phys:fv_nphys+nh_phys,1:nhr_phys,:),&
         254           0 :                    fld_tmp)
         255           0 :               fld(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth,k,itr,ie) = fld_tmp(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth)
         256             :             end do
         257             :           end do
         258             :         end do
         259           0 :         deallocate(fld_tmp)
         260           0 :       else if (nphys==nc) then
         261           0 :         if (ndepth>nhr) &
         262           0 :              call endrun("fill_halo_and_extend_panel: ndepth>nhr")
         263           0 :         nhe_phys= 0
         264           0 :         nht_phys= nhe_phys+nhr
         265           0 :         nh_phys = nhr
         266           0 :         allocate(fld_tmp(1-nht_phys:nphys+nht_phys,1-nht_phys:nphys+nht_phys))
         267           0 :         do ie=nets,nete
         268           0 :           do itr=1,num_flds
         269           0 :             do k=1,numlev
         270             :               call extend_panel_interpolate(nc,nhc,nhr,nht_phys,ns,nh_phys,&
         271           0 :                    fld(:,:,k,itr,ie),fvm(ie)%cubeboundary,&
         272           0 :                    fvm(ie)%halo_interp_weight(1:ns,1-nh_phys:nc+nh_phys,1:nhr,:),&
         273             :                    fvm(ie)%ibase(1-nh_phys:nc+nh_phys,1:nhr,:),&
         274           0 :                    fld_tmp)
         275           0 :               fld(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth,k,itr,ie) = fld_tmp(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth)
         276             :             end do
         277             :           end do
         278             :         end do
         279           0 :         deallocate(fld_tmp)
         280             :       else
         281           0 :         call endrun("fill_halo_and_extend_panel: resolution not supported")
         282             :       end if
         283             :     end if
         284           0 :   end subroutine fill_halo_and_extend_panel
         285             :
         286             :
         287             :   ! initialize global buffers shared by all threads
         288        1536 :   subroutine fvm_init1(par,elem)
         289           0 :     use parallel_mod,           only: parallel_t
         290             :     use cam_abortutils,         only: endrun
         291             :     use cam_logfile,            only: iulog
         292             :     use control_mod,            only: rsplit
         293             :     use dimensions_mod,         only: qsize, qsize_d
         294             :     use dimensions_mod,         only: fvm_supercycling, fvm_supercycling_jet
         295             :     use dimensions_mod,         only: nc,nhe, nhc, nlev,ntrac, ntrac_d,ns, nhr, use_cslam
         296             :     use dimensions_mod,         only: large_Courant_incr
         297             :     use dimensions_mod,         only: kmin_jet,kmax_jet
         298             :
         299             :     type (parallel_t) :: par
         300             :     type (element_t),intent(inout)            :: elem(:)
         301             :     !
         302        1536 :     if (use_cslam) then
         303        1536 :       if (par%masterproc) then
         304           2 :         write(iulog,*) "                                           "
         305           2 :         write(iulog,*) "|-----------------------------------------|"
         306           2 :         write(iulog,*) "| FVM tracer transport scheme information |"
         307           2 :         write(iulog,*) "|-----------------------------------------|"
         308           2 :         write(iulog,*) "                                           "
         309             :       end if
         310        1536 :       if (use_cslam) then
         311        1536 :         if (par%masterproc) then
         312           2 :           write(iulog,*) "Running consistent SE-CSLAM, Lauritzen et al. (2017, MWR)."
         313           2 :           write(iulog,*) "CSLAM = Conservative Semi-LAgrangian Multi-tracer scheme"
         314           2 :           write(iulog,*) "Lauritzen et al., (2010), J. Comput. Phys."
         315           2 :           write(iulog,*) "  "
         316             :         end if
         317             :       end if
         318             :       !
         319             :       ! PARAMETER ERROR CHECKING
         320             :       !
         321        1536 :       if (kmin_jet>kmax_jet) &
         322           0 :            call endrun("PARAMETER ERROR for fvm: kmin_jet must be < kmax_jet")
         323        1536 :       if (ntrac>ntrac_d) &
         324           0 :            call endrun("PARAMETER ERROR for fvm: ntrac > ntrac_d")
         325             :
         326        1536 :       if (qsize>0.and.mod(rsplit,fvm_supercycling).ne.0) then
         327           0 :         if (par%masterproc) then
         328           0 :           write(iulog,*)'cannot supercycle fvm tracers with respect to se tracers'
         329           0 :           write(iulog,*)'with this choice of rsplit =',rsplit
         330           0 :           write(iulog,*)'rsplit must be a multiple of fvm_supercycling=',fvm_supercycling
         331             :         end if
         332           0 :         call endrun("PARAMETER ERROR for fvm: mod(rsplit,fvm_supercycling)<>0")
         333             :       endif
         334             :
         335        1536 :       if (qsize>0.and.mod(rsplit,fvm_supercycling_jet).ne.0) then
         336           0 :         if (par%masterproc) then
         337           0 :           write(iulog,*)'cannot supercycle fvm tracers with respect to se tracers'
         338           0 :           write(iulog,*)'with this choice of rsplit =',rsplit
         339           0 :           write(iulog,*)'rsplit must be a multiple of fvm_supercycling_jet=',fvm_supercycling_jet
         340             :         end if
         341           0 :         call endrun("PARAMETER ERROR for fvm: mod(rsplit,fvm_supercycling_jet)<>0")
         342             :       endif
         343             :
         344        1536 :       if (large_Courant_incr.and.(fvm_supercycling.ne.fvm_supercycling_jet)) then
         345           0 :         if (par%masterproc) then
         346           0 :           write(iulog,*)'Large Courant number increment requires no level dependent supercycling'
         347           0 :           write(iulog,*)'i.e. fvm_supercycling must be equal to fvm_supercycling_jet'
         348             :         end if
         349           0 :         call endrun("PARAMETER ERROR for fvm: large_courant_incr requires fvm_supercycling=fvm_supercycling_jet")
         350             :       endif
         351             :
         352        1536 :       if (par%masterproc) then
         353           2 :         write(iulog,*) "                                            "
         354           2 :         write(iulog,*) "Done Tracer transport scheme information    "
         355           2 :         write(iulog,*) "                                            "
         356             :       end if
         357             :
         358             :
         359        1536 :       if (par%masterproc) write(iulog,*) "fvm resolution is nc*nc in each element: nc = ",nc
         360        1536 :       if (par%masterproc) write(iulog,*)'ntrac,ntrac_d=',ntrac,ntrac_d
         361        1536 :       if (par%masterproc) write(iulog,*)'qsize,qsize_d=',qsize,qsize_d
         362             :
         363             :       if (nc.ne.3) then
         364             :         if (par%masterproc) then
         365             :           write(iulog,*) "Only nc==3 is supported for CSLAM"
         366             :         endif
         367             :         call endrun("PARAMETER ERRROR for fvm: only nc=3 supported for CSLAM")
         368             :       end if
         369             :
         370        1536 :       if (par%masterproc) then
         371           2 :         write(iulog,*) "  "
         372             :         if (ns==1) then
         373             :           write(iulog,*) "ns==1: using no interpolation for mapping cell averages values across edges"
         374             :           write(iulog,*) "Note: this is not a recommended setting - large errors at panel edges!"
         375             :         else if (ns==2) then
         376             :           write(iulog,*) "ns==2: using linear interpolation for mapping cell averages values across edges"
         377             :           write(iulog,*) "Note that ns=4 is default CSLAM setting used in Lauritzen et al. (2010)"
         378             :           write(iulog,*) "so this option is slightly less accurate (but the stencil is smaller near panel edges!)"
         379             :
         380             :         else if (ns==3) then
         381           2 :           write(iulog,*) "ns==3: using quadratic interpolation for mapping cell averages values across edges"
         382           2 :           write(iulog,*) "Note that ns=4 is default CSLAM setting used in Lauritzen et al. (2010)"
         383           2 :           write(iulog,*) "so this option is slightly less accurate (but the stencil is smaller near panel edges!)"
         384             :         else if (ns==4) then
         385             :           write(iulog,*) "ns==4: using cubic interpolation for mapping cell averages values across edges"
         386             :           write(iulog,*) "This is default CSLAM setting used in Lauritzen et al. (2010)"
         387             :         else
         388             :           write(iulog,*) "Not a tested value for ns but it should work! You choose ns = ",ns
         389             :         end if
         390             :
         391             :         !       if (ns.NE.3) then
         392             :         !         write(*,*) "In fvm_reconstruction_mod function matmul_w has been hard-coded for ns=3 for performance"
         393             :         !         write(*,*) "Revert to general code - outcommented above"
         394             :         !         call endrun("stopping")
         395             :         !       end if
         396             :       end if
         397             :
         398             :       if (MOD(ns,2)==0.and.nhr+(nhe-1)+ns/2>nc+nc) then
         399             :         write(iulog,*) "to run this combination of ns and nhr you need to increase nc to ",nhr+ns/2+nhe-1
         400             :         write(iulog,*) "You choose (ns,nhr,nc,nhe)=",ns,nhr,nc,nhe
         401             :         call endrun("stopping")
         402             :       end if
         403             :       if (MOD(ns,2)==1.and.nhr+(ns-1)/2+(nhe-1)>nc+nc) then
         404             :         write(iulog,*) "to run this combination of ns and nhr you need to increase nc to ",nhr+(ns-1)/2+nhe-1
         405             :         write(iulog,*) "You choose (ns,nhr,nc,nhe)=",ns,nhr,nc,nhe
         406             :         call endrun("stopping")
         407             :       end if
         408             :
         409             :       if (nc==3.and.ns.ne.3) then
         410             :         if (par%masterproc) then
         411             :           write(iulog,*) "Recommended setting for nc=3 is ns=3 (linear interpolation in halo)"
         412             :           write(iulog,*) "You choose ns=",ns
         413             :           write(iulog,*) "Goto dimensions_mod to change value of ns"
         414             :           write(iulog,*) "or outcomment call haltmop below (i.e. you know what you are doing!)"
         415             :         endif
         416             :         call endrun("stopping")
         417             :       end if
         418             :     end if
         419             :
         420             :     if (nc==4.and.ns.ne.4) then
         421             :       if (par%masterproc) then
         422             :         write(iulog,*) "Recommended setting for nc=4 is ns=4 (cubic interpolation in halo)"
         423             :         write(iulog,*) "You choose ns=",ns
         424             :         write(iulog,*) "Goto dimensions_mod to change value of ns"
         425             :         write(iulog,*) "or outcomment call haltmop below (i.e. you know what you are doing!)"
         426             :       endif
         427             :       call endrun("stopping")
         428             :     end if
         429             :
         430             :     if (nhe .ne. 1) then
         431             :       if (par%masterproc) then
         432             :         write(iulog,*) "PARAMETER ERROR for fvm: Number of halo zone for the extended"
         433             :         write(iulog,*) "element nhe has to be 1, only this is available now! STOP!"
         434             :       endif
         435             :       call endrun("stopping")
         436             :     end if
         437        1536 :   end subroutine fvm_init1
         438             :
         439             :
         440             :
         441             :
         442             :
         443             :   ! initialization that can be done in threaded regions
         444        1536 :   subroutine fvm_init2(elem,fvm,hybrid,nets,nete)
         445             :     use fvm_control_volume_mod, only: fvm_mesh,fvm_set_cubeboundary
         446             :     use bndry_mod,              only: compute_ghost_corner_orientation
         447             :     use dimensions_mod,         only: nlev, nc, nhc, nhe, ntrac, ntrac_d, np
         448             :     use dimensions_mod,         only: nhc_phys, fv_nphys
         449             :     use dimensions_mod,         only: fvm_supercycling, fvm_supercycling_jet
         450             :     use dimensions_mod,         only: kmin_jet,kmax_jet
         451             :     use hycoef,                 only: hyai, hybi, ps0
         452             :     use derivative_mod,         only: subcell_integration
         453             :     use air_composition,        only: thermodynamic_active_species_num
         454             :
         455             :     type (fvm_struct) :: fvm(:)
         456             :     type (element_t)  :: elem(:)
         457             :     type (hybrid_t)   :: hybrid
         458             :     integer           :: ie,nets,nete,k,klev
         459             :     real(kind=r8)     :: one(np,np)
         460             :
         461       32256 :     one = 1.0_r8
         462       12336 :     do ie=nets,nete
         463      293136 :       do k = 1, nlev
         464      280800 :         fvm(ie)%dp_ref(k)         = ( hyai(k+1) - hyai(k) )*ps0 + ( hybi(k+1) - hybi(k) )*ps0
         465      291600 :         fvm(ie)%dp_ref_inverse(k) = 1.0_r8/fvm(ie)%dp_ref(k)
         466             :       end do
         467             :     end do
         468             :
         469        1536 :     call compute_ghost_corner_orientation(hybrid,elem,nets,nete)
         470             :     ! run some tests:
         471             :     !    call test_ghost(hybrid,elem,nets,nete)
         472             :
         473       12336 :     do ie=nets,nete
         474       10800 :       call fvm_set_cubeboundary(elem(ie),fvm(ie))
         475             :       call fvm_mesh(elem(ie),fvm(ie))
         476      140400 :       fvm(ie)%inv_area_sphere    = 1.0_r8/fvm(ie)%area_sphere
         477             :       !
         478             :       ! compute CSLAM areas consistent with SE area (at 1 degree they can be up to
         479             :       ! 1E-6 different than the correct spherical areas used in CSLAM)
         480             :       !
         481       10800 :       call subcell_integration(one, np, nc, elem(ie)%metdet,fvm(ie)%inv_se_area_sphere)
         482      140400 :       fvm(ie)%inv_se_area_sphere = 1.0_r8/fvm(ie)%inv_se_area_sphere
         483             :
         484    36622800 :       fvm(ie)%fc(:,:,:,:) = 0.0_r8
         485    51397200 :       fvm(ie)%fm(:,:,:,:) = 0.0_r8
         486    25575936 :       fvm(ie)%ft(:,:,:  ) = 0.0_r8
         487             :     enddo
         488             :     ! Need to allocate ghostBufQnhc after compute_ghost_corner_orientation because it
         489             :     ! changes the values for reverse
         490             :
         491        1536 :     call initghostbuffer(hybrid%par,ghostBufQnhc_s,elem,nlev*(ntrac+1),nhc,nc,nthreads=1)
         492        1536 :     call initghostbuffer(hybrid%par,ghostBufQnhc_t1,elem,nlev, nhc,nc,nthreads=1)
         493        1536 :     call initghostbuffer(hybrid%par,ghostBufQnhc_h,elem,nlev*(ntrac+1),nhc,nc,nthreads=horz_num_threads)
         494        1536 :     call initghostbuffer(hybrid%par,ghostBufQnhc_vh,elem,nlev*(ntrac+1),nhc,nc,nthreads=vert_num_threads*horz_num_threads)
         495        1536 :     klev = kmax_jet-kmin_jet+1
         496        1536 :     call initghostbuffer(hybrid%par,ghostBufQ1_h,elem,klev*(ntrac+1),1,nc,nthreads=horz_num_threads)
         497        1536 :     call initghostbuffer(hybrid%par,ghostBufQ1_vh,elem,klev*(ntrac+1),1,nc,nthreads=vert_num_threads*horz_num_threads)
         498             : !    call initghostbuffer(hybrid%par,ghostBufFlux_h,elem,4*nlev,nhe,nc,nthreads=horz_num_threads)
         499        1536 :     call initghostbuffer(hybrid%par,ghostBufFlux_vh,elem,4*nlev,nhe,nc,nthreads=vert_num_threads*horz_num_threads)
         500        1536 :     call initghostbuffer(hybrid%par,ghostBuf_cslam2gll,elem,nlev*thermodynamic_active_species_num,nhc,nc,nthreads=1)
         501             :     !
         502             :     ! preallocate buffers for physics-dynamics coupling
         503             :     !
         504        1536 :     if (fv_nphys.ne.nc) then
         505           0 :        call initghostbuffer(hybrid%par,ghostBufPG_s,elem,nlev*(4+ntrac),nhc_phys,fv_nphys,nthreads=1)
         506             :     else
         507        1536 :        call initghostbuffer(hybrid%par,ghostBufPG_s,elem,nlev*3,nhc_phys,fv_nphys,nthreads=1)
         508             :     end if
         509             :
         510        1536 :     if (fvm_supercycling.ne.fvm_supercycling_jet) then
         511             :       !
         512             :       ! buffers for running different fvm time-steps in the jet region
         513             :       !
         514           0 :       klev = kmax_jet-kmin_jet+1
         515           0 :       call initghostbuffer(hybrid%par,ghostBufQnhcJet_h,elem,klev*(ntrac+1),nhc,nc,nthreads=horz_num_threads)
         516           0 :       call initghostbuffer(hybrid%par,ghostBufFluxJet_h,elem,4*klev,nhe,nc,nthreads=horz_num_threads)
         517             :     end if
         518        1536 :   end subroutine fvm_init2
         519             :
         520             :
         521        1536 :   subroutine fvm_init3(elem,fvm,hybrid,nets,nete,irecons)
         522        1536 :     use control_mod     ,       only: neast, nwest, seast, swest
         523             :     use fvm_analytic_mod,       only: compute_reconstruct_matrix
         524             :     use dimensions_mod  ,       only: fv_nphys, use_cslam
         525             :     use dimensions_mod,         only: nlev, nc, nhe, nlev, nhc
         526             :     use coordinate_systems_mod, only: cartesian2D_t,cartesian3D_t
         527             :     use coordinate_systems_mod, only: cubedsphere2cart, cart2cubedsphere
         528             :     implicit none
         529             :     type (element_t) ,intent(inout)  :: elem(:)
         530             :     type (fvm_struct),intent(inout)  :: fvm(:)
         531             :     type (hybrid_t)  ,intent(in)     :: hybrid
         532             :     integer          ,intent(in)     :: nets,nete,irecons
         533             :     !
         534        1536 :     type (edgeBuffer_t)     :: cellghostbuf
         535             :     integer                 :: ie, ixy, ivertex, i, j,istart,itot,ishft,imin,imax
         536             :     integer, dimension(2,4) :: unit_vec
         537             :     integer                 :: rot90_matrix(2,2), iside
         538             :
         539             :     type (cartesian2D_t)                :: tmpgnom
         540             :     type (cartesian2D_t)                :: gnom
         541             :     type(cartesian3D_t)                 :: tmpcart3d
         542             :
         543        1536 :     if (use_cslam.and.nc.ne.fv_nphys) then
         544             :       !
         545             :       ! fill the fvm halo for mapping in d_p_coupling if
         546             :       ! physics grid resolution is different than fvm resolution
         547             :       !
         548           0 :       call fill_halo_fvm(elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
         549             :     end if
         550             :
         551             :
         552        1536 :     imin=1-nhc
         553        1536 :     imax=nc+nhc
         554             :     !
         555             :     ! fill halo start
         556             :     !
         557        1536 :     itot=9+irecons-1+2
         558        1536 :     call initghostbuffer(hybrid%par,cellghostbuf,elem,itot,nhc,nc)
         559       12336 :     do ie=nets,nete
         560       10800 :       istart = 0
         561      982800 :       call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord(1,:,:),1,istart,ie)
         562       10800 :       istart = istart+1
         563      982800 :       call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord(2,:,:),1,istart,ie)
         564       10800 :       istart = istart+1
         565       32400 :       do ixy=1,2
         566      118800 :         do ivertex=1,4
         567     7862400 :           call ghostpack(cellghostbuf, fvm(ie)%vtx_cart(ivertex,ixy,:,:) ,1,istart,ie)
         568      108000 :           istart = istart+1
         569             :         end do
         570             :       end do
         571      982800 :       call ghostpack(cellghostbuf, fvm(ie)%flux_orient(1,:,:) ,1,istart,ie)
         572       66336 :       do ixy=1,irecons-1
         573       54000 :         istart=istart+1
         574     4924800 :         call ghostpack(cellghostbuf, fvm(ie)%spherecentroid(ixy,:,:) ,1,istart,ie)
         575             :       end do
         576             :     end do
         577        1536 :     call ghost_exchange(hybrid,cellghostbuf,location='fvm_init3')
         578       12336 :     do ie=nets,nete
         579       10800 :       istart = 0
         580     1954800 :       call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord(1,:,:),1,istart,ie)
         581       10800 :       istart = istart+1
         582     1954800 :       call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord(2,:,:),1,istart,ie)
         583       10800 :       istart = istart+1
         584       32400 :       do ixy=1,2
         585      118800 :         do ivertex=1,4
         586    15638400 :           call ghostunpack(cellghostbuf, fvm(ie)%vtx_cart(ivertex,ixy,:,:) ,1,istart,ie)
         587      108000 :           istart = istart+1
         588             :         end do
         589             :       end do
         590     1954800 :       call ghostunpack(cellghostbuf, fvm(ie)%flux_orient(1,:,:) ,1,istart,ie)
         591       66336 :       do ixy=1,irecons-1
         592       54000 :         istart=istart+1
         593     9784800 :         call ghostunpack(cellghostbuf, fvm(ie)%spherecentroid(ixy,:,:) ,1,istart,ie)
         594             :       end do
         595             :     enddo
         596        1536 :     call freeghostbuffer(cellghostbuf)
         597             :     !
         598             :     ! indicator for non-existing cells
         599             :     ! set vtx_cart to corner value in non-existent cells
         600             :     !
         601       12336 :     do ie=nets,nete
         602       12336 :        if (fvm(ie)%cubeboundary==nwest) then
         603         372 :          fvm(ie)%flux_orient     (:  ,1-nhc      :0     ,nc      +1 :nc      +nhc      ) = -1
         604         696 :          fvm(ie)%spherecentroid  (:,    1-nhc      :0     ,nc      +1 :nc      +nhc    ) = -1e5_r8
         605         588 :          fvm(ie)%vtx_cart(:,1,1-nhc:0     ,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(4,1,1,nc)
         606         588 :          fvm(ie)%vtx_cart(:,2,1-nhc:0     ,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(4,2,1,nc)
         607       10788 :        else if (fvm(ie)%cubeboundary==swest) then
         608         372 :          fvm(ie)%flux_orient     (:,1-nhc      :0     ,1-nhc      :0   ) = -1
         609         696 :          fvm(ie)%spherecentroid  (:,1-nhc      :0     ,1-nhc      :0   ) = -1e5_r8
         610         588 :          fvm(ie)%vtx_cart(:,1,1-nhc:0     ,1-nhc:0     ) = fvm(ie)%vtx_cart(1,1,1,1)
         611         588 :          fvm(ie)%vtx_cart(:,2,1-nhc:0     ,1-nhc:0     ) = fvm(ie)%vtx_cart(1,2,1,1)
         612       10776 :        else if (fvm(ie)%cubeboundary==neast) then
         613         372 :          fvm(ie)%flux_orient     (:,nc      +1 :nc      +nhc      ,nc      +1 :nc      +nhc    ) = -1
         614         696 :          fvm(ie)%spherecentroid  (:,nc      +1 :nc      +nhc      ,nc      +1 :nc      +nhc    ) = -1e5_r8
         615         588 :          fvm(ie)%vtx_cart(:,1,nc+1 :nc+nhc,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(3,1,nc,nc)
         616         588 :          fvm(ie)%vtx_cart(:,2,nc+1 :nc+nhc,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(3,2,nc,nc)
         617       10764 :        else if (fvm(ie)%cubeboundary==seast) then
         618         372 :          fvm(ie)%flux_orient     (:,nc      +1 :nc      +nhc      ,1-nhc      :0   ) = -1
         619         696 :          fvm(ie)%spherecentroid  (:,nc      +1 :nc      +nhc      ,1-nhc      :0   ) = -1e5_r8
         620         588 :          fvm(ie)%vtx_cart(:,1,nc+1 :nc+nhc,1-nhc:0     ) = fvm(ie)%vtx_cart(2,1,nc,1)
         621         588 :          fvm(ie)%vtx_cart(:,2,nc+1 :nc+nhc,1-nhc:0     ) = fvm(ie)%vtx_cart(2,2,nc,1)
         622             :        end if
         623             :      end do
         624             :
         625             :      !
         626             :      ! set vectors for perpendicular flux vector
         627             :      !
         628        1536 :      rot90_matrix(1,1) = 0; rot90_matrix(2,1) =  1 !counter-clockwise rotation matrix
         629        1536 :      rot90_matrix(1,2) =-1; rot90_matrix(2,2) =  0 !counter-clockwise rotation matrix
         630             :
         631        1536 :      iside = 1
         632        1536 :      unit_vec(1,iside) = 0 !x-component of displacement vector for side 1
         633        1536 :      unit_vec(2,iside) = 1 !y-component of displacement vector for side 1
         634             :
         635        6144 :      do iside=2,4
         636        6144 :        unit_vec(:,iside) = MATMUL(rot90_matrix(:,:),unit_vec(:,iside-1))
         637             :      end do
         638             :
         639             :      !
         640             :      ! fill halo done
         641             :      !
         642             :      !-------------------------------
         643             :
         644       12336 :      do ie=nets,nete
         645     3942000 :        fvm(ie)%displ_max = 0.0_r8
         646      109536 :        do j=imin,imax
         647      982800 :          do i=imin,imax
         648             :            !
         649             :            ! rotate gnomonic coordinate vector
         650             :            !
         651             :            !           fvm(ie)%norm_elem_coord(:,i,j) = MATMUL(fvm(ie)%rot_matrix(:,:,i,j),fvm(ie)%norm_elem_coord(:,i,j))
         652             :            !
         653      874800 :            ishft = NINT(fvm(ie)%flux_orient(2,i,j))
         654     2721600 :            do ixy=1,2
         655             :              !
         656             :              ! rotate coordinates if needed through permutation
         657             :              !
         658     8748000 :              fvm(ie)%vtx_cart(1:4,ixy,i,j) = cshift(fvm(ie)%vtx_cart(1:4,ixy,i,j),shift=ishft)
         659     8748000 :              fvm(ie)%flux_vec        (ixy,i,j,1:4) = cshift(unit_vec                (ixy,1:4    ),shift=ishft)
         660             :              !
         661             :              ! set flux vector to zero in non-existent cells (corner halo)
         662             :              !
         663     8748000 :              fvm(ie)%flux_vec        (ixy,i,j,1:4) = fvm(ie)%ifct(i,j)*fvm(ie)%flux_vec(ixy,i,j,1:4)
         664             :
         665     1749600 :              iside=1
         666             :              fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
         667     1749600 :                   ABS(fvm(ie)%vtx_cart(4,ixy,i,j)-fvm(ie)%vtx_cart(1,ixy,i,j))
         668     1749600 :              iside=2
         669             :              fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
         670     1749600 :                   ABS(fvm(ie)%vtx_cart(1,ixy,i,j)-fvm(ie)%vtx_cart(2,ixy,i,j))
         671     1749600 :              iside=3
         672             :              fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
         673     1749600 :                   ABS(fvm(ie)%vtx_cart(2,ixy,i,j)-fvm(ie)%vtx_cart(3,ixy,i,j))
         674     1749600 :              iside=4
         675             :              fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
         676     2624400 :                   ABS(fvm(ie)%vtx_cart(2,ixy,i,j)-fvm(ie)%vtx_cart(1,ixy,i,j))
         677             :            end do
         678             :          end do
         679             :        end do
         680             :      end do
         681             :      !
         682             :      ! pre-compute derived metric terms used for integration, polynomial
         683             :      ! evaluation at fvm cell vertices, etc.
         684             :      !
         685       12336 :      do ie=nets,nete
         686       10800 :        call compute_reconstruct_matrix(nc,nhe,nhc,irecons,fvm(ie)%dalpha,fvm(ie)%dbeta,&
         687             :            fvm(ie)%spherecentroid,fvm(ie)%vtx_cart,fvm(ie)%centroid_stretch,&
         688       23136 :            fvm(ie)%vertex_recons_weights,fvm(ie)%recons_metrics,fvm(ie)%recons_metrics_integral)
         689             :      end do
         690             :       !
         691             :       ! create a normalized element coordinate system with a halo
         692             :       !
         693       12336 :      do ie=nets,nete
         694      109536 :        do j=1-nhc,nc+nhc
         695      982800 :          do i=1-nhc,nc+nhc
         696             :            !
         697             :            ! only compute for physically existent cells
         698             :            !
         699      972000 :            if (fvm(ie)%ifct(i,j)>0) then
         700      874368 :              gnom%x = fvm(ie)%norm_elem_coord(1,i,j)
         701      874368 :              gnom%y = fvm(ie)%norm_elem_coord(2,i,j)
         702             :              !
         703             :              ! coordinate transform only necessary for points on another panel
         704             :              !
         705      874368 :              if (NINT(fvm(ie)%flux_orient(1,1,1)).NE.NINT(fvm(ie)%flux_orient(1,i,j))) then
         706       38016 :                tmpcart3d=cubedsphere2cart(gnom,NINT(fvm(ie)%flux_orient(1,i,j)))
         707       38016 :                tmpgnom=cart2cubedsphere(tmpcart3d,NINT(fvm(ie)%flux_orient(1,1,1)))
         708             :              else
         709      836352 :                tmpgnom%x = fvm(ie)%norm_elem_coord(1,i,j)
         710      836352 :                tmpgnom%y = fvm(ie)%norm_elem_coord(2,i,j)
         711             :              end if
         712             :              !
         713             :              ! convert to element normalized coordinates
         714             :              !
         715      874368 :              fvm(ie)%norm_elem_coord(1,i,j) =(tmpgnom%x-elem(ie)%corners(1)%x)/&
         716      874368 :                   (0.5_r8*dble(nc)*fvm(ie)%dalpha)-1.0_r8
         717             :              fvm(ie)%norm_elem_coord(2,i,j) =(tmpgnom%y-elem(ie)%corners(1)%y)/&
         718      874368 :                   (0.5_r8*dble(nc)*fvm(ie)%dalpha)-1.0_r8
         719             :            else
         720         432 :              fvm(ie)%norm_elem_coord(1,i,j) = 1D9
         721         432 :              fvm(ie)%norm_elem_coord(2,i,j) = 1D9
         722             :            end if
         723             :          end do
         724             :        end do
         725             :      end do
         726             :
         727        1536 :    end subroutine fvm_init3
         728             :
         729             :
         730        1536 :   subroutine fvm_pg_init(elem, fvm, hybrid, nets, nete,irecons)
         731             :     use coordinate_systems_mod, only : cartesian2D_t,cartesian3D_t
         732             :     use control_mod, only : neast, nwest, seast, swest
         733             :     use coordinate_systems_mod, only : cubedsphere2cart, cart2cubedsphere
         734             :     use dimensions_mod, only: fv_nphys, nhe_phys,nhc_phys
         735             :     use cube_mod       ,only: dmap
         736             :     use control_mod    ,only: cubed_sphere_map
         737             :     use fvm_analytic_mod, only: compute_reconstruct_matrix
         738             :
         739             :     type (element_t) , intent(in)    :: elem(:)
         740             :     type (fvm_struct), intent(inout) :: fvm(:)
         741             :     type (hybrid_t)  , intent(in)    :: hybrid
         742             :
         743             :     type (cartesian2D_t)                :: gnom
         744             :     type(cartesian3D_t)                 :: tmpcart3d
         745             :     type (cartesian2D_t)                :: tmpgnom
         746             :
         747             :
         748             :     integer, intent(in)                     :: nets  ! starting thread element number (private)
         749             :     integer, intent(in)                     :: nete,irecons  ! ending thread element number   (private)
         750             :
         751             :     ! ==================================
         752             :     ! Local variables
         753             :     ! ==================================
         754             :
         755             :     integer                 :: ie, ixy, ivertex, i, j,istart,itot,ishft,imin,imax
         756             :     integer, dimension(2,4) :: unit_vec
         757             :     integer                 :: rot90_matrix(2,2), iside
         758             :
         759        1536 :     type (edgeBuffer_t)                     :: cellghostbuf
         760             :
         761             :     ! D is derivative of gnomonic mapping
         762        3072 :     real (kind=r8)              :: D(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,2,2)
         763             :     real (kind=r8)              :: detD,x1,x2
         764             :
         765        1536 :     if (fv_nphys>0) then
         766             :       !
         767             :       ! do the same as fvm_init3 for the metric terms of physgrid
         768             :       !
         769        1536 :       imin=1-nhc_phys
         770        1536 :       imax=fv_nphys+nhc_phys
         771             :       !
         772             :       ! fill halo start
         773             :       !
         774        1536 :       itot=9+irecons-1+2
         775        1536 :       call initghostbuffer(hybrid%par,cellghostbuf,elem,itot,nhc_phys,fv_nphys)
         776       12336 :       do ie=nets,nete
         777       10800 :         istart = 0
         778      982800 :         call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(1,:,:),1,istart,ie)
         779       10800 :         istart = istart+1
         780      982800 :         call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(2,:,:),1,istart,ie)
         781       10800 :         istart = istart+1
         782       32400 :         do ixy=1,2
         783      118800 :           do ivertex=1,4
         784     7862400 :             call ghostpack(cellghostbuf, fvm(ie)%vtx_cart_physgrid(ivertex,ixy,:,:) ,1,istart,ie)
         785      108000 :             istart = istart+1
         786             :           end do
         787             :         end do
         788      982800 :         call ghostpack(cellghostbuf, fvm(ie)%flux_orient_physgrid(1,:,:) ,1,istart,ie)
         789       66336 :         do ixy=1,irecons-1
         790       54000 :           istart=istart+1
         791     4924800 :           call ghostpack(cellghostbuf, fvm(ie)%spherecentroid_physgrid(ixy,:,:) ,1,istart,ie)
         792             :         end do
         793             :       end do
         794        1536 :       call ghost_exchange(hybrid,cellghostbuf,location='fvm_pg_init')
         795       12336 :       do ie=nets,nete
         796       10800 :         istart = 0
         797     1954800 :         call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(1,:,:),1,istart,ie)
         798       10800 :         istart = istart+1
         799     1954800 :         call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(2,:,:),1,istart,ie)
         800       10800 :         istart = istart+1
         801       32400 :         do ixy=1,2
         802      118800 :           do ivertex=1,4
         803    15638400 :             call ghostunpack(cellghostbuf, fvm(ie)%vtx_cart_physgrid(ivertex,ixy,:,:) ,1,istart,ie)
         804      108000 :             istart = istart+1
         805             :           end do
         806             :         end do
         807     1954800 :         call ghostunpack(cellghostbuf, fvm(ie)%flux_orient_physgrid(1,:,:) ,1,istart,ie)
         808       66336 :         do ixy=1,irecons-1
         809       54000 :           istart=istart+1
         810     9784800 :           call ghostunpack(cellghostbuf, fvm(ie)%spherecentroid_physgrid(ixy,:,:) ,1,istart,ie)
         811             :         end do
         812             :       enddo
         813        1536 :       call freeghostbuffer(cellghostbuf)
         814             :       !
         815             :       ! indicator for non-existing cells
         816             :       ! set vtx_cart to corner value in non-existent cells
         817             :       !
         818       12336 :       do ie=nets,nete
         819       12336 :         if (fvm(ie)%cubeboundary==nwest) then
         820         372 :           fvm(ie)%flux_orient_physgrid   (:  ,1-nhc_phys      :0     ,fv_nphys      +1 :fv_nphys      +nhc_phys    ) = -1
         821         696 :           fvm(ie)%spherecentroid_physgrid(:,  1-nhc_phys      :0     ,fv_nphys      +1 :fv_nphys      +nhc_phys    ) = -1e5_r8
         822           0 :           fvm(ie)%vtx_cart_physgrid(:,1,1-nhc_phys:0     ,fv_nphys+1 :fv_nphys+nhc_phys) = &
         823         588 :                fvm(ie)%vtx_cart_physgrid(4,1,1,fv_nphys)
         824           0 :           fvm(ie)%vtx_cart_physgrid(:,2,1-nhc_phys:0     ,fv_nphys+1 :fv_nphys+nhc_phys) = &
         825         588 :                fvm(ie)%vtx_cart_physgrid(4,2,1,fv_nphys)
         826       10788 :         else if (fvm(ie)%cubeboundary==swest) then
         827         372 :           fvm(ie)%flux_orient_physgrid   (:,1-nhc_phys      :0     ,1-nhc_phys      :0   ) = -1
         828         696 :           fvm(ie)%spherecentroid_physgrid(:,1-nhc_phys      :0     ,1-nhc_phys      :0   ) = -1e5_r8
         829         588 :           fvm(ie)%vtx_cart_physgrid(:,1,1-nhc_phys:0     ,1-nhc_phys:0     ) = fvm(ie)%vtx_cart_physgrid(1,1,1,1)
         830         588 :           fvm(ie)%vtx_cart_physgrid(:,2,1-nhc_phys:0     ,1-nhc_phys:0     ) = fvm(ie)%vtx_cart_physgrid(1,2,1,1)
         831       10776 :         else if (fvm(ie)%cubeboundary==neast) then
         832           0 :           fvm(ie)%flux_orient_physgrid   (:,fv_nphys      +1 :fv_nphys      +nhc_phys      , &
         833         372 :                fv_nphys      +1 :fv_nphys      +nhc_phys      ) = -1
         834           0 :           fvm(ie)%spherecentroid_physgrid(:,fv_nphys      +1 :fv_nphys      +nhc_phys      , &
         835         696 :                fv_nphys      +1 :fv_nphys      +nhc_phys      ) = -1e5_r8
         836           0 :           fvm(ie)%vtx_cart_physgrid(:,1,fv_nphys+1 :fv_nphys+nhc_phys,fv_nphys+1 :fv_nphys+nhc_phys) = &
         837         588 :                fvm(ie)%vtx_cart_physgrid(3,1,fv_nphys,fv_nphys)
         838           0 :           fvm(ie)%vtx_cart_physgrid(:,2,fv_nphys+1 :fv_nphys+nhc_phys,fv_nphys+1 :fv_nphys+nhc_phys) = &
         839         588 :                fvm(ie)%vtx_cart_physgrid(3,2,fv_nphys,fv_nphys)
         840       10764 :         else if (fvm(ie)%cubeboundary==seast) then
         841         372 :           fvm(ie)%flux_orient_physgrid   (:,fv_nphys      +1 :fv_nphys      +nhc_phys      ,1-nhc_phys      :0   ) = -1
         842         696 :           fvm(ie)%spherecentroid_physgrid(:,fv_nphys      +1 :fv_nphys      +nhc_phys      ,1-nhc_phys      :0   ) = -1e5_r8
         843           0 :           fvm(ie)%vtx_cart_physgrid(:,1,fv_nphys+1 :fv_nphys+nhc_phys,1-nhc_phys:0     ) = &
         844         588 :                fvm(ie)%vtx_cart_physgrid(2,1,fv_nphys,1)
         845           0 :           fvm(ie)%vtx_cart_physgrid(:,2,fv_nphys+1 :fv_nphys+nhc_phys,1-nhc_phys:0     ) = &
         846         588 :                fvm(ie)%vtx_cart_physgrid(2,2,fv_nphys,1)
         847             :         end if
         848             :       end do
         849             :
         850             :       !
         851             :       ! set vectors for perpendicular flux vector
         852             :       !
         853        1536 :       rot90_matrix(1,1) = 0; rot90_matrix(2,1) =  1 !counter-clockwise rotation matrix
         854        1536 :       rot90_matrix(1,2) =-1; rot90_matrix(2,2) =  0 !counter-clockwise rotation matrix
         855             :
         856        1536 :       iside = 1
         857        1536 :       unit_vec(1,iside) = 0 !x-component of displacement vector for side 1
         858        1536 :       unit_vec(2,iside) = 1 !y-component of displacement vector for side 1
         859             :
         860        6144 :       do iside=2,4
         861        6144 :         unit_vec(:,iside) = MATMUL(rot90_matrix(:,:),unit_vec(:,iside-1))
         862             :       end do
         863             :
         864             :       !
         865             :       ! fill halo done
         866             :       !
         867             :       !-------------------------------
         868             :
         869       12336 :       do ie=nets,nete
         870      109536 :         do j=imin,imax
         871      982800 :           do i=imin,imax
         872             :             !
         873             :             ! rotate gnomonic coordinate vector
         874             :             !
         875      874800 :             ishft = NINT(fvm(ie)%flux_orient_physgrid(2,i,j))
         876     2721600 :             do ixy=1,2
         877             :               !
         878             :               ! rotate coordinates if needed through permutation
         879             :               !
         880     9622800 :               fvm(ie)%vtx_cart_physgrid(1:4,ixy,i,j) = cshift(fvm(ie)%vtx_cart_physgrid(1:4,ixy,i,j),shift=ishft)
         881             :             end do
         882             :           end do
         883             :         end do
         884             :       end do
         885             :       !
         886             :       ! pre-compute derived metric terms used for integration, polynomial
         887             :       ! evaluation at fvm cell vertices, etc.
         888             :       !
         889       12336 :       do ie=nets,nete
         890       10800 :         call compute_reconstruct_matrix(fv_nphys,nhe_phys,nhc_phys,irecons,fvm(ie)%dalpha_physgrid,fvm(ie)%dbeta_physgrid,&
         891             :              fvm(ie)%spherecentroid_physgrid,fvm(ie)%vtx_cart_physgrid,fvm(ie)%centroid_stretch_physgrid,&
         892       23136 :              fvm(ie)%vertex_recons_weights_physgrid,fvm(ie)%recons_metrics_physgrid,fvm(ie)%recons_metrics_integral_physgrid)
         893             :       end do
         894             :       !
         895             :       ! code specific for physgrid
         896             :       !
         897             :       !
         898             :       ! create a normalized element coordinate system with a halo
         899             :       !
         900       12336 :       do ie=nets,nete
         901      109536 :         do j=1-nhc_phys,fv_nphys+nhc_phys
         902      982800 :           do i=1-nhc_phys,fv_nphys+nhc_phys
         903             :             !
         904             :             ! only compute for physically existent cells
         905             :             !
         906      972000 :             if (fvm(ie)%ifct_physgrid(i,j)>0) then
         907      874368 :               gnom%x = fvm(ie)%norm_elem_coord_physgrid(1,i,j)
         908      874368 :               gnom%y = fvm(ie)%norm_elem_coord_physgrid(2,i,j)
         909             :               !
         910             :               ! coordinate transform only necessary for points on another panel
         911             :               !
         912      874368 :               if (NINT(fvm(ie)%flux_orient_physgrid(1,1,1)).NE.NINT(fvm(ie)%flux_orient_physgrid(1,i,j))) then
         913       38016 :                 tmpcart3d=cubedsphere2cart(gnom,NINT(fvm(ie)%flux_orient_physgrid(1,i,j)))
         914       38016 :                 tmpgnom=cart2cubedsphere(tmpcart3d,NINT(fvm(ie)%flux_orient_physgrid(1,1,1)))
         915             :               else
         916      836352 :                 tmpgnom%x = fvm(ie)%norm_elem_coord_physgrid(1,i,j)
         917      836352 :                 tmpgnom%y = fvm(ie)%norm_elem_coord_physgrid(2,i,j)
         918             :               end if
         919             :               !
         920             :               ! convert to element normalized coordinates
         921             :               !
         922     1748736 :               fvm(ie)%norm_elem_coord_physgrid(1,i,j) =(tmpgnom%x-elem(ie)%corners(1)%x)/&
         923     1748736 :                    (0.5_r8*dble(fv_nphys)*fvm(ie)%dalpha_physgrid)-1.0_r8
         924           0 :               fvm(ie)%norm_elem_coord_physgrid(2,i,j) =(tmpgnom%y-elem(ie)%corners(1)%y)/&
         925      874368 :                    (0.5_r8*dble(fv_nphys)*fvm(ie)%dalpha_physgrid)-1.0_r8
         926             :             else
         927         432 :               fvm(ie)%norm_elem_coord_physgrid(1,i,j) = 1D9
         928         432 :               fvm(ie)%norm_elem_coord_physgrid(2,i,j) = 1D9
         929             :             end if
         930             :           end do
         931             :         end do
         932             :       end do
         933             :       !
         934             :       ! compute Dinv
         935             :       !
         936       12336 :       do ie=nets,nete
         937      109536 :         do j=1-nhc_phys,fv_nphys+nhc_phys
         938      982800 :           do i=1-nhc_phys,fv_nphys+nhc_phys
         939      874800 :             x1 = fvm(ie)%norm_elem_coord_physgrid(1,i,j)
         940      874800 :             x2 = fvm(ie)%norm_elem_coord_physgrid(2,i,j)
         941     6123600 :             call Dmap(D(i,j,:,:),x1,x2,elem(ie)%corners3D,cubed_sphere_map,elem(ie)%corners,elem(ie)%u2qmap,elem(ie)%facenum)
         942      874800 :             detD = D(i,j,1,1)*D(i,j,2,2) - D(i,j,1,2)*D(i,j,2,1)
         943             :
         944      874800 :             fvm(ie)%Dinv_physgrid(i,j,1,1) =  D(i,j,2,2)/detD
         945      874800 :             fvm(ie)%Dinv_physgrid(i,j,1,2) = -D(i,j,1,2)/detD
         946      874800 :             fvm(ie)%Dinv_physgrid(i,j,2,1) = -D(i,j,2,1)/detD
         947      972000 :             fvm(ie)%Dinv_physgrid(i,j,2,2) =  D(i,j,1,1)/detD
         948             :           end do
         949             :         end do
         950             :       end do
         951             :     end if
         952             :
         953        1536 :   end subroutine fvm_pg_init
         954             :
         955             :
         956             : end module fvm_mod

Generated by: LCOV version 1.14