Merge branch 'time-averaging' of github.com:chowland/AFiD-MuRPhFi int…

…o time-averaging

Makefile

@@ -28,15 +28,15 @@ endif
 
 ifeq ($(MACHINE),PC)
 # GNU Debug Flags
-	FC += -O0 -g -fbacktrace -Wall -Wextra
+	# FC += -O0 -g -fbacktrace -Wall -Wextra
 	# FC += -Wpedantic
 	# FC += -Warray-temporaries
 	# FC += -fcheck=all -finit-real=snan -ffpe-trap=invalid #-std=f2018
 # FC += -O0 -pg -fbacktrace -fbounds-check
 # Intel Debug Flags
 # FC += -O0 -g -traceback -check bounds
 	ifeq ($(FLAVOUR),GNU)
-		LDFLAGS = -L$(HOME)/fftw-install/lib -lfftw3 -llapack -ldl
+		LDFLAGS = -lfftw3 -llapack -ldl
 	else
 		LDFLAGS = -lfftw3 -qmkl=sequential
 	endif

src/flow_solver/CheckDivergence.F90

@@ -13,7 +13,7 @@ subroutine CheckDivergence(qmax,qmaxr)
     use local_arrays, only: vy,vx,vz
     use afid_salinity, only: vyr,vxr,vzr
     use mpih
-    use decomp_2d, only: xstart,xend,xstartr,xendr,nrank
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,nrank
     implicit none
     real,intent(out) :: qmax,qmaxr
     integer :: jc,kc,kp,jp,ic,ip

src/flow_solver/CreateGrid.F90

@@ -14,10 +14,7 @@ subroutine CreateGrid
     use GridModule
     implicit none
 
-    real :: x1,x2,x3
-
-    integer :: i, j, kc, km, kp
-    logical :: fexist
+    integer :: kc
 
     do kc=1,nxm
         kmv(kc)=kc-1

src/flow_solver/CreateInitialConditions.F90

@@ -16,11 +16,10 @@ subroutine CreateInitialConditions
     use afid_salinity, only: RayS
     use afid_phasefield, only: pf_eps, read_phase_field_params
     implicit none
-    integer :: j,k,i,kmid, io
+    integer :: j,k,i,kmid
     real :: xxx,yyy,zzz,eps,varptb,amp
     real :: h0,t0,Lambda,r, x0, A, B, alpha
     real, dimension(11) :: yh, zh
-    logical :: exists
 
     call random_seed()
 
@@ -37,15 +36,29 @@ subroutine CreateInitialConditions
     end if
 
     if (gAxis == 1) then
-        if ((RayT < 0) .and. (RayS <0)) then
-            !CJH: Stratified shear layer initial condition
-            do i=xstart(3),xend(3)
-                do j=xstart(2),xend(2)
-                    do k=1,nxm
-                        vy(k,j,i) = tanh(xm(k) - alx3/2.0)
+        if (RayT < 0) then
+            if (RayS <0) then
+                !CJH: Stratified shear layer initial condition
+                do i=xstart(3),xend(3)
+                    do j=xstart(2),xend(2)
+                        do k=1,nxm
+                            vy(k,j,i) = tanh(xm(k) - alx3/2.0)
+                            ! vz(k,j,i) = 1.0/cosh(xm(k) - alx3/2.0)
+                        end do
                     end do
                 end do
-            end do
+            else
+                !CJH: Salt-fingering initial condition
+                do i=xstart(3),xend(3)
+                    do j=xstart(2),xend(2)
+                        do k=1,nxm
+                            vx(k,j,i) = 0.0
+                            vy(k,j,i) = 0.0
+                            vz(k,j,i) = 0.0
+                        end do
+                    end do
+                end do
+            end if
         else
             !CJH: RBC initial condition as used in AFiD 1.0
             eps = 0.01
@@ -178,8 +191,8 @@ subroutine CreateInitialConditions
                 do k=1,nxm
                     call random_number(varptb)
                     if (abs(xm(k)-0.5) + eps > 0.5) then
