From 233309b77fc45157eb6dc13d7ed868ee7c563ac6 Mon Sep 17 00:00:00 2001
From: s2006749 <s2006749@ed.ac.uk>
Date: Thu, 24 Oct 2024 13:10:20 +0100
Subject: [PATCH] added extra bit to deal with final part of ypraf in periodic
 case

---
 src/BC-Ellipsoid.f90 |  2 +-
 src/genepsi3d.f90    | 27 ++++++++++++++++++---------
 src/variables.f90    |  2 +-
 3 files changed, 20 insertions(+), 11 deletions(-)

diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index 0abe7501..474331e1 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -84,7 +84,7 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,dz,remp)
         if (ym /= ym) then
             write(*,*) "ym = ", ym, " should be ", ((real(j-1,mytype))*dy), ", as j = ", j,  " (or maybe it should be ", ((real(j-1,mytype))*dy)/real(nraf,mytype), ")"
             if (j.lt.(nyraf)) then 
-                write(*,*) "yp(j-1) = ", yp(j-1), " yp(j+1,..) = ", yp(j+1:)
+                write(*,*) "yp(j-1) = ", yp(j-1), " yp(j+1) = ", yp(j+1)
             endif
         endif
 
diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90
index 6bd7f4b8..2ecd5d40 100644
--- a/src/genepsi3d.f90
+++ b/src/genepsi3d.f90
@@ -55,11 +55,11 @@ contains
 
     IMPLICIT NONE
 
-    INTEGER :: nxi,nxf,ny,nyi,nyf,nzi,nzf
+    INTEGER, intent(in) :: nxi,nxf,ny,nyi,nyf,nzi,nzf
     REAL(mytype),DIMENSION(nxi:nxf,nyi:nyf,nzi:nzf) :: epsi
-    REAL(mytype)               :: dx,dz
-    REAL(mytype),DIMENSION(ny) :: yp
-    REAL(mytype)               :: remp
+    REAL(mytype), intent(in)               :: dx,dz
+    REAL(mytype),DIMENSION(ny), intent(in) :: yp
+    REAL(mytype), intent(in)               :: remp
 
     IF (itype.EQ.itype_cyl) THEN
 
@@ -276,7 +276,7 @@ contains
     real(mytype),dimension(ysize(1),nyraf,ysize(3))    :: yepsi
     real(mytype),dimension(zsize(1),zsize(2),nzraf)    :: zepsi
     real(mytype),dimension(ny)                         :: yp
-    real(mytype),dimension(nyraf+20)                      :: ypraf
+    real(mytype),dimension(nyraf)                      :: ypraf
     real(mytype)                     :: dxraf,dyraf,dzraf
     integer                          :: i,j,k
     integer                          :: ii,jj,kk
@@ -314,16 +314,25 @@ contains
     else
        dyraf =yly/real(nyraf-1, mytype)
     endif
+    write(*,*) ny, size(yp), size(ypraf), nraf
     do j=1,ny-1
        do jraf=1,nraf
           ypraf(jraf+nraf*(j-1))=yp(j)+real(jraf-1, mytype)*(yp(j+1)-yp(j))/real(nraf, mytype)
-          if (ypraf(jraf+nraf*(j-1)) /= ypraf(jraf+nraf*(j-1))) then 
-            write(*,*) "At j = ", j, ", jraf = ", jraf, "ypraf = ", yp(j)+real(jraf-1, mytype)*(yp(j+1)-yp(j))/real(nraf, mytype)
-          endif
+         !  if (ypraf(jraf+nraf*(j-1)) /= ypraf(jraf+nraf*(j-1))) then 
+         !    write(*,*) "At j = ", j, ", jraf = ", jraf, "ypraf = ", yp(j)+real(jraf-1, mytype)*(yp(j+1)-yp(j))/real(nraf, mytype)
+         !  endif
        enddo
     enddo
+    if (ncly) then 
+      do jraf = 1,nraf
+         ypraf(jraf+nraf*(ny-1))=yp(ny)+real(jraf-1,mytype)*(yly-yp(ny))/real(nraf,mytype)
+      enddo
+    endif
+   !  write(*,*) yp
     if(.not.ncly)ypraf(nyraf)=yp(ny)
-    if(.not.ncly)write(*,*) "Changed ypraf (", nyraf, "). To ", yp(ny)
+   !  if(.not.ncly)write(*,*) "Changed ypraf (", nyraf, "). To ", yp(ny)
+   !  write(*,*) ypraf
+
     yepsi=zero
     call geomcomplex(yepsi,ystart(1),yend(1),nyraf,1,nyraf,ystart(3),yend(3),dx,ypraf,dz,one)
     ! if (nrank==0) print*,'    step 3'
diff --git a/src/variables.f90 b/src/variables.f90
index f67446a4..bb391771 100644
--- a/src/variables.f90
+++ b/src/variables.f90
@@ -591,7 +591,7 @@ contains
        !complex_geometry
        nxraf=(nxm)*nraf
        if (.not.nclx) nxraf=nxraf+1
-       nyraf=(nym)*nraf+20
+       nyraf=(nym)*nraf
        if (.not.ncly) nyraf=nyraf+1
        nzraf=(nzm)*nraf
        if (.not.nclz) nzraf=nzraf+1
-- 
GitLab