Juan 10/03/2023:ZSOLVER/discretisation.f90, Cray OPENACC Opt, pass ZTAB by arg...

Juan 10/03/2023:ZSOLVER/discretisation.f90, Cray OPENACC Opt, pass ZTAB by arg + dims in apply_tridiag_solve_mnh_allT

src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90

@@ -2003,7 +2003,7 @@ end subroutine line_Jacobi_mnh
 
     real, dimension(:,:,:) , pointer, contiguous :: zSr_st , zb_st , zSu_in_st , zSu_out_st
     real, dimension(:)     , pointer, contiguous :: za_k, zb_k, zc_k, zd_k 
-    integer :: ii,ij
+!!$    integer :: ii,ij
 
     integer , parameter :: max_lev = 128
     logical , save                                 :: Lfirst_call_tridiag_mnhallT = .true.
@@ -2020,6 +2020,7 @@ end subroutine line_Jacobi_mnh
 
     real(kind=rl) :: Tijs
 
+    integer :: nib,nie,njb,nje,nzb,nze
     
     if (LUseT ) then
 
@@ -2036,6 +2037,10 @@ end subroutine line_Jacobi_mnh
       zc_k  => vert_coeff%c
       zd_k  => vert_coeff%d
 
+      nib = Lbound(zb_st,1) ; nie = Ubound(zb_st,1)
+      njb = Lbound(zb_st,2) ; nje = Ubound(zb_st,2)
+      nzb = Lbound(zb_st,3) ; nze = Ubound(zb_st,3)
+
       if ( Lfirst_call_level_tridiag_mnhallT(level) ) then
          Lfirst_call_level_tridiag_mnhallT(level) = .false.
 
@@ -2081,48 +2086,82 @@ end subroutine line_Jacobi_mnh
       
       tmp_k => Ttridiag_mnh(level)%tmp_k
       c_k => Ttridiag_mnh(level)%c_k
-         
-      ! acc kernels present(zSr_st,zSu_out_st)
-      !$acc parallel present(zSr_st,zSu_out_st)
-      !$mnh_do_concurrent ( ii=iib:iie , ij=ijb:ije )   
-         iz=1 
-         zSr_st(ii,ij,iz) = zb_st(ii,ij,iz)
-         !$acc loop seq
-         do iz=2,nz-1
-            zSr_st(ii,ij,iz) = zb_st(ii,ij,iz) - zd_k(iz) * ( &
-                 zSu_in_st(ii+1,ij,iz) + &
-                 zSu_in_st(ii-1,ij,iz) + &
-                 zSu_in_st(ii,ij+1,iz) + &
-                 zSu_in_st(ii,ij-1,iz) )
-         end do
-         iz=nz
-         zSr_st(ii,ij,iz) = zb_st(ii,ij,iz)
-      !
-      ! Thomas algorithm
-      !
-         iz = 1     
-         zSu_out_st(ii,ij,iz) = zSr_st(ii,ij,iz) / (tmp_k(iz)*Tijs*zd_k(iz))
-         !$acc loop seq
-         do iz=2,nz-1
-            zSu_out_st(ii,ij,iz) = (zSr_st(ii,ij,iz) / (Tijs*zd_k(iz)) & 
-                 - zSu_out_st(ii,ij,iz-1)*zc_k(iz)) / tmp_k(iz)
-         end do
-         iz=nz
-         zSu_out_st(ii,ij,iz) = (zSr_st(ii,ij,iz) / (Tijs*zd_k(iz)) &
-              -  zSu_out_st(ii,ij,iz-1)*zc_k(iz)) / tmp_k(iz)         
-      !     
-      ! STEP 2: back substitution
-         !$acc loop seq
-         do iz=nz-1,1,-1             
-            zSu_out_st(ii,ij,iz) = zSu_out_st(ii,ij,iz) & 
-                 - c_k(iz) * zSu_out_st(ii,ij,iz+1)
-         end do
-      !end do; end do
-      !$mnh_end_do()   
-      !$acc end parallel
-      ! acc end kernels
+
+!!$      print*,"apply_tridiag::level,Sr%ix_min/%ix_max/%icompx_min/%icompx_max",&
+!!$           Sr%ix_min,Sr%ix_max,Sr%icompx_min,Sr%icompx_max
+!!$      print*,"apply_tridiag::level,iib,iie,ijb,ije,nz",level,iib,iie,ijb,ije,nz
+!!$      print*,"apply_tridiag::level,nib,nie,njb,nje,nzb,nze",level,&
+!!$                                   nib,nie,njb,nje,nzb,nze
+!!$      print*,"apply_tridiag::level, zSr_st=L/U/S",level,Lbound(zSr_st),Ubound(zSr_st),Shape(zSr_st)
+!!$      print*,"apply_tridiag::lvel, zc_k  =L/U/S",level,Lbound(zc_k),Ubound(zc_k),Shape(zc_k)
+
+      call  apply_tridiag_solve_mnh_allT_dim(&
+           zSr_st,zSu_out_st,zb_st,zSu_in_st,&
+           zc_k,zd_k,tmp_k,c_k &
+      )
        
-    end if
+      end if
+
+    contains
+      
+      subroutine  apply_tridiag_solve_mnh_allT_dim(&
+           pzSr_st,pzSu_out_st,pzb_st,pzSu_in_st,&
+           pzc_k,pzd_k,ptmp_k,pc_k &
+        )
+
+        implicit none
+
+        real :: pzSr_st    (nib:nie,njb:nje,nzb:nze) ,&
+                pzSu_out_st(nib:nie,njb:nje,nzb:nze) ,&
+                pzb_st     (nib:nie,njb:nje,nzb:nze) ,&
+                pzSu_in_st (nib:nie,njb:nje,nzb:nze) ,&
+                pzc_k(nz),pzd_k(nz),ptmp_k(nz),pc_k(nz)
+        !
+        ! local var
+        !
+        integer :: ii,ij,iz
+                   
+        ! acc kernels present(pzSr_st,pzSu_out_st,pzb_st,pzSu_in_st)
+        !$acc parallel present(pzSr_st,pzSu_out_st,pzb_st,pzSu_in_st)
+        !$acc &        present(pzc_k,pzd_k,ptmp_k,pc_k)
+        !$mnh_do_concurrent ( ii=iib:iie , ij=ijb:ije )   
+        iz=1 
+        pzSr_st(ii,ij,iz) = pzb_st(ii,ij,iz)
+        !$acc loop seq
+        do iz=2,nz-1
+           pzSr_st(ii,ij,iz) = pzb_st(ii,ij,iz) - pzd_k(iz) * ( &
+                pzSu_in_st(ii+1,ij,iz) + &
+                pzSu_in_st(ii-1,ij,iz) + &
+                pzSu_in_st(ii,ij+1,iz) + &
+                pzSu_in_st(ii,ij-1,iz) )
+        end do
+        iz=nz
+        pzSr_st(ii,ij,iz) = pzb_st(ii,ij,iz)
+        !
+        ! Thomas algorithm
+        !
+        iz = 1     
+        pzSu_out_st(ii,ij,iz) = pzSr_st(ii,ij,iz) / (ptmp_k(iz)*Tijs*pzd_k(iz))
+        !$acc loop seq
+        do iz=2,nz-1
+           pzSu_out_st(ii,ij,iz) = (pzSr_st(ii,ij,iz) / (Tijs*pzd_k(iz)) & 
+                - pzSu_out_st(ii,ij,iz-1)*pzc_k(iz)) / ptmp_k(iz)
+        end do
+        iz=nz
+        pzSu_out_st(ii,ij,iz) = (pzSr_st(ii,ij,iz) / (Tijs*pzd_k(iz)) &
+             -  pzSu_out_st(ii,ij,iz-1)*pzc_k(iz)) / ptmp_k(iz)         
+        !     
+        ! STEP 2: back substitution
+        !$acc loop seq
+        do iz=nz-1,1,-1             
+           pzSu_out_st(ii,ij,iz) = pzSu_out_st(ii,ij,iz) & 
+                - pc_k(iz) * pzSu_out_st(ii,ij,iz+1)
+        end do
+        !$mnh_end_do()   
+        !$acc end parallel
+        ! acc end kernels
+         
+      end subroutine apply_tridiag_solve_mnh_allT_dim
 
   end subroutine apply_tridiag_solve_mnh_allT
   !==================================================================