-                    amp = 0.5 - abs(xm(k)-0.5) ! CJH Prevent values of |T| exceeding 0.5
-                    temp(k,j,i) = temp(k,j,i) + amp*(2.d0*varptb - 1.d0)
+                        amp = 0.5 - abs(xm(k)-0.5) ! CJH Prevent values of |T| exceeding 0.5
+                        temp(k,j,i) = temp(k,j,i) + amp*(2.d0*varptb - 1.d0)
                     else
                     temp(k,j,i) = temp(k,j,i) + eps*(2.d0*varptb - 1.d0)
                     end if

src/flow_solver/ExplicitTermsTemp.F90

@@ -14,7 +14,7 @@ subroutine ExplicitTermsTemp
     use decomp_2d, only: xstart,xend
     implicit none
     integer :: jc,kc,ic
-    integer :: km,kp,jm,jp,im,ip
+    integer :: jm,jp,im,ip
     real    :: htx,hty,htz,udy,udz
     real    :: udzq,udyq
     real    :: dyyt,dzzt

src/flow_solver/ImplicitAndUpdateTemp.F90

@@ -17,9 +17,7 @@ subroutine ImplicitAndUpdateTemp
     use ibm_param
     implicit none
     integer :: jc,kc,ic
-    integer :: km,kp,ke
     real    :: alpec,dxxt
-    real    :: app,acc,amm
 
     alpec=al/pect
 

src/flow_solver/ImplicitAndUpdateVX.F90

@@ -18,7 +18,7 @@ subroutine ImplicitAndUpdateVX
     use ibm_param
     implicit none
     integer :: jc,kc
-    integer :: km,kp,ic,ke
+    integer :: km,kp,ic
     real    :: alre,udx3
     real    :: amm,acc,app,dxp,dxxvx
 

src/flow_solver/ImplicitAndUpdateVY.F90

@@ -18,7 +18,7 @@ subroutine ImplicitAndUpdateVY
     use ibm_param
     implicit none
     integer :: kc,jmm,jc,ic
-    integer :: kpp,kmm,ke
+    integer :: kpp,kmm
     real    :: alre,udy
     real    :: amm,acc,app
     real    :: dyp,dxxvy

src/flow_solver/ImplicitAndUpdateVZ.F90

@@ -18,7 +18,7 @@ subroutine ImplicitAndUpdateVZ
     use ibm_param
     implicit none
     integer :: kc,jc,ic,imm
-    integer :: kmm,kpp,ke
+    integer :: kmm,kpp
     real    :: alre,amm,acc,app,udz
     real    :: dxxvz,dzp
 

src/flow_solver/ReadInputFile.F90

@@ -13,10 +13,9 @@ subroutine ReadInputFile
     use afid_phasefield, only: pf_A, pf_D, pf_eps, pf_Lambda, pf_S, pf_Tm
     implicit none
     character(len=4) :: dummy
-    integer flagstat,flagbal,flagmelt
+    integer :: flagmelt
     integer :: flagMR, flagsal, flagPF
-    integer :: FFscaleS!, pf_IBM
-    logical fexist
+    integer :: FFscaleS
 
     open(unit=15,file='bou.in',status='old')
         read(15,301) dummy

src/grid.F90

@@ -227,7 +227,7 @@ subroutine scallop_grid(c_grd, m_grd, Nm, grd_len, Retau, dw)
         real, intent(in) :: grd_len, Retau, dw
 
         integer :: k, ks
-        real :: alpha, kb, dyw, sig, dxlo, dxup, dxsmooth
+        real :: alpha, kb, sig, dxlo, dxup, dxsmooth
 
         ! Index of roughness height
         kb = 0.2*Retau/dw

src/h5tools/MakeMovieXCut.F90

@@ -10,13 +10,13 @@
 subroutine Mkmov_xcut
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
     use local_arrays, only: vz,vy,vx,temp
-    use afid_salinity, only: sal
-    use afid_phasefield, only: phi
+    use afid_salinity, only: sal, RayS
+    ! use afid_phasefield, only: phi
     use afid_moisture, only: humid
     use h5_tools
-    use param, only: nxm, nxmr, IBM
+    use param, only: nxm, nxmr, IBM, RayT
     implicit none
     character(70) :: filename
     character(4) :: varname
@@ -28,6 +28,7 @@ subroutine Mkmov_xcut
 
     ! Select plane - plane next to lower wall
     ic = 1
+    if (RayS < 0 .and. RayT < 0) ic = nxm/2
     if (IBM) ic = nxm/2
 
     ! Record filename as string
@@ -64,6 +65,7 @@ subroutine Mkmov_xcut
         !! Repeat on refined grid to save salinity
         ! Select wall plane index for refined grid
         ic = 1
+        if (RayS < 0 .and. RayT < 0) ic = nxmr/2
         if (IBM) ic = nxmr/2
 
         call h5_open_or_create(file_id, filename, comm, fexist)

src/h5tools/MakeMovieYCut.F90

@@ -11,7 +11,7 @@ subroutine Mkmov_ycut
     use param, only: nym, nymr
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
     use local_arrays, only: vz,vy,vx,temp
     use afid_salinity, only: sal
     use afid_phasefield, only: phi

src/h5tools/MakeMovieZCut.F90

@@ -11,11 +11,10 @@ subroutine Mkmov_zcut
     use param, only: nzm, nzmr
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
-    use local_arrays, only: vz,vy,vx,temp, pr
-    ! use mgrd_arrays, only: sal, phi, phic, tempr
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
+    use local_arrays, only: vz,vy,vx,temp!, pr
     use afid_salinity, only: sal
-    use afid_phasefield, only: phi, phic, tempr
+    use afid_phasefield, only: phi!, phic, tempr
     use afid_moisture, only: humid
     use h5_tools
     implicit none

src/h5tools/h5_tools.F90

@@ -71,7 +71,7 @@ subroutine write_H5_plane(file_id, varname, var, axis)
     character, intent(in) :: axis
 
     integer, parameter :: ndims = 2
-    integer(HID_T) :: filespace, slabspace, memspace, dset
+    integer(HID_T) :: filespace, memspace, dset
     integer(HSIZE_T), dimension(2) :: dims, data_count, data_offset
     character(len=5) :: frame
     character(len=10) :: dsetname

src/h5tools/mean_zplane.F90

@@ -3,7 +3,7 @@ subroutine mean_zplane
     use param, only: nx
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
     use local_arrays, only: vz,vy,vx,temp
     use afid_salinity, only: sal
     use afid_phasefield, only: phi
@@ -132,7 +132,7 @@ subroutine mean_yplane
     use param, only: nx
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
     use local_arrays, only: vz,vy,vx,temp
     use afid_salinity, only: sal
     use afid_phasefield, only: phi

src/ibm/IBMTools.F90

@@ -12,7 +12,7 @@ subroutine calc_interface_height(ph, h)
     !< input variable 
     real, intent(out) :: h(xstartr(2)-lvlhalo:,xstartr(3)-lvlhalo:)
     real :: dxx(nxmr)
-    integer :: i,j,k,kp
+    integer :: i,j,k
 
     h(:,:) = 0.0
 

src/main.F90

@@ -26,7 +26,7 @@ program AFiD
     integer :: prow=0,pcol=0
     integer :: lfactor,lfactor2
     character(100) :: arg
-    logical :: nanexist, write_mean_planes=.true.
+    logical :: write_mean_planes=.true.!, nanexist
     ! real,allocatable,dimension(:,:) :: dummy,dscan,dbot
     ! integer :: comm,ierror,row_id,row_coords(2),ic,jc,kc
 

src/moisture.F90

@@ -150,7 +150,7 @@ end subroutine SetHumidityBCs
 subroutine CreateInitialHumidity
     use mpih
     integer :: ic, jc, kc
-    real :: rnum, r, r2, amp
+    real :: rnum, r2, amp
     real :: bz(nxm), qz(nxm)
     logical :: exists
     character(len=30) :: dsetname, filename