LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - fvm_analytic_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 380 | 485 | 78.4 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 17 | 23 | 73.9 %

              Line data    Source code
           1             : !MODULE FVM_ANALYTIC_MOD--------------------------------------------CE-for FVM!
           2             : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011                                    !
           3             : ! This module contains all analytical terms for fvm                           !
           4             : !-----------------------------------------------------------------------------!
           5             : module fvm_analytic_mod
           6             :   use shr_kind_mod,   only: r8=>shr_kind_r8
           7             :   use control_mod,    only : north, south, east, west, neast, nwest, seast, swest
           8             :   use cam_abortutils, only: endrun
           9             :
          10             :   implicit none
          11             :   private
          12             :
          13             :   public :: get_high_order_weights_over_areas,  compute_reconstruct_matrix
          14             :   public :: compute_halo_vars, init_flux_orient
          15             :   public :: I_00, I_10, I_01, I_20, I_02, I_11, gauss_points
          16             :   public :: F_00, F_10, F_01, F_20, F_02, F_11
          17             :   public :: create_interpolation_points, compute_basic_coordinate_vars
          18             :
          19             : CONTAINS
          20             :
          21       21600 :   subroutine compute_basic_coordinate_vars(elem,&
          22       21600 :        nc,irecons,dalpha,dbeta,vtx_cart,center_cart,area_sphere,spherecentroid)
          23             :     use coordinate_systems_mod, only: cart2spherical
          24             :     use element_mod,            only: element_t
          25             :     use coordinate_systems_mod, only: spherical_polar_t
          26             :
          27             :     type (element_t),         intent(in ) :: elem
          28             :     integer,                  intent(in)  :: nc,irecons
          29             :
          30             :     real (kind=r8),           intent(out) :: dalpha, dbeta
          31             :     real (kind=r8),           intent(out) :: vtx_cart   (4,2,nc,nc)
          32             :     real (kind=r8),           intent(out) :: area_sphere(nc,nc)
          33             :     real (kind=r8),           intent(out) :: spherecentroid(irecons-1,nc,nc)
          34             :     type (spherical_polar_t), intent(out) :: center_cart(nc,nc) ! Spherical coordinates of fvm grid
          35             :
          36             :     integer        :: i,j
          37             :     real (kind=r8) :: centerx,centery
          38       43200 :     real (kind=r8) :: acartx(nc+1), acarty(nc+1)
          39             :
          40       21600 :     dalpha=abs(elem%corners(1)%x-elem%corners(2)%x)/nc
          41       21600 :     dbeta =abs(elem%corners(1)%y-elem%corners(4)%y)/nc
          42             :
          43      108000 :     do i=1,nc+1
          44       86400 :       acartx(i) = tan(elem%corners(1)%x+(i-1)*dalpha)
          45      108000 :       acarty(i) = tan(elem%corners(1)%y+(i-1)*dbeta)
          46             :     end do
          47             :
          48       86400 :     do j=1,nc
          49      280800 :       do i=1,nc
          50      194400 :         centerx = tan(elem%corners(1)%x+(i-0.5_r8)*dalpha)
          51      194400 :         centery = tan(elem%corners(1)%y+(j-0.5_r8)*dbeta)
          52      259200 :         center_cart(i,j) = cart2spherical(centerx,centery,elem%FaceNum)
          53             :       enddo
          54             :     enddo
          55             :
          56     2224800 :     vtx_cart = -9D9
          57       86400 :     do j=1,nc
          58      280800 :       do i=1,nc
          59      194400 :         vtx_cart(1,1,i,j) = acartx(i  )
          60      194400 :         vtx_cart(1,2,i,j) = acarty(j  )
          61             :
          62      194400 :         vtx_cart(2,1,i,j) = acartx(i+1)
          63      194400 :         vtx_cart(2,2,i,j) = acarty(j  )
          64             :
          65      194400 :         vtx_cart(3,1,i,j) = acartx(i+1)
          66      194400 :         vtx_cart(3,2,i,j) = acarty(j+1)
          67             :
          68      194400 :         vtx_cart(4,1,i,j) = acartx(i  )
          69      259200 :         vtx_cart(4,2,i,j) = acarty(j+1)
          70             :       end do
          71             :     end do
          72             :     ! compute area and centroid for the interior and halo zone of interior elements
          73       21600 :     call moment_onsphere(nc,irecons,area_sphere,vtx_cart,.true.,spherecentroid)
          74       21600 :   end subroutine compute_basic_coordinate_vars
          75             :
          76       21600 :   subroutine compute_halo_vars(faceno,cubeboundary,nc,nhc,nhe,&
          77       21600 :        jx_min,jx_max,jy_min,jy_max,flux_orient, ifct, rot_matrix)
          78             :     use control_mod, only : north, south, east, west, neast, nwest, seast, swest
          79             :
          80             :     integer,          intent(in)  :: faceno,nc,nhc,nhe,cubeboundary
          81             :
          82             :     integer,          intent(out) :: jx_min(3),jx_max(3),jy_min(3),jy_max(3)
          83             :     real (kind=r8),   intent(out) :: flux_orient(2, 1-nhc:nc+nhc,1-nhc:nc+nhc)
          84             :     integer,          intent(out) :: ifct          (1-nhc:nc+nhc,1-nhc:nc+nhc)
          85             :     integer,          intent(out) :: rot_matrix(2,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
          86             :
          87             :     integer :: i,j
          88             :     integer :: rot90_matrix(2,2)
          89             :     integer :: ishft
          90             :
          91             :
          92       21600 :     jx_min(2) = 0; jx_max(2) = -1; jy_min(2) = 0; jy_max(2) = -1
          93       21600 :     jx_min(3) = 0; jx_max(3) = -1; jy_min(3) = 0; jy_max(3) = -1
          94             :
          95       40416 :     select case (cubeboundary)
          96             :     case (0)
          97       18816 :       jx_min(1)=1-nhe; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1+nhe
          98             :     case (west)
          99         672 :       jx_min(1)=1    ; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1+nhe
         100         672 :       jx_min(2)=1-nhe; jx_max(2)=1       ; jy_min(2)=1-nhe; jy_max(2)=nc+1+nhe
         101             :     case(east)
         102         672 :       jx_min(1)=1-nhe; jx_max(1)=nc+1    ; jy_min(1)=1-nhe; jy_max(1)=nc+1+nhe
         103         672 :       jx_min(2)=nc+1 ; jx_max(2)=nc+1+nhe; jy_min(2)=1-nhe; jy_max(2)=nc+1+nhe
         104             :     case(north)
         105         672 :       jx_min(1)=1-nhe; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1
         106         672 :       jx_min(2)=1-nhe; jx_max(2)=nc+1+nhe; jy_min(2)=nc+1 ; jy_max(2)=nc+1+nhe
         107             :     case(south)
         108         672 :       jx_min(1)=1-nhe; jx_max(1)=nc+1+nhe; jy_min(1)=1    ; jy_max(1)=nc+1+nhe
         109         672 :       jx_min(2)=1-nhe; jx_max(2)=nc+1+nhe; jy_min(2)=1-nhe; jy_max(2)=1
         110             :     case(swest)
         111          24 :       jx_min(1)=1    ; jx_max(1)=nc+1+nhe; jy_min(1)=1    ; jy_max(1)=nc+1+nhe
         112          24 :       jx_min(2)=1    ; jx_max(2)=nc+1+nhe; jy_min(2)=1-nhe; jy_max(2)=1
         113          24 :       jx_min(3)=1-nhe; jx_max(3)=1       ; jy_min(3)=1    ; jy_max(3)=nc+1+nhe
         114             :     case(seast)
         115          24 :       jx_min(1)=1-nhe; jx_max(1)=nc+1    ; jy_min(1)=1    ; jy_max(1)=nc+1+nhe
         116          24 :       jx_min(2)=1-nhe; jx_max(2)=nc+1    ; jy_min(2)=1-nhe; jy_max(2)=1
         117          24 :       jx_min(3)=nc+1 ; jx_max(3)=nc+1+nhe; jy_min(3)=1    ; jy_max(3)=nc+1+nhe
         118             :     case(neast)
         119          24 :       jx_min(1)=1-nhe; jx_max(1)=nc+1    ; jy_min(1)=1-nhe; jy_max(1)=nc+1
         120          24 :       jx_min(2)=1-nhe; jx_max(2)=nc+1    ; jy_min(2)=nc+1 ; jy_max(2)=nc+1+nhe
         121          24 :       jx_min(3)=nc+1 ; jx_max(3)=nc+1+nhe; jy_min(3)=1-nhe; jy_max(3)=nc+1
         122             :     case(nwest)
         123          24 :       jx_min(1)=1    ; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1
         124          24 :       jx_min(2)=1    ; jx_max(2)=nc+1+nhe; jy_min(2)=nc+1 ; jy_max(2)=nc+1+nhe
         125          24 :       jx_min(3)=1-nhe; jx_max(3)=1       ; jy_min(3)=1-nhe; jy_max(3)=nc+1
         126             :
         127             :     case default
         128           0 :       print *, 'Fatal Error in fvm_line_integrals_mod.F90.'
         129       21600 :       call endrun('Selected case for cubeboundary does not exists!')
         130             :     end select
         131             :     !
         132             :     ! init location of flux-sides
         133             :     !
         134       21600 :     call init_flux_orient(flux_orient,ifct,nc,nhc,cubeboundary,faceno)
         135     3909600 :     rot_matrix(1,1,:,:) = 1; rot_matrix(1,2,:,:) = 0;
         136     3909600 :     rot_matrix(2,1,:,:) = 0; rot_matrix(2,2,:,:) = 1;
         137             :
         138       21600 :     if (cubeboundary>0) then
         139             :       !
         140             :       ! clockwise 90 rotation of vectors
         141             :       !
         142        2784 :       rot90_matrix(1,1) = 0; rot90_matrix(2,1) = -1;
         143        2784 :       rot90_matrix(1,2) = 1; rot90_matrix(2,2) =  0;
         144       27840 :       do j=1-nhc,nc+nhc
         145      253344 :         do i=1-nhc,nc+nhc
         146     1076544 :           do ishft=1,4-nint(flux_orient(2,i,j))
         147     6007392 :             rot_matrix(:,:,i,j) = MATMUL(rot90_matrix,rot_matrix(:,:,i,j))
         148             :           end do
         149             :         enddo
         150             :       enddo
         151             :     end if
         152       21600 :   end subroutine compute_halo_vars
         153             :
         154             :
         155             :   ! ----------------------------------------------------------------------------------!
         156             :   !SUBROUTINE MOMENT_ONSPHERE-----------------------------------------------CE-for FVM!
         157             :   ! AUTHOR: CHRISTOPH ERATH, 20.July 2011                                             !
         158             :   ! DESCRIPTION: Compute area and centroids/moments via line integrals                !
         159             :   !                                                                                   !
         160             :   ! INPUT:  x  ...  x cartesian coordinats of the arrival grid on the cube            !
         161             :   !         y  ...  y cartesian coordinats of the arrival grid on the cube            !
         162             :   !            ... cell boundaries in x and y directions                              !
         163             :   ! INPUT/OUTPUT:                                                                     !
         164             :   !         area      ... area of cells on the sphere                                 !
         165             :   !         centroid  ... x,y,x^2,y^2,xy                                              !
         166             :   !-----------------------------------------------------------------------------------!
         167       21600 :   subroutine moment_onsphere(nc,irecons,area,vtx_cart,lanalytic,spherecentroid)
         168             :     use dimensions_mod, only: ngpc
         169             :
         170             :     integer, intent(in)  :: nc,irecons
         171             :     real (kind=r8), dimension(nc,nc)                              , intent(out) :: area
         172             :     real (kind=r8), dimension(irecons-1,nc,nc), intent(out) :: spherecentroid
         173             :     real (kind=r8), dimension(4,2,nc,nc)      , intent(in)  :: vtx_cart
         174             :     logical, optional, intent(in) :: lanalytic
         175             :     integer                                              :: i,j
         176             :     !
         177             :     ! variables for call to get_high_order_weights_over_areas
         178             :     !
         179             :     integer, parameter :: num_area=1, num_seg_max=2
         180             :     REAL(KIND=r8), dimension(2,num_seg_max,num_area) :: xx, dxx
         181             :     integer             , dimension(num_area               ), parameter :: num_seg=2
         182       43200 :     REAL(KIND=r8), dimension(irecons,num_area):: weights
         183       43200 :     real (kind=r8), dimension(nc+1) :: x, y
         184             :
         185             :
         186             :     real (kind=r8), dimension(ngpc):: gsweights, gspts
         187             :     !
         188             :     ! initialize quadrature weights for get_high_order_weights_over_areas
         189             :     !
         190       21600 :     call gauss_points(ngpc,gsweights,gspts) !set gauss points/weights
         191       86400 :     gspts = 0.5_r8*(gspts+1.0_r8) !shift location so in [0:1] instead of [-1:1]
         192             :
         193       86400 :     x(1:nc) = vtx_cart(1,1,1:nc,1   )
         194       86400 :     y(1:nc) = vtx_cart(1,2,1   ,1:nc)
         195       21600 :     x(nc+1) = vtx_cart(2,1,  nc,1   )
         196       21600 :     y(nc+1) = vtx_cart(3,2,1   ,nc  )
         197             :
         198           0 :     select case (irecons)
         199             :     case(1)
         200           0 :       if (present(lanalytic)) then
         201           0 :         do j=1,nc
         202           0 :           do i=1,nc
         203           0 :             area(i,j) = (I_00(x(i+1),y(j+1)) - I_00(x(i),y(j+1)) + &
         204           0 :                  I_00(x(i),y(j)) - I_00(x(i+1),y(j)))
         205             :           end do
         206             :         end do
         207             :       else
         208           0 :         call endrun("non-analytic moments not coded for irecons=1")
         209             :       end if
         210             :
         211             :     case(3)
         212           0 :       if (present(lanalytic)) then
         213           0 :         do j=1,nc
         214           0 :           do i=1,nc
         215           0 :             area(i,j) = (I_00(x(i+1),y(j+1)) - I_00(x(i),y(j+1)) + &
         216           0 :                  I_00(x(i),y(j)) - I_00(x(i+1),y(j)))
         217             :             ! Compute centroids via line integrals
         218           0 :             spherecentroid(1,i,j) = (I_10(x(i+1),y(j+1)) - I_10(x(i),y(j+1)) + &
         219           0 :                  I_10(x(i),y(j)) - I_10(x(i+1),y(j))) / area(i,j)
         220           0 :             spherecentroid(2,i,j) = (I_01(x(i+1),y(j+1)) - I_01(x(i),y(j+1)) + &
         221           0 :                  I_01(x(i),y(j)) - I_01(x(i+1),y(j))) / area(i,j)
         222             :           end do
         223             :         end do
         224             :       else
         225           0 :         call endrun("non-analytic moments not coded for irecons=3")
         226             :       end if
         227             :
         228             :
         229             :     case(6)
         230       21600 :       if (present(lanalytic)) then
         231       86400 :         do j=1,nc
         232      280800 :           do i=1,nc
         233             :             !       area(i,j) = surfareaxy(x(i),x(i+1),y(j),y(j+1))
         234      777600 :             area(i,j) = (I_00(x(i+1),y(j+1)) - I_00(x(i),y(j+1)) + &
         235      972000 :                  I_00(x(i),y(j)) - I_00(x(i+1),y(j)))
         236             :             ! Compute centroids via line integrals
         237      194400 :             spherecentroid(1,i,j) = (I_10(x(i+1),y(j+1)) - I_10(x(i),y(j+1)) + &
         238      194400 :                  I_10(x(i),y(j)) - I_10(x(i+1),y(j))) / area(i,j)
         239      194400 :             spherecentroid(2,i,j) = (I_01(x(i+1),y(j+1)) - I_01(x(i),y(j+1)) + &
         240      194400 :                  I_01(x(i),y(j)) - I_01(x(i+1),y(j))) / area(i,j)
         241             :             ! TAN(alpha)^2 component
         242      194400 :             spherecentroid(3,i,j) = (I_20(x(i+1),y(j+1)) - I_20(x(i),y(j+1)) + &
         243      194400 :                  I_20(x(i),y(j)) - I_20(x(i+1),y(j))) / area(i,j)
         244             :             ! TAN(beta)^2 component
         245      194400 :             spherecentroid(4,i,j) = (I_02(x(i+1),y(j+1)) - I_02(x(i),y(j+1)) + &
         246      194400 :                  I_02(x(i),y(j)) - I_02(x(i+1),y(j))) / area(i,j)
         247             :             ! TAN(alpha) TAN(beta) component
         248      194400 :             spherecentroid(5,i,j) = (I_11(x(i+1),y(j+1)) - I_11(x(i),y(j+1)) + &
         249      259200 :                  I_11(x(i),y(j)) - I_11(x(i+1),y(j))) / area(i,j)
         250             :           end do
         251             :         end do
         252             :       else
         253           0 :         do j=1,nc
         254           0 :           do i=1,nc
         255             :
         256           0 :             xx (1,1,1) = x(i)       ; xx (2,1,1) = y(j+1);
         257           0 :             dxx(1,1,1) = x(i+1)-x(i); dxx(2,1,1) = 0.0_r8 ;
         258             :
         259           0 :             xx (1,2,1) = x(i+1)     ; xx (2,2,1) = y(j)  ;
         260           0 :             dxx(1,2,1) = x(i)-x(i+1); dxx(2,2,1) = 0.0_r8 ;
         261             :
         262           0 :             call get_high_order_weights_over_areas(xx,dxx,num_seg,num_seg_max,num_area,weights,ngpc,gsweights,gspts,irecons)
         263             :
         264           0 :             area(i,j)         = weights(1,1)
         265             :
         266           0 :             spherecentroid(1:5,i,j) = weights(2:6,1)/area(i,j)
         267             :           end do
         268             :         end do
         269             :       end if
         270             :     case default
         271       21600 :       call endrun('SUBROUTINE moment_on_sphere: irecons out of range')
         272             :     end select
         273       21600 :   end subroutine moment_onsphere
         274             :
         275             :
         276             :   ! ----------------------------------------------------------------------------------!
         277             :   !SUBROUTINES I_00, I_01, I_20, I_02, I11----------------------------------CE-for FVM!
         278             :   ! AUTHOR: CHRISTOPH ERATH, 17.October 2011                                          !
         279             :   ! DESCRIPTION: calculates the exact integrals                                       !
         280             :   !                                                                                   !
         281             :   ! CALLS: none                                                                       !
         282             :   ! INPUT: x    ... x coordinate of the evaluation point (Cartesian on the cube)      !
         283             :   !        y    ... y coordinate of the evaluation point (Cartesian on the cube)      !
         284             :   ! OUTPUT: I_00, I_01, I_20, I_02, I11                                               !
         285             :   !-----------------------------------------------------------------------------------!
         286     1555200 :   function I_00(x,y)
         287             :     implicit none
         288             :     real (kind=r8)                 :: I_00
         289             :     real (kind=r8), intent(in)     :: x,y
         290             :
         291     1555200 :     I_00 = ATAN(x*y/SQRT(1.0_r8+x*x+y*y))
         292     1555200 :   end function I_00
         293             :
         294     1555200 :   function I_10(x,y)
         295             :     implicit none
         296             :     real (kind=r8)                 :: I_10
         297             :     real (kind=r8), intent(in)     :: x,y
         298             :     real (kind=r8)                 :: tmp
         299             :
         300             :     !   tmp = ATAN(x)
         301             :     !   I_10 = -ASINH(y*COS(tmp))
         302     1555200 :     tmp = y*COS(ATAN(x))
         303     1555200 :     I_10 = -log(tmp+sqrt(tmp*tmp+1))
         304     1555200 :   end function I_10
         305             :
         306             :
         307     1555200 :   function I_01(x,y)
         308             :     implicit none
         309             :     real (kind=r8)                 :: I_01
         310             :     real (kind=r8), intent(in)     :: x,y
         311             :     real (kind=r8)                 :: tmp
         312             :
         313             :     !   I_01 = -ASINH(x/SQRT(1+y*y))
         314     1555200 :     tmp=x/SQRT(1+y*y)
         315     1555200 :     I_01 = -log(tmp+sqrt(tmp*tmp+1))
         316     1555200 :   end function I_01
         317             :
         318     1555200 :   function I_20(x,y)
         319             :     implicit none
         320             :     real (kind=r8)                 :: I_20
         321             :     real (kind=r8), intent(in)     :: x,y
         322             :     real (kind=r8)                 :: tmp,tmp1
         323             :
         324     1555200 :     tmp = 1.0_r8+y*y
         325     1555200 :     tmp1=x/SQRT(tmp)
         326     1555200 :     I_20 = y*log(tmp1+sqrt(tmp1*tmp1+1))+ACOS(x*y/(SQRT((1.0_r8+x*x)*tmp)))
         327     1555200 :   end function I_20
         328             :
         329     1555200 :   function I_02(x,y)
         330             :     implicit none
         331             :     real (kind=r8)                  :: I_02
         332             :     real (kind=r8), intent(in)      :: x,y
         333             :     real (kind=r8)                  :: tmp,tmp1
         334             :
         335             :     !   tmp=1.0_r8+x*x
         336             :     !   I_02 = x*ASINH(y/SQRT(tmp))+ACOS(x*y/SQRT(tmp*(1+y*y)))
         337     1555200 :     tmp=1.0_r8+x*x
         338     1555200 :     tmp1=y/SQRT(tmp)
         339             :
         340     1555200 :     I_02 = x*log(tmp1+sqrt(tmp1*tmp1+1))+ACOS(x*y/SQRT(tmp*(1+y*y)))
         341             :
         342     1555200 :   end function I_02
         343             :
         344     1555200 :   function I_11(x,y)
         345             :     implicit none
         346             :     real (kind=r8)                   :: I_11
         347             :     real (kind=r8), intent(in)       :: x,y
         348             :
         349     1555200 :     I_11 = -SQRT(1+x*x+y*y)
         350     1555200 :   end function I_11
         351             :   !END SUBROUTINES I_00, I_01, I_20, I_02, I11------------------------------CE-for FVM!
         352             :
         353             :
         354           0 :   real (kind=r8) function F_00(x_in,y_in)
         355             :     implicit none
         356             :     real (kind=r8), intent(in) :: x_in,y_in
         357             :     real (kind=r8) :: x,y
         358             :     !
         359           0 :     x = x_in
         360           0 :     y = y_in
         361           0 :     F_00 =y/((1.0_r8+x*x)*SQRT(1.0_r8+x*x+y*y))
         362           0 :   end function F_00
         363             :
         364           0 :   real (kind=r8) function F_10(x_in,y_in)
         365             :     implicit none
         366             :     real (kind=r8), intent(in) :: x_in,y_in
         367             :     real (kind=r8) :: x,y
         368             :
         369           0 :     x = x_in
         370           0 :     y = y_in
         371             :
         372           0 :     F_10 =x*y/((1.0_r8+x*x)*SQRT(1.0_r8+x*x+y*y))
         373           0 :   end function F_10
         374             :
         375           0 :   real (kind=r8) function F_01(x_in,y_in)
         376             :     implicit none
         377             :     real (kind=r8), intent(in) :: x_in,y_in
         378             :     real (kind=r8) :: x,y
         379             :
         380           0 :     x = x_in
         381           0 :     y = y_in
         382             :
         383           0 :     F_01 =-1.0_r8/(SQRT(1.0_r8+x*x+y*y))
         384           0 :   end function F_01
         385             :
         386           0 :   real (kind=r8) function F_20(x_in,y_in)
         387             :     implicit none
         388             :     real (kind=r8), intent(in) :: x_in,y_in
         389             :     real (kind=r8) :: x,y
         390             :
         391           0 :     x = x_in
         392           0 :     y = y_in
         393             :
         394           0 :     F_20 =x*x*y/((1.0_r8+x*x)*SQRT(1.0_r8+x*x+y*y))
         395           0 :   end function F_20
         396             :
         397           0 :   real (kind=r8) function F_02(x_in,y_in)
         398             :     implicit none
         399             :     real (kind=r8), intent(in) :: x_in,y_in
         400             :     real (kind=r8) :: x,y,alpha,tmp
         401             :
         402           0 :     x = x_in
         403           0 :     y = y_in
         404             :
         405             :     alpha = ATAN(x)
         406             : !     F_02 =-y/SQRT(1.0_r8+x*x+y*y)+ASINH(y*COS(alpha))
         407           0 :     tmp=y*COS(alpha)
         408           0 :     F_02 =-y/SQRT(1.0_r8+x*x+y*y)+log(tmp+sqrt(tmp*tmp+1))
         409             :
         410             :     !
         411             :     ! cos(alpha) = 1/sqrt(1+x*x)
         412             :     !
         413           0 :   end function F_02
         414             :
         415           0 :   real (kind=r8) function F_11(x_in,y_in)
         416             :     implicit none
         417             :     real (kind=r8), intent(in) :: x_in,y_in
         418             :     real (kind=r8) :: x,y
         419             :
         420           0 :     x = x_in
         421           0 :     y = y_in
         422             :
         423           0 :     F_11 =-x/(SQRT(1.0_r8+x*x+y*y))
         424           0 :   end function F_11
         425             :
         426             :
         427             :
         428             :   !
         429             :   ! matrix version of reconstruct_cubic_onface
         430             :   !
         431       21600 :   subroutine compute_reconstruct_matrix(nc,nhe,nhc,irecons,dalpha,dbeta,spherecentroid,vtx_cart,&
         432       21600 :        centroid_stretch,vertex_recons_weights,recons_metrics,recons_metrics_integral)
         433             :     implicit none
         434             :     integer              , intent(in) :: nc,nhe,irecons,nhc
         435             :     real (kind=r8), intent(in) :: dalpha,dbeta
         436             :     real (kind=r8), dimension(irecons-1,1-nhc:nc+nhc,1-nhc:nc+nhc), intent(in) :: spherecentroid
         437             :     real (kind=r8), dimension(4,2,1-nhc:nc+nhc,1-nhc:nc+nhc)      , intent(in) :: vtx_cart
         438             :
         439             :     real (kind=r8), dimension(7,1-nhe:nc+nhe,1-nhe:nc+nhe)        , intent(out):: centroid_stretch
         440             :     real (kind=r8), dimension(4,1:irecons-1,1-nhe:nc+nhe,1-nhe:nc+nhe), intent(out):: vertex_recons_weights
         441             :     real (kind=r8), dimension(3,1-nhe:nc+nhe,1-nhe:nc+nhe)        , intent(out):: recons_metrics
         442             :     real (kind=r8), dimension(3,1-nhe:nc+nhe,1-nhe:nc+nhe)        , intent(out):: recons_metrics_integral
         443             :
         444             :     !
         445             :     integer  :: i, j, count, m, n
         446             :     real (kind=r8) :: coef,tmp,cartx,carty
         447             :     !
         448             :     ! pre-compute variables for reconstruction
         449             :     !
         450             :     select case (irecons)
         451             :     case(3)
         452           0 :       do j= 1-nhe,nc+nhe
         453           0 :         do i=1-nhe,nc+nhe
         454           0 :           count = 1
         455           0 :           do n = j, j+1
         456           0 :             do m = i, i+1
         457           0 :               cartx = vtx_cart(count,1,i,j); carty = vtx_cart(count,2,i,j);
         458             :
         459           0 :               vertex_recons_weights(count,1,i,j) = cartx - spherecentroid(1,i,j)
         460           0 :               vertex_recons_weights(count,2,i,j) = carty - spherecentroid(2,i,j)
         461             :
         462           0 :               count=count+1
         463             :             end do
         464             :           enddo
         465             :         end do
         466             :       end do
         467           0 :       call endrun("recons_metrics and recons_metrics_integral not initialize")
         468             :       !
         469             :       ! for reconstruction
         470             :       !
         471           0 :       do j= 1-nhe,nc+nhe
         472           0 :         do i=1-nhe,nc+nhe
         473             :           !
         474             :           !***************
         475             :           !*    dfdx     *
         476             :           !***************
         477             :           !
         478           0 :           coef = 1.0_r8/(12.0_r8 * dalpha)                   !finite difference coefficient
         479           0 :           coef = coef /( 1.0_r8 + spherecentroid(1,i,j)**2) !stretching coefficient
         480             :
         481           0 :           centroid_stretch(1,i,j) = coef
         482             :           !
         483             :           !***************
         484             :           !*    dfdy     *
         485             :           !***************
         486             :           !
         487           0 :           coef = 1.0_r8/(12.0_r8 * dbeta)                    !finite difference coefficient
         488           0 :           coef = coef /( 1.0_r8 + spherecentroid(2,i,j)**2) !stretching coefficient
         489             :
         490           0 :           centroid_stretch(2,i,j) = coef
         491             :         end do
         492             :       end do
         493             :     case(6)
         494      108000 :       do j= 1-nhe,nc+nhe
         495      475200 :         do i=1-nhe,nc+nhe
         496     1922400 :           do count=1,4
         497     1468800 :             cartx = vtx_cart(count,1,i,j); carty = vtx_cart(count,2,i,j);
         498             :
         499     1468800 :             vertex_recons_weights(count,1,i,j) = cartx - spherecentroid(1,i,j)
         500     1468800 :             vertex_recons_weights(count,2,i,j) = carty - spherecentroid(2,i,j)
         501             :
         502     1468800 :             vertex_recons_weights(count,3,i,j) = (spherecentroid(1,i,j)**2 - &
         503             :                  spherecentroid(3,i,j))   + &
         504     1468800 :                  (cartx - spherecentroid(1,i,j))**2
         505     1468800 :             vertex_recons_weights(count,4,i,j) = (spherecentroid(2,i,j)**2 - &
         506             :                  spherecentroid(4,i,j)) + &
         507     1468800 :                  (carty - spherecentroid(2,i,j))**2
         508             :
         509     1468800 :             vertex_recons_weights(count,5,i,j) = (cartx - spherecentroid(1,i,j))*     &
         510             :                  (carty - spherecentroid(2,i,j))+     &
         511             :                  (spherecentroid(1,i,j) *             &
         512             :                  spherecentroid(2,i,j) - &
         513     1836000 :                  spherecentroid(5,i,j))
         514             :           end do
         515             :         end do
         516             :       end do
         517             :
         518      108000 :       do j= 1-nhe,nc+nhe
         519      475200 :         do i=1-nhe,nc+nhe
         520      367200 :           recons_metrics(1,i,j) = spherecentroid(1,i,j)**2 -spherecentroid(3,i,j)
         521      367200 :           recons_metrics(2,i,j) = spherecentroid(2,i,j)**2 -spherecentroid(4,i,j)
         522             :           recons_metrics(3,i,j) = spherecentroid(1,i,j)*spherecentroid(2,i,j)-&
         523      367200 :                spherecentroid(5,i,j)
         524             :
         525             :           recons_metrics_integral(1,i,j) = &
         526      367200 :                2.0_r8*spherecentroid(1,i,j)**2 -spherecentroid(3,i,j)
         527             :           recons_metrics_integral(2,i,j) = &
         528      367200 :                2.0_r8*spherecentroid(2,i,j)**2 -spherecentroid(4,i,j)
         529             :           recons_metrics_integral(3,i,j) = &
         530             :                2.0_r8*spherecentroid(1,i,j)*spherecentroid(2,i,j)-&
         531      453600 :                spherecentroid(5,i,j)
         532             :         end do
         533             :       end do
         534             :
         535             :
         536             :
         537             :       !
         538             :       ! pre-compute variables for reconstruction
         539             :       !
         540      108000 :       do j= 1-nhe,nc+nhe
         541      475200 :         do i=1-nhe,nc+nhe
         542             :           !
         543             :           !***************
         544             :           !*    dfdx     *
         545             :           !***************
         546             :           !
         547      367200 :           coef = 1.0_r8/(12.0_r8 * dalpha)                   !finite difference coefficient
         548      367200 :           coef = coef /( 1.0_r8 + spherecentroid(1,i,j)**2) !stretching coefficient
         549             :
         550      367200 :           centroid_stretch(1,i,j) = coef
         551             :           !
         552             :           !***************
         553             :           !*    dfdy     *
         554             :           !***************
         555             :           !
         556      367200 :           coef = 1.0_r8/(12.0_r8 * dbeta)                    !finite difference coefficient
         557      367200 :           coef = coef /( 1.0_r8 + spherecentroid(2,i,j)**2) !stretching coefficient
         558             :
         559      367200 :           centroid_stretch(2,i,j) = coef
         560             :
         561             :           !*****************
         562             :           !*    d2fdx2     *
         563             :           !*****************
         564             :           !
         565      367200 :           coef = 1.0_r8 / (12.0_r8 * dalpha**2)                  !finite difference coefficient
         566             :           !
         567             :           ! stretching coefficient part 2
         568             :           !      recons(3,i,j) = (a * recons(1,i,j)+ recons(3,i,j))*b
         569             :           !
         570      367200 :           tmp = 0.5_r8/((1.0_r8 + spherecentroid(1,i,j)**2)**2)
         571             :
         572      367200 :           centroid_stretch(3,i,j) = coef*tmp
         573      367200 :           centroid_stretch(6,i,j) = -spherecentroid(1,i,j)/(1.0_r8 + spherecentroid(1,i,j)**2)
         574             :
         575             :           !
         576             :           !*****************
         577             :           !*    d2fdy2     *
         578             :           !*****************
         579             :           !
         580             :           !
         581      367200 :           coef = 1.0_r8 / (12.0_r8 * dbeta**2)                     !finite difference coefficient
         582             :           !
         583             :           ! stretching coefficient part 2
         584             :           !
         585             :           !      recons(4,i,j) = (a * recons(1,i,j)+ recons(4,i,j))*b
         586             :           !
         587      367200 :           tmp =0.5_r8/((1.0_r8 + spherecentroid(2,i,j)**2)**2)
         588             :
         589      367200 :           centroid_stretch(4,i,j) = coef*tmp
         590      367200 :           centroid_stretch(7,i,j) = -spherecentroid(2,i,j)/(1.0_r8 + spherecentroid(2,i,j)**2)
         591             :           !
         592             :           !*****************
         593             :           !*    d2fdxdy    *
         594             :           !*****************
         595             :           !
         596             :           !
         597      367200 :           coef = 1.0_r8 / (4.0_r8 * dalpha * dbeta)            !finite difference coefficient
         598             :           coef = coef  / ((1.0_r8 + spherecentroid(1,i,j)**2) * &
         599      367200 :                (1.0_r8 + spherecentroid(2,i,j)**2))    !stretching coefficient
         600             :
         601      453600 :           centroid_stretch(5,i,j) = coef
         602             :         enddo
         603             :       enddo
         604             :     case default
         605       21600 :       call endrun('SUBROUTINE compute_reconstruct_matrix: irecons out of range')
         606             :     end select
         607       21600 :   end subroutine compute_reconstruct_matrix
         608             :
         609             :
         610 11594793592 :   subroutine get_high_order_weights_over_areas(x,dx,num_seg,num_seg_max,num_area,weights,ngpc,gsweights, gspts,irecons)
         611             :     implicit none
         612             :     integer                                                 , intent(in)    :: num_area, num_seg_max, irecons
         613             :     REAL(KIND=r8), dimension(2,num_seg_max,num_area ), intent(inout) :: x, dx
         614             :     integer                                                 , intent(in)    :: ngpc
         615             :     integer             , dimension(num_area               ), intent(in)    :: num_seg
         616             :     REAL(KIND=r8), dimension(irecons,num_area), intent(out)   :: weights
         617             :
         618 23189587184 :     real (kind=r8), dimension(ngpc,num_seg_max               ) :: xq,yq        !quadrature points along line segments
         619 23189587184 :     real (kind=r8), dimension(ngpc,num_seg_max,irecons) :: F            !potentials
         620 23189587184 :     real (kind=r8), dimension(                 irecons) :: weights_area
         621 23189587184 :     real (kind=r8), dimension(ngpc,num_seg_max) :: xq2, yrh, rho, tmp !intermediate variables for optimization
         622 11594793592 :     REAL(KIND=r8) , dimension(ngpc,num_seg_max) :: xq2ir, xq2i, rhoi  !intermediate variables for optimization
         623             :
         624             :     integer :: iseg,iarea,i,j,k
         625             :
         626             :     real (kind=r8), dimension(ngpc) :: gsweights, gspts
         627             :
         628 >41741*10^7 :     weights(1:irecons,1:num_area) = 0.0_r8 !may not be necessary dbgxxx
         629 69568761552 :     do iarea=1,num_area
         630 >14061*10^7 :       do iseg=1,num_seg(iarea)
         631 >33058*10^7 :         xq(:,iseg) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(:)
         632 >38855*10^7 :         yq(:,iseg) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(:)
         633             :       end do
         634             :       !
         635             :       ! potentials (equation's 23-28 in CSLAM paper; Lauritzen et al., 2010):
         636             :       !
         637             :       ! (Rory Kelly optimization)
         638             :       !
         639 >14061*10^7 :       do j=1,num_seg(iarea)
         640             : !DIR$ SIMD
         641 >33058*10^7 :         do i=1,ngpc
         642 >24793*10^7 :           xq2(i,j)   =  xq(i,j)*xq(i,j)
         643 >24793*10^7 :           xq2i(i,j)  =  1.0_r8/(1.0_r8+xq2(i,j))
         644 >24793*10^7 :           xq2ir(i,j) =  SQRT(xq2i(i,j))
         645 >24793*10^7 :           rho(i,j)   =  SQRT(1.0_r8+xq2(i,j)+yq(i,j)*yq(i,j))
         646 >24793*10^7 :           rhoi(i,j)  =  1.0_r8/rho(i,j)
         647 >24793*10^7 :           yrh(i,j)   =  yq(i,j)*rhoi(i,j)
         648 >24793*10^7 :           tmp(i,j)   =  yq(i,j)*xq2ir(i,j)
         649 >24793*10^7 :           F(i,j,1)   =  yrh(i,j)*xq2i(i,j)                 !F_00 !F_00
         650 >24793*10^7 :           F(i,j,2)   =  xq(i,j)*yrh(i,j)*xq2i(i,j)         !F_10 !F_10
         651 >24793*10^7 :           F(i,j,3)   = -1.0_r8*rhoi(i,j)                    !F_01 !F_01
         652 >24793*10^7 :           F(i,j,4)   =  xq2(i,j)*yrh(i,j)*xq2i(i,j)        !F_20 !F_20
         653 >33058*10^7 :           F(i,j,6)   = -xq(i,j)*rhoi(i,j)                  !F_11 !F_11
         654             :         enddo
         655             :         !
         656             :         ! take F(i,j,5) out of loop above since it prevents vectorization
         657             :         !
         658 >38855*10^7 :         do i=1,ngpc
         659 >33058*10^7 :           F(i,j,5)   = -yq(i,j)*rhoi(i,j)+log(tmp(i,j)+rho(i,j)*xq2ir(i,j))  !F_02 !F_02
         660             :         end do
         661             :       enddo
         662 >40581*10^7 :       weights_area = 0.0_r8
         663 >40581*10^7 :       do k=1,irecons
         664 >90169*10^7 :         do iseg=1,num_seg(iarea)
         665 >23313*10^8 :           weights_area(k) = weights_area(k) + sum(gsweights(:)*F(:,iseg,k))*0.5_r8*dx(1,iseg,iarea)
         666             :         end do
         667             :       end do
         668 >41741*10^7 :       weights(1:irecons,iarea) = weights_area(1:irecons)
         669             :     end do
         670 11594793592 :   end subroutine get_high_order_weights_over_areas
         671             :
         672             :
         673             :   !********************************************************************************
         674             :   !
         675             :   ! Gauss-Legendre quadrature
         676             :   !
         677             :   ! Tabulated values
         678             :   !
         679             :   !********************************************************************************
         680      760416 :   subroutine gauss_points(n,weights,points)
         681             :     implicit none
         682             :     integer, intent(in ) :: n
         683             :     real (kind=r8), dimension(:), intent(out) :: weights, points !dimension(n)
         684             :
         685      760416 :     select case (n)
         686             :       !    CASE(1)
         687             :       !       abscissae(1) = 0.0_r8
         688             :       !       weights(1)   = 2.0_r8
         689             :     case(2)
         690           0 :       points(1)    = -sqrt(1.0_r8/3.0_r8)
         691           0 :       points(2)    =  sqrt(1.0_r8/3.0_r8)
         692           0 :       weights(1)   =  1.0_r8
         693           0 :       weights(2)   =  1.0_r8
         694             :     case(3)
         695      760416 :       points(1)    = -0.774596669241483377035853079956_r8
         696      760416 :       points(2)    =  0.0_r8
         697      760416 :       points(3)    =  0.774596669241483377035853079956_r8
         698      760416 :       weights(1)   =  0.555555555555555555555555555556_r8
         699      760416 :       weights(2)   =  0.888888888888888888888888888889_r8
         700      760416 :       weights(3)   =  0.555555555555555555555555555556_r8
         701             :     case(4)
         702           0 :       points(1)    = -0.861136311594052575223946488893_r8
         703           0 :       points(2)    = -0.339981043584856264802665659103_r8
         704           0 :       points(3)    =  0.339981043584856264802665659103_r8
         705           0 :       points(4)    =  0.861136311594052575223946488893_r8
         706           0 :       weights(1)   =  0.347854845137453857373063949222_r8
         707           0 :       weights(2)   =  0.652145154862546142626936050778_r8
         708           0 :       weights(3)   =  0.652145154862546142626936050778_r8
         709           0 :       weights(4)   =  0.347854845137453857373063949222_r8
         710             :     case(5)
         711           0 :       points(1)    = -(1.0_r8/3.0_r8)*sqrt(5.0_r8+2.0_r8*sqrt(10.0_r8/7.0_r8))
         712           0 :       points(2)    = -(1.0_r8/3.0_r8)*sqrt(5.0_r8-2.0_r8*sqrt(10.0_r8/7.0_r8))
         713           0 :       points(3)    =  0.0_r8
         714           0 :       points(4)    =  (1.0_r8/3.0_r8)*sqrt(5.0_r8-2.0_r8*sqrt(10.0_r8/7.0_r8))
         715           0 :       points(5)    =  (1.0_r8/3.0_r8)*sqrt(5.0_r8+2.0_r8*sqrt(10.0_r8/7.0_r8))
         716           0 :       weights(1)   = (322.0_r8-13.0_r8*sqrt(70.0_r8))/900.0_r8
         717           0 :       weights(2)   = (322.0_r8+13.0_r8*sqrt(70.0_r8))/900.0_r8
         718           0 :       weights(3)   = 128.0_r8/225.0_r8
         719           0 :       weights(4)   = (322.0_r8+13.0_r8*sqrt(70.0_r8))/900.0_r8
         720           0 :       weights(5)   = (322.0_r8-13.0_r8*sqrt(70.0_r8))/900.0_r8
         721             :     case default
         722      760416 :       call endrun('SUBROUTINE gauss_points: n out of range in (0<n<5)')
         723             :     end select
         724             :
         725      760416 :   end subroutine gauss_points
         726             :
         727             :
         728             :
         729       21600 : subroutine init_flux_orient(flux_orient,ifct,nc,nhc,cubeboundary,faceno)
         730             :   implicit none
         731             :   integer              , intent(in)  :: cubeboundary,faceno,nc,nhc
         732             :   real (kind=r8), intent(out) :: flux_orient(2  ,1-nhc:nc+nhc,1-nhc:nc+nhc)
         733             :   integer              , intent(out) :: ifct           (1-nhc:nc+nhc,1-nhc:nc+nhc)
         734             :
         735             :   integer :: ib
         736             :
         737     5464800 :   flux_orient      = 99.9D9 !for debugging
         738             :   !
         739             :   ! halo of flux_orient will be filled through boundary exchange
         740             :   !
         741      280800 :   flux_orient (1,1:nc,1:nc) = dble(faceno)
         742     1965600 :   flux_orient (2,:,:) = 0.0_r8
         743     1965600 :   ifct(:,:) = 1
         744       21600 :   if (cubeboundary>0) then
         745             :
         746             :      !
         747             :      ! cshift (permute) value needed to be applied to vertex number so that they match orientation
         748             :      ! of the interior of the panel
         749             :      !
         750             :      !
         751        2784 :      ib = cubeboundary
         752        2784 :      if (faceno==2) then
         753        4184 :         if (ib==north.or.ib==nwest.or.ib==neast) flux_orient(2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 1
         754        4184 :         if (ib==south.or.ib==swest.or.ib==seast) flux_orient(2,1-nhc:nc+nhc,1-nhc:0     ) = 3
         755             :      end if
         756        2784 :      if (faceno==3) then
         757        4184 :         if (ib==north.or.ib==nwest.or.ib==neast) flux_orient (2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 2
         758        4184 :         if (ib==south.or.ib==swest.or.ib==seast) flux_orient (2,1-nhc:nc+nhc,1-nhc:0     ) = 2
         759             :      end if
         760        2784 :      if (faceno==4) then
         761        4184 :         if (ib==north.or.ib==nwest.or.ib==neast) flux_orient (2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 3
         762        4184 :         if (ib==south.or.ib==swest.or.ib==seast) flux_orient (2,1-nhc:nc+nhc,1-nhc:0     ) = 1
         763             :      end if
         764        2784 :      if (faceno==5) then
         765        4184 :         if (ib==south.or.ib==swest.or.ib==seast) flux_orient (2,1-nhc:nc+nhc,1-nhc:0     ) = 2
         766        4904 :         if (ib== west.or.ib==swest.or.ib==nwest) flux_orient (2,1-nhc:0     ,1-nhc:nc+nhc) = 3
         767        4904 :         if (ib== east.or.ib==seast.or.ib==neast) flux_orient (2, nc+1:nc+nhc,1-nhc:nc+nhc) = 1
         768             :      end if
         769             :
         770        2784 :      if (faceno==6) then
         771        4184 :         if (ib==north.or.ib==nwest.or.ib==neast ) flux_orient (2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 2
         772        4904 :         if (ib==west .or.ib==swest.or.ib==nwest ) flux_orient (2,1-nhc:0    ,1-nhc:nc+nhc)  = 1
         773        4904 :         if (ib==east .or.ib==seast.or.ib==neast ) flux_orient (2,nc+1:nc+nhc,1-nhc:nc+nhc)  = 3
         774             :      end if
         775             :      !
         776             :      ! non-existent cells in physical space
         777             :      !
         778        2784 :      if (cubeboundary==nwest) then
         779         336 :         flux_orient(2,1-nhc:0     ,nc+1 :nc+nhc) = 0
         780         312 :         ifct       (  1-nhc:0     ,nc+1 :nc+nhc) = 0
         781        2760 :      else if (cubeboundary==swest) then
         782         336 :         flux_orient (2,1-nhc:0     ,1-nhc:0     ) = 0
         783         312 :         ifct        (  1-nhc:0     ,1-nhc:0     ) = 0
         784        2736 :      else if (cubeboundary==neast) then
         785         336 :         flux_orient (2,nc+1 :nc+nhc,nc+1 :nc+nhc) = 0
         786         312 :         ifct        (  nc+1 :nc+nhc,nc+1 :nc+nhc) = 0
         787        2712 :      else if (cubeboundary==seast) then
         788         336 :         flux_orient (2,nc+1 :nc+nhc,1-nhc:0     ) = 0
         789         312 :         ifct        (  nc+1 :nc+nhc,1-nhc:0     ) = 0
         790             :      end if
         791             :    end if
         792             :
         793       21600 :  end subroutine init_flux_orient
         794             :
         795             : !
         796             : !
         797             : !
         798             :
         799             : ! ----------------------------------------------------------------------------------!
         800             : !SUBROUTINE CREATE_INTERPOLATIION_POINTS----------------------------------CE-for FVM!
         801             : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011                                          !
         802             : ! DESCRIPTION: for elements, which share a cube edge, we have to do some            !
         803             : !        interpolation on different cubic faces, also in the halo region:           !
         804             : !        because we also need the reconstruction coefficients in the halo zone,     !
         805             : !        which is basically calculated twice, on the original cell of an element    !
         806             : !        on face A and on a cell in the halo region of an element of face B         !
         807             : !        The crux is, that the interpolation has to be the same to ensure           !
         808             : !        conservation of the scheme                                                 !
         809             : !        SYMMETRY of the CUBE is used for calucaltion the interpolation_point       !
         810             : !                                                                                   !
         811             : ! CALLS: interpolation_point                                                        !
         812             : ! INPUT/OUTPUT:                                                                     !
         813             : !        elem     ...  element structure from HOMME                                 !
         814             : !        fvm   ...  structure                                                       !
         815             : !-----------------------------------------------------------------------------------!
         816       21600 :   subroutine create_interpolation_points(elem,&
         817             :        nc,nhc,nhr,ns,nh,cubeboundary,&
         818       21600 :        dalpha,dbeta,ibase,halo_interp_weight)
         819             :     use element_mod           , only: element_t
         820             :     use coordinate_systems_mod, only: cartesian2D_t
         821             :     use control_mod           , only: north, south, east, west, neast, nwest, seast, swest
         822             :     use cube_mod              , only: cube_xstart, cube_xend, cube_ystart, cube_yend
         823             :
         824             :     implicit none
         825             :     type (element_t), intent(in)     :: elem
         826             :
         827             :     integer              , intent(in)   :: nc,nhc,nhr,ns,nh,cubeboundary
         828             :     integer              , intent(out) :: ibase(1-nh:nc+nh,1:nhr,2)
         829             :     real (kind=r8), intent(out) :: halo_interp_weight(1:ns,1-nh:nc+nh,1:nhr,2)
         830             :     !
         831             :     ! pre-compute weight/index matrices
         832             :     !
         833             :     integer :: imin,imax,jmin,jmax,iinterp
         834             :     real (kind=r8), intent(in)   :: dalpha,dbeta
         835             :
         836       43200 :     real    (kind=r8), dimension(1-nhc:nc+nhc) :: gnomxstart, gnomxend, gnomystart, gnomyend
         837             :     integer                                         :: i, halo, ida, ide, iref1
         838             :     type (cartesian2D_t)                            :: tmpgnom
         839       43200 :     real (kind=r8)  :: interp(1-nh:nc+nh,1:nhr,2)
         840             :     integer                              ::ibaseref
         841       21600 :     integer                             :: ibase_tmp(1-nh:nc+nh,1:nhr,2)
         842             :
         843      756000 :   ibase = 99999 !dbg
         844     2570400 :   halo_interp_weight(:,:,:,:) = 9.99E9_r8 !dbg
         845             :
         846             :   ! element is not on a corner, but shares a cube edge (call of subroutine)
         847       21600 :   if(cubeboundary <= 4) then
         848       21504 :     gnomxstart(1-nhc)=elem%corners(1)%x-(nhc-0.5_r8)*dalpha
         849       21504 :     gnomystart(1-nhc)=elem%corners(1)%y-(nhc-0.5_r8)*dbeta
         850      193536 :     do i=2-nhc,nc+nhc
         851      172032 :       gnomxstart(i)=gnomxstart(i-1)+dalpha
         852      193536 :       gnomystart(i)=gnomystart(i-1)+dbeta
         853             :     end do
         854       21504 :     ida=1-nhc  !lower bound
         855       21504 :     ide=nc+nhc !upper bound
         856             :     select case (cubeboundary)
         857             :       !INTERIOR element
         858             :       case(0)
         859             :         ! nothing to do!
         860             :       !CASE WEST
         861             :       case(west)
         862        2016 :         do halo=1,nhr
         863             : !          iref1=ida
         864        1344 :           tmpgnom%x=cube_xstart-(halo-0.5_r8)*dalpha
         865       10080 :           do i=halo-nh,nc+nh-(halo-1) !see fvm_reconstruction to understand these boundaries
         866        8064 :             iref1=ida
         867        8064 :             tmpgnom%y=gnomystart(i)
         868        8064 :             call interpolation_point(nc,ns,tmpgnom,gnomystart,1,4,1,interp(i,halo,1),&
         869       17472 :                                      ida,ide,iref1,ibase_tmp(i,halo,1))
         870             :           end do
         871             :         end do
         872             :
         873             :       !CASE EAST
         874             :       case(east)
         875             :         ! east zone
         876        2016 :         do halo=1,nhr
         877        1344 :           iref1=ida
         878        1344 :           tmpgnom%x=cube_xend+(halo-0.5_r8)*dalpha
         879       10080 :           do i=halo-nh,nc+nh-(halo-1)
         880        8064 :             tmpgnom%y=gnomystart(i)
         881        8064 :             call interpolation_point(nc,ns,tmpgnom,gnomystart,1,2,1,interp(i,halo,1),&
         882       17472 :                                      ida,ide,iref1,ibase_tmp(i,halo,1))
         883             :           end do
         884             :         end do
         885             :
         886             :       !CASE NORTH
         887             :       case(north)
         888             :         ! north zone
         889        2016 :         do halo=1,nhr
         890        1344 :           tmpgnom%y=cube_yend+(halo-0.5_r8)*dbeta
         891        1344 :           iref1=ida
         892       10080 :           do i=halo-nh,nc+nh-(halo-1)
         893        8064 :             tmpgnom%x=gnomxstart(i)
         894             :             !
         895             :             ! dbg - change to interp(i,halo,1) instead of interp(i,halo,2)
         896             :             !       so that I can get rid of iinterp = 1 in fvm_reconstruction_mod
         897             :             !
         898        8064 :             call interpolation_point(nc,ns,tmpgnom,gnomxstart,1,6,0,interp(i,halo,2),&
         899       17472 :                                      ida,ide,iref1,ibase_tmp(i,halo,2))
         900             :           end do
         901             :         end do
         902             :       !CASE SOUTH
         903             :       case(south)
         904             :        !south zone
         905        2016 :        do halo=1,nhr
         906        1344 :           iref1=ida
         907        1344 :           tmpgnom%y=cube_ystart-(halo-0.5_r8)*dbeta
         908       10080 :           do i=halo-nh,nc+nh-(halo-1)
         909        8064 :             tmpgnom%x=gnomxstart(i)
         910        8064 :             call interpolation_point(nc,ns,tmpgnom,gnomxstart,1,5,0,interp(i,halo,2),&
         911       17472 :                                      ida,ide,iref1,ibase_tmp(i,halo,2))
         912             :           end do
         913             :         end do
         914             :
         915             :         !
         916             :         !THIS CASE SHOULD NOT HAPPEN!
         917             :      case default
         918           0 :            print *,'Fatal Error in first select statement:'
         919       21504 :            call endrun('fvm_reconstruction_mod.F90 subroutine fillhalo_cubic!' )
         920             :     end select
         921             :   !CORNER TREATMENT
         922             :   else
         923          96 :     gnomxstart(1-nhc)=cube_xstart-(nhc-0.5_r8)*dalpha
         924          96 :     gnomxend(nc+nhc)=cube_xend+(nhc-0.5_r8)*dalpha
         925          96 :     gnomystart(1-nhc)=cube_ystart-(nhc-0.5_r8)*dbeta
         926          96 :     gnomyend(nc+nhc)=cube_yend+(nhc-0.5_r8)*dbeta
         927         864 :     do i=2-nhc,nc+nhc
         928         768 :       gnomxstart(i)=gnomxstart(i-1)+dalpha
         929         768 :       gnomxend(nc+1-i)=gnomxend(nc+2-i)-dalpha
         930         768 :       gnomystart(i)=gnomystart(i-1)+dbeta
         931         864 :       gnomyend(nc+1-i)=gnomyend(nc+2-i)-dbeta
         932             :     end do
         933             :
         934             :     select case (cubeboundary)
         935             :       !CASE SOUTH WEST
         936             :       case(swest)
         937             :         ! west zone
         938          72 :         do halo=1,nhr
         939          48 :           tmpgnom%x=cube_xstart-(halo-0.5_r8)*dalpha
         940          48 :           ida=1
         941          48 :           ide=nc+nc
         942          48 :           iref1=ida
         943         336 :           do i=0,nc+nh-(halo-1)
         944         264 :             tmpgnom%y=gnomystart(i)
         945         264 :             call interpolation_point(nc,ns,tmpgnom,gnomystart,1,4,1,interp(i,halo,1),&
         946         576 :                                      ida,ide,iref1,ibase_tmp(i,halo,1))
         947             :           end do
         948             :        end do
         949             :       !CASE SOUTH EAST
         950             :       case(seast)
         951             :         ! east zone
         952          72 :         do halo=1,nhr
         953          48 :           tmpgnom%x=cube_xend+(halo-0.5_r8)*dalpha
         954          48 :           ida=1
         955          48 :           ide=nc+nc
         956          48 :           iref1=ida
         957         336 :           do i=0,nc+nh-(halo-1)
         958         264 :             tmpgnom%y=gnomystart(i)
         959         264 :             call interpolation_point(nc,ns,tmpgnom,gnomystart,1,2,1, interp(i,halo,1),&
         960         576 :                                      ida,ide,iref1,ibase_tmp(i,halo,1))
         961             :           end do
         962             :        end do
         963             :       !CASE NORTH EAST
         964             :       case(neast)
         965             :         ! east zone
         966          72 :         do halo=1,nhr
         967          48 :           tmpgnom%x=cube_xend+(halo-0.5_r8)*dalpha
         968          48 :           ida=1-nc
         969          48 :           ide=nc
         970          48 :           iref1=ida
         971         336 :           do i=halo-nh,nc+1
         972         264 :             tmpgnom%y=gnomyend(i)
         973         264 :             call interpolation_point(nc,ns,tmpgnom,gnomyend,1,2,1, interp(i,halo,1),&
         974         576 :                                      ida,ide,iref1,ibase_tmp(i,halo,1))
         975             :           end do
         976             :        end do
         977             :       !CASE NORTH WEST
         978             :       case(nwest)
         979             :         ! west zone
         980          72 :         do halo=1,2
         981          48 :           tmpgnom%x=cube_xstart-(halo-0.5_r8)*dalpha
         982          48 :           ida=1-nc
         983          48 :           ide=nc
         984          48 :           iref1=ida
         985         336 :           do i=halo-nh,nc+1
         986         264 :             tmpgnom%y=gnomyend(i)
         987         264 :             call interpolation_point(nc,ns,tmpgnom,gnomyend,1,4,1, interp(i,halo,1),&
         988         576 :                                      ida,ide,iref1,ibase_tmp(i,halo,1))
         989             :           end do
         990             :        end do
         991             :         !THIS CASE SHOULD NOT HAPPEN!
         992             :           case default
         993           0 :             print *,'Fatal Error in second select statement:'
         994          96 :             call endrun('fvm_reconstruction_mod.F90 subroutine create_interpolationpoint!')
         995             :          end select
         996             :   endif
         997             :
         998             :   !**************************
         999             :   !
        1000             :   ! compute haloe weights and indices
        1001             :   !
        1002       21600 :   if (cubeboundary>0) then
        1003        2784 :      if (cubeboundary<5) then
        1004             :         !
        1005             :         ! element is located at a panel side but is not a corner element
        1006             :         ! (west,east,south,north) = (1,2,3,4)
        1007             :         !
        1008        2688 :         if (cubeboundary==west .or.cubeboundary==east ) then
        1009             :            iinterp = 1
        1010             :         end if
        1011        2688 :         if (cubeboundary==north.or.cubeboundary==south) iinterp = 2
        1012        8064 :         do halo=1,nhr
        1013       40320 :            do i=halo-nh,nc+nh-(halo-1)
        1014       32256 :               ibaseref=ibase_tmp(i,halo,iinterp)
        1015       32256 :               ibase(i,halo,1) = ibaseref
        1016             :               call get_equispace_weights(dbeta, interp(i,halo,iinterp),&
        1017       37632 :                    halo_interp_weight(:,i,halo,1),ns)
        1018             :            end do
        1019             :         end do
        1020             :      else
        1021             :         !
        1022             :         ! element is located at a cube corner
        1023             :         ! (swest,seast,nwest,neast)=(5,6,7,8)
        1024             :         !
        1025         288 :         do halo=1,nhr
        1026         192 :            if (cubeboundary==swest .or.cubeboundary==seast) then
        1027          96 :               imin = 0      ; imax = nc+nh-(halo-1);
        1028          96 :               jmin = halo-nh; jmax = nc+1;
        1029             :            else
        1030          96 :               jmin = 0      ; jmax = nc+nh-(halo-1);
        1031          96 :               imin = halo-nh; imax = nc+1;
        1032             :            end if
        1033        1248 :            do i=imin,imax
        1034        1056 :               ibaseref=ibase_tmp(i,halo,1)
        1035        1056 :               ibase(i,halo,1) = ibaseref
        1036        1248 :               call get_equispace_weights(dbeta, interp(i,halo,1),halo_interp_weight(:,i,halo,1),ns)
        1037             :            end do
        1038             :            !
        1039             :            ! reverse weights/indices for fotherpanel (see details on reconstruct_matrix)
        1040             :            !
        1041        4416 :            halo_interp_weight(1:ns,jmin:jmax,halo,2) = halo_interp_weight(ns:1:-1,imax:imin:-1,halo,1)
        1042        1344 :            ibase       (jmin:jmax,halo     ,2) = nc+1-(ns-1)-ibase(imax:imin:-1,halo        ,1)
        1043             :         end do
        1044             :      end if
        1045             :
        1046             :   end if
        1047             :
        1048             :
        1049       21600 : end subroutine create_interpolation_points
        1050             :
        1051             :
        1052             :
        1053             :
        1054             : !END SUBROUTINE CREATE_INTERPOLATION_POINTS-------------------------------CE-for FVM!
        1055             :
        1056             :
        1057             :
        1058             : ! ----------------------------------------------------------------------------------!
        1059             : !SUBROUTINE INTERPOLATION_POINT-------------------------------------------CE-for FVM!
        1060             : ! AUTHOR: CHRISTOPH ERATH, 14.November 2011                                         !
        1061             : ! DESCRIPTION: calculates the interpolation point on from face 1 in face 2 in       !
        1062             : !        alpha/beta coordinates, only 1D                                            !
        1063             : !                                                                                   !
        1064             : ! CALLS: cubedsphere2cart, cart2cubedsphere                                         !
        1065             : ! INPUT: gnom... 1D coordinates                                                     !
        1066             : !        gnom1d... 1d coordinates                                                   !
        1067             : !        face1... orginal face                                                      !
        1068             : !        face2... target face (where the interpolation has to be done)              !
        1069             : !        xy ... 0 for alpha coordinate, any other for beta                          !
        1070             : !        except.which type, interior, left edge (-1), right edge (1)                !
        1071             : !        point... interpolation point                                               !
        1072             : !        ida  ... begin of interpval                                                !
        1073             : !        ide  ... end of interpval                                                  !
        1074             :
        1075             :
        1076             : ! INPUT/OUTPUT/RETURN:                                                              !
        1077             : !        iref ... where we start the search, is also an OUTPUT, so we know for the  !
        1078             : !                 next point where to start                                         !
        1079             : !-----------------------------------------------------------------------------------!
        1080             :                                       !
        1081             : ! DESCRIPTION: searchs where the interpolation point has to be (iref), two values   !
        1082             : !        of interpval on the left and on the right, except if we are out of range   !
        1083             : !        which is indicated through ia and ie, respectively                         !
        1084             : !        It is a 1D interpolation, use alpha/beta coordinates!!!                    !
        1085             : !                                                                                   !
        1086             : ! CALLS: cubic_equispace_interp                                                     !
        1087             : ! INPUT: iref ... where we start the search, is also an OUTPUT, so we know for the  !
        1088             : !                 next point where to start                                         !
        1089             : !        ibaseref ... startindex of the four tracer value for the reconstruction    !
        1090             : !        point    ... provides the difference of the interpolation point to use it  !
        1091             : !                     directly in CUBIC_EQUISPACE_INTERP                            !
        1092             : !-----------------------------------------------------------------------------------!
        1093       33312 : function get_gno_point(gnom,face1,face2,xy) result(point)
        1094             :   use coordinate_systems_mod, only : cubedsphere2cart, cart2cubedsphere, &
        1095             :                                      cartesian2D_t,cartesian3D_t
        1096             :   implicit none
        1097             :   type (cartesian2D_t), intent(in)                     :: gnom
        1098             :   integer, intent(in)                                  :: face1, face2, xy
        1099             :   real (kind=r8)                                :: point
        1100             :
        1101             :   type(cartesian3D_t)                 :: tmpcart3d
        1102             :   type (cartesian2D_t)                :: tmpgnom
        1103             :
        1104       33312 :   tmpcart3d=cubedsphere2cart(gnom,face1)
        1105       33312 :   tmpgnom=cart2cubedsphere(tmpcart3d,face2)
        1106       33312 :   if(xy==0) then
        1107       16128 :     point=tmpgnom%x
        1108             :   else
        1109       17184 :     point=tmpgnom%y
        1110             :   end if
        1111       33312 : end function get_gno_point
        1112             :
        1113       33312 : subroutine interpolation_point(nc,ns,gnom,gnom1d,face1,face2,xy,point,ida,ide,iref,ibaseref)
        1114             :   use coordinate_systems_mod, only : cartesian2D_t
        1115             :   implicit none
        1116             :   integer                                , intent(in) :: nc,ns
        1117             :   type (cartesian2D_t), intent(in)                    :: gnom
        1118             :   real (kind=r8), dimension(1-nc:), intent(in) :: gnom1d  !dimension(1-nhc:nc+nhc)
        1119             :   integer, intent(in)                                 :: face1, face2, xy
        1120             :   integer,intent(in)                                  :: ida, ide
        1121             :   integer,intent(inout)                               :: iref,ibaseref
        1122             :   real (kind=r8), intent(inout)                :: point
        1123             :
        1124             : !  type(cartesian3D_t)                 :: tmpcart3d
        1125             : !  type (cartesian2D_t)                :: tmpgnom
        1126             :
        1127       33312 :   point = get_gno_point(gnom,face1,face2,xy)
        1128             :
        1129             : !  tmpcart3d=cubedsphere2cart(gnom,face1)
        1130             : !  tmpgnom=cart2cubedsphere(tmpcart3d,face2)
        1131             : !  if(xy==0) then
        1132             : !    point=tmpgnom%x
        1133             : !  else
        1134             : !    point=tmpgnom%y
        1135             : !  end if
        1136             :   !
        1137             :   ! in which cell is interpolation point located? gno(iref) is location of point to the right that is closest
        1138             :   !
        1139             :   ! |----------|---------|------x---|----------|------|------
        1140             :   !                 gno(iref-1)  gno(iref)
        1141             :   !
        1142       33312 :   iref=ida
        1143      181584 :   do while (point>gnom1d(iref))
        1144      148320 :     iref = iref + 1
        1145      148320 :     if (iref>ide+1) then
        1146           0 :       call endrun("error in search - ABORT; probably invalid ns-nc combination")
        1147             :     end if
        1148      181584 :     if (iref>ide) then
        1149             : !       write(*,*) "extrapolation in interpolation_point",iref,ide
        1150          48 :        iref=ide
        1151          48 :        exit
        1152             :     endif
        1153             :   end do
        1154             :   !
        1155             :   ! this routine works for ns=1 and ns even
        1156             :   !
        1157       33312 :   if (MOD(ns,2)==1) then
        1158       33312 :      iref = max(iref,ida+1)!make sure gnom1d does not go out of bounds for extrapolation
        1159       33312 :      if (gnom1d(iref)-point>point-gnom1d(iref-1)) iref=iref-1
        1160       33312 :      iref=iref-((ns-1)/2)
        1161       33312 :      ibaseref = min(max(iref,ida),ide-(ns-1))!extrapolation
        1162       33312 :      point=point-gnom1d(ibaseref)
        1163           0 :   else if (MOD(ns, 2)==0) then
        1164             :      !
        1165             :      ! this code is only coded for ns even
        1166             :      !
        1167             :      ! ibaseref is the left most index used for 1D interpolation
        1168             :      ! (hence iref = iref-ns/2 except near corners)
        1169             :      !
        1170           0 :      iref = iref-ns/2
        1171           0 :      ibaseref = min(max(iref,ida),ide-(ns-1))
        1172           0 :      point=point-gnom1d(ibaseref)
        1173             :   end if
        1174       33312 : end subroutine interpolation_point
        1175             : !END SUBROUTINE INTERPOLATION_POINT---------------------------------------CE-for FVM!
        1176             : ! ---------------------------------------------------------------------!
        1177             : !                                                                      !
        1178             : ! Precompute weights for Lagrange interpolation                        !
        1179             : ! for equi-distant source grid values                                  !
        1180             : !                                                                      !
        1181             : !----------------------------------------------------------------------!
        1182             :
        1183       33312 : subroutine get_equispace_weights(dx, x, w,ns)
        1184             :   !
        1185             :   ! Coordinate system for Lagrange interpolation:
        1186             :   !
        1187             :   ! |------|------|------|------|
        1188             :   ! 0     dx    2*dx   3*dx   ns*dx
        1189             :   !
        1190             :   implicit none
        1191             :   real (kind=r8),intent(in)                  :: dx  ! spacing of points, alpha/beta
        1192             :   real (kind=r8),intent(in)                  :: x   ! X coordinate where interpolation is to be applied
        1193             :   real (kind=r8),dimension(:),intent(out)    :: w   ! dimension(ns)
        1194             :   integer              ,intent(in)                  :: ns
        1195             :   !
        1196             :   integer :: j,k
        1197             :   !
        1198             :   ! use Lagrange interpolation formulae, e.g.,:
        1199             :   !
        1200             :   !                http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html
        1201             :   !
        1202      133248 :   w = 1.0_r8
        1203       33312 :   if (ns.ne.1) then
        1204      133248 :      do j=1,ns
        1205      433056 :         do k=1,ns
        1206      399744 :            if (k.ne.j) then
        1207      199872 :               w(j)=w(j)*(x-dble(k-1)*dx)/(dble(j-1)*dx-dble(k-1)*dx)
        1208             :            end if
        1209             :         end do
        1210             :      end do
        1211             :   end if
        1212       33312 : end subroutine get_equispace_weights
        1213             :
        1214             : end module fvm_analytic_mod

Generated by: LCOV version 1.14