LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - fvm_overlap_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 139 | 330 | 42.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 8 | 11 | 72.7 %

              Line data    Source code
           1             : module fvm_overlap_mod
           2             :   use shr_kind_mod,   only: r8=>shr_kind_r8
           3             :
           4             :   real (kind=r8),parameter, private  :: bignum = 1.0e20_r8
           5             :   real (kind=r8),parameter, private  :: tiny   = 1.0e-12_r8
           6             :   real (kind=r8),parameter, private  :: fuzzy_width = 10.0_r8*tiny
           7             :
           8             :   public:: compute_weights_cell
           9             :
          10             :   private
          11             :   integer, parameter :: max_cross = 10
          12             : contains
          13      194400 :   subroutine compute_weights_cell(nvertex,lexact_horizontal_line_integrals,&
          14      194400 :        xcell_in,ycell_in,jx,jy,nreconstruction,xgno,ygno,igno_min,igno_max,&
          15             :        jx_min, jx_max, jy_min, jy_max,&
          16      194400 :        ngauss,gauss_weights,abscissae,weights,weights_eul_index,jcollect,jmax_segments)
          17             :
          18             :     implicit none
          19             :     integer , intent(in) :: nvertex
          20             :     logical, intent(in) :: lexact_horizontal_line_integrals
          21             :     integer , intent(in):: nreconstruction, jx,jy,ngauss,jmax_segments
          22             :     !
          23             :     ! dimension(nvertex)
          24             :     !
          25             :     real (kind=r8)   ,  dimension(4), intent(in):: xcell_in,ycell_in
          26             :     !
          27             :     integer , intent(in)               :: jx_min, jy_min, jx_max, jy_max,igno_min,igno_max
          28             :     !
          29             :     ! dimension(-ihalo:nc+2+ihalo)
          30             :     !
          31             :     real (kind=r8), dimension(igno_min:igno_max), intent(in) :: xgno, ygno
          32             :     !
          33             :     ! for Gaussian quadrature
          34             :     !
          35             :     real (kind=r8), dimension(:), intent(in) :: gauss_weights, abscissae !dimension(ngauss)
          36             :     !
          37             :     ! Number of Eulerian sub-cell integrals for the cell in question
          38             :     !
          39             :     integer , intent(out)       :: jcollect
          40             :     !
          41             :     ! local workspace
          42             :     !
          43             :     !
          44             :     ! max number of line segments is:
          45             :     !
          46             :     ! (number of longitudes)*(max average number of crossings per line segment = 3)*ncube*2
          47             :     !
          48             :     real (kind=r8)   ,  &
          49             :          dimension(jmax_segments,nreconstruction), intent(out) :: weights
          50             :     integer ,  &
          51             :          dimension(jmax_segments,2), intent(out)      :: weights_eul_index
          52             :
          53             :     integer :: jsegment
          54             :     !
          55             :     ! variables for registering crossings with Eulerian latitudes and longitudes
          56             :     !
          57             :     integer :: jcross_lat
          58             :     !
          59             :     ! max. crossings per side is 2*ihalo
          60             :     !
          61             :     real (kind=r8), &
          62             :          dimension(max_cross,2) :: r_cross_lat
          63             :     integer , &
          64             :          dimension(max_cross,2) :: cross_lat_eul_index
          65      388800 :     real (kind=r8)   ,  dimension(nvertex) :: xcell,ycell
          66             :
          67     1166400 :     xcell = xcell_in(1:nvertex)
          68      972000 :     ycell = ycell_in(1:nvertex)
          69             :
          70      194400 :     jsegment   = 0
          71    59680800 :     weights    = 0.0_r8
          72      194400 :     jcross_lat = 0
          73             :
          74             :     call side_integral(lexact_horizontal_line_integrals,xcell,ycell,nvertex,jsegment,jmax_segments,&
          75             :          weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,igno_min,igno_max,jx_min, jx_max, jy_min, jy_max,&
          76             :          ngauss,gauss_weights,abscissae,&
          77      194400 :          jcross_lat,r_cross_lat,cross_lat_eul_index)
          78             :     !
          79             :     !**********************
          80             :     !
          81             :     ! Do inner integrals
          82             :     !
          83             :     !**********************
          84             :     !
          85             :     call compute_inner_line_integrals_lat(lexact_horizontal_line_integrals,&
          86             :          r_cross_lat,cross_lat_eul_index,&
          87             :          jcross_lat,jsegment,xgno,igno_min,igno_max,jx_min, jx_max, jy_min, jy_max,&
          88             :          weights,weights_eul_index,&
          89      194400 :          nreconstruction,ngauss,gauss_weights,abscissae)
          90             :
          91      194400 :     IF (ABS((jcross_lat/2)-DBLE(jcross_lat)/2.0_r8)>tiny) then
          92           0 :       WRITE(*,*) "number of latitude crossings are not even: ABORT",jcross_lat,jx,jy
          93           0 :       STOP
          94             :     END IF
          95             :
          96             :     !
          97             :     ! collect line-segment that reside in the same Eulerian cell
          98             :     !
          99      194400 :     if (jsegment>0) then
         100      194400 :       call collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments)
         101             :     else
         102           0 :       jcollect = 0
         103             :     end if
         104      194400 :   end subroutine compute_weights_cell
         105             :   !
         106             :   !****************************************************************************
         107             :   !
         108             :   ! organize data and store it
         109             :   !
         110             :   !****************************************************************************
         111             :   !
         112      194400 :   subroutine collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments)
         113             :     implicit none
         114             :     integer ,                                  INTENT(IN   ) :: jsegment,jmax_segments
         115             :     integer , intent(in)    :: nreconstruction
         116             :     !
         117             :     real (kind=r8)  , dimension(:,:), intent(inout) :: weights !dimension(jmax_segments,nreconstruction)
         118             :     integer , dimension(:,:), intent(inout) :: weights_eul_index !dimension(jmax_segments,2)
         119             :     integer ,                                  INTENT(OUT  ) :: jcollect
         120             :     !
         121             :     ! local workspace
         122             :     !
         123             :     integer :: imin, imax, jmin, jmax, i,j,k,h
         124             :     logical                 :: ltmp
         125             :
         126      388800 :     real (kind=r8)  , dimension(jmax_segments,nreconstruction) :: weights_out
         127      388800 :     integer , dimension(jmax_segments,2     ) :: weights_eul_index_out
         128             :
         129    59680800 :     weights_out           = 0.0_r8
         130    20023200 :     weights_eul_index_out = -100
         131             :
         132     1166400 :     imin = MINVAL(weights_eul_index(1:jsegment,1))
         133      972000 :     imax = MAXVAL(weights_eul_index(1:jsegment,1))
         134      972000 :     jmin = MINVAL(weights_eul_index(1:jsegment,2))
         135      972000 :     jmax = MAXVAL(weights_eul_index(1:jsegment,2))
         136             :
         137      194400 :     ltmp = .FALSE.
         138             :
         139      194400 :     jcollect = 1
         140             :
         141      388800 :     do j=jmin,jmax
         142      583200 :        do i=imin,imax
         143      972000 :           do k=1,jsegment
         144      972000 :              if (weights_eul_index(k,1)==i.AND.weights_eul_index(k,2)==j) then
         145             :                 weights_out(jcollect,1:nreconstruction) = &
         146     5443200 :                 weights_out(jcollect,1:nreconstruction) + weights(k,1:nreconstruction)
         147             :                 ltmp = .TRUE.
         148             :                 h = k
         149             :              endif
         150             :           enddo
         151      194400 :           if (ltmp) then
         152      583200 :              weights_eul_index_out(jcollect,:) = weights_eul_index(h,:)
         153      194400 :              jcollect = jcollect+1
         154             :           endif
         155      388800 :           ltmp = .FALSE.
         156             :        enddo
         157             :     enddo
         158      194400 :     jcollect = jcollect-1
         159    59875200 :     weights           = weights_out
         160    20217600 :     weights_eul_index = weights_eul_index_out
         161      194400 :   end subroutine collect
         162             :   !
         163             :   !*****************************************************************************************
         164             :   !
         165             :   ! compute crossings with Eulerian latitudes and longitudes
         166             :   !
         167             :   !*****************************************************************************************
         168             :   !
         169      583200 :   subroutine compute_inner_line_integrals_lat(lexact_horizontal_line_integrals,r_cross_lat,&
         170      194400 :        cross_lat_eul_index,&
         171      194400 :        jcross_lat,jsegment,xgno,igno_min,igno_max,jx_min,jx_max,jy_min, jy_max,weights,weights_eul_index,&
         172      194400 :        nreconstruction,ngauss,gauss_weights,abscissae)
         173             :     implicit none
         174             :     logical, intent(in) :: lexact_horizontal_line_integrals
         175             :     !
         176             :     ! variables for registering crossings with Eulerian latitudes and longitudes
         177             :     !
         178             :     integer ,         intent(in):: jcross_lat, nreconstruction,ngauss,igno_min,igno_max
         179             :     integer ,         intent(inout):: jsegment
         180             :     !
         181             :     ! for Gaussian quadrature
         182             :     !
         183             :     real (kind=r8), dimension(ngauss), intent(in) :: gauss_weights, abscissae
         184             :     !
         185             :     ! max. crossings per side is 2*ihalo
         186             :     !
         187             :
         188             :     real (kind=r8)  , dimension(:,:), intent(in):: r_cross_lat ! dimension(8*ihalo,2)
         189             :     integer , dimension(:,:), intent(in):: cross_lat_eul_index ! ! dimension(8*ihalo,2)
         190             :     integer , intent(in):: jx_min, jx_max, jy_min, jy_max
         191             :
         192             :     real (kind=r8), dimension(igno_min:igno_max), intent(in)  :: xgno !dimension(-ihalo:nc+2+ihalo)
         193             :     !
         194             :     ! dimension(jmax_segments,nreconstruction)
         195             :     !
         196             :     real (kind=r8), dimension(:,:), intent(inout) :: weights
         197             :     !
         198             :     ! dimension(jmax_segments,2)
         199             :     !
         200             :     integer , dimension(:,:), intent(inout)               :: weights_eul_index
         201             :
         202      388800 :     real (kind=r8)   , dimension(nreconstruction)         :: weights_tmp
         203             :     integer :: imin,imax,i,j,k,h
         204             :     real (kind=r8), dimension(2)  :: rstart,rend,rend_tmp
         205             :     real (kind=r8), dimension(2)  :: xseg, yseg
         206             :     5   FORMAT(10e14.6)
         207      194400 :     if (jcross_lat>0) then
         208           0 :        do i=MINVAL(cross_lat_eul_index(1:jcross_lat,2)),MAXVAL(cross_lat_eul_index(1:jcross_lat,2))
         209             :           !
         210             :           ! find "first" crossing with Eulerian cell i
         211             :           !
         212           0 :           do k=1,jcross_lat
         213           0 :              if (cross_lat_eul_index(k,2)==i) exit
         214             :           enddo
         215           0 :           do j=k+1,jcross_lat
         216             :              !
         217             :              ! find "second" crossing with Eulerian cell i
         218             :              !
         219           0 :              if (cross_lat_eul_index(j,2)==i) then
         220           0 :                 if (r_cross_lat(k,1)<r_cross_lat(j,1)) then
         221           0 :                    rstart = r_cross_lat(k,1:2)
         222           0 :                    rend   = r_cross_lat(j,1:2)
         223           0 :                    imin   = cross_lat_eul_index(k,1)
         224           0 :                    imax   = cross_lat_eul_index(j,1)
         225             :                 else
         226           0 :                    rstart = r_cross_lat(j,1:2)
         227           0 :                    rend   = r_cross_lat(k,1:2)
         228           0 :                    imin   = cross_lat_eul_index(j,1)
         229           0 :                    imax   = cross_lat_eul_index(k,1)
         230             :                 endif
         231           0 :                 do h=imin,imax
         232           0 :                    if (h==imax) then
         233           0 :                       rend_tmp = rend
         234             :                    else
         235           0 :                       rend_tmp(1) = xgno(h+1)
         236           0 :                       rend_tmp(2) = r_cross_lat(k,2)
         237             :                    endif
         238           0 :                    xseg(1) = rstart(1)
         239           0 :                    xseg(2) = rend_tmp(1)
         240           0 :                    yseg(1) = rstart(2)
         241           0 :                    yseg(2) = rend_tmp(2)
         242             :                    call get_weights_exact(lexact_horizontal_line_integrals, weights_tmp,xseg,yseg,&
         243           0 :                         nreconstruction,ngauss,gauss_weights,abscissae)
         244             :
         245           0 :                    if (i.LE.jy_max-1.AND.i.GE.jy_min.AND.h.LE.jx_max-1.AND.h.GE.jx_min) then
         246           0 :                       jsegment=jsegment+1
         247           0 :                       weights_eul_index(jsegment,1) = h
         248           0 :                       weights_eul_index(jsegment,2) = i
         249           0 :                       weights(jsegment,1:nreconstruction) = -weights_tmp
         250             :                    endif
         251             :                    !
         252             :                    ! subtract the same weights on the west side of the line
         253             :                    !
         254           0 :                    if (i.LE.jy_max.AND.i.GE.jy_min+1.AND.h.LE.jx_max-1.AND.h.GE.jx_min) then
         255           0 :                       jsegment = jsegment+1
         256           0 :                       weights_eul_index(jsegment,1) = h
         257           0 :                       weights_eul_index(jsegment,2) = i-1
         258           0 :                       weights(jsegment,1:nreconstruction) = weights_tmp
         259             :                    endif
         260             :                    !
         261             :                    ! prepare for next iteration
         262             :                    !
         263           0 :                    rstart = rend_tmp
         264             :                 enddo
         265             :              endif
         266             :           enddo
         267             :        enddo
         268             :     endif
         269      194400 :   end subroutine compute_inner_line_integrals_lat
         270             :
         271             :   !
         272             :   ! line integral from (a1_in,a2_in) to (b1_in,b2_in)
         273             :   ! If line is coniciding with an Eulerian longitude or latitude the routine
         274             :   ! needs to know where an adjacent side is located to determine which
         275             :   ! reconstruction must be used. therefore (c1,c2) is passed to the routine
         276             :   !
         277             :   !
         278             :
         279      194400 :   subroutine side_integral(lexact_horizontal_line_integrals,&
         280      194400 :        x_in,y_in,nvertex,jsegment,jmax_segments,&
         281      194400 :        weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,igno_min,igno_max,&
         282             :        jx_min,jx_max,jy_min,jy_max,&
         283      194400 :        ngauss,gauss_weights,abscissae,&!)!phl add jx_min etc.
         284             :        jcross_lat,r_cross_lat,cross_lat_eul_index)
         285             :     implicit none
         286             :
         287             :
         288             :     logical, intent(in) :: lexact_horizontal_line_integrals
         289             :     integer ,            intent(in)    :: nreconstruction,jx,jy,jmax_segments,ngauss
         290             :     integer , intent(in)               :: nvertex,igno_min,igno_max
         291             :     !
         292             :     ! for Gaussian quadrature
         293             :     !
         294             :     real (kind=r8), dimension(:), intent(in) :: gauss_weights, abscissae !dimension(ngauss)
         295             :     real (kind=r8), dimension(:)        , intent(in)    :: x_in,y_in !dimension(1:nvertex)
         296             :
         297             :     integer , intent(in)               :: jx_min, jy_min, jx_max, jy_max
         298             :     real (kind=r8), dimension(igno_min:igno_max), intent(in) :: xgno, ygno !dimension(-ihalo:nc+2+ihalo)
         299             :     integer ,            intent(inout) :: jsegment
         300             : !    integer ,dimension(0:2),intent(in)    :: jx_eul_in, jy_eul_in
         301             :     real (kind=r8)   , dimension(:,:), intent(out) :: weights !dimension(jmax_segments,nreconstruction)
         302             :     integer ,  &
         303             :          dimension(jmax_segments,2), intent(out) :: weights_eul_index
         304             :
         305             :     !
         306             :     ! variables for registering crossings with Eulerian latitudes and longitudes
         307             :     !
         308             :     integer ,         intent(inout):: jcross_lat
         309             :     !
         310             :     ! max. crossings per side is 2*ihalo
         311             :     !
         312             :     real (kind=r8), &
         313             :          dimension(max_cross,2), intent(inout):: r_cross_lat
         314             :     integer , &
         315             :          dimension(max_cross,2), intent(inout):: cross_lat_eul_index
         316             :     !
         317             :     ! local variables
         318             :     !
         319             :     real (kind=r8), dimension(2) :: xseg,yseg
         320             :     real (kind=r8), dimension(0:3) :: x,y
         321             :     real (kind=r8)               :: xeul,yeul,xcross,ycross,slope
         322             :     integer ::    jx_eul_tmp,jy_eul_tmp
         323             :     integer :: xsgn1,ysgn1,xsgn2,ysgn2
         324             :     integer :: iter
         325             :     logical :: lcontinue, lsame_cell_x, lsame_cell_y
         326             :
         327             :     integer :: jx_eul, jy_eul, side_count
         328      388800 :     real (kind=r8), dimension(0:nvertex+2)  :: xcell,ycell
         329             :
         330             :
         331             : 5   FORMAT(10e14.6)
         332             :     !
         333             :     !***********************************************
         334             :     !
         335             :     ! find jx_eul and jy_eul for (x(1),y(1))
         336             :     !
         337             :     !***********************************************
         338             :     !
         339      194400 :     jx_eul = jx; jy_eul = jy
         340     1749600 :     xcell(1:nvertex)=x_in; ycell(1:nvertex)=y_in
         341      972000 :     DO iter=1,nvertex
         342      777600 :       CALL truncate_vertex(xcell(iter),jx_eul,xgno,igno_min,igno_max)
         343      972000 :       CALL truncate_vertex(ycell(iter),jy_eul,ygno,igno_min,igno_max)
         344             :     END DO
         345      194400 :     xcell(0) = xcell(nvertex); xcell(nvertex+1)=xcell(1); xcell(nvertex+2)=xcell(2);
         346      194400 :     ycell(0) = ycell(nvertex); ycell(nvertex+1)=ycell(1); ycell(nvertex+2)=ycell(2);
         347             :
         348             :
         349     7192800 :     IF ((&
         350      388800 :          MAXVAL(xcell).LE.xgno(jx_min).OR.MINVAL(xcell).GE.xgno(jx_max).OR.&
         351      388800 :          MAXVAL(ycell).LE.ygno(jy_min).OR.MINVAL(ycell).GE.ygno(jy_max))) THEN
         352             :       !
         353             :       ! entire cell off panel
         354             :       !
         355             :     ELSE
         356      194400 :       jx_eul = MIN(MAX(jx,jx_min),jx_max-1)
         357      194400 :       jy_eul = MIN(MAX(jy,jy_min),jy_max-1)
         358      194400 :       CALL which_eul_cell(xcell(1:3),jx_eul,xgno,igno_min,igno_max)
         359      194400 :       CALL which_eul_cell(ycell(1:3),jy_eul,ygno,igno_min,igno_max)
         360             :
         361      194400 :       side_count = 1
         362      972000 :       DO WHILE (side_count<nvertex+1)
         363      777600 :         iter = 0
         364      777600 :         lcontinue = .TRUE.
         365     6998400 :         x(0:3) = xcell(side_count-1:side_count+2); y(0:3) = ycell(side_count-1:side_count+2);
         366     1555200 :         DO while (lcontinue)
         367      777600 :           iter = iter+1
         368      777600 :           IF (iter>10) THEN
         369           0 :             WRITE(*,*) "search not converging",iter
         370           0 :             STOP
         371             :           END IF
         372      777600 :           lsame_cell_x = (x(2).GE.xgno(jx_eul).AND.x(2).LE.xgno(jx_eul+1))
         373      777600 :           lsame_cell_y = (y(2).GE.ygno(jy_eul).AND.y(2).LE.ygno(jy_eul+1))
         374      777600 :           IF (lsame_cell_x.AND.lsame_cell_y) THEN
         375             :             !
         376             :             !****************************
         377             :             !
         378             :             ! same cell integral
         379             :             !
         380             :             !****************************
         381             :             !
         382      777600 :             xseg(1) = x(1); yseg(1) = y(1); xseg(2) = x(2); yseg(2) = y(2)
         383      777600 :             jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
         384      777600 :             lcontinue = .FALSE.
         385             :             !
         386             :             ! prepare for next side if (x(2),y(2)) is on a grid line
         387             :             !
         388      777600 :             IF (x(2).EQ.xgno(jx_eul+1).AND.x(3)>xgno(jx_eul+1)) THEN
         389             :               !
         390             :               ! cross longitude jx_eul+1
         391             :               !
         392           0 :               jx_eul=jx_eul+1
         393      777600 :             ELSE IF (x(2).EQ.xgno(jx_eul  ).AND.x(3)<xgno(jx_eul)) THEN
         394             :               !
         395             :               ! cross longitude jx_eul
         396             :               !
         397           0 :               jx_eul=jx_eul-1
         398             :             END IF
         399      777600 :             IF (y(2).EQ.ygno(jy_eul+1).AND.y(3)>ygno(jy_eul+1)) THEN
         400             :               !
         401             :               ! register crossing with latitude: line-segments point Northward
         402             :               !
         403           0 :               jcross_lat = jcross_lat + 1
         404           0 :               jy_eul     = jy_eul     + 1
         405           0 :               cross_lat_eul_index(jcross_lat,1) = jx_eul
         406           0 :               cross_lat_eul_index(jcross_lat,2) = jy_eul
         407           0 :               r_cross_lat(jcross_lat,1) = x(2)
         408           0 :               r_cross_lat(jcross_lat,2) = y(2)
         409             : !              write(*,*) "A register crossing with latitude",x(2),y(2),jx_eul,jy_eul
         410      777600 :             ELSE IF (y(2).EQ.ygno(jy_eul  ).AND.y(3)<ygno(jy_eul)) THEN
         411             :               !
         412             :               ! register crossing with latitude: line-segments point Southward
         413             :               !
         414           0 :               jcross_lat = jcross_lat+1
         415           0 :               cross_lat_eul_index(jcross_lat,1) = jx_eul
         416           0 :               cross_lat_eul_index(jcross_lat,2) = jy_eul
         417           0 :               r_cross_lat(jcross_lat,1) = x(2)
         418           0 :               r_cross_lat(jcross_lat,2) = y(2)
         419             : !              write(*,*) "B register crossing with latitude",x(2),y(2),jx_eul,jy_eul
         420             :               !
         421           0 :               jy_eul=jy_eul-1
         422             :             END IF
         423             :             lcontinue=.FALSE.
         424             :           ELSE
         425             :             !
         426             :             !****************************
         427             :             !
         428             :             ! not same cell integral
         429             :             !
         430             :             !****************************
         431             :             !
         432           0 :             IF (lsame_cell_x) THEN
         433           0 :               ysgn1 = (1+INT(SIGN(1.0_r8,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0"
         434           0 :               ysgn2 = INT(SIGN(1.0_r8,y(2)-y(1)))       !"1" if y(2)>y(1) else "-1"
         435             :               !
         436             :               !*******************************************************************************
         437             :               !
         438             :               ! there is at least one crossing with latitudes but no crossing with longitudes
         439             :               !
         440             :               !*******************************************************************************
         441             :               !
         442           0 :               yeul   = ygno(jy_eul+ysgn1)
         443           0 :               IF (x(1).EQ.x(2)) THEN
         444             :                 !
         445             :                 ! line segment is parallel to longitude (infinite slope)
         446             :                 !
         447             :                 xcross = x(1)
         448             :               ELSE
         449           0 :                 slope  = (y(2)-y(1))/(x(2)-x(1))
         450           0 :                 xcross = x_cross_eul_lat(x(1),y(1),yeul,slope)
         451             :                 !
         452             :                 ! constrain crossing to be "physically" possible
         453             :                 !
         454           0 :                 xcross = MIN(MAX(xcross,xgno(jx_eul)),xgno(jx_eul+1))
         455             :                 !
         456             :                 ! debugging
         457             :                 !
         458           0 :                 IF (xcross.GT.xgno(jx_eul+1).OR.xcross.LT.xgno(jx_eul)) THEN
         459           0 :                   WRITE(*,*) "xcross is out of range",jx,jy
         460           0 :                   WRITE(*,*) "xcross-xgno(jx_eul+1), xcross-xgno(jx_eul))",&
         461           0 :                        xcross-xgno(jx_eul+1), xcross-ygno(jx_eul)
         462           0 :                   STOP
         463             :                 END IF
         464             :               END IF
         465           0 :               xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul
         466           0 :               jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
         467             :               !
         468             :               ! prepare for next iteration
         469             :               !
         470           0 :               x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2
         471             :               !
         472             :               ! register crossing with latitude
         473             :               !
         474           0 :               jcross_lat = jcross_lat+1
         475           0 :               cross_lat_eul_index(jcross_lat,1) = jx_eul
         476           0 :               if (ysgn2>0) then
         477           0 :                 cross_lat_eul_index(jcross_lat,2) = jy_eul
         478             :               else
         479           0 :                 cross_lat_eul_index(jcross_lat,2) = jy_eul+1
         480             :               end if
         481           0 :               r_cross_lat(jcross_lat,1) = xcross
         482           0 :               r_cross_lat(jcross_lat,2) = yeul
         483           0 :             ELSE IF (lsame_cell_y) THEN
         484             :               !
         485             :               !*******************************************************************************
         486             :               !
         487             :               ! there is at least one crossing with longitudes but no crossing with latitudes
         488             :               !
         489             :               !*******************************************************************************
         490             :               !
         491           0 :               xsgn1 = (1+INT(SIGN(1.0_r8,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0"
         492           0 :               xsgn2 = INT(SIGN(1.0_r8,x(2)-x(1))) !"1" if x(2)>x(1) else "-1"
         493           0 :               xeul   = xgno(jx_eul+xsgn1)
         494           0 :               IF (ABS(x(2)-x(1))<fuzzy_width) THEN
         495             :                 ! fuzzy crossing
         496           0 :                 ycross = 0.5_r8*(y(2)-y(1))
         497             :               ELSE
         498           0 :                 slope  = (y(2)-y(1))/(x(2)-x(1))
         499           0 :                 ycross = y_cross_eul_lon(x(1),y(1),xeul,slope)
         500             :               END IF
         501             :               !
         502             :               ! constrain crossing to be "physically" possible
         503             :               !
         504           0 :               ycross = MIN(MAX(ycross,ygno(jy_eul)),ygno(jy_eul+1))
         505             :               !
         506             :               ! debugging
         507             :               !
         508           0 :               IF (ycross.GT.ygno(jy_eul+1).OR.ycross.LT.ygno(jy_eul)) THEN
         509           0 :                 WRITE(*,*) "ycross is out of range"
         510           0 :                 WRITE(*,*) "jx,jy,jx_eul,jy_eul",jx,jy,jx_eul,jy_eul
         511           0 :                 WRITE(*,*) "ycross-ygno(jy_eul+1), ycross-ygno(jy_eul))",&
         512           0 :                      ycross-ygno(jy_eul+1), ycross-ygno(jy_eul)
         513           0 :                 STOP
         514             :               END IF
         515           0 :               xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross
         516           0 :               jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
         517             :               !
         518             :               ! prepare for next iteration
         519             :               !
         520           0 :               x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2
         521             :             ELSE
         522             :               !
         523             :               !*******************************************************************************
         524             :               !
         525             :               ! there are crossings with longitude(s) and latitude(s)
         526             :               !
         527             :               !*******************************************************************************
         528             :               !
         529           0 :               xsgn1 = (1+INT(SIGN(1.0_r8,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0"
         530           0 :               xsgn2 = (INT(SIGN(1.0_r8,x(2)-x(1)))) !"1" if x(2)>x(1) else "0"
         531           0 :               xeul   = xgno(jx_eul+xsgn1)
         532           0 :               ysgn1 = (1+INT(SIGN(1.0_r8,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0"
         533           0 :               ysgn2 = INT(SIGN(1.0_r8,y(2)-y(1)))       !"1" if y(2)>y(1) else "-1"
         534           0 :               yeul   = ygno(jy_eul+ysgn1)
         535             :
         536           0 :               slope  = (y(2)-y(1))/(x(2)-x(1))
         537           0 :               IF (ABS(x(2)-x(1))<fuzzy_width) THEN
         538           0 :                 ycross = 0.5_r8*(y(2)-y(1))
         539             :               ELSE
         540           0 :                 ycross = y_cross_eul_lon(x(1),y(1),xeul,slope)
         541             :               END IF
         542           0 :               xcross = x_cross_eul_lat(x(1),y(1),yeul,slope)
         543             :
         544             :
         545           0 :               IF ((xsgn2>0.AND.xcross.LE.xeul).OR.(xsgn2<0.AND.xcross.GE.xeul)) THEN
         546             :                 !
         547             :                 ! cross latitude
         548             :                 !
         549           0 :                 xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul
         550           0 :                 jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
         551             :                 !
         552             :                 ! prepare for next iteration
         553             :                 !
         554           0 :                 x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2
         555             :                 !
         556             :                 ! register crossing with latitude
         557             :                 !
         558           0 :                 jcross_lat = jcross_lat+1
         559           0 :                 cross_lat_eul_index(jcross_lat,1) = jx_eul
         560           0 :                 if (ysgn2>0) then
         561           0 :                   cross_lat_eul_index(jcross_lat,2) = jy_eul
         562             :                 else
         563           0 :                   cross_lat_eul_index(jcross_lat,2) = jy_eul+1
         564             :                 end if
         565           0 :                 r_cross_lat(jcross_lat,1) = xcross
         566           0 :                 r_cross_lat(jcross_lat,2) = yeul
         567             : !              write(*,*) "D register crossing with latitude",xcross,yeul,jx_eul,cross_lat_eul_index(jcross_lat,2)
         568             :               ELSE
         569             :                 !
         570             :                 ! cross longitude
         571             :                 !
         572           0 :                 xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross
         573           0 :                 jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
         574             :                 !
         575             :                 ! prepare for next iteration
         576             :                 !
         577           0 :                 x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2
         578             :               END IF
         579             :
         580             :             END IF
         581             :           END IF
         582             :           !
         583             :           ! register line-segment (don't register line-segment if outside of panel)
         584             :           !
         585             :           if (jx_eul_tmp>=jx_min.AND.jy_eul_tmp>=jy_min.AND.&
         586     1555200 :                jx_eul_tmp<=jx_max-1.AND.jy_eul_tmp<=jy_max-1) then
         587      777600 :             jsegment=jsegment+1
         588      777600 :             weights_eul_index(jsegment,1) = jx_eul_tmp
         589      777600 :             weights_eul_index(jsegment,2) = jy_eul_tmp
         590             :
         591             :             call get_weights_exact(lexact_horizontal_line_integrals.AND.ABS(yseg(2)-yseg(1))<tiny,&
         592             :                  weights(jsegment,:),&
         593      777600 :                  xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
         594             : !old            call get_weights_gauss(weights(jsegment,1:nreconstruction),&
         595             : !old                 xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
         596             :           ELSE
         597             :             !
         598             :             ! segment outside of panel
         599             :             !
         600             :           END IF
         601             :
         602             :         END DO
         603      777600 :         side_count = side_count+1
         604             :       END DO
         605             :     END IF
         606      194400 :   end subroutine side_integral
         607             :
         608             :
         609           0 :   real (kind=r8) function y_cross_eul_lon(x,y,xeul,slope)
         610             :     implicit none
         611             :     real (kind=r8), intent(in) :: x,y
         612             :     real (kind=r8)              , intent(in) :: xeul,slope
         613             :     !
         614             :     ! line: y=a*x+b
         615             :     !
         616             :     real (kind=r8) :: b
         617             :
         618           0 :     b = y-slope*x
         619           0 :     y_cross_eul_lon = slope*xeul+b
         620           0 :   end function y_cross_eul_lon
         621             :
         622           0 :   real (kind=r8) function x_cross_eul_lat(x,y,yeul,slope)
         623             :     implicit none
         624             :     real (kind=r8), intent(in) :: x,y
         625             :     real (kind=r8)              , intent(in) :: yeul,slope
         626             :
         627           0 :     if (fuzzy(ABS(slope),fuzzy_width)>0) THEN
         628           0 :         x_cross_eul_lat = x+(yeul-y)/slope
         629             :     ELSE
         630             :       x_cross_eul_lat = bignum
         631             :     END IF
         632           0 :   end function x_cross_eul_lat
         633             :
         634      777600 :   subroutine get_weights_exact(lexact_horizontal_line_integrals,weights,xseg,yseg,nreconstruction,&
         635      777600 :        ngauss,gauss_weights,abscissae)
         636             :     use fvm_analytic_mod, only: I_00, I_10, I_01, I_20, I_02, I_11
         637             :     implicit none
         638             :     logical, intent(in) :: lexact_horizontal_line_integrals
         639             :     integer , intent(in) :: nreconstruction, ngauss
         640             :     real (kind=r8), intent(out) :: weights(:)
         641             :     real (kind=r8), dimension(:), intent(in) :: gauss_weights, abscissae !dimension(ngauss)
         642             :
         643             :
         644             :     real (kind=r8), dimension(:), intent(in) :: xseg,yseg !dimension(2)
         645             :     !
         646             :     ! compute weights
         647             :     !
         648      777600 :     if(lexact_horizontal_line_integrals) then
         649      388800 :       weights(1) = ((I_00(xseg(2),yseg(2))-I_00(xseg(1),yseg(1))))
         650      388800 :       if (ABS(weights(1))>1.0_r8) THEN
         651           0 :         WRITE(*,*) "1 exact weights(jsegment)",weights(1),xseg,yseg
         652           0 :         stop
         653             :       end if
         654      388800 :       if (nreconstruction>1) then
         655      388800 :          weights(2) = ((I_10(xseg(2),yseg(2))-I_10(xseg(1),yseg(1))))
         656      388800 :          weights(3) = ((I_01(xseg(2),yseg(2))-I_01(xseg(1),yseg(1))))
         657             :       endif
         658      388800 :       if (nreconstruction>3) then
         659      388800 :          weights(4) = ((I_20(xseg(2),yseg(2))-I_20(xseg(1),yseg(1))))
         660      388800 :          weights(5) = ((I_02(xseg(2),yseg(2))-I_02(xseg(1),yseg(1))))
         661      388800 :          weights(6) = ((I_11(xseg(2),yseg(2))-I_11(xseg(1),yseg(1))))
         662             :       endif
         663             :     else
         664      388800 :       call get_weights_gauss(weights,xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
         665             :     endif
         666      777600 :   end subroutine get_weights_exact
         667             :
         668             :
         669             :
         670      388800 :   subroutine get_weights_gauss(weights,xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
         671             :     use fvm_analytic_mod, only: F_00, F_10, F_01, F_20, F_02, F_11
         672             :     implicit none
         673             :     integer , intent(in) :: nreconstruction,ngauss
         674             :     real (kind=r8), intent(out) :: weights(:)
         675             :     real (kind=r8), dimension(2     ), intent(in) :: xseg,yseg
         676             :     real (kind=r8) :: slope
         677             :     !
         678             :     ! compute weights
         679             :     !
         680             :     !
         681             :     ! for Gaussian quadrature
         682             :     !
         683             :     real (kind=r8), dimension(ngauss), intent(in) :: gauss_weights, abscissae
         684             :
         685             :     ! if line-segment parallel to x or y use exact formulaes else use qudrature
         686             :     !
         687             :     real (kind=r8) :: b,integral,dx2,xc,x,y
         688             :     integer :: i
         689             :
         690             : !    if (fuzzy(abs(xseg(1) -xseg(2)),fuzzy_width)==0)then
         691      388800 :     if (xseg(1).EQ.xseg(2))then
         692     2721600 :       weights = 0.0_r8
         693             :     else
         694           0 :       slope    = (yseg(2)-yseg(1))/(xseg(2)-xseg(1))
         695           0 :       b        = yseg(1)-slope*xseg(1)
         696           0 :       dx2      = 0.5_r8*(xseg(2)-xseg(1))
         697           0 :       xc       = 0.5_r8*(xseg(1)+xseg(2))
         698           0 :       integral = 0.0_r8
         699           0 :       do i=1,ngauss
         700           0 :         x        = xc+abscissae(i)*dx2
         701           0 :         y        = slope*x+b
         702           0 :         integral = integral+gauss_weights(i)*F_00(x,y)
         703             :       enddo
         704           0 :       weights(1) = integral*dx2
         705           0 :       if (nreconstruction>1) then
         706             :         integral = 0.0_r8
         707           0 :         do i=1,ngauss
         708           0 :           x        = xc+abscissae(i)*dx2
         709           0 :           y        = slope*x+b
         710           0 :           integral = integral+gauss_weights(i)*F_10(x,y)
         711             :         enddo
         712           0 :         weights(2) = integral*dx2
         713           0 :         integral = 0.0_r8
         714           0 :         do i=1,ngauss
         715           0 :           x        = xc+abscissae(i)*dx2
         716           0 :           y        = slope*x+b
         717           0 :           integral = integral+gauss_weights(i)*F_01(x,y)
         718             :         enddo
         719           0 :         weights(3) = integral*dx2
         720             :       endif
         721           0 :       if (nreconstruction>3) then
         722             :         integral = 0.0_r8
         723           0 :         do i=1,ngauss
         724           0 :           x        = xc+abscissae(i)*dx2
         725           0 :           y        = slope*x+b
         726           0 :           integral = integral+gauss_weights(i)*F_20(x,y)
         727             :         enddo
         728           0 :         weights(4) = integral*dx2
         729           0 :         integral = 0.0_r8
         730           0 :         do i=1,ngauss
         731           0 :           x        = xc+abscissae(i)*dx2
         732           0 :           y        = slope*x+b
         733           0 :           integral = integral+gauss_weights(i)*F_02(x,y)
         734             :         enddo
         735           0 :         weights(5) = integral*dx2
         736           0 :         integral = 0.0_r8
         737           0 :         do i=1,ngauss
         738           0 :           x        = xc+abscissae(i)*dx2
         739           0 :           y        = slope*x+b
         740           0 :           integral = integral+gauss_weights(i)*F_11(x,y)
         741             :         enddo
         742           0 :         weights(6) = integral*dx2
         743             :       endif
         744             :     end if
         745      388800 :   end subroutine get_weights_gauss
         746             :
         747     1555200 :   subroutine truncate_vertex(x,j_eul,gno,igno_min,igno_max)
         748             :     implicit none
         749             :     integer , intent(inout) :: j_eul
         750             :     integer , intent(in)    :: igno_min,igno_max
         751             :
         752             :     real (kind=r8)                    , intent(inout)    :: x
         753             :     real (kind=r8), dimension(igno_min:igno_max), intent(in)    :: gno !dimension(-ihalo:nc+2+ihalo)
         754             : !    real (kind=r8), intent(in)    :: eps
         755             :
         756             :     logical                 :: lcontinue
         757             :     integer :: iter, xsgn
         758             :     real (kind=r8) :: dist,dist_new,tmp
         759             :
         760     1555200 :     lcontinue = .TRUE.
         761     1555200 :     iter = 0
         762     1555200 :     dist = bignum
         763             :
         764     1555200 :     xsgn     = INT(SIGN(1.0_r8,x-gno(j_eul)))
         765             :
         766     3693600 :     DO WHILE (lcontinue)
         767     2138400 :       if ((j_eul<igno_min) .or. (j_eul>igno_max)) then
         768           0 :         write(*,*) 'something is wrong', j_eul,igno_min,igno_max, iter
         769           0 :         stop
         770             :       endif
         771     2138400 :       iter     = iter+1
         772     2138400 :       tmp      = x-gno(j_eul)
         773     2138400 :       dist_new = ABS(tmp)
         774     2138400 :       IF (dist_new>dist) THEN
         775             :         lcontinue = .FALSE.
         776     2138400 :       ELSE IF (ABS(tmp)<1.0E-9_r8) THEN
         777     1555200 :         x = gno(j_eul)
         778     1555200 :         lcontinue = .FALSE.
         779             :       ELSE
         780      583200 :         j_eul = j_eul+xsgn
         781      583200 :         dist = dist_new
         782             :       END IF
         783     3693600 :       IF (iter>100) THEN
         784           0 :         WRITE(*,*) "truncate vertex not converging"
         785           0 :         STOP
         786             :       END IF
         787             :     END DO
         788     1555200 :   END subroutine truncate_vertex
         789             :
         790      388800 :   subroutine which_eul_cell(x,j_eul,gno,igno_min,igno_max)
         791             :     implicit none
         792             :     integer , intent(inout) :: j_eul
         793             :     integer , intent(in)    :: igno_min,igno_max
         794             :     real (kind=r8), dimension(:)                    , intent(in):: x !dimension(3)
         795             :     real (kind=r8), dimension(igno_min:igno_max), intent(in)    :: gno ! dimension(-ihalo:nc+2+ihalo)
         796             :
         797             :     logical :: lcontinue
         798             :     integer :: iter
         799             :
         800      388800 :     lcontinue = .TRUE.
         801      388800 :     iter = 0
         802             :
         803      777600 :     DO WHILE (lcontinue)
         804      388800 :       iter = iter+1
         805      388800 :       IF (x(1).GE.gno(j_eul).AND.x(1).LT.gno(j_eul+1)) THEN
         806      388800 :         lcontinue = .FALSE.
         807             :         !
         808             :         ! special case when x(1) is on top of grid line
         809             :         !
         810      388800 :         IF (x(1).EQ.gno(j_eul)) THEN
         811             :           !
         812             :           ! x(1) is on top of gno(J_eul)
         813             :           !
         814      388800 :           IF (x(2).GT.gno(j_eul)) THEN
         815             :             j_eul = j_eul
         816      194400 :           ELSE IF (x(2).LT.gno(j_eul)) THEN
         817           0 :             j_eul = j_eul-1
         818             :           ELSE
         819             :             !
         820             :             ! x(2) is on gno(j_eul) grid line; need x(3) to determine Eulerian cell
         821             :             !
         822      194400 :             IF (x(3).GT.gno(j_eul)) THEN
         823             :               !
         824             :               ! x(3) to the right
         825             :               !
         826             :               j_eul = j_eul
         827           0 :             ELSE IF (x(3).LT.gno(j_eul)) THEN
         828             :               !
         829             :               ! x(3) to the left
         830             :               !
         831           0 :               j_eul = j_eul-1
         832             :             ELSE
         833           0 :               WRITE(*,*) "inconsistent cell: x(1)=x(2)=x(3)",x(1),x(2),x(3)
         834           0 :               STOP
         835             :             END IF
         836             :           END IF
         837             :         END IF
         838             :       ELSE
         839             :         !
         840             :         ! searching - prepare for next iteration
         841             :         !
         842           0 :         IF (x(1).GE.gno(j_eul+1)) THEN
         843           0 :           j_eul = j_eul + 1
         844             :         ELSE
         845             :           !
         846             :           ! x(1).LT.gno(j_eul)
         847             :           !
         848           0 :           j_eul = j_eul - 1
         849             :         END IF
         850             :       END IF
         851      777600 :       IF (iter>1000.OR.j_eul<igno_min.OR.j_eul>igno_max) THEN
         852           0 :         WRITE(*,*) "search is which_eul_cell not converging!", iter, j_eul,igno_min,igno_max
         853           0 :         WRITE(*,*) "gno", gno(igno_min), gno(igno_max)
         854           0 :         write(*,*) gno
         855           0 :         STOP
         856             :       END IF
         857             :     END DO
         858      388800 :   END subroutine which_eul_cell
         859             :
         860             :
         861           0 :   function fuzzy(x,epsilon)
         862             :     implicit none
         863             :
         864             :     integer :: fuzzy
         865             :     real (kind=r8), intent(in) :: epsilon
         866             :     real (kind=r8) :: x
         867             :
         868           0 :     IF (ABS(x)<epsilon) THEN
         869             :       fuzzy = 0
         870           0 :     ELSE IF (x >epsilon) THEN
         871             :       fuzzy = 1
         872             :     ELSE !IF (x < fuzzy_width) THEN
         873           0 :       fuzzy = -1
         874             :     ENDIF
         875           0 :   end function
         876             :
         877             : end module fvm_overlap_mod

Generated by: LCOV version 1.14