diff --git a/Makefile b/Makefile index 48550c73..cdd741e6 100644 --- a/Makefile +++ b/Makefile @@ -1,12 +1,17 @@ # Choose the machine being used -# Options: PC_GNU, PC_INTEL, (i)SNELLIUS, IRENE(_SKL/_ROME), MARENOSTRUM, SUPERMUC -MACHINE=PC_GNU +# Options: PC, SNELLIUS, IRENE, MARENOSTRUM, SUPERMUC, DISCOVERER +MACHINE=PC +FLAVOUR=GNU # Modules required for each HPC system as follows: -# SNELLIUS: 2022 foss/2022a HDF5/1.12.2-gompi-2022a -# iSNELLIUS: 2022 intel/2022a FFTW/3.3.10-GCC-11.3.0 HDF5/1.12.2-iimpi-2021a -# IRENE: flavor/hdf5/parallel hdf5 fftw3/gnu -# MARENOSTRUM: fabric intel mkl impi hdf5 fftw szip -# SUPERMUC: fftw hdf5 +# SNELLIUS: +# GNU: 2022 foss/2022a HDF5/1.12.2-gompi-2022a +# Intel: 2022 intel/2022a FFTW/3.3.10-GCC-11.3.0 HDF5/1.12.2-iimpi-2021a +# IRENE (Intel): flavor/hdf5/parallel hdf5 fftw3/gnu +# MARENOSTRUM (Intel): fabric intel mkl impi hdf5 fftw szip +# SUPERMUC (Intel): fftw hdf5 +# DISCOVERER: +# GNU: hdf5/1/1.14/latest-gcc-openmpi fftw/3/latest-gcc-openmpi lapack +# Intel: hdf5/1/1.14/latest-intel-openmpi fftw/3/latest-gcc-openmpi mkl #======================================================================= # Compiler options @@ -15,39 +20,49 @@ MACHINE=PC_GNU # Object and module directory: OBJDIR=obj -ifeq ($(MACHINE),PC_GNU) - FC = h5pfc -cpp -fdefault-real-8 -fdefault-double-8 -O3 - # FC += -pg -fbacktrace -fbounds-check - LDFLAGS = -lfftw3 -llapack -lblas -ldl +ifeq ($(FLAVOUR),GNU) + FC = h5pfc -cpp -fdefault-real-8 -fdefault-double-8 -fallow-argument-mismatch +else + FC = h5pfc -fpp -r8 endif -ifeq ($(MACHINE),PC_INTEL) - FC = h5pfc -fpp -r8 -O3 -## Traceback / Debug -# FC += -r8 -O0 -g -traceback -check bounds -# FC += -DSHM -DSHM_DEBUG - LDFLAGS = -lfftw3 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lhdf5_fortran -lhdf5 -lsz -lz -ldl -lm + +ifeq ($(MACHINE),PC) +# GNU Debug Flags + # 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 + else + LDFLAGS = -lfftw3 -qmkl=sequential + endif endif -ifeq ($(MACHINE),iSNELLIUS) - FC = h5pfc -fpp -r8 -O3 -align array64byte -fma -ftz -fomit-frame-pointer - LDFLAGS = -lfftw3 -qmkl=sequential +ifeq ($(MACHINE),DISCOVERER) + ifeq ($(FLAVOUR),GNU) + LDFLAGS += -lfftw3 -llapack -ldl + else + LDFLAGS += -lfftw3 -qmkl=sequential + endif endif ifeq ($(MACHINE),SNELLIUS) - FC = h5pfc -cpp -fdefault-real-8 -fdefault-double-8 -w -fallow-argument-mismatch - FC += -O2 -march=znver1 -mtune=znver1 -mfma -mavx2 -m3dnow -fomit-frame-pointer - BLAS_LIBS = -lscalapack -lopenblas -ldl - LDFLAGS = -lfftw3 $(BLAS_LIBS) + ifeq ($(FLAVOUR),GNU) + FC += -O2 -march=znver1 -mtune=znver1 -mfma -mavx2 -m3dnow -fomit-frame-pointer + LDFLAGS = -lfftw3 -lopenblas -ldl + else + FC += -align array64byte -fma -ftz -fomit-frame-pointer + LDFLAGS = -lfftw3 -qmkl=sequential + endif endif -ifeq ($(MACHINE),IRENE_SKL) - FC = h5pfc -fpp -r8 -O3 -mtune=skylake -xCORE-AVX512 -m64 -fPIC $(FFTW3_FFLAGS) +ifeq ($(MACHINE),IRENE) + FC += -mtune=skylake -xCORE-AVX512 -m64 -fPIC $(FFTW3_FFLAGS) LDFLAGS = $(FFTW3_LDFLAGS) $(MKL_LDFLAGS) -ldl endif -ifeq ($(MACHINE),IRENE_ROME) - FC = h5pfc -fpp -r8 -O3 -mavx2 $(FFTW3_FFLAGS) - HDF5_LIBS = -lhdf5_fortran -lhdf5 -lz -ldl -lm - LDFLAGS = $(FFTW3_LDFLAGS) $(MKL_LDFLAGS) $(HDF5_LIBS) -endif ifeq ($(MACHINE),MARENOSTRUM) - FC = h5pfc -fpp -r8 -O3 -mtune=skylake -xCORE-AVX512 -m64 -fPIC $(FFTW_FFLAGS) + FC += -mtune=skylake -xCORE-AVX512 -m64 -fPIC $(FFTW_FFLAGS) LDFLAGS = $(FFTW_LIBS) -mkl=sequential endif ifeq ($(MACHINE),SUPERMUC) @@ -55,14 +70,10 @@ ifeq ($(MACHINE),SUPERMUC) LDFLAGS = $(FFTW_LIB) $(HDF5_F90_SHLIB) $(HDF5_SHLIB) -qmkl=sequential endif -ifeq ($(MACHINE),SNELLIUS) +ifeq ($(FLAVOUR),GNU) FC += -J $(OBJDIR) else - ifeq ($(MACHINE),PC_GNU) - FC += -J $(OBJDIR) - else - FC += -module $(OBJDIR) - endif + FC += -module $(OBJDIR) endif #======================================================================= @@ -73,42 +84,31 @@ EXTRA_DIST = transpose_z_to_x.F90 transpose_x_to_z.F90 transpose_x_to_y.F90\ factor.F90 halo.F90 fft_common.F90 alloc.F90 halo_common.F90 # Object files associated with standard flow solver -OBJS = obj/main.o obj/CalcLocalDivergence.o obj/CalcMaxCFL.o \ - obj/CalcMeanProfiles.o obj/CheckDivergence.o obj/CorrectPressure.o \ +OBJS = obj/main.o obj/CalcMaxCFL.o \ + obj/CalcMeanProfiles.o obj/CheckDivergence.o \ obj/CorrectVelocity.o obj/CreateGrid.o obj/CreateInitialConditions.o \ obj/DeallocateVariables.o obj/DebugRoutines.o obj/ExplicitTermsTemp.o \ obj/ExplicitTermsVX.o obj/ExplicitTermsVY.o obj/ExplicitTermsVZ.o \ obj/factorize.o obj/HdfReadContinua.o obj/HdfRoutines.o \ obj/ImplicitAndUpdateTemp.o obj/ImplicitAndUpdateVX.o obj/ImplicitAndUpdateVY.o \ - obj/ImplicitAndUpdateVZ.o obj/InitPressureSolver.o obj/InitTimeMarchScheme.o \ + obj/ImplicitAndUpdateVZ.o obj/InitTimeMarchScheme.o \ obj/InitVariables.o obj/LocateLargeDivergence.o obj/MakeMovieXCut.o \ obj/MakeMovieYCut.o obj/MakeMovieZCut.o obj/MpiAuxRoutines.o \ obj/QuitRoutine.o obj/ReadInputFile.o obj/ResetLogs.o \ obj/SetTempBCs.o obj/SolveImpEqnUpdate_Temp.o obj/SolveImpEqnUpdate_X.o \ - obj/SolveImpEqnUpdate_YZ.o obj/SolvePressureCorrection.o obj/SpecRoutines.o \ + obj/SolveImpEqnUpdate_YZ.o obj/SpecRoutines.o \ obj/TimeMarcher.o obj/WriteFlowField.o obj/WriteGridInfo.o \ obj/CalcWriteQ.o obj/GlobalQuantities.o obj/ReadFlowInterp.o # Object files associated with multiple resolution grids OBJS += obj/CreateMgrdGrid.o obj/InitMgrdVariables.o \ - obj/DeallocateMgrdVariables.o obj/CreateMgrdStencil.o# obj/CreateMgrdStencil.o + obj/DeallocateMgrdVariables.o obj/CreateMgrdStencil.o # Object files associated with initial condition interpolation OBJS += obj/CreateNewInputStencil.o obj/CreateOldGrid.o obj/CreateNewSalStencil.o \ - obj/InterpInputSal.o obj/InterpInputVel.o obj/InterpSalMgrd.o \ + obj/InterpInputSal.o obj/InterpInputVel.o \ obj/InterpVelMgrd.o obj/InitInputVars.o obj/DeallocateInputVars.o \ - obj/InterpInputPhi.o# obj/CreateInputStencil.o obj/CreateSalStencil.o - -# Object files associated with the salinity field -OBJS += obj/ExplicitTermsSal.o obj/ImplicitAndUpdateSal.o obj/SolveImpEqnUpdate_Sal.o \ - obj/UpdateScalarBCs.o obj/CreateICSal.o obj/InitSalVariables.o \ - obj/DeallocateSalVariables.o obj/SetSalBCs.o - -# Object files associated with the phase-field method -OBJS += obj/AddLatentHeat.o obj/DeallocatePFVariables.o obj/ExplicitTermsPhi.o \ - obj/ImplicitAndUpdatePhi.o obj/InitPFVariables.o obj/InterpPhiMgrd.o \ - obj/InterpTempMgrd.o obj/SolveImpEqnUpdate_Phi.o obj/CreateICPF.o \ - obj/ImmersedBoundary.o + obj/InterpInputPhi.o # # Object files associated with the immersed boundary method OBJS += obj/SolveImpEqnUpdate_Temp_ibm.o obj/SolveImpEqnUpdate_X_ibm.o \ @@ -120,13 +120,15 @@ OBJS += obj/mean_zplane.o # Module object files MOBJS = obj/param.o obj/decomp_2d.o obj/AuxiliaryRoutines.o obj/decomp_2d_fft.o \ - obj/HermiteInterpolations.o obj/GridModule.o obj/h5_tools.o obj/means.o obj/ibm_param.o + obj/pressure.o obj/HermiteInterpolations.o obj/grid.o obj/h5_tools.o obj/means.o \ + obj/ibm_param.o obj/IBMTools.o obj/moisture.o obj/salinity.o obj/phasefield.o #======================================================================= # Files that create modules: #======================================================================= MFILES = param.F90 decomp_2d.F90 AuxiliaryRoutines.F90 decomp_2d_fft.F90 \ - HermiteInterpolations.F90 GridModule.F90 ibm_param.F90 + pressure.F90 HermiteInterpolations.F90 grid.F90 ibm_param.F90 IBMTools.F90 \ + moisture.F90 salinity.F90 phasefield.F90 #============================================================================ # make PROGRAM @@ -151,14 +153,24 @@ $(OBJDIR)/decomp_2d_fft.o: src/flow_solver/2decomp/decomp_2d_fft.F90 $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/ibm_param.o: src/ibm/ibm_param.F90 $(FC) -c -o $@ $< $(LDFLAGS) -$(OBJDIR)/GridModule.o: src/flow_solver/GridModule.F90 +$(OBJDIR)/grid.o: src/grid.F90 $(FC) -c -o $@ $< $(LDFLAGS) -$(OBJDIR)/HermiteInterpolations.o: src/multires/HermiteInterpolations.F90 +$(OBJDIR)/pressure.o: src/pressure.F90 + $(FC) -c -o $@ $< $(LDFLAGS) +$(OBJDIR)/HermiteInterpolations.o: src/multires/HermiteInterpolations.F90 obj/ibm_param.o $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/h5_tools.o: src/h5tools/h5_tools.F90 $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/means.o: src/h5tools/means.F90 obj/ibm_param.o $(FC) -c -o $@ $< $(LDFLAGS) +$(OBJDIR)/IBMTools.o: src/ibm/IBMTools.F90 + $(FC) -c -o $@ $< $(LDFLAGS) +$(OBJDIR)/salinity.o: src/salinity.F90 + $(FC) -c -o $@ $< $(LDFLAGS) +$(OBJDIR)/phasefield.o: src/phasefield.F90 obj/salinity.o + $(FC) -c -o $@ $< $(LDFLAGS) +$(OBJDIR)/moisture.o: src/moisture.F90 + $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/%.o: src/%.F90 $(MOBJS) $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/%.o: src/flow_solver/%.F90 $(MOBJS) @@ -169,10 +181,6 @@ $(OBJDIR)/%.o: src/multires/%.F90 $(MOBJS) $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/%.o: src/multires/IC_interpolation/%.F90 $(MOBJS) $(FC) -c -o $@ $< $(LDFLAGS) -$(OBJDIR)/%.o: src/multires/phase-field/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) -$(OBJDIR)/%.o: src/multires/salinity/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) $(OBJDIR)/%.o: src/ibm/%.F90 $(MOBJS) $(FC) -c -o $@ $< $(LDFLAGS) diff --git a/examples/RainyBenard/bou.in b/examples/RainyBenard/bou.in new file mode 100644 index 00000000..8e5f3443 --- /dev/null +++ b/examples/RainyBenard/bou.in @@ -0,0 +1,55 @@ +Base grid resolution / number of substeps +NXM NYM NZM NSST(>1=RK,else AB) +96 1024 1 3 + +Refined grid: on/off and grid size +MULTIRES NXMR NYMR NZMR +0 128 128 1 + +Salinity and Phase-field (on/off) +FLAGSAL FLAGPF +0 0 + +Read flow / Reset logs +NREAD IRESET +0 0 + +Time steps / time limit (wall/flow) +NTST WALLTIMEMAX TMAX +1000000 120 2500.0 + +Time interval for saving stats/movie frames +TOUT TFRAME SAVE_3D +1.0 5.0 100.0 + +Domain size (keep ALX3 = 1.0) +ALX3 YLEN ZLEN +1.0 20.0 0.01 + +Grid stretching parameters +ISTR3 STR3 ISTR3R +6 12.0 3 + +Physical flow parameters / Free-fall timescale +RAYT PRAT RAYS PRAS FFscaleS +2e5 1.0 1e7 1.0 0 + +Time stepping parameters +IDTV DT RESID CFLMAX DTMIN DTMAX +1 1e-1 1.e-3 1.0 1e-8 1e-1 + +Dirichlet/Neumann BC flag on upper(N) and lower(S) walls +inslwS inslwN TfixS TfixN SfixS SfixN +1 1 1 1 1 1 + +Active/passive scalar flags / Gravitational axis (x=1,y=2,z=3) +active_T active_S gAxis +1 0 1 + +Wall velocity / Pressure gradient / Melt boundary condition +xplusU xminusU dPdy dPdz MELT +0.0 0.0 0.0 0.0 0 + +Phase-field parameters +pf_D(*Pe) pf_A pf_S pf_Tm IBM pf_IC +0.1 10.0 10.0 0.0 0 1 diff --git a/examples/RainyBenard/humid.in b/examples/RainyBenard/humid.in new file mode 100644 index 00000000..074da1a6 --- /dev/null +++ b/examples/RainyBenard/humid.in @@ -0,0 +1,24 @@ +An input file for parameters used in moist convection +--- + +alpha (Dimensionless saturation parameter) +3.0 + +beta (Dimensionless lapse rate) +1.2 + +gamma (Dimensioness latent heat parameter) +0.19 + +tau (Condensation time scale) +0.375 + +Diffusive scale for tau? (.true. / .false.) +.false. + +Fixed humidity boundary conditions (1 / 0) +top bottom +1 1 + +Sm (Humid-to-thermal diffusivity ratio) +1.0 \ No newline at end of file diff --git a/src/ReadFlowInterp.F90 b/src/ReadFlowInterp.F90 index 0c3065c0..b83a5b20 100644 --- a/src/ReadFlowInterp.F90 +++ b/src/ReadFlowInterp.F90 @@ -4,7 +4,9 @@ subroutine ReadFlowInterp(prow,pcol) use local_arrays use param use input_grids - use mgrd_arrays + use afid_salinity + use afid_phasefield + use afid_moisture, only: humid use AuxiliaryRoutines implicit none @@ -13,8 +15,6 @@ subroutine ReadFlowInterp(prow,pcol) character*70 :: filnam,dsetname - integer :: i, j, k, ic, jc, kc - real :: yleno, zleno logical :: fexist @@ -160,7 +160,11 @@ subroutine ReadFlowInterp(prow,pcol) inquire(file=filnam, exist=fexist) if (fexist) then call HdfReadContinua(nz, ny, nx, xstart(2), xend(2), & - xstart(3), xend(3), 1, pr) + xstart(3), xend(3), 7, pr) + end if + if (moist) then + call HdfReadContinua(nz, ny, nx, xstart(2), xend(2), & + xstart(3), xend(3), 8, humid) end if end if diff --git a/src/flow_solver/CalcMaxCFL.F90 b/src/flow_solver/CalcMaxCFL.F90 index 5cb122ae..a3c3108f 100644 --- a/src/flow_solver/CalcMaxCFL.F90 +++ b/src/flow_solver/CalcMaxCFL.F90 @@ -8,66 +8,66 @@ ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine CalcMaxCFL(cflm,cflmr) - use param - use local_arrays, only: vx,vy,vz - use mgrd_arrays, only: vxr,vyr,vzr - use decomp_2d - use mpih - implicit none - real,intent(out) :: cflm, cflmr - integer :: j,k,jp,kp,i,ip - real :: qcf - - cflm=0.00000001d0 -! -!$OMP PARALLEL DO & -!$OMP DEFAULT(none) & -!$OMP SHARED(xstart,xend,nxm,vz,vy,vx) & -!$OMP SHARED(dz,dy,udx3m) & -!$OMP PRIVATE(i,j,k,ip,jp,kp,qcf) & -!$OMP REDUCTION(max:cflm) - do i=xstart(3),xend(3) +subroutine CalcMaxCFL(cflm,cflmr) + use param + use local_arrays, only: vx,vy,vz + use afid_salinity, only: vxr,vyr,vzr + use decomp_2d + use mpih + implicit none + real,intent(out) :: cflm, cflmr + integer :: j,k,jp,kp,i,ip + real :: qcf + + cflm=0.00000001d0 + ! + !$OMP PARALLEL DO & + !$OMP DEFAULT(none) & + !$OMP SHARED(xstart,xend,nxm,vz,vy,vx) & + !$OMP SHARED(dz,dy,udx3m) & + !$OMP PRIVATE(i,j,k,ip,jp,kp,qcf) & + !$OMP REDUCTION(max:cflm) + do i=xstart(3),xend(3) ip=i+1 do j=xstart(2),xend(2) - jp=j+1 - do k=1,nxm - kp=k+1 - qcf=( abs((vz(k,j,i)+vz(k,j,ip))*0.5d0*dz) & - +abs((vy(k,j,i)+vy(k,jp,i))*0.5d0*dy) & - +abs((vx(k,j,i)+vx(kp,j,i))*0.5d0*udx3m(k))) - - cflm = max(cflm,qcf) - enddo - enddo - enddo -!$OMP END PARALLEL DO - - call MpiAllMaxRealScalar(cflm) - - if (salinity) then + jp=j+1 + do k=1,nxm + kp=k+1 + qcf=( abs((vz(k,j,i)+vz(k,j,ip))*0.5d0*dz) & + +abs((vy(k,j,i)+vy(k,jp,i))*0.5d0*dy) & + +abs((vx(k,j,i)+vx(kp,j,i))*0.5d0*udx3m(k))) + + cflm = max(cflm,qcf) + end do + end do + end do + !$OMP END PARALLEL DO + + call MpiAllMaxRealScalar(cflm) + + if (salinity) then ! Refined mesh - + cflmr = 1.d-8 do i=xstartr(3),xendr(3) - ip = i + 1 - do j=xstartr(2),xendr(2) - jp = j + 1 - do k=1,nxmr - kp = k + 1 - qcf=( abs((vzr(k,j,i) + vzr(k,j,ip))*0.5d0*dzr) & - +abs((vyr(k,j,i) + vyr(k,jp,i))*0.5d0*dyr) & - +abs((vxr(k,j,i) + vxr(kp,j,i))*0.5d0*udx3mr(k))) - - cflmr = max(cflmr,qcf) + ip = i + 1 + do j=xstartr(2),xendr(2) + jp = j + 1 + do k=1,nxmr + kp = k + 1 + qcf=( abs((vzr(k,j,i) + vzr(k,j,ip))*0.5d0*dzr) & + +abs((vyr(k,j,i) + vyr(k,jp,i))*0.5d0*dyr) & + +abs((vxr(k,j,i) + vxr(kp,j,i))*0.5d0*udx3mr(k))) + + cflmr = max(cflmr,qcf) + end do end do - end do end do - + call MpiAllMaxRealScalar(cflmr) - else + else cflmr = cflm - end if - - return - end + end if + + return +end \ No newline at end of file diff --git a/src/flow_solver/CalcMeanProfiles.F90 b/src/flow_solver/CalcMeanProfiles.F90 index 7ef1245c..b39cb6e0 100644 --- a/src/flow_solver/CalcMeanProfiles.F90 +++ b/src/flow_solver/CalcMeanProfiles.F90 @@ -12,25 +12,24 @@ subroutine CalcMeanProfiles use param use local_arrays, only: vx,vy,vz,temp - use mgrd_arrays, only: vxr,vyr,vzr,sal,phi use mpih - use decomp_2d, only: xstart,xend,xstartr,xendr - use ibm_param, only: solidr + use decomp_2d, only: xstart,xend + use afid_salinity, only: CalcSalStats, CreateSalinityH5Groups + use afid_moisture, only: CalcMoistStats, CreateMoistH5Groups + use afid_phasefield, only: CalcPhiStats, CreatePhaseH5Groups implicit none integer :: i, j, k, im, ip, jm, jp, km, kp - real :: inym, inzm, inymr, inzmr, tdx, tdxr + real :: inym, inzm, tdx real :: tii(2) - real, dimension(nxm) :: Tbar, Trms, chiT, chiT2 + real, dimension(nxm) :: Tbar, Trms, chiT!, chiT2 real, dimension(nxm) :: vxT, vyT, vzT real, dimension(nxm) :: vxrms, vyrms, vzrms real, dimension(nxm) :: vybar, vzbar real, dimension(nxm) :: epsilon real, dimension(nxm) :: vxvy, vxvz, vyvz - real, dimension(nxmr) :: Sbar, Srms, chiS - real, dimension(nxmr) :: vxS, vyS, vzS - real, dimension(nxmr) :: phibar, phirms + character(5) :: nstat character(30) :: dsetname,filename logical :: fexist @@ -43,11 +42,6 @@ subroutine CalcMeanProfiles Trms(:) =0.0; vxvy(:) =0.0; vxvz(:) =0.0 chiT(:) =0.0; epsilon(:)=0.0; vyvz(:) =0.0 - Sbar(:)=0.0; Srms(:)=0.0; chiS(:)=0.0 - vxS(:) =0.0; vyS(:) =0.0; vzS(:) =0.0 - - phibar(:) = 0.0; phirms(:) = 0.0; - inym = 1.d0/nym inzm = 1.d0/nzm @@ -212,167 +206,6 @@ subroutine CalcMeanProfiles epsilon(k) = epsilon(k)*inym*inzm/ren end do - if (salinity) then - inymr = 1.d0/nymr - inzmr = 1.d0/nzmr - - if (IBM) then - do i=xstartr(3),xendr(3) - do j=xstartr(2),xendr(2) - do k=1,nxmr - ! Only record data from fluid phase - if (.not. solidr(k,j,i)) then - Sbar(k) = Sbar(k) + sal(k,j,i)*inymr*inzmr - end if - ! NB This is zero in solid anyway, so no need for an if - vxS(k) = vxS(k) + 0.5*(vxr(k,j,i)+vxr(k+1,j,i))*sal(k,j,i)*inymr*inzmr - vyS(k) = vyS(k) + 0.5*(vyr(k,j,i)+vyr(k,j+1,i))*sal(k,j,i)*inymr*inzmr - vzS(k) = vzS(k) + 0.5*(vzr(k,j,i)+vzr(k,j,i+1))*sal(k,j,i)*inymr*inzmr - end do - end do - end do - else - do i=xstartr(3),xendr(3) - do j=xstartr(2),xendr(2) - do k=1,nxmr - Sbar(k) = Sbar(k) + sal(k,j,i)*inymr*inzmr - - vxS(k) = vxS(k) + 0.5*(vxr(k,j,i)+vxr(k+1,j,i))*sal(k,j,i)*inymr*inzmr - vyS(k) = vyS(k) + 0.5*(vyr(k,j,i)+vyr(k,j+1,i))*sal(k,j,i)*inymr*inzmr - vzS(k) = vzS(k) + 0.5*(vzr(k,j,i)+vzr(k,j,i+1))*sal(k,j,i)*inymr*inzmr - end do - end do - end do - end if - - if (IBM) then - do i=xstartr(3),xendr(3) - do j=xstartr(2),xendr(2) - do k=1,nxmr - if (.not. solidr(k,j,i)) then - Srms(k) = Srms(k) + sal(k,j,i)**2*inymr*inzmr - end if - end do - end do - end do - else - do i=xstartr(3),xendr(3) - do j=xstartr(2),xendr(2) - do k=1,nxmr - Srms(k) = Srms(k) + sal(k,j,i)**2*inymr*inzmr - end do - end do - end do - end if - - call MpiSumReal1D(Sbar,nxmr) - call MpiSumReal1D(vxS,nxmr) - call MpiSumReal1D(vyS,nxmr) - call MpiSumReal1D(vzS,nxmr) - call MpiSumReal1D(Srms,nxmr) - - do k=1,nxmr - Srms(k) = sqrt(Srms(k)) - end do - - if (IBM) then - do i=xstartr(3),xendr(3) - ip = i + 1 - im = i - 1 - do j=xstartr(2),xendr(2) - jp = j + 1 - jm = j - 1 - do k=1,nxmr - if (.not. solidr(k,j,i)) then - if (solidr(k,j,ip)) then - chiS(k) = chiS(k) + ((sal(k,j,i)-sal(k,j,im))*0.5*dzr)**2*inymr*inzmr - elseif (solidr(k,j,im)) then - chiS(k) = chiS(k) + ((sal(k,j,ip)-sal(k,j,i))*0.5*dzr)**2*inymr*inzmr - else - chiS(k) = chiS(k) + ((sal(k,j,ip)-sal(k,j,im))*0.5*dzr)**2*inymr*inzmr - end if - if (solidr(k,jp,i)) then - chiS(k) = chiS(k) + ((sal(k,j,i)-sal(k,jm,i))*0.5*dyr)**2*inymr*inzmr - elseif (solidr(k,jm,i)) then - chiS(k) = chiS(k) + ((sal(k,jp,i)-sal(k,j,i))*0.5*dyr)**2*inymr*inzmr - else - chiS(k) = chiS(k) + ((sal(k,jp,i)-sal(k,jm,i))*0.5*dyr)**2*inymr*inzmr - end if - end if - end do - if (.not. solidr(1,j,i)) then - tdxr = 0.5*dxr/g3rmr(1) - chiS(1) = chiS(1) + ((sal(2,j,i)-sal(1,j,i)+2.0*SfixS*(sal(1,j,i)-salbp(1,j,i)))& - & *tdxr)**2*inymr*inzmr - end if - do k=2,nxmr-1 - if (.not. solidr(k,j,i)) then - tdxr = 0.5*dxr/g3rmr(k) - chiS(k) = chiS(k) + ((sal(k+1,j,i)-sal(k-1,j,i)) * tdxr)**2*inymr*inzmr - end if - end do - if (.not. solidr(nxmr,j,i)) then - tdxr = 0.5*dxr/g3rmr(nxmr) - chiS(nxmr) = chiS(nxmr) + ((sal(nxmr,j,i)-sal(nxmr-1,j,i)+2.0*SfixN*(saltp(1,j,i)-sal(nxmr,j,i)))& - & *tdxr)**2*inymr*inzmr - end if - end do - end do - else - do i=xstartr(3),xendr(3) - ip = i + 1 - im = i - 1 - do j=xstartr(2),xendr(2) - jp = j + 1 - jm = j - 1 - do k=1,nxmr - chiS(k) = chiS(k) + ((sal(k,j,ip)-sal(k,j,im))*0.5*dzr)**2*inymr*inzmr - chiS(k) = chiS(k) + ((sal(k,jp,i)-sal(k,jm,i))*0.5*dyr)**2*inymr*inzmr - end do - tdxr = 0.5*dxr/g3rmr(1) - chiS(1) = chiS(1) + ((sal(2,j,i)-sal(1,j,i)+2.0*SfixS*(sal(1,j,i)-salbp(1,j,i)))& - & *tdxr)**2*inymr*inzmr - do k=2,nxmr-1 - tdxr = 0.5*dxr/g3rmr(k) - chiS(k) = chiS(k) + ((sal(k+1,j,i)-sal(k-1,j,i)) * tdxr)**2*inymr*inzmr - end do - tdxr = 0.5*dxr/g3rmr(nxmr) - chiS(nxmr) = chiS(nxmr) + ((sal(nxmr,j,i)-sal(nxmr-1,j,i)+2.0*SfixN*(saltp(1,j,i)-sal(nxmr,j,i)))& - & *tdxr)**2*inymr*inzmr - end do - end do - end if - - - call MpiSumReal1D(chiS,nxmr) - - do k=1,nxmr - chiS(k) = chiS(k)/pecs - end do - end if - - if (phasefield) then - inymr = 1.d0/nymr - inzmr = 1.d0/nzmr - - do i=xstartr(3),xendr(3) - do j=xstartr(2),xendr(2) - do k=1,nxmr - phibar(k) = phibar(k) + phi(k,j,i)*inymr*inzmr - - phirms(k) = phirms(k) + phi(k,j,i)**2*inymr*inzmr - end do - end do - end do - - call MpiSumReal1D(phibar,nxmr) - call MpiSumReal1D(phirms,nxmr) - - do k=1,nxmr - phirms(k) = sqrt(phirms(k)) - end do - end if - write(nstat,"(i5.5)")nint(time/tout) filename = trim("outputdir/means.h5") @@ -380,6 +213,9 @@ subroutine CalcMeanProfiles if (.not.fexist) then if (ismaster) then call HdfCreateMeansFile(filename) + if (salinity) call CreateSalinityH5Groups(filename) + if (phasefield) call CreatePhaseH5Groups(filename) + if (moist) call CreateMoistH5Groups(filename) end if end if @@ -431,35 +267,14 @@ subroutine CalcMeanProfiles dsetname = trim("epsilon/"//nstat) call HdfSerialWriteReal1D(dsetname,filename,epsilon,nxm) - - if (salinity) then - dsetname = trim("Sbar/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,Sbar,nxmr) - dsetname = trim("vxS/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,vxS,nxmr) - - dsetname = trim("vyS/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,vyS,nxmr) + end if - dsetname = trim("vzS/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,vzS,nxmr) + if (salinity) call CalcSalStats - dsetname = trim("Srms/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,Srms,nxmr) + if (phasefield) call CalcPhiStats - dsetname = trim("chiS/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,chiS,nxmr) - end if - - if (phasefield) then - dsetname = trim("phibar/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,phibar,nxmr) - - dsetname = trim("phirms/"//nstat) - call HdfSerialWriteReal1D(dsetname,filename,phirms,nxmr) - end if - end if + if (moist) call CalcMoistStats call MpiBarrier @@ -514,26 +329,6 @@ subroutine HdfCreateMeansFile(filename) ! call h5gclose_f(group_id,hdf_error) call h5gcreate_f(file_id,"epsilon",group_id,hdf_error) call h5gclose_f(group_id,hdf_error) - if (salinity) then - call h5gcreate_f(file_id,"Sbar",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - call h5gcreate_f(file_id,"vxS",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - call h5gcreate_f(file_id,"vyS",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - call h5gcreate_f(file_id,"vzS",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - call h5gcreate_f(file_id,"Srms",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - call h5gcreate_f(file_id,"chiS",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - end if - if (phasefield) then - call h5gcreate_f(file_id,"phibar",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - call h5gcreate_f(file_id,"phirms",group_id,hdf_error) - call h5gclose_f(group_id,hdf_error) - end if - + call h5fclose_f(file_id, hdf_error) end subroutine HdfCreateMeansFile \ No newline at end of file diff --git a/src/flow_solver/CheckDivergence.F90 b/src/flow_solver/CheckDivergence.F90 index 8fc406a7..91458532 100644 --- a/src/flow_solver/CheckDivergence.F90 +++ b/src/flow_solver/CheckDivergence.F90 @@ -8,81 +8,81 @@ ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine CheckDivergence(qmax,qmaxr) - use param - use local_arrays, only: vy,vx,vz - use mgrd_arrays, only: vyr,vxr,vzr - use mpih - use decomp_2d, only: xstart,xend,xstartr,xendr,nrank - implicit none - real,intent(out) :: qmax,qmaxr - integer :: jc,kc,kp,jp,ic,ip - real :: dqcap - - qmax =-huge(0.0) - +subroutine CheckDivergence(qmax,qmaxr) + use param + 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 + implicit none + real,intent(out) :: qmax,qmaxr + integer :: jc,kc,kp,jp,ic,ip + real :: dqcap + + qmax =-huge(0.0) + !$OMP PARALLEL DO & !$OMP DEFAULT(none) & !$OMP SHARED(xstart,xend,nxm,vz,vy,vx,dz,dy,udx3m) & !$OMP PRIVATE(ic,jc,kc,ip,jp,kp) & !$OMP PRIVATE(dqcap) & !$OMP REDUCTION(max:qmax) - do ic=xstart(3),xend(3) + do ic=xstart(3),xend(3) ip=ic+1 do jc=xstart(2),xend(2) - jp=jc+1 + jp=jc+1 do kc=1,nxm - kp=kc+1 - dqcap= (vz(kc,jc,ip)-vz(kc,jc,ic))*dz & - +(vy(kc,jp,ic)-vy(kc,jc,ic))*dy & - +(vx(kp,jc,ic)-vx(kc,jc,ic))*udx3m(kc) - qmax = max(abs(dqcap),qmax) - enddo - enddo - enddo -!$OMP END PARALLEL DO - - call MpiMaxRealScalar(qmax) - - if (salinity) then -!------------------------------------ + kp=kc+1 + dqcap= (vz(kc,jc,ip) - vz(kc,jc,ic))*dz & + +(vy(kc,jp,ic) - vy(kc,jc,ic))*dy & + +(vx(kp,jc,ic) - vx(kc,jc,ic))*udx3m(kc) + qmax = max(abs(dqcap),qmax) + end do + end do + end do + !$OMP END PARALLEL DO + + call MpiMaxRealScalar(qmax) + + if (salinity) then + !------------------------------------ qmaxr =-huge(0.0) - ! if(nrank.eq.0) write(*,*) "I J K RANK" - ! !$OMP PARALLEL DO & - ! !$OMP DEFAULT(none) & - ! !$OMP SHARED(xstartr,xendr,nxmr,vzr,vyr,vxr,dzr,dyr,udx3mr) & - ! !$OMP PRIVATE(ic,jc,kc,ip,jp,kp) & - ! !$OMP PRIVATE(dqcap) & - ! !$OMP REDUCTION(max:qmaxr) + ! if(nrank.eq.0) write(*,*) "I J K RANK" +! !$OMP PARALLEL DO & +! !$OMP DEFAULT(none) & +! !$OMP SHARED(xstartr,xendr,nxmr,vzr,vyr,vxr,dzr,dyr,udx3mr) & +! !$OMP PRIVATE(ic,jc,kc,ip,jp,kp) & +! !$OMP PRIVATE(dqcap) & +! !$OMP REDUCTION(max:qmaxr) do ic=xstartr(3),xendr(3) - ip=ic+1 - do jc=xstartr(2),xendr(2) - jp=jc+1 - do kc=1,nxmr - kp=kc+1 - dqcap= (vzr(kc,jc,ip)-vzr(kc,jc,ic))*dzr & - +(vyr(kc,jp,ic)-vyr(kc,jc,ic))*dyr & - +(vxr(kp,jc,ic)-vxr(kc,jc,ic))*udx3mr(kc) - !if (abs(dqcap).gt.resid) then - !write(*,*) ic,jc,kc,nrank - !write(*,*) "vz",(vz(kc,jc,ip)-vz(kc,jc,ic))*dz - !write(*,*) "vy",(vy(kc,jp,ic)-vy(kc,jc,ic))*dy - !write(*,*) "vx",(vx(kp,jc,ic)-vx(kc,jc,ic))*udx3m(kc) - !write(*,*) "vym",ic,jc,kc,vy(kc,jc,ic) - !write(*,*) "vyp",ic,jp,kc,vy(kc,jp,ic) - !endif - qmaxr = max(abs(dqcap),qmaxr) - enddo - enddo - enddo - ! !$OMP END PARALLEL DO - + ip=ic+1 + do jc=xstartr(2),xendr(2) + jp=jc+1 + do kc=1,nxmr + kp=kc+1 + dqcap= (vzr(kc,jc,ip) - vzr(kc,jc,ic))*dzr & + +(vyr(kc,jp,ic) - vyr(kc,jc,ic))*dyr & + +(vxr(kp,jc,ic) - vxr(kc,jc,ic))*udx3mr(kc) + !if (abs(dqcap).gt.resid) then + !write(*,*) ic,jc,kc,nrank + !write(*,*) "vz",(vz(kc,jc,ip)-vz(kc,jc,ic))*dz + !write(*,*) "vy",(vy(kc,jp,ic)-vy(kc,jc,ic))*dy + !write(*,*) "vx",(vx(kp,jc,ic)-vx(kc,jc,ic))*udx3m(kc) + !write(*,*) "vym",ic,jc,kc,vy(kc,jc,ic) + !write(*,*) "vyp",ic,jp,kc,vy(kc,jp,ic) + !endif + qmaxr = max(abs(dqcap),qmaxr) + end do + end do + end do + ! !$OMP END PARALLEL DO + call MpiMaxRealScalar(qmaxr) - - else + + else qmaxr = qmax - end if + end if + - - return - end + return +end \ No newline at end of file diff --git a/src/flow_solver/CorrectVelocity.F90 b/src/flow_solver/CorrectVelocity.F90 index 4534558a..bfc6e385 100644 --- a/src/flow_solver/CorrectVelocity.F90 +++ b/src/flow_solver/CorrectVelocity.F90 @@ -10,44 +10,18 @@ subroutine CorrectVelocity use param - use local_arrays, only: vy,vx,dphhalo,vz,temp - use mgrd_arrays, only: sal + use local_arrays, only: vy,vx,vz,temp + use afid_salinity, only: sal use decomp_2d, only: xstart,xend,xstartr,xendr + use afid_pressure, only: RemoveDivergence use mpih implicit none - integer :: jc,jm,kc,km,ic,im,kmid - real :: usukm,udy,udz,locdph + integer :: jc,kc,ic,kmid real, dimension(nxm) :: vxbar, Gshape real, dimension(nxmr):: GshapeR real :: vybulk, vzbulk, Tbulk, Sbulk, idx, vz_target, lam, lfac - udy = al*dt*dy - udz = al*dt*dz - -!$OMP PARALLEL DO & -!$OMP DEFAULT(none) & -!$OMP SHARED(vz,vy,vx,dphhalo,udz,udy,udx3c) & -!$OMP SHARED(xstart,xend,nxm,kmv,dt,al) & -!$OMP PRIVATE(ic,jc,kc) & -!$OMP PRIVATE(im,jm,km,usukm,locdph) - do ic=xstart(3),xend(3) - im=ic-1 - do jc=xstart(2),xend(2) - jm=jc-1 - do kc=1,nxm - km=kmv(kc) - usukm = al*dt*udx3c(kc) - locdph=dphhalo(kc,jc,ic) - vx(kc,jc,ic)=vx(kc,jc,ic)- & - (locdph-dphhalo(km,jc,ic))*usukm - vy(kc,jc,ic)=vy(kc,jc,ic)- & - (locdph-dphhalo(kc,jm,ic))*udy - vz(kc,jc,ic)=vz(kc,jc,ic)- & - (locdph-dphhalo(kc,jc,im))*udz - enddo - enddo - enddo -!$OMP END PARALLEL DO + call RemoveDivergence !CJH Prescribe mean volume flux !Treat dPdz input as a desired Re_b diff --git a/src/flow_solver/CreateInitialConditions.F90 b/src/flow_solver/CreateInitialConditions.F90 index dadd2256..e0582e8e 100644 --- a/src/flow_solver/CreateInitialConditions.F90 +++ b/src/flow_solver/CreateInitialConditions.F90 @@ -13,6 +13,8 @@ subroutine CreateInitialConditions use local_arrays, only: vy,vx,temp,vz use decomp_2d, only: xstart,xend use mpih + use afid_salinity, only: RayS + use afid_phasefield, only: pf_eps, read_phase_field_params implicit none integer :: j,k,i,kmid, io real :: xxx,yyy,zzz,eps,varptb,amp @@ -133,7 +135,7 @@ subroutine CreateInitialConditions end if ! Set velocity to zero if we are using the phase-field method - if (melt .or. phasefield .or. IBM) then + if (melt .or. phasefield .or. IBM .or. moist) then do i=xstart(3),xend(3) do j=xstart(2),xend(2) do k=1,nxm @@ -224,6 +226,16 @@ subroutine CreateInitialConditions end do end if + if (moist) then + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + temp(k,j,i) = 0.0 + end do + end do + end do + end if + if (phasefield) then if (pf_IC == 1) then ! 1D moving interface example @@ -391,6 +403,21 @@ subroutine CreateInitialConditions end if end if + ! if (IBM) then + ! do i=xstart(3),xend(3) + ! do j=xstart(2),xend(2) + ! h0 = 0.25 + (ym(j) - 0.5)**2 + ! do k=1,nxm + ! if (xm(k) < h0) then + ! temp(k,j,i) = 0.0 + ! else + ! temp(k,j,i) = 1.0 + ! end if + ! end do + ! end do + ! end do + ! end if + end if if (melt) then diff --git a/src/flow_solver/DeallocateVariables.F90 b/src/flow_solver/DeallocateVariables.F90 index 10dabdc7..b9e5abe3 100644 --- a/src/flow_solver/DeallocateVariables.F90 +++ b/src/flow_solver/DeallocateVariables.F90 @@ -1,82 +1,73 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! ! +! ! ! FILE: DeallocateVariables.F90 ! ! CONTAINS: subroutine DeallocateVariables ! -! ! +! ! ! PURPOSE: Finalization routine. Deallocates all ! ! variables used in the code ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine DeallocateVariables - use param - use local_arrays - use mgrd_arrays - use stat_arrays - use AuxiliaryRoutines - implicit none - - call DestroyReal1DArray(zc) - call DestroyReal1DArray(zm) - call DestroyReal1DArray(ak1) - call DestroyReal1DArray(ao) - - call DestroyReal1DArray(yc) - call DestroyReal1DArray(ym) - call DestroyReal1DArray(ak2) - call DestroyReal1DArray(ap) - - call DestroyReal1DArray(xc) - call DestroyReal1DArray(xm) - call DestroyReal1DArray(g3rc) - call DestroyReal1DArray(g3rm) - - call DestroyReal1DArray(udx3c) - call DestroyReal1DArray(udx3m) - - call DestroyReal1DArray(d3xc) - call DestroyReal1DArray(d3xm) - - call DestroyReal1DArray(ap3ck) - call DestroyReal1DArray(ac3ck) - call DestroyReal1DArray(am3ck) - - call DestroyReal1DArray(ap3sk) - call DestroyReal1DArray(ac3sk) - call DestroyReal1DArray(am3sk) - - call DestroyReal1DArray(ap3ssk) - call DestroyReal1DArray(ac3ssk) - call DestroyReal1DArray(am3ssk) - - call DestroyReal1DArray(amphk) - call DestroyReal1DArray(acphk) - call DestroyReal1DArray(apphk) - - call DestroyInt1dArray(kmc) - call DestroyInt1dArray(kpc) - call DestroyInt1dArray(kmv) - call DestroyInt1dArray(kpv) - - call DestroyReal3DArray(tempbp) - call DestroyReal3DArray(temptp) - - call DestroyReal3DArray(vx) - call DestroyReal3DArray(vy) - call DestroyReal3DArray(vz) - call DestroyReal3DArray(temp) - - call DestroyReal3DArray(pr) - call DestroyReal3DArray(rhs) - - call DestroyReal3DArray(dph) - call DestroyReal3DArray(dphhalo) - - call DestroyReal3DArray(rux) - call DestroyReal3DArray(ruy) - call DestroyReal3DArray(ruz) - call DestroyReal3DArray(rutemp) - - return - end +subroutine DeallocateVariables + use param + use local_arrays + use mgrd_arrays + use stat_arrays + use AuxiliaryRoutines + implicit none + + call DestroyReal1DArray(zc) + call DestroyReal1DArray(zm) + + call DestroyReal1DArray(yc) + call DestroyReal1DArray(ym) + + call DestroyReal1DArray(xc) + call DestroyReal1DArray(xm) + call DestroyReal1DArray(g3rc) + call DestroyReal1DArray(g3rm) + + call DestroyReal1DArray(udx3c) + call DestroyReal1DArray(udx3m) + + call DestroyReal1DArray(d3xc) + call DestroyReal1DArray(d3xm) + + call DestroyReal1DArray(ap3ck) + call DestroyReal1DArray(ac3ck) + call DestroyReal1DArray(am3ck) + + call DestroyReal1DArray(ap3sk) + call DestroyReal1DArray(ac3sk) + call DestroyReal1DArray(am3sk) + + call DestroyReal1DArray(ap3ssk) + call DestroyReal1DArray(ac3ssk) + call DestroyReal1DArray(am3ssk) + + call DestroyInt1dArray(kmc) + call DestroyInt1dArray(kpc) + call DestroyInt1dArray(kmv) + call DestroyInt1dArray(kpv) + + call DestroyReal3DArray(tempbp) + call DestroyReal3DArray(temptp) + + call DestroyReal3DArray(vx) + call DestroyReal3DArray(vy) + call DestroyReal3DArray(vz) + call DestroyReal3DArray(temp) + + call DestroyReal3DArray(pr) + call DestroyReal3DArray(rhs) + + call DestroyReal3DArray(dph) + call DestroyReal3DArray(dphhalo) + + call DestroyReal3DArray(rux) + call DestroyReal3DArray(ruy) + call DestroyReal3DArray(ruz) + call DestroyReal3DArray(rutemp) + return +end \ No newline at end of file diff --git a/src/flow_solver/ExplicitTermsTemp.F90 b/src/flow_solver/ExplicitTermsTemp.F90 index 17ed356e..e0de8c4f 100644 --- a/src/flow_solver/ExplicitTermsTemp.F90 +++ b/src/flow_solver/ExplicitTermsTemp.F90 @@ -8,128 +8,79 @@ ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine ExplicitTermsTemp - use param - use local_arrays, only: vy,vx,temp,vz,hro - use decomp_2d, only: xstart,xend - implicit none - integer :: jc,kc,ic - integer :: km,kp,jm,jp,im,ip - real :: htx,hty,htz,udy,udz - real :: udzq,udyq - real :: dyyt,dzzt +subroutine ExplicitTermsTemp + use param + use local_arrays, only: vy,vx,temp,vz,hro + use decomp_2d, only: xstart,xend + implicit none + integer :: jc,kc,ic + integer :: km,kp,jm,jp,im,ip + real :: htx,hty,htz,udy,udz + real :: udzq,udyq + real :: dyyt,dzzt + + udz=dz*0.5d0 + udy=dy*0.5d0 + udzq=dzq/pect + udyq=dyq/pect + + !$OMP PARALLEL DO & + !$OMP DEFAULT(none) & + !$OMP SHARED(xstart,xend,vz,vy,vx,nxm) & + !$OMP SHARED(kmv,kpv,am3sk,ac3sk,ap3sk,udz) & + !$OMP SHARED(udy,udzq,udyq,udx3c,temp,hro) & + !$OMP PRIVATE(ic,jc,kc,im,ip,km,kp,jm,jp) & + !$OMP PRIVATE(htx,hty,htz,dyyt,dzzt) + do ic=xstart(3),xend(3) + im=ic-1 + ip=ic+1 + do jc=xstart(2),xend(2) + jm=jc-1 + jp=jc+1 + do kc=1,nxm - udz=dz*0.5d0 - udy=dy*0.5d0 - udzq=dzq/pect - udyq=dyq/pect + ! ! x-advection d/dx(vx * T) + if (kc==1) then + htx = ( & + vx(kc+1,jc,ic)*(temp(kc+1,jc,ic) + temp(kc,jc,ic)) & + - vx(kc ,jc,ic)*2.0*tempbp(1,jc,ic) & + )*0.5*udx3m(kc) + elseif (kc==nxm) then + htx = ( & + vx(kc+1,jc,ic)*2.0*temptp(1,jc,ic) & + - vx(kc ,jc,ic)*(temp(kc,jc,ic) + temp(kc-1,jc,ic)) & + )*0.5*udx3m(kc) + else + htx = ( & + vx(kc+1,jc,ic)*(temp(kc+1,jc,ic) + temp(kc ,jc,ic)) & + - vx(kc ,jc,ic)*(temp(kc ,jc,ic) + temp(kc-1,jc,ic)) & + )*0.5*udx3m(kc) + end if -!$OMP PARALLEL DO & -!$OMP DEFAULT(none) & -!$OMP SHARED(xstart,xend,vz,vy,vx,nxm) & -!$OMP SHARED(kmv,kpv,am3sk,ac3sk,ap3sk,udz) & -!$OMP SHARED(udy,udzq,udyq,udx3c,temp,hro) & -!$OMP PRIVATE(ic,jc,kc,im,ip,km,kp,jm,jp) & -!$OMP PRIVATE(htx,hty,htz,dyyt,dzzt) - do ic=xstart(3),xend(3) - im=ic-1 - ip=ic+1 - do jc=xstart(2),xend(2) - jm=jc-1 - jp=jc+1 + ! ! z-advection d/dx(vz * T) + htz = ( & + vz(kc,jc,ip)*(temp(kc,jc,ip)+temp(kc,jc,ic)) & + - vz(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jc,im)) & + )*udz + + ! ! y-advection d/dx(vy * T) + hty=( & + vy(kc,jp,ic)*(temp(kc,jp,ic)+temp(kc,jc,ic)) & + - vy(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jm,ic)) & + )*udy - kc = 1 - kp = 2 - htx=(vx(kp,jc,ic)*(temp(kp,jc,ic) + temp(kc,jc,ic)) - & - vx(kc,jc,ic)*2.d0*tempbp(1,jc,ic))*udx3m(kc)*0.5d0 - htz=(vz(kc,jc,ip)*(temp(kc,jc,ip)+temp(kc,jc,ic))- & - vz(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jc,im)) & - )*udz - hty=(vy(kc,jp,ic)*(temp(kc,jp,ic)+temp(kc,jc,ic))- & - vy(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jm,ic)) & - )*udy - dzzt=(temp(kc,jc,ip) - 2.0*temp(kc,jc,ic) + temp(kc,jc,im))*udzq - dyyt=(temp(kc,jp,ic) - 2.0*temp(kc,jc,ic) + temp(kc,jm,ic))*udyq - hro(kc,jc,ic) = -(htx+hty+htz)+dyyt+dzzt - ! hro(kc,jc,ic) = dyyt+dzzt + ! zz second derivatives of temp + dzzt=(temp(kc,jc,ip) - 2.0*temp(kc,jc,ic) + temp(kc,jc,im))*udzq + ! yy second derivatives of temp + dyyt=(temp(kc,jp,ic) - 2.0*temp(kc,jc,ic) + temp(kc,jm,ic))*udyq - do kc=2,nxm-1 - km=kc-1 - kp=kc+1 - ! - ! rho vx term - ! - ! - ! d rho q_x - ! ----------- - ! d x - ! - htx=(vx(kp,jc,ic)*(temp(kp,jc,ic)+temp(kc,jc,ic))- & - vx(kc,jc,ic)*(temp(kc,jc,ic)+temp(km,jc,ic)) & - )*udx3m(kc)*0.5d0 -! -! -! rho vz term -! -! -! d rho q_z -! ----------- -! d z -! - htz=(vz(kc,jc,ip)*(temp(kc,jc,ip)+temp(kc,jc,ic))- & - vz(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jc,im)) & - )*udz -! -! -! rho vy term -! -! -! d rho q_y -! ----------- -! d y -! - hty=(vy(kc,jp,ic)*(temp(kc,jp,ic)+temp(kc,jc,ic))- & - vy(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jm,ic)) & - )*udy -! -! -! zz second derivatives of temp -! - dzzt=(temp(kc,jc,ip) & - -2.0*temp(kc,jc,ic) & - +temp(kc,jc,im))*udzq - -! -! yy second derivatives of temp -! - dyyt=(temp(kc,jp,ic) & - -2.0*temp(kc,jc,ic) & - +temp(kc,jm,ic))*udyq -! - hro(kc,jc,ic) = -(htx+hty+htz)+dyyt+dzzt - ! hro(kc,jc,ic) = dyyt+dzzt - enddo - - kc = nxm - kp = nx - km = nxm - 1 - htx=(vx(kp,jc,ic)*2.d0*temptp(1,jc,ic) - & - vx(kc,jc,ic)*(temp(kc,jc,ic)+temp(km,jc,ic)) & - )*udx3m(kc)*0.5d0 - htz=(vz(kc,jc,ip)*(temp(kc,jc,ip)+temp(kc,jc,ic))- & - vz(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jc,im)) & - )*udz - hty=(vy(kc,jp,ic)*(temp(kc,jp,ic)+temp(kc,jc,ic))- & - vy(kc,jc,ic)*(temp(kc,jc,ic)+temp(kc,jm,ic)) & - )*udy - dzzt=(temp(kc,jc,ip) - 2.0*temp(kc,jc,ic) + temp(kc,jc,im))*udzq - dyyt=(temp(kc,jp,ic) - 2.0*temp(kc,jc,ic) + temp(kc,jm,ic))*udyq - hro(kc,jc,ic) = -(htx+hty+htz)+dyyt+dzzt - ! hro(kc,jc,ic) = dyyt+dzzt - - enddo - enddo -!$OMP END PARALLEL DO - - return - end + hro(kc,jc,ic) = -(htx+hty+htz)+dyyt+dzzt + ! hro(kc,jc,ic) = dyyt+dzzt + end do + + end do + end do + !$OMP END PARALLEL DO + + return +end \ No newline at end of file diff --git a/src/flow_solver/ExplicitTermsVX.F90 b/src/flow_solver/ExplicitTermsVX.F90 index bc40c370..3ce121cb 100644 --- a/src/flow_solver/ExplicitTermsVX.F90 +++ b/src/flow_solver/ExplicitTermsVX.F90 @@ -11,19 +11,16 @@ subroutine ExplicitTermsVX use param use local_arrays, only: vz,vy,vx,temp,qcap - use mgrd_arrays, only: salc,phic use decomp_2d, only: xstart,xend implicit none integer :: jc,kc integer :: km,kp,jm,jp,ic,im,ip real :: hxx,hxy,hxz - real :: udz,udy,tempit,salit,volpen + real :: udz,udy,tempit real :: udzq,udyq - real :: dzzvx,dyyvx,pf_eta + real :: dzzvx,dyyvx real, dimension(nxm) :: idx - pf_eta = ren*(1.51044385*pf_eps)**2 - udy=dy*0.25 udz=dz*0.25 @@ -127,44 +124,6 @@ subroutine ExplicitTermsVX end do end do end do - - !CJH Add salinity component of buoyancy if used - if (salinity) then - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - do kc=2,nxm - salit = active_S*salc(kc,jc,ic) - qcap(kc,jc,ic) = qcap(kc,jc,ic) - bycs*salit - end do - end do - end do - - if (melt) then - !CJH For melting boundary, remove far-field buoyancy - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - do kc=2,nxm - qcap(kc,jc,ic) = qcap(kc,jc,ic) + bycs - byct - end do - end do - end do - end if - end if - end if - - if (phasefield) then - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - do kc=2,nxm - km = kc - 1 - volpen = 0.5d0*(phic(kc,jc,ic) + phic(km,jc,ic))* & - vx(kc,jc,ic)/pf_eta - ! volpen = 0.5*(1.0+tanh((xc(kc)-0.5)/2/pf_eps))* & - ! vx(kc,jc,ic)/pf_gamma - qcap(kc,jc,ic) = qcap(kc,jc,ic) - volpen - end do - end do - end do end if !$OMP END PARALLEL DO diff --git a/src/flow_solver/ExplicitTermsVY.F90 b/src/flow_solver/ExplicitTermsVY.F90 index 3e023af2..f074412c 100644 --- a/src/flow_solver/ExplicitTermsVY.F90 +++ b/src/flow_solver/ExplicitTermsVY.F90 @@ -11,17 +11,14 @@ subroutine ExplicitTermsVY use param use local_arrays, only: vx,vy,vz,dph,temp - use mgrd_arrays, only: salc,phic use decomp_2d, only: xstart,xend implicit none integer :: kc,kp,jpp,jmm,jc,ic,imm,ipp integer :: kpp,kmm real :: udzq,udyq real :: udy,udz,hyx,hyy,hyz - real :: dyyvy, dzzvy, pf_eta - real :: tempit, salit, volpen, Gy - - pf_eta = ren*(1.51044385*pf_eps)**2 + real :: dyyvy, dzzvy + real :: tempit, Gy udyq=dyq/ren udzq=dzq/ren @@ -122,44 +119,8 @@ subroutine ExplicitTermsVY end do end do end do - - !CJH Add salinity component of buoyancy if used - if (salinity) then - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - jmm=jc-1 - do kc=1,nxm - salit =active_S*salc(kc,jc,ic) - dph(kc,jc,ic) = dph(kc,jc,ic) - bycs*salit - end do - end do - end do - - if (melt) then - !CJH For melting boundary, remove far-field buoyancy - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - do kc=1,nxm - dph(kc,jc,ic) = dph(kc,jc,ic) + bycs - byct - end do - end do - end do - end if - end if end if - if (phasefield) then - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - jmm=jc-1 - do kc=1,nxm - volpen = 0.5d0*(phic(kc,jc,ic) + phic(kc,jmm,ic))* & - vy(kc,jc,ic)/pf_eta - dph(kc,jc,ic) = dph(kc,jc,ic) - volpen - end do - end do - end do - end if !$OMP END PARALLEL DO return diff --git a/src/flow_solver/ExplicitTermsVZ.F90 b/src/flow_solver/ExplicitTermsVZ.F90 index 968973d0..0d6b3d18 100644 --- a/src/flow_solver/ExplicitTermsVZ.F90 +++ b/src/flow_solver/ExplicitTermsVZ.F90 @@ -11,17 +11,14 @@ subroutine ExplicitTermsVZ use param use local_arrays, only: vx,vy,vz,dq,temp - use mgrd_arrays, only: salc,phic use decomp_2d, only: xstart,xend implicit none integer :: kc,kp,jpp,jmm,jc,ic,imm,ipp integer :: kmm,kpp real :: hzx,hzy,hzz,udy,udz real :: udyq,udzq - real :: dzzvz,dyyvz,pf_eta - real :: tempit, salit,volpen!, Gz - - pf_eta = ren*(1.51044385*pf_eps)**2 + real :: dzzvz,dyyvz + real :: tempit!, Gz udyq=dyq/ren udzq=dzq/ren @@ -118,45 +115,8 @@ subroutine ExplicitTermsVZ end do end do end do - - !CJH Add salinity component of buoyancy if used - if (salinity) then - do ic=xstart(3),xend(3) - imm=ic-1 - do jc=xstart(2),xend(2) - do kc=1,nxm - salit =active_S*salc(kc,jc,ic) - dq(kc,jc,ic) = dq(kc,jc,ic) - bycs*salit - end do - end do - end do - - if (melt) then - !CJH For melting boundary, remove far-field buoyancy - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - do kc=1,nxm - dq(kc,jc,ic) = dq(kc,jc,ic) + bycs - byct - end do - end do - end do - end if - end if end if !$OMP END PARALLEL DO - if (phasefield) then - do ic=xstart(3),xend(3) - imm=ic-1 - do jc=xstart(2),xend(2) - do kc=1,nxm - volpen = 0.5d0*(phic(kc,jc,ic) + phic(kc,jc,imm))* & - vz(kc,jc,ic)/pf_eta - dq(kc,jc,ic) = dq(kc,jc,ic) - volpen - end do - end do - end do - end if - return end subroutine ExplicitTermsVZ \ No newline at end of file diff --git a/src/flow_solver/ImplicitAndUpdateTemp.F90 b/src/flow_solver/ImplicitAndUpdateTemp.F90 index 0d737b8e..80400bb0 100644 --- a/src/flow_solver/ImplicitAndUpdateTemp.F90 +++ b/src/flow_solver/ImplicitAndUpdateTemp.F90 @@ -68,7 +68,7 @@ subroutine ImplicitAndUpdateTemp enddo !$OMP END PARALLEL DO - if (IBM) then + if (IBM .and. .not. phasefield) then call SolveImpEqnUpdate_Temp_ibm else ! Solve equation and update temperature diff --git a/src/flow_solver/InitVariables.F90 b/src/flow_solver/InitVariables.F90 index 6426fb9d..87dbf3c5 100644 --- a/src/flow_solver/InitVariables.F90 +++ b/src/flow_solver/InitVariables.F90 @@ -1,234 +1,95 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! ! +! ! ! FILE: InitVariables.F90 ! ! CONTAINS: subroutine InitVariables ! -! ! +! ! ! PURPOSE: Initialization routine. Sets to zero all ! ! variables used in the code ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine InitVariables - use param - use local_arrays - ! use mgrd_arrays - use stat_arrays - use decomp_2d - use AuxiliaryRoutines - implicit none - -!------------------------------------------------- -! Arrays for grid making -!------------------------------------------------- - - call AllocateReal1DArray(zc,1,nz) - call AllocateReal1DArray(zm,0,nz+1) - call AllocateReal1DArray(ak1,1,nz) - call AllocateReal1DArray(ao,1,nz) - - call AllocateReal1DArray(yc,1,ny) - call AllocateReal1DArray(ym,0,ny+1) - call AllocateReal1DArray(ak2,1,ny) - call AllocateReal1DArray(ap,1,ny) - - call AllocateReal1DArray(xc,1,nx) - call AllocateReal1DArray(xm,1,nx) - call AllocateReal1DArray(g3rc,1,nx) - call AllocateReal1DArray(g3rm,1,nx) - - !CJH New terms for modified derivatives - call AllocateReal1DArray(d3xc,1,nx) - call AllocateReal1DArray(d3xm,1,nx) - - call AllocateReal1DArray(udx3c,1,nx) - call AllocateReal1DArray(udx3m,1,nx) - - call AllocateReal1DArray(ap3ck,1,nx) - call AllocateReal1DArray(ac3ck,1,nx) - call AllocateReal1DArray(am3ck,1,nx) - - call AllocateReal1DArray(ap3sk,1,nx) - call AllocateReal1DArray(ac3sk,1,nx) - call AllocateReal1DArray(am3sk,1,nx) - - call AllocateReal1DArray(ap3ssk,1,nx) - call AllocateReal1DArray(ac3ssk,1,nx) - call AllocateReal1DArray(am3ssk,1,nx) - - call AllocateReal1DArray(amphk,1,nx) - call AllocateReal1DArray(acphk,1,nx) - call AllocateReal1DArray(apphk,1,nx) - - call AllocateInt1dArray(kmc,1,nx) - call AllocateInt1dArray(kpc,1,nx) - call AllocateInt1dArray(kmv,1,nx) - call AllocateInt1dArray(kpv,1,nx) +subroutine InitVariables + use param + use local_arrays + use stat_arrays + use decomp_2d + use AuxiliaryRoutines + implicit none !------------------------------------------------- -! Arrays for mgrd grid making +! Arrays for grid making !------------------------------------------------- - ! call AllocateReal1DArray(zcr,1,nzr) - ! call AllocateReal1DArray(zmr,0,nzr+1) !CS mgrd - - ! call AllocateReal1DArray(ycr,1,nyr) - ! call AllocateReal1DArray(ymr,0,nyr+1) !CS mgrd - - ! call AllocateReal1DArray(xcr,1,nxr) - ! call AllocateReal1DArray(xmr,0,nxr) !CS mgrd - - ! call AllocateReal1DArray(g3rcr,1,nxr) - ! call AllocateReal1DArray(g3rmr,1,nxr) + call AllocateReal1DArray(zc,1,nz) + call AllocateReal1DArray(zm,0,nz+1) - ! !CJH New terms for modified derivatives - ! call AllocateReal1DArray(d3xcr,1,nxr) - ! call AllocateReal1DArray(d3xmr,1,nxr) + call AllocateReal1DArray(yc,1,ny) + call AllocateReal1DArray(ym,0,ny+1) - ! call AllocateReal1DArray(udx3cr,1,nxr) - ! call AllocateReal1DArray(udx3mr,1,nxr) + call AllocateReal1DArray(xc,1,nx) + call AllocateReal1DArray(xm,1,nx) + call AllocateReal1DArray(g3rc,1,nx) + call AllocateReal1DArray(g3rm,1,nx) - ! call AllocateReal1DArray(ap3ckr,1,nxr) - ! call AllocateReal1DArray(ac3ckr,1,nxr) - ! call AllocateReal1DArray(am3ckr,1,nxr) + !CJH New terms for modified derivatives + call AllocateReal1DArray(d3xc,1,nx) + call AllocateReal1DArray(d3xm,1,nx) - ! call AllocateReal1DArray(ap3sskr,1,nxr) - ! call AllocateReal1DArray(ac3sskr,1,nxr) - ! call AllocateReal1DArray(am3sskr,1,nxr) + call AllocateReal1DArray(udx3c,1,nx) + call AllocateReal1DArray(udx3m,1,nx) - ! call AllocateInt1dArray(kmcr,1,nxr) - ! call AllocateInt1dArray(kpcr,1,nxr) - ! call AllocateInt1dArray(kmvr,1,nxr) - ! call AllocateInt1dArray(kpvr,1,nxr) + call AllocateReal1DArray(ap3ck,1,nx) + call AllocateReal1DArray(ac3ck,1,nx) + call AllocateReal1DArray(am3ck,1,nx) - ! call AllocateInt1dArray(irangs,0,nx) - ! call AllocateInt1dArray(jrangs,0,ny) - ! call AllocateInt1dArray(krangs,0,nz) - - ! call AllocateInt1dArray(irangc,0,nx) - ! call AllocateInt1dArray(jrangc,0,ny) - ! call AllocateInt1dArray(krangc,0,nz) - - ! call AllocateReal2DArray(cxvx,1,4,0,nxr) - ! call AllocateReal2DArray(cxvy,1,4,0,nxr) - ! call AllocateReal2DArray(cxvz,1,4,0,nxr) - - ! call AllocateReal2DArray(cyvx,1,4,0,nyr) - ! call AllocateReal2DArray(cyvy,1,4,0,nyr) - ! call AllocateReal2DArray(cyvz,1,4,0,nyr) - - ! call AllocateReal2DArray(czvx,1,4,0,nzr) - ! call AllocateReal2DArray(czvy,1,4,0,nzr) - ! call AllocateReal2DArray(czvz,1,4,0,nzr) - - ! call AllocateReal2DArray(cxrs,1,4,0,nxr) - ! call AllocateReal2DArray(cyrs,1,4,0,nyr) - ! call AllocateReal2DArray(czrs,1,4,0,nzr) - - ! call AllocateReal2DArray(cxsalc,1,4,0,nx) - ! call AllocateReal2DArray(cysalc,1,4,0,ny) - ! call AllocateReal2DArray(czsalc,1,4,0,nz) - -!------------------------------------------------- -! Arrays for temperature boundary conditions -!------------------------------------------------- + call AllocateReal1DArray(ap3sk,1,nx) + call AllocateReal1DArray(ac3sk,1,nx) + call AllocateReal1DArray(am3sk,1,nx) - call AllocateReal3DArray(tempbp,1,1,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - call AllocateReal3DArray(temptp,1,1,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal1DArray(ap3ssk,1,nx) + call AllocateReal1DArray(ac3ssk,1,nx) + call AllocateReal1DArray(am3ssk,1,nx) - !-- For salinity - ! call AllocateReal3DArray(salbp,1,1,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) - ! call AllocateReal3DArray(saltp,1,1,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + call AllocateInt1dArray(kmc,1,nx) + call AllocateInt1dArray(kpc,1,nx) + call AllocateInt1dArray(kmv,1,nx) + call AllocateInt1dArray(kpv,1,nx) !------------------------------------------------- -! Arrays for statistics +! Arrays for temperature boundary conditions !------------------------------------------------- - ! if (statcal) then - ! call AllocateReal1DArray(vx_me,1,nxm) - ! call AllocateReal1DArray(vy_me,1,nxm) - ! call AllocateReal1DArray(vz_me,1,nxm) - - ! call AllocateReal1DArray(vx_rms,1,nxm) - ! call AllocateReal1DArray(vy_rms,1,nxm) - ! call AllocateReal1DArray(vz_rms,1,nxm) - - ! call AllocateReal1DArray(vx_me_buf,1,nxm) - ! call AllocateReal1DArray(vy_me_buf,1,nxm) - ! call AllocateReal1DArray(vz_me_buf,1,nxm) - - ! call AllocateReal1DArray(vx_msq_buf,1,nxm) - ! call AllocateReal1DArray(vy_msq_buf,1,nxm) - ! call AllocateReal1DArray(vz_msq_buf,1,nxm) - - ! call AllocateReal1DArray(temp_me,1,nxm) - ! call AllocateReal1DArray(temp_rms,1,nxm) - ! call AllocateReal1DArray(tempvx_me,1,nxm) - - ! call AllocateReal1DArray(vxvy_corr,1,nxm) - ! if (disscal) then - ! call AllocateReal1DArray(disste,1,nxm) - ! call AllocateReal1DArray(dissth,1,nxm) - ! end if - ! end if - - !------------------------------------------------- - ! Arrays with ghost cells - !------------------------------------------------- - call AllocateReal3DArray(vy,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - call AllocateReal3DArray(vz,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - call AllocateReal3DArray(vx,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - call AllocateReal3DArray(pr,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - call AllocateReal3DArray(temp,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - call AllocateReal3DArray(dphhalo,1,nxm,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - - !-- For salinity - ! call AllocateReal3DArray(sal,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) - - !-- For Q-criteria - call AllocateReal3DArray(qtens,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) - - !----------------------------------------------- - ! Arrays without ghost cells - !----------------------------------------------- - call AllocateReal3DArray(rhs,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(dph,1,nxm,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(dq,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(qcap,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(rux,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(ruy,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(ruz,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(hro,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - call AllocateReal3DArray(rutemp,1,nx,xstart(2),xend(2),xstart(3),xend(3)) - - !-- For salinity - ! call AllocateReal3DArray(rhsr,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) - ! call AllocateReal3DArray(rusal,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) - ! call AllocateReal3DArray(hsal,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) - - return - end - -! -------------------------------- -! subroutine InitMgrdVariables -! use param -! use mgrd_arrays -! use decomp_2d -! use AuxiliaryRoutines -! implicit none - -! call AllocateReal3DArray(vyr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) -! call AllocateReal3DArray(vzr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) -! call AllocateReal3DArray(vxr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) - -! call AllocateReal3DArray(tpdv,-1,nx+1,xstart(2)-2,xend(2)+2,xstart(3)-2,xend(3)+2) -! !CS For tpdvr, larger array needed to prevent memory overflow in InterpVelMgrd -! call AllocateReal3DArray(tpdvr,1,nxmr,xstartr(2)-2,xendr(2)+2,xstartr(3)-2,xendr(3)+2) - -! call AllocateReal3DArray(salc,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) -! !CJH Needed for melt boundary condition -! if (melt) then -! call AllocateReal3DArray(tempr,1,1,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) -! end if - -! return -! end + call AllocateReal3DArray(tempbp,1,1,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal3DArray(temptp,1,1,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + + !------------------------------------------------- + ! Arrays with ghost cells + !------------------------------------------------- + call AllocateReal3DArray(vy,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal3DArray(vz,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal3DArray(vx,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal3DArray(pr,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal3DArray(temp,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + call AllocateReal3DArray(dphhalo,1,nxm,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + + !-- For salinity + ! call AllocateReal3DArray(sal,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + + !-- For Q-criteria + call AllocateReal3DArray(qtens,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + + !----------------------------------------------- + ! Arrays without ghost cells + !----------------------------------------------- + call AllocateReal3DArray(rhs,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(dph,1,nxm,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(dq,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(qcap,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(rux,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(ruy,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(ruz,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(hro,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(rutemp,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + + return +end \ No newline at end of file diff --git a/src/flow_solver/QuitRoutine.F90 b/src/flow_solver/QuitRoutine.F90 index 4cb7a40d..1dccc9ab 100644 --- a/src/flow_solver/QuitRoutine.F90 +++ b/src/flow_solver/QuitRoutine.F90 @@ -13,6 +13,10 @@ subroutine QuitRoutine(tin,normalexit,errorcode) use param use decomp_2d, only: nrank, decomp_2d_finalize use decomp_2d_fft + use afid_pressure, only: DeallocatePressureVars + use afid_moisture, only: DeallocateMoistVariables + use afid_phasefield, only: DeallocatePFVariables + use afid_salinity, only: DeallocateSalVariables implicit none logical, intent(in) :: normalexit integer :: errorcode @@ -34,10 +38,12 @@ subroutine QuitRoutine(tin,normalexit,errorcode) endif call DeallocateVariables + call DeallocatePressureVars if (multires) call DeallocateMgrdVariables if (salinity) call DeallocateSalVariables if (phasefield) call DeallocatePFVariables if (IBM) call DeallocateIBMVariables + if (moist) call DeallocateMoistVariables call HdfClose call decomp_2d_fft_finalize call decomp_2d_finalize diff --git a/src/flow_solver/ReadInputFile.F90 b/src/flow_solver/ReadInputFile.F90 index cf21a8f3..470a8716 100644 --- a/src/flow_solver/ReadInputFile.F90 +++ b/src/flow_solver/ReadInputFile.F90 @@ -9,6 +9,8 @@ subroutine ReadInputFile use param use ibm_param, only: solidtype + use afid_salinity, only: RayS, PraS, bycs, PecS, SfixN, SfixS + 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 diff --git a/src/flow_solver/SetTempBCs.F90 b/src/flow_solver/SetTempBCs.F90 index 9aa5fb0a..a49a30f4 100644 --- a/src/flow_solver/SetTempBCs.F90 +++ b/src/flow_solver/SetTempBCs.F90 @@ -11,6 +11,7 @@ subroutine SetTempBCs use param use decomp_2d + use afid_moisture, only: beta_q implicit none integer :: ic,jc @@ -56,6 +57,15 @@ subroutine SetTempBCs end do end if end if + + if (moist) then + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + tempbp(1,jc,ic) = 0.0 + temptp(1,jc,ic) = beta_q - 1.0 + end do + end do + end if call update_halo(temptp,lvlhalo) call update_halo(tempbp,lvlhalo) diff --git a/src/flow_solver/SpecRoutines.F90 b/src/flow_solver/SpecRoutines.F90 index ff010360..efbb8a0e 100644 --- a/src/flow_solver/SpecRoutines.F90 +++ b/src/flow_solver/SpecRoutines.F90 @@ -2,6 +2,7 @@ subroutine CalcFourierCoef(var,fouvar) use, intrinsic :: iso_c_binding use param use fftw_params +use afid_pressure use decomp_2d use decomp_2d_fft use mpih @@ -19,7 +20,7 @@ subroutine CalcFourierCoef(var,fouvar) !RO Allocate variables for FFT transform -call CreateFFTTmpArrays +! call CreateFFTTmpArrays nymh=nym/2+1 @@ -27,62 +28,7 @@ subroutine CalcFourierCoef(var,fouvar) !RO Plan FFT transforms if not planned previously -if (.not.planned) then - iodim(1)%n=nzm - iodim(1)%is=(sp%zen(1)-sp%zst(1)+1)*(sp%zen(2)-sp%zst(2)+1) - iodim(1)%os=(sp%zen(1)-sp%zst(1)+1)*(sp%zen(2)-sp%zst(2)+1) - iodim_howmany(1)%n=(sp%zen(1)-sp%zst(1)+1) - iodim_howmany(1)%is=1 - iodim_howmany(1)%os=1 - iodim_howmany(2)%n=(sp%zen(2)-sp%zst(2)+1) - iodim_howmany(2)%is=(sp%zen(1)-sp%zst(1)+1) - iodim_howmany(2)%os=(sp%zen(1)-sp%zst(1)+1) - fwd_guruplan_z=fftw_plan_guru_dft(1,iodim, & - 2,iodim_howmany,cz1,cz1, & - FFTW_FORWARD,FFTW_ESTIMATE) - iodim(1)%n=nzm - bwd_guruplan_z=fftw_plan_guru_dft(1,iodim, & - 2,iodim_howmany,cz1,cz1, & - FFTW_BACKWARD,FFTW_ESTIMATE) - - if (.not.c_associated(bwd_guruplan_z)) then - if (ismaster) print*,'Failed to create guru plan. You should' - if (ismaster) print*,'link with FFTW3 before MKL' - if (ismaster) print*,'Please check linking order.' - call MPI_Abort(MPI_COMM_WORLD,1,info) - endif - - iodim(1)%n=nym - iodim(1)%is=ph%yen(1)-ph%yst(1)+1 - iodim(1)%os=sp%yen(1)-sp%yst(1)+1 - iodim_howmany(1)%n=(ph%yen(1)-ph%yst(1)+1) - iodim_howmany(1)%is=1 - iodim_howmany(1)%os=1 - iodim_howmany(2)%n=(ph%yen(3)-ph%yst(3)+1) - iodim_howmany(2)%is=(ph%yen(1)-ph%yst(1)+1) & - *(ph%yen(2)-ph%yst(2)+1) - iodim_howmany(2)%os=(sp%yen(1)-sp%yst(1)+1) & - *(sp%yen(2)-sp%yst(2)+1) - fwd_guruplan_y=fftw_plan_guru_dft_r2c(1,iodim, & - 2,iodim_howmany,ry1,cy1, & - FFTW_ESTIMATE) - - iodim(1)%n=nym - iodim(1)%is=sp%yen(1)-sp%yst(1)+1 - iodim(1)%os=ph%yen(1)-ph%yst(1)+1 - iodim_howmany(1)%n=(sp%yen(1)-sp%yst(1)+1) - iodim_howmany(1)%is=1 - iodim_howmany(1)%os=1 - iodim_howmany(2)%n=(sp%yen(3)-sp%yst(3)+1) - iodim_howmany(2)%is=(sp%yen(1)-sp%yst(1)+1) & - *(sp%yen(2)-sp%yst(2)+1) - iodim_howmany(2)%os=(ph%yen(1)-ph%yst(1)+1) & - *(ph%yen(2)-ph%yst(2)+1) - bwd_guruplan_y=fftw_plan_guru_dft_c2r(1,iodim, & - 2,iodim_howmany,cy1,ry1, & - FFTW_ESTIMATE) - planned=.true. -endif +if (.not.planned) call PlanFourierTransform call dfftw_execute_dft_r2c(fwd_guruplan_y,ry1,cy1) @@ -95,7 +41,7 @@ subroutine CalcFourierCoef(var,fouvar) call transpose_z_to_x(cz1,fouvar,sp) -call DestroyFFTTmpArrays +! call DestroyFFTTmpArrays return end subroutine CalcFourierCoef @@ -293,7 +239,7 @@ end subroutine CalcPowerSpec subroutine WritePowerSpec use param use local_arrays, only: vz,vy,vx,temp - use mgrd_arrays, only: salc + use afid_salinity, only: salc use decomp_2d, only: xstart,xend,xstartr,xendr use stat_arrays use mpih diff --git a/src/flow_solver/TimeMarcher.F90 b/src/flow_solver/TimeMarcher.F90 index 95daee6b..5a5cc11e 100644 --- a/src/flow_solver/TimeMarcher.F90 +++ b/src/flow_solver/TimeMarcher.F90 @@ -13,10 +13,14 @@ subroutine TimeMarcher use param use local_arrays - use mgrd_arrays, only: vxr,vyr,vzr,salc,sal,phi,phic,tempr + ! use mgrd_arrays, only: vxr,vyr,vzr,salc,sal,phi,phic,tempr + use afid_pressure + use afid_salinity + use afid_phasefield use mpih use decomp_2d use ibm_param, only: aldto + use afid_moisture implicit none integer :: ns integer :: j,k,i @@ -35,7 +39,7 @@ subroutine TimeMarcher ga = gam(ns) ro = rom(ns) - if (melt) call UpdateBCs + ! if (melt) call UpdateBCs ! iF(ANY(IsNaN(phi))) write(*,*)nrank,'NaN in PHI pre-explicit' call ExplicitTermsVX @@ -43,26 +47,55 @@ subroutine TimeMarcher call ExplicitTermsVZ call ExplicitTermsTemp - if (phasefield) call ExplicitTermsPhi - ! iF(ANY(IsNaN(phi))) write(*,*)nrank,'NaN in PHI pre-implicit' - if (phasefield) call ImplicitAndUpdatePhi - ! iF(ANY(IsNaN(phi))) write(*,*)nrank,'NaN in PHI post-implicit' + if (salinity) then + call ExplicitSalinity + + ! If using salinity as an active scalar, add its buoyancy contribution + ! to the relevant component of the momentum equation + if (active_S==1) then + if (gAxis==1) then + call AddSalBuoyancy(qcap) + elseif (gAxis==2) then + call AddSalBuoyancy(dph) + elseif (gAxis==3) then + call AddSalBuoyancy(dq) + end if + end if + end if + + if (phasefield) then + call ExplicitPhase + call AddVolumePenalty + if (salinity) then + call AddSaltFluxInterface + call AdjustMeltPoint + end if + call ImplicitPhase + ! Add the latent heat and salt terms *after* computing the implicit solve for phi + call AddLatentHeat + if (salinity) call AddLatentSalt + end if + + if (moist) call ExplicitHumidity - !CJH: Phi must be updated before computing S explicit terms and latent heat - ! varaible rhsr used to store d(phi)/dt for the following subroutines - if (salinity) call ExplicitTermsSal !Refined - if (phasefield) call AddLatentHeat + if (moist) call AddCondensation ! iF(ANY(IsNaN(vx))) write(*,*)nrank,'NaN in VX pre-implicit',ns ! iF(ANY(IsNaN(vy))) write(*,*)nrank,'NaN in VY pre-implicit',ns ! iF(ANY(IsNaN(vz))) write(*,*)nrank,'NaN in VZ pre-implicit',ns ! iF(ANY(IsNaN(temp))) write(*,*)nrank,'NaN in TEMP pre-implicit',ns + if (phasefield) call update_halo(phi,lvlhalo) + ! if (phasefield .and. IBM) call UpdateIBMLocation + call ImplicitAndUpdateVX call ImplicitAndUpdateVY call ImplicitAndUpdateVZ call ImplicitAndUpdateTemp - if (salinity) call ImplicitAndUpdateSal !Refined + + if (salinity) call ImplicitSalinity + + if (moist) call ImplicitHumidity ! iF(ANY(IsNaN(vx))) write(*,*)nrank,'NaN in VX post-implicit',ns ! iF(ANY(IsNaN(vy))) write(*,*)nrank,'NaN in VY post-implicit',ns @@ -70,6 +103,7 @@ subroutine TimeMarcher ! iF(ANY(IsNaN(temp))) write(*,*)nrank,'NaN in TEMP post-implicit',ns ! if (phasefield .and. IBM) call ImmersedBoundary + if (phasefield .and. (pf_direct_force > 0)) call ForceIceVelZero call update_halo(vy,lvlhalo) call update_halo(vz,lvlhalo) @@ -103,23 +137,26 @@ subroutine TimeMarcher call update_halo(temp,lvlhalo) if (salinity) call update_halo(sal,lvlhalo) if (phasefield) call update_halo(phi,lvlhalo) + if (moist) call update_halo(humid,lvlhalo) if (salinity) then call InterpVelMgrd !Vel from base mesh to refined mesh call update_halo(vxr,lvlhalo) call update_halo(vyr,lvlhalo) call update_halo(vzr,lvlhalo) - call InterpSalMgrd !Sal from refined mesh to base mesh + call InterpSalMultigrid !Sal from refined mesh to base mesh call update_halo(salc,lvlhalo) end if if (phasefield) then - call InterpTempMgrd + call InterpTempMultigrid call update_halo(tempr,lvlhalo) - call InterpPhiMgrd + call InterpPhiMultigrid call update_halo(phic,lvlhalo) end if + if (moist) call UpdateSaturation + ! iF(ANY(IsNaN(vx))) write(*,*)nrank,'NaN in VX',ns ! iF(ANY(IsNaN(vy))) write(*,*)nrank,'NaN in VY',ns ! iF(ANY(IsNaN(vz))) write(*,*)nrank,'NaN in VZ',ns diff --git a/src/flow_solver/WriteFlowField.F90 b/src/flow_solver/WriteFlowField.F90 index 1abeb8dc..3047ab50 100644 --- a/src/flow_solver/WriteFlowField.F90 +++ b/src/flow_solver/WriteFlowField.F90 @@ -11,7 +11,9 @@ subroutine WriteFlowField(final) use param use local_arrays, only: vz,vy,vx,temp,pr - use mgrd_arrays, only: sal,phi + use afid_salinity, only: sal + use afid_phasefield, only: phi + use afid_moisture, only: humid implicit none logical, intent(in) :: final character(30) :: filnam1,dsetname,basename @@ -30,6 +32,10 @@ subroutine WriteFlowField(final) call HdfWriteRealHalo3D(filnam1,vy) filnam1 = trim(basename)//'_vz.h5' call HdfWriteRealHalo3D(filnam1,vz) + if (moist) then + filnam1 = trim(basename)//'_humid.h5' + call HdfWriteRealHalo3D(filnam1,humid) + end if if (salinity) then filnam1 = trim(basename)//'_sal.h5' call HdfWriteRealHalo3DR(filnam1,sal) diff --git a/src/flow_solver/param.F90 b/src/flow_solver/param.F90 index c60748a4..ff7def15 100644 --- a/src/flow_solver/param.F90 +++ b/src/flow_solver/param.F90 @@ -13,13 +13,12 @@ module param real :: alx3,str3 integer :: istr3 real :: ylen,zlen - real :: rayt,prat,rays,pras,dt,resid + real :: rayt,prat,dt,resid integer :: inslws,inslwn - integer :: TfixS,TfixN,SfixS,SfixN !CJH option for fixed T/S BCs + integer :: TfixS,TfixN !CJH option for fixed T/S BCs integer :: gAxis !CJH option to choose gravity axis real :: xminusU,xplusU,dPdz,dPdy - real :: pf_A, pf_eps, pf_D, pf_S, pf_Tm, pf_Lambda - logical :: IBM + logical :: IBM, moist ! integer :: starea,tsta real :: dtmin,dtmax,limitCFL integer :: nson,idtv @@ -54,23 +53,16 @@ module param !=========================================================== !******* Metric coefficients ******************************* real, allocatable, dimension(:) :: ap3ck,ac3ck,am3ck - real, allocatable, dimension(:) :: ap3ckr,ac3ckr,am3ckr !CS mgrd real, allocatable, dimension(:) :: ap3sk,ac3sk,am3sk real, allocatable, dimension(:) :: ap3ssk,ac3ssk,am3ssk - real, allocatable, dimension(:) :: ap3sskr,ac3sskr,am3sskr !CS mgrd - real, allocatable, dimension(:) :: ap3spkr,ac3spkr,am3spkr !CJH phase-field !============================================================ - !******* Variables for FFTW and Poisson solver**************** - real, allocatable, dimension(:) :: ak2,ap - real, allocatable, dimension(:) :: ak1,ao - real, allocatable, dimension(:) :: amphk,acphk,apphk !=========================================================== !******* Other variables *********************************** integer :: nxm, nym, nzm integer :: nxmr, nymr, nzmr !CS mgrd - real :: ren, pect, pecs - real :: Rrho, byct, bycs + real :: ren, pect!, pecs + real :: Rrho, byct!, bycs real :: pi real :: al,ga,ro real :: beta @@ -82,7 +74,6 @@ module param real, dimension(1:ndv) :: vmax real, dimension(1:3) :: gam,rom,alm real, allocatable, dimension(:,:,:) :: tempbp,temptp !CJH make BCs 3D arrays so we can use update_halo - real, allocatable, dimension(:,:,:) :: salbp,saltp integer, dimension(5) :: spec_idx logical :: dumpslabs=.false. @@ -114,7 +105,6 @@ module local_arrays real,allocatable,dimension(:,:,:) :: rux,ruy,ruz,rutemp real,allocatable,dimension(:,:,:) :: dph,qcap,dq,hro,dphhalo real,allocatable,dimension(:,:,:) :: qtens - ! real,allocatable,dimension(:,:,:) :: sal,rhsr,rusal,hsal end module local_arrays module mgrd_arrays @@ -124,15 +114,15 @@ module mgrd_arrays integer,allocatable,dimension(:) :: irangc,jrangc,krangc integer,allocatable,dimension(:) :: irangb,jrangb,krangb integer,allocatable,dimension(:) :: irangr,jrangr,krangr - real,allocatable,dimension(:,:,:) :: sal,rhsr,rusal,hsal + integer,allocatable,dimension(:) :: yc_to_ymr, zc_to_zmr + real,allocatable,dimension(:,:,:) :: rhsr real,allocatable,dimension(:,:) :: cxvx, cxvy, cxvz, cxrs, cxsalc, cxphic real,allocatable,dimension(:,:) :: cyvx, cyvy, cyvz, cyrs, cysalc, cyphic real,allocatable,dimension(:,:) :: czvx, czvy, czvz, czrs, czsalc, czphic - real,allocatable,dimension(:,:,:) :: vxr,vyr,vzr !CS mgrd + real,allocatable,dimension(:,:) :: cych, czch real,allocatable,dimension(:,:,:) :: tpdv,tpdvr !CS mgrd - real,allocatable,dimension(:,:,:) :: salc - real,allocatable,dimension(:,:,:) :: tempr,Tplaner - real,allocatable,dimension(:,:,:) :: phi,phic,ruphi,hphi + real,allocatable,dimension(:,:,:) :: Tplaner + real,allocatable,dimension(:,:) :: solid_height, height_vx, height_vy, height_vz end module mgrd_arrays !=============================================================== module stat_arrays @@ -178,6 +168,9 @@ module mpih integer :: ierr integer, parameter :: master=0 integer :: MDP = MPI_DOUBLE_PRECISION + integer :: comm_yz !> MPI communicator across y and z (2D decomposition) + integer :: comm_xy !> MPI communicator across x and y (1D decomposition) + integer :: comm_xz !> MPI communicator across x and z (1D decomposition) end module mpih !==================================================== module fftw_params diff --git a/src/flow_solver/GridModule.F90 b/src/grid.F90 similarity index 100% rename from src/flow_solver/GridModule.F90 rename to src/grid.F90 diff --git a/src/h5tools/HdfReadContinua.F90 b/src/h5tools/HdfReadContinua.F90 index 06ec824f..4659c0c6 100644 --- a/src/h5tools/HdfReadContinua.F90 +++ b/src/h5tools/HdfReadContinua.F90 @@ -53,6 +53,8 @@ subroutine HdfReadContinua(n1o,n2o,n3o,xs2,xe2,xs3,xe3,intvar,qua) filnam1 = trim('outputdir/continua_phi.h5') case (7) filnam1 = trim('outputdir/continua_pr.h5') + case (8) + filnam1 = trim('outputdir/continua_humid.h5') end select !RO Set offsets and element counts diff --git a/src/h5tools/HdfRoutines.F90 b/src/h5tools/HdfRoutines.F90 index 3f102461..aac3df7e 100644 --- a/src/h5tools/HdfRoutines.F90 +++ b/src/h5tools/HdfRoutines.F90 @@ -492,8 +492,11 @@ subroutine HdfWriteRealHalo3D(filnam1,qua) call h5pclose_f(plist_id, hdf_error) + ! call h5pcreate_f(H5P_DATASET_CREATE_F, plist_id, hdf_error) + ! call h5pset_chunk_f(plist_id, ndims, data_count, hdf_error) call h5dcreate_f(file_id, 'var', H5T_NATIVE_DOUBLE, & & filespace, dset, hdf_error) + ! call h5pclose_f(plist_id, hdf_error) call h5screate_simple_f(ndims, data_count, memspace, hdf_error) diff --git a/src/h5tools/MakeMovieXCut.F90 b/src/h5tools/MakeMovieXCut.F90 index ea3038c1..de690e8a 100644 --- a/src/h5tools/MakeMovieXCut.F90 +++ b/src/h5tools/MakeMovieXCut.F90 @@ -12,7 +12,9 @@ subroutine Mkmov_xcut use hdf5 use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X use local_arrays, only: vz,vy,vx,temp - use mgrd_arrays, only: sal,phi + use afid_salinity, only: sal + use afid_phasefield, only: phi + use afid_moisture, only: humid use h5_tools use param, only: nxm, nxmr, IBM implicit none @@ -32,7 +34,8 @@ subroutine Mkmov_xcut filename = trim("outputdir/flowmov/movie_xcut.h5") ! MPI - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true.,.true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true.,.true./), comm, ierr) + comm = comm_yz info = MPI_INFO_NULL call h5_open_or_create(file_id, filename, comm, fexist) @@ -50,6 +53,11 @@ subroutine Mkmov_xcut varname = 'temp' call write_H5_plane(file_id, varname, temp(ic, xstart(2):xend(2), xstart(3):xend(3)), 'x') + if (moist) then + varname = 'qhum' + call write_H5_plane(file_id, varname, humid(ic, xstart(2):xend(2), xstart(3):xend(3)), 'x') + end if + call h5fclose_f(file_id, hdf_error) if (salinity) then diff --git a/src/h5tools/MakeMovieYCut.F90 b/src/h5tools/MakeMovieYCut.F90 index 9ff28e1d..8596f5af 100644 --- a/src/h5tools/MakeMovieYCut.F90 +++ b/src/h5tools/MakeMovieYCut.F90 @@ -13,7 +13,9 @@ subroutine Mkmov_ycut use hdf5 use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X use local_arrays, only: vz,vy,vx,temp - use mgrd_arrays, only: sal,phi + use afid_salinity, only: sal + use afid_phasefield, only: phi + use afid_moisture, only: humid use h5_tools implicit none character(70) :: filename @@ -31,7 +33,8 @@ subroutine Mkmov_ycut filename = trim("outputdir/flowmov/movie_ycut.h5") ! MPI - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false.,.true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false.,.true./), comm, ierr) + comm = comm_xz info = MPI_INFO_NULL if (ic.le.xend(2) .and. ic.ge.xstart(2)) then @@ -51,6 +54,11 @@ subroutine Mkmov_ycut varname = 'temp' call write_H5_plane(file_id, varname, temp(1:nxm, ic, xstart(3):xend(3)), 'y') + if (moist) then + varname = 'qhum' + call write_H5_plane(file_id, varname, humid(1:nxm, ic, xstart(3):xend(3)), 'y') + end if + call h5fclose_f(file_id, hdf_error) end if diff --git a/src/h5tools/MakeMovieZCut.F90 b/src/h5tools/MakeMovieZCut.F90 index 6241befe..36e7a33f 100644 --- a/src/h5tools/MakeMovieZCut.F90 +++ b/src/h5tools/MakeMovieZCut.F90 @@ -13,7 +13,10 @@ subroutine Mkmov_zcut 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 mgrd_arrays, only: sal, phi, phic, tempr + use afid_salinity, only: sal + use afid_phasefield, only: phi, phic, tempr + use afid_moisture, only: humid use h5_tools implicit none character(70) :: filename @@ -31,7 +34,8 @@ subroutine Mkmov_zcut filename = trim("outputdir/flowmov/movie_zcut.h5") ! MPI - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true.,.false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true.,.false./), comm, ierr) + comm = comm_xy info = MPI_INFO_NULL if(ic.le.xend(3) .and. ic.ge.xstart(3)) then @@ -51,6 +55,11 @@ subroutine Mkmov_zcut varname = 'temp' call write_H5_plane(file_id, varname, temp(1:nxm, xstart(2):xend(2), ic), 'z') + if (moist) then + varname = 'qhum' + call write_H5_plane(file_id, varname, humid(1:nxm, xstart(2):xend(2), ic), 'z') + end if + call h5fclose_f(file_id, hdf_error) end if diff --git a/src/h5tools/h5_tools.F90 b/src/h5tools/h5_tools.F90 index 7348d366..33466dd5 100644 --- a/src/h5tools/h5_tools.F90 +++ b/src/h5tools/h5_tools.F90 @@ -1,7 +1,7 @@ module h5_tools use HDF5 use mpih, only: MPI_INFO_NULL - use param, only: salinity, phasefield, nx, nxm, nym, nzm, time, tframe, nxmr, nymr, nzmr, IBM + use param, only: salinity, phasefield, nx, nxm, nym, nzm, time, tframe, nxmr, nymr, nzmr, IBM, moist use decomp_2d, only: xstart, xstartr implicit none @@ -57,6 +57,11 @@ subroutine h5_add_slice_groups(file_id) call h5gcreate_f(file_id, "phi", group_id, hdf_error) call h5gclose_f(group_id, hdf_error) end if + + if (moist) then + call h5gcreate_f(file_id, "qhum", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + end if end subroutine subroutine write_H5_plane(file_id, varname, var, axis) @@ -130,4 +135,19 @@ subroutine write_H5_plane(file_id, varname, var, axis) end subroutine +!> Initialise the MPI communicators used for writing 2D slices +subroutine InitSliceCommunicators + use mpih, only: comm_xy, comm_xz, comm_yz + use decomp_2d, only: DECOMP_2D_COMM_CART_X + integer :: ierr + + !! yz-plane (cut of constant x) + call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .true./), comm_yz, ierr) + !! xy-plane (cut of constant z) + call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true.,.false./), comm_xy, ierr) + !! xz-plane (cut of constant y) + call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false.,.true./), comm_xz, ierr) + +end subroutine InitSliceCommunicators + end module h5_tools \ No newline at end of file diff --git a/src/h5tools/mean_zplane.F90 b/src/h5tools/mean_zplane.F90 index 4300b00f..3d329317 100644 --- a/src/h5tools/mean_zplane.F90 +++ b/src/h5tools/mean_zplane.F90 @@ -1,19 +1,21 @@ subroutine mean_zplane - use param, only: nx, nzm + use param, only: nx use mpih use hdf5 - use decomp_2d, only: xstart,xend,xsize,xstartr,xendr,DECOMP_2D_COMM_CART_X,nrank + use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X use local_arrays, only: vz,vy,vx,temp - use mgrd_arrays, only: sal, phi + use afid_salinity, only: sal + use afid_phasefield, only: phi + use afid_moisture, only: humid use h5_tools use means implicit none character(70) :: filename character(4) :: varname - real, allocatable :: uplane(:,:), vplane(:,:), wplane(:,:), Tplane(:,:) - integer :: comm, zrank, i, j, bufsize + real, allocatable :: uplane(:,:), vplane(:,:), wplane(:,:), Tplane(:,:), qplane(:,:) + integer :: comm, zrank integer(HID_T) :: file_id logical :: fexist @@ -23,16 +25,20 @@ subroutine mean_zplane allocate(vplane(1:nxm,xstart(2):xend(2))) allocate(wplane(1:nxm,xstart(2):xend(2))) allocate(Tplane(1:nxm,xstart(2):xend(2))) + if (moist) allocate(qplane(1:nxm,xstart(2):xend(2))) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + comm = comm_xz call MPI_COMM_RANK(comm, zrank, ierr) call zmean(vx(1:nx ,xstart(2):xend(2),xstart(3):xend(3)), uplane, comm, zrank) call zmean(vy(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), vplane, comm, zrank) call zmean(vz(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), wplane, comm, zrank) call zmean(temp(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), Tplane, comm, zrank) + if (moist) call zmean(humid(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), qplane, comm, zrank) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + comm = comm_xy if (zrank==0) then @@ -49,6 +55,10 @@ subroutine mean_zplane call write_H5_plane(file_id, varname, wplane, 'z') varname='temp' call write_H5_plane(file_id, varname, Tplane, 'z') + if (moist) then + varname='qhum' + call write_H5_plane(file_id, varname, qplane, 'z') + end if ! Close HDF5 file call h5fclose_f(file_id, hdf_error) @@ -58,16 +68,19 @@ subroutine mean_zplane if (allocated(vplane)) deallocate(vplane) if (allocated(wplane)) deallocate(wplane) if (allocated(Tplane)) deallocate(Tplane) + if (allocated(qplane)) deallocate(qplane) if (salinity) then allocate(Tplane(1:nxmr,xstartr(2):xendr(2))) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + comm = comm_xz call MPI_COMM_RANK(comm, zrank, ierr) call zmeanr(sal(1:nxmr,xstartr(2):xendr(2),xstartr(3):xendr(3)), Tplane, comm, zrank) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + comm = comm_xy if (zrank==0) then ! Open the movie file (create if it doesn't exist) @@ -88,12 +101,14 @@ subroutine mean_zplane if (phasefield) then allocate(Tplane(1:nxmr,xstartr(2):xendr(2))) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + comm = comm_xz call MPI_COMM_RANK(comm, zrank, ierr) call zmeanr(phi(1:nxmr,xstartr(2):xendr(2),xstartr(3):xendr(3)), Tplane, comm, zrank) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + comm = comm_xy if (zrank==0) then ! Open the movie file (create if it doesn't exist) @@ -114,20 +129,22 @@ subroutine mean_zplane end subroutine mean_zplane subroutine mean_yplane - use param, only: nx, nym + use param, only: nx use mpih use hdf5 - use decomp_2d, only: xstart,xend,xsize,xstartr,xendr,DECOMP_2D_COMM_CART_X,nrank + use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X use local_arrays, only: vz,vy,vx,temp - use mgrd_arrays, only: sal, phi + use afid_salinity, only: sal + use afid_phasefield, only: phi + use afid_moisture, only: humid use h5_tools use means implicit none character(70) :: filename character(4) :: varname - real, allocatable :: uplane(:,:), vplane(:,:), wplane(:,:), Tplane(:,:) - integer :: comm, yrank, i, j, bufsize + real, allocatable :: uplane(:,:), vplane(:,:), wplane(:,:), Tplane(:,:), qplane(:,:) + integer :: comm, yrank integer(HID_T) :: file_id logical :: fexist @@ -137,16 +154,20 @@ subroutine mean_yplane allocate(vplane(1:nxm,xstart(3):xend(3))) allocate(wplane(1:nxm,xstart(3):xend(3))) allocate(Tplane(1:nxm,xstart(3):xend(3))) + if (moist) allocate(qplane(1:nxm,xstart(3):xend(3))) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + comm = comm_xy call MPI_COMM_RANK(comm, yrank, ierr) call ymean(vx(1:nx ,xstart(2):xend(2),xstart(3):xend(3)), uplane, comm, yrank) call ymean(vy(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), vplane, comm, yrank) call ymean(vz(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), wplane, comm, yrank) call ymean(temp(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), Tplane, comm, yrank) + if (moist) call ymean(humid(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), qplane, comm, yrank) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + comm = comm_xz if (yrank==0) then @@ -163,6 +184,10 @@ subroutine mean_yplane call write_H5_plane(file_id, varname, wplane, 'y') varname='temp' call write_H5_plane(file_id, varname, Tplane, 'y') + if (moist) then + varname='qhum' + call write_H5_plane(file_id, varname, qplane, 'y') + end if ! Close HDF5 file call h5fclose_f(file_id, hdf_error) @@ -172,16 +197,19 @@ subroutine mean_yplane if (allocated(vplane)) deallocate(vplane) if (allocated(wplane)) deallocate(wplane) if (allocated(Tplane)) deallocate(Tplane) + if (allocated(qplane)) deallocate(qplane) if (salinity) then allocate(Tplane(1:nxmr,xstartr(3):xendr(3))) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + comm = comm_xy call MPI_COMM_RANK(comm, yrank, ierr) call ymeanr(sal(1:nxmr,xstartr(2):xendr(2),xstartr(3):xendr(3)), Tplane, comm, yrank) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + comm = comm_xz if (yrank==0) then ! Open the movie file (create if it doesn't exist) @@ -202,12 +230,14 @@ subroutine mean_yplane if (phasefield) then allocate(Tplane(1:nxmr,xstartr(3):xendr(3))) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.true., .false./), comm, ierr) + comm = comm_xy call MPI_COMM_RANK(comm, yrank, ierr) call ymeanr(phi(1:nxmr,xstartr(2):xendr(2),xstartr(3):xendr(3)), Tplane, comm, yrank) - call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + ! call MPI_CART_SUB(DECOMP_2D_COMM_CART_X, (/.false., .true./), comm, ierr) + comm = comm_xz if (yrank==0) then ! Open the movie file (create if it doesn't exist) diff --git a/src/ibm/IBMTools.F90 b/src/ibm/IBMTools.F90 new file mode 100644 index 00000000..53ad72e8 --- /dev/null +++ b/src/ibm/IBMTools.F90 @@ -0,0 +1,185 @@ +module IBMTools + use decomp_2d + use param + use mgrd_arrays + +contains + +!! May have to make h(1,:,:) so that we can update halo +!! Or alternatively, include halo in loop, ensuring that phi has an up-to-date halo +subroutine calc_interface_height(ph, h) + real, intent(in) :: ph(1:,xstartr(2)-lvlhalo:,xstartr(3)-lvlhalo:) + !< input variable + real, intent(out) :: h(xstartr(2)-lvlhalo:,xstartr(3)-lvlhalo:) + real :: dxx(nxmr) + integer :: i,j,k,kp + + h(:,:) = 0.0 + + dxx = d3xcr(1:nxmr)/dxr + do i=xstartr(3)-lvlhalo,xendr(3)+lvlhalo + do j=xstartr(2)-lvlhalo,xendr(2)+lvlhalo + do k=1,nxmr + h(j,i) = h(j,i) + ph(k,j,i)*dxx(k) + ! kp = k + 1 + ! if (ph(k,j,i) <= 0.5 .and. ph(kp,j,i) > 0.5) then + ! h(j,i) = xmr(k) + (xmr(kp) - xmr(k))*(0.5 - ph(k,j,i))/(ph(kp,j,i) - ph(k,j,i)) + ! end if + end do + end do + end do + +end subroutine calc_interface_height + +subroutine interp_height_to_vel_grid(hr, hcx, hcy, hcz) + real, intent(in) :: hr(xstartr(2)-lvlhalo:,xstartr(3)-lvlhalo:) + real, intent(out) :: hcx(xstart(2):,xstart(3):) + real, intent(out) :: hcy(xstart(2):,xstart(3):) + real, intent(out) :: hcz(xstart(2):,xstart(3):) + + real :: hv2(4,4) + real :: hv1(4) + integer :: ic, jc, icr, jcr, icc, jcc + + do ic=xstart(3),xend(3) + icr = krangs(ic) + icc = zc_to_zmr(ic) + do jc=xstart(2),xend(2) + jcr = jrangs(jc) + jcc = yc_to_ymr(jc) + + hv2 = hr(jcr-2:jcr+1,icr-2:icr+1) + hv1 = czphic(1,ic)*hv2(:,1) + czphic(2,ic)*hv2(:,2) & + + czphic(3,ic)*hv2(:,3) + czphic(4,ic)*hv2(:,4) + hcx(jc,ic) = sum(cyphic(1:4,jc)*hv1) + + hcy(jc,ic) = sum(cych(1:4,jc)*hv1) + + hv1 = czch(1,ic)*hv2(:,1) + czch(2,ic)*hv2(:,2) & + + czch(3,ic)*hv2(:,3) + czch(4,ic)*hv2(:,4) + hcz(jc,ic) = sum(cyphic(1:4,jc)*hv1) + end do + end do + +end subroutine + + +!> Takes an input height profile `h` (function of y and z) +!> and outputs a 3D integer array `ibmask` equal to 2 in liquid +!> and equal to 0 in the solid phase +!> grd should be 'x', 'y', or 'z' to denote whether the profile +!> is on the grid for vx, vy, or vz +!> ASSUMES SOLID BELOW LIQUID +subroutine mask_below_height(h, ibmask, grd) + real, intent(in) :: h(xstart(2):,xstart(3):) + integer, intent(out) :: ibmask(1:,xstart(2):,xstart(3):) + character, intent(in) :: grd + + real :: ze, ye, xe + integer :: i,j,k + + ibmask(:,:,:) = 2 + + do i=xstart(3),xend(3) + ze = zm(i) + if (grd=='z') ze = zc(i) + do j=xstart(2),xend(2) + ye = ym(j) + if (grd=='y') ye = yc(j) + do k=1,nxm + xe = xm(k) + if (grd=='x') xe = xc(k) + if (xe < h(j,i)) then + ibmask(k,j,i) = 0 + end if + end do + end do + end do + +end subroutine mask_below_height + +!> Takes an input height profile `h` (function of y and z) +!> and outputs a 3D integer array `ibmask` equal to 2 in liquid +!> and equal to 0 in the solid phase +!> grd should be 'x', 'y', or 'z' to denote whether the profile +!> is on the grid for vx, vy, or vz +!> ASSUMES LIQUID BELOW SOLID +subroutine mask_above_height(h, ibmask, grd) + real, intent(in) :: h(xstart(2):,xstart(3):) + integer, intent(out) :: ibmask(1:,xstart(2):,xstart(3):) + character, intent(in) :: grd + + real :: ze, ye, xe + integer :: i,j,k + + ibmask(:,:,:) = 2 + + do i=xstart(3),xend(3) + ze = zm(i) + if (grd=='z') ze = zc(i) + do j=xstart(2),xend(2) + ye = ym(j) + if (grd=='y') ye = yc(j) + do k=1,nxm + xe = alx3 - xm(k) + if (grd=='x') xe = alx3 - xc(k) + if (xe < h(j,i)) then + ibmask(k,j,i) = 0 + end if + end do + end do + end do + +end subroutine mask_above_height + +subroutine calc_IBM_interpolation(h, ibmask, dist, grd) + real, intent(in) :: h(xstart(2):,xstart(3):) + integer, intent(inout) :: ibmask(1:,xstart(2):,xstart(3):) + real, intent(out) :: dist(:) + character, intent(in) :: grd + + real :: xe,ye,ze + integer :: i,j,k,km,kp,n + + n = 0 + + do i=xstart(3),xend(3) + ze = zm(i) + if (grd=='z') ze = zc(i) + do j=xstart(2),xend(2) + ye = ym(j) + if (grd=='y') ye = yc(j) + do k=1,nxm + km = kmv(k) + kp = kpv(k) + xe = xm(k) + xem = xm(km) + xep = xm(kp) + if (grd=='x') then + xe = xc(k) + xem = xc(km) + xep = xc(kp) + end if + + ! Liquid over solid + if ((ibmask(k,j,i)==2) .and. (ibmask(km,j,i)==0)) then + n = n + 1 + d1x = xep - xe + d2x = xe - h(j,i) + dist(n) = d2x/(d1x + d2x) + ibmask(k,j,i) = 1 + + ! Liquid under solid + elseif ((ibmask(k,j,i)==2) .and. (ibmask(kp,j,i)==0)) then + n = n + 1 + d1x = xe - xem + d2x = alx3 - h(j,i) - xe + dist(n) = d2x/(d1x + d2x) + ibmask(k,j,i) = -1 + end if + end do + end do + end do +end subroutine calc_IBM_interpolation + +end module IBMTools \ No newline at end of file diff --git a/src/ibm/SolveImpEqnUpdate_Sal_ibm.F90 b/src/ibm/SolveImpEqnUpdate_Sal_ibm.F90 index 7837f166..9c568731 100644 --- a/src/ibm/SolveImpEqnUpdate_Sal_ibm.F90 +++ b/src/ibm/SolveImpEqnUpdate_Sal_ibm.F90 @@ -10,12 +10,13 @@ subroutine SolveImpEqnUpdate_Sal_ibm use param - use mgrd_arrays, only: sal,rhsr + use mgrd_arrays, only: rhsr + use afid_salinity, only: sal, PecS, ap3sskr, ac3sskr, am3sskr use decomp_2d, only: xstartr,xendr - use ibm_param, only: ibmaskr, distbr + use ibm_param, only: ibmaskr implicit none real, dimension(nxr) :: amkl,apkl,ackl,fkl - integer :: jc,kc,info,ipkv(nxmr),ic,nrhs,km,kp,n + integer :: jc,kc,info,ipkv(nxmr),ic,km,kp,n real :: betadx,ackl_b real :: amkT(nxmr-1),ackT(nxmr),apkT(nxmr-1),appk(nxmr-2) diff --git a/src/ibm/SolveImpEqnUpdate_Temp_ibm.F90 b/src/ibm/SolveImpEqnUpdate_Temp_ibm.F90 index 03a6b584..ac7f4b27 100644 --- a/src/ibm/SolveImpEqnUpdate_Temp_ibm.F90 +++ b/src/ibm/SolveImpEqnUpdate_Temp_ibm.F90 @@ -10,13 +10,13 @@ subroutine SolveImpEqnUpdate_Temp_ibm use param - use local_arrays, only : temp,rhs,hro + use local_arrays, only : temp,rhs use decomp_2d, only: xstart,xend use ibm_param, only: distbt, ibmaskt ! use param_particle ! SL implicit none real, dimension(nx) :: amkl,apkl,ackl,fkl - integer :: jc,kc,info,ipkv(nxm),ic,nrhs,km,kp,n + integer :: jc,kc,info,ipkv(nxm),ic,km,kp,n real :: betadx,ackl_b, Tb real :: amkT(nxm-1),ackT(nxm),apkT(nxm-1),appk(nxm-2) ! real :: kcpp,kcpm ! SL diff --git a/src/ibm/ibm_param.F90 b/src/ibm/ibm_param.F90 index 8fb5ff96..53a29d7a 100644 --- a/src/ibm/ibm_param.F90 +++ b/src/ibm/ibm_param.F90 @@ -13,52 +13,4 @@ module ibm_param real,dimension(mpun) :: distx, disty, distz real,dimension(mpun) :: distbt real,dimension(mpunr) :: distbr - ! real,dimension(mpun) :: temb,salfix - ! real,dimension(mpun) :: q1bo,q2bo,q3bo,densb,salb - end module ibm_param - -! SL ======================================================= -! module param_particle -! implicit none -! integer :: Nparticle,NL,Ne,npp,nll,nee,partshape,thermal_couple,iangle -! integer :: imonitor,MonitorNumber,iJointPDF,HSectPosi -! integer, allocatable, dimension(:,:) :: MonitorPosition -! integer, allocatable, dimension(:,:,:) :: allsdind -! real :: radius,diameter,kparticle,friction,friction2,TOUTHSect,HSectCoor,rhocpparticle -! real :: dVL,rx,ry,Am,delta,angular_velocity -! real, allocatable, dimension(:) :: vec_angle,part_surf_vel_x,part_surf_vel_y -! real, allocatable, dimension(:) :: HSectTemp,HSectVx,HSectVy -! real, allocatable, dimension(:,:) :: ParticleCenter -! real, allocatable, dimension(:,:,:) :: phaindx,phaindy,phaindz,kcp,insidex,insidey,insidez,rhocpcp -! real, allocatable, dimension(:,:,:) :: vxcp,vycp,vzcp,vxcp_mean,vycp_mean -! real, allocatable, dimension(:,:,:) :: vxc_ibm,vyc_ibm,vzc_ibm -! real, allocatable, dimension(:,:) :: lagnodx,lagnody,ffx,ffy -! real, allocatable, dimension(:,:,:) :: mls_transfer -! real, allocatable, dimension(:,:) :: mls_h,mls_c -! ! real, allocatable, dimension(:,:,:) :: delta_vx_index,delta_vy_index,delta_vz_index -! end module param_particle - - -! module param_tracer -! implicit none -! real, allocatable, dimension(:,:,:) :: vxc,vyc,vzc,tempc -! real, allocatable, dimension(:,:,:) :: str_11,str_22,str_33 -! real, allocatable, dimension(:,:,:) :: str_12,str_13,str_23 -! real, allocatable, dimension(:,:,:) :: vor_y,vor_z,vor_x - -! real, allocatable, dimension(:,:) :: xp,xpo -! real, allocatable, dimension(:) :: vxp,vyp,vzp,Tp -! real, allocatable, dimension(:) :: vxop,vyop,vzop - -! real, allocatable, dimension(:) :: str11p,str22p,str33p -! real, allocatable, dimension(:) :: str12p,str13p,str23p -! real, allocatable, dimension(:) :: voryp,vorzp,vorxp - -! integer iftracer,Ntracer,ONtracer,tracer_write,tracer_read,if_tracer_pair -! real tracer_init_sep -! real x_i,y_i,z_i,x_f,y_f,z_f,z_c -! real timeONtracer,TOUT_tracer -! integer dimen -! end module param_tracer -! ========================================================== diff --git a/src/ibm/invtr1_ibm.F90 b/src/ibm/invtr1_ibm.F90 deleted file mode 100644 index fb0c42ab..00000000 --- a/src/ibm/invtr1_ibm.F90 +++ /dev/null @@ -1,22 +0,0 @@ -!RO Extra instructions in the case of Immersed boundary for invtr1 - forclo=1.d0 - - if(infig.gt.0) then - usaldto = 1./aldto - do n=1,npunz - ic=indgeo(1,n,1) - jc=indgeo(1,n,2) - kc=indgeo(1,n,3) - forclo(kc,jc,ic)=0.d0 -! ie=indgeoe(1,n,1) -! je=indgeoe(1,n,2) - ke=indgeoe(1,n,3) - q1e=((al*dt+aldto)*vz(ke,jc,ic)-al*dt*q1bo(n))*usaldto - rhs(kc,jc,ic) = -vz(kc,jc,ic) + q1e*distb(1,n) - q1bo(n)= vz(ke,jc,ic) - end do - endif - - - call SolveImpEqnUpdate_YZ_ibm(vz,rhs,forclo) - diff --git a/src/ibm/invtr2_ibm.F90 b/src/ibm/invtr2_ibm.F90 deleted file mode 100644 index 81e22198..00000000 --- a/src/ibm/invtr2_ibm.F90 +++ /dev/null @@ -1,21 +0,0 @@ -!RO Extra instructions in the case of Immersed boundary for invtr2 - - forclo=1.d0 - if(infig.gt.0) then - usaldto = 1./aldto - do n=1,npuny - ic=indgeo(2,n,1) - jc=indgeo(2,n,2) - kc=indgeo(2,n,3) - forclo(kc,jc,ic)=0.d0 -! ie=indgeoe(2,n,1) -! je=indgeoe(2,n,2) - ke=indgeoe(2,n,3) - q2e=((al*dt+aldto)*vy(ke,jc,ic)-al*dt*q2bo(n))*usaldto - rhs(kc,jc,ic) = -vy(kc,jc,ic) + q2e*distb(2,n) - q2bo(n)= vy(ke,jc,ic) - end do - endif - - call SolveImpEqnUpdate_YZ_ibm(vy,rhs,forclo) - diff --git a/src/ibm/invtr3_ibm.F90 b/src/ibm/invtr3_ibm.F90 deleted file mode 100644 index 8c102d0e..00000000 --- a/src/ibm/invtr3_ibm.F90 +++ /dev/null @@ -1,18 +0,0 @@ - forclo = 1.d0 - if(infig.gt.0) then - usaldto = 1./aldto - do n=1,npunx - ic=indgeo(3,n,1) - jc=indgeo(3,n,2) - kc=indgeo(3,n,3) - forclo(kc,jc,ic)=0.d0 - ! ie=indgeoe(3,n,1) - ! je=indgeoe(3,n,2) - ke=indgeoe(3,n,3) - q3e=((al*dt+aldto)*vx(ke,jc,ic)-al*dt*q3bo(n))*usaldto - rhs(kc,jc,ic) = -vx(kc,jc,ic) + q3e*distb(3,n) - q3bo(n)= vx(ke,jc,ic) - end do - end if - - call SolveImpEqnUpdate_X_ibm diff --git a/src/ibm/invtrro_ibm.F90 b/src/ibm/invtrro_ibm.F90 deleted file mode 100644 index 3a289754..00000000 --- a/src/ibm/invtrro_ibm.F90 +++ /dev/null @@ -1,37 +0,0 @@ - forclo=1.d0 - if(infig.gt.0) then - usaldto = 1./aldto - do n=1,npunte - ic=indgeot(n,1) - jc=indgeot(n,2) - kc=indgeot(n,3) - forclo(kc,jc,ic)=0.d0 -! ie=indgeoet(n,1) -! je=indgeoet(n,2) - ke=indgeoet(n,3) - dense=((al*dt+aldto)*temp(ke,jc,ic)-al*dt*densb(n))*usaldto - rhs(kc,jc,ic) = -temp(kc,jc,ic) + dense*distbt(n) & - +(1.-distbt(n))*temb(n) - densb(n)= temp(ke,jc,ic) - - if(ifnoslipy.eq.1) then ! SL - if (jc==1) then - rhs(kc,jc,ic)=temp(kc,jc+1,ic)-temp(kc,jc,ic) - elseif (jc==nym) then - rhs(kc,jc,ic)=temp(kc,jc-1,ic)-temp(kc,jc,ic) - endif - endif - - if(ifnoslipz.eq.1) then ! SL - if (ic==1) then - rhs(kc,jc,ic)=temp(kc,jc,ic+1)-temp(kc,jc,ic) - elseif (ic==nzm) then - rhs(kc,jc,ic)=temp(kc,jc,ic-1)-temp(kc,jc,ic) - endif - endif - - end do - end if - - - call SolveImpEqnUpdate_Temp_ibm diff --git a/src/ibm/topogr_ibm.F90 b/src/ibm/topogr_ibm.F90 index 31f81d92..825bf1ea 100644 --- a/src/ibm/topogr_ibm.F90 +++ b/src/ibm/topogr_ibm.F90 @@ -6,26 +6,23 @@ subroutine topogr use param use decomp_2d, only: xstart,xend,xstartr,xendr use ibm_param - use mgrd_arrays, only: sal + use afid_salinity, only: RayS use mpih implicit none - integer :: i,j,k,l,kstartp, nc, Npart, ncz - integer :: km,kp,jm,jp,im,ip,mm!,kup,klo + integer :: i,j,k,l, nc, Npart, ncz + integer :: km,kp - real :: xe, xem, xep!, xlo, xup - real :: ye, yem, yep - real :: ze, zem, zep + real :: xe, xem, xep + real :: ye, yem + real :: ze, zem real :: delta1x, delta2x, r2, Lhex, radius, porosity real :: solid_temp, rp, tp, amp, rx, ry, rz - integer,allocatable :: ind1(:), ind2(:) real,allocatable :: xpart(:), ypart(:), zpart(:) integer :: ibmask(1:nx,xstart(2):xend(2),xstart(3):xend(3)) !! Variables for writing details of solid centres to file character(len=30) :: dsetname, filename - ! Flag for boundary interpolation for velocity - logical, parameter :: velBCinterp = .true. logical :: fexist @@ -628,6 +625,7 @@ subroutine topogr ! n = 0 ibmaskr(:,:,:) = 2 + solidr(:,:,:) = .false. do i=xstartr(3)-1,xendr(3)+1 ze = zmr(i) diff --git a/src/main.F90 b/src/main.F90 index 45b386bb..063be8cf 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -2,30 +2,40 @@ program AFiD use mpih use param use local_arrays, only: vx,vy,vz,temp,pr - use mgrd_arrays, only: vxr,vyr,vzr,salc,sal,phic,tempr,phi + ! use mgrd_arrays, only: phic,tempr,phi!,vxr,vyr,vzr,salc,sal use AuxiliaryRoutines use hdf5 use decomp_2d use decomp_2d_fft + use afid_pressure + use afid_moisture + use afid_salinity + use afid_phasefield + use h5_tools, only: InitSliceCommunicators ! use stat_arrays, only: nstatsamples,vx_global,vy_global,vz_global !$ use omp_lib implicit none - integer :: errorcode, nthreads, i, j, k + integer :: errorcode!, nthreads, i, j, k real :: instCFL,CFLmr,dmax,dmaxr real :: ti(2), tin(3), minwtdt - real :: ts, varptb,chksum + real :: ts!, varptb,chksum + real :: td(2) !< debugging time measure integer :: prow=0,pcol=0 integer :: lfactor,lfactor2 character(100) :: arg logical :: nanexist, write_mean_planes=.true. - real,allocatable,dimension(:,:) :: dummy,dscan,dbot - integer :: comm,ierror,row_id,row_coords(2),ic,jc,kc + ! real,allocatable,dimension(:,:) :: dummy,dscan,dbot + ! integer :: comm,ierror,row_id,row_coords(2),ic,jc,kc !******************************************************* !******* Read input file bou.in by all processes******** !******************************************************* ! + + ! Set `moist` to .true. if humid.in exists + inquire(file="humid.in", exist=moist) + call ReadInputFile if (nzm==1 .or. nym==1) write_mean_planes = .false. @@ -85,9 +95,13 @@ program AFiD call InitTimeMarchScheme call InitVariables + call InitPressureVars if (multires) call InitMgrdVariables !CS mgrd if (salinity) call InitSalVariables if (phasefield) call InitPFVariables + if (moist) call InitMoistVariables + + call InitSliceCommunicators call CreateGrid if (multires) call CreateMgrdGrid !CS mgrd @@ -136,6 +150,7 @@ program AFiD call InitPressureSolver call SetTempBCs if (salinity) call SetSalBCs + if (moist) call SetHumidityBCs if(readflow) then @@ -152,12 +167,22 @@ program AFiD instCFL=0.d0 call CreateInitialConditions - if (salinity) call CreateICSal - if (phasefield) call CreateICPF + if (salinity) call CreateInitialSalinity + if (phasefield) call CreateInitialPhase + if (moist) call CreateInitialHumidity endif - if (IBM) call topogr +!CS Create multigrid stencil for interpolation + if (multires) call CreateMgrdStencil + if (phasefield .and. IBM) call CreatePFStencil + + if (phasefield) call update_halo(phi,lvlhalo) + + if (IBM) then + call topogr + ! if (phasefield) call UpdateIBMLocation + end if !EP Update all relevant halos call update_halo(vx,lvlhalo) @@ -165,20 +190,19 @@ program AFiD call update_halo(vz,lvlhalo) call update_halo(temp,lvlhalo) if (salinity) call update_halo(sal,lvlhalo) - if (phasefield) call update_halo(phi,lvlhalo) call update_halo(pr,lvlhalo) + if (moist) call update_halo(humid,lvlhalo) + if (moist) call UpdateSaturation -!CS Create multigrid stencil for interpolation - if (multires) call CreateMgrdStencil !CS Interpolate initial values if (salinity) then call InterpVelMgrd - call InterpSalMgrd + call InterpSalMultigrid end if if (phasefield) then - call InterpTempMgrd - call InterpPhiMgrd + call InterpTempMultigrid + call InterpPhiMultigrid end if !EP Update all relevant halos @@ -205,7 +229,12 @@ program AFiD end if if (ismaster) write(*,*) "Writing 3D fields" + call MpiBarrier + td(1) = MPI_WTIME() call WriteFlowField(.false.) + call MpiBarrier + td(2) = MPI_WTIME() + if (ismaster) write(*,*) "Flow field writing took: ",td(2)-td(1) !EP Check divergence. Should be reduced to machine precision after the first !phcalc. Here it can still be high. diff --git a/src/moisture.F90 b/src/moisture.F90 new file mode 100644 index 00000000..c5af3885 --- /dev/null +++ b/src/moisture.F90 @@ -0,0 +1,537 @@ +!> Module adding evolution of a specific humidity field to AFiD +!! following Vallis et al. (2019) J. Fluid Mech. +module afid_moisture + use param + use decomp_2d, only: xstart, xend, nrank, update_halo + use AuxiliaryRoutines + use local_arrays, only: temp, vx, vy, vz, rhs, hro + implicit none + + real, allocatable, dimension(:,:,:) :: humid !! Specific humidity field + real, allocatable, dimension(:,:,:) :: qsat !! Saturation vapour field + real, allocatable, dimension(:,:,:) :: rkhumid !! RK storage array for humidity + real, allocatable, dimension(:,:,:) :: hkhumid !! RK storage array for humidity (previous substep) + + real :: gamma_q !! Coefficient for latent heat effect of condensing vapour + real :: tau_q !! Time-scale for condensation + real :: kappa_q !! Diffusivity of humidity field + real :: alpha_q !! Saturation coefficient + real :: beta_q !! Lapse rate coefficient + real :: pecq !! Peclet number for humidity (inverse dimensionless diffusivity) + + integer :: qfixS !! Flag for whether humidity takes fixed value at lower plate + integer :: qfixN !! Flag for whether humidity takes fixed value at upper plate + + real, allocatable, dimension(:,:,:) :: humbp !! Humidity boundary value (bottom plate) + real, allocatable, dimension(:,:,:) :: humtp !! Humidity boundary value (top plate) + +contains + +!> Subroutine to allocate memory for all moisture-related variables +subroutine InitMoistVariables + ! Allocate array for 3D specific humidity field (needs to be nx size in x for writing out continua file) + call AllocateReal3DArray(humid,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + ! Allocate array for 3D saturation values + call AllocateReal3DArray(qsat,1,nxm,xstart(2),xend(2),xstart(3),xend(3)) + + ! Runge-Kutta storage arrays + call AllocateReal3DArray(rkhumid,1,nxm,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(hkhumid,1,nxm,xstart(2),xend(2),xstart(3),xend(3)) + + call AllocateReal3DArray(humbp,1,1,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(humtp,1,1,xstart(2),xend(2),xstart(3),xend(3)) + + ! Default parameter values to match Vallis et al (2019) examples + + beta_q = 1.2 + alpha_q = 3.0 + pecq = pect ! Set humid Peclet equal to thermal + kappa_q = 1.0/pecq + gamma_q = 0.19 + tau_q = 5e-5*pect ! Vallis et al use diffusive scaling for nondimensionalisation + qfixN = 1 + qfixS = 1 + + ! Check if there is an extra input file, and update parameters + call ReadMoistParameters + + dtmax = min(dtmax, 0.1*tau_q) + if (ismaster) write(*,*) 'tau, dtmax: ', tau_q, dtmax + +end subroutine InitMoistVariables + +!> Read parameters from the humid.in input file +subroutine ReadMoistParameters + logical :: exists, ifdiff + integer :: io + real :: Sm !! humid-to-thermal diffusivity ratio + + inquire(file="humid.in", exist=exists) + if (exists) then + open(file="humid.in", newunit=io, status="old", action="read") + read(io,*) + read(io,*) + read(io,*) + read(io,*) + read(io,*) alpha_q + read(io,*) + read(io,*) + read(io,*) beta_q + read(io,*) + read(io,*) + read(io,*) gamma_q + read(io,*) + read(io,*) + read(io,*) tau_q + read(io,*) + read(io,*) + read(io,*) ifdiff + read(io,*) + read(io,*) + read(io,*) + read(io,*) qfixN, qfixS + read(io,*) + read(io,*) + read(io,*) Sm + close(io) + pecq = pect/Sm + kappa_q = 1.0/pecq + if (ifdiff) tau_q = pect*tau_q + end if + + ! If the input file has a negative gamma, set the value of gamma + ! such that delta m is equal to 1 + if (gamma_q < 0) then + if (qfixN == 2) then ! (unsaturated top boundary) + gamma_q = beta_q + else + gamma_q = beta_q/(1.0 - exp(-alpha_q)) + end if + end if + + if (ismaster) then + write(*,*) 'al, be, ga, tau' + write(*,*) alpha_q, beta_q, gamma_q, tau_q + end if + +end subroutine + +!> Deallocate the variables used for Rainy-Benard model +subroutine DeallocateMoistVariables + + call DestroyReal3DArray(humid) + call DestroyReal3DArray(qsat) + call DestroyReal3DArray(rkhumid) + call DestroyReal3DArray(hkhumid) + +end subroutine DeallocateMoistVariables + +subroutine SetHumidityBCs + integer :: ic, jc + + ! NOTE: tempbp = 0 => humbp = 1 + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + humbp(1,jc,ic) = exp(alpha_q*tempbp(1,jc,ic)) + humtp(1,jc,ic) = exp(alpha_q*(temptp(1,jc,ic) - beta_q)) + end do + end do + ! Use the input qfixN==2 to set an unsaturated top boundary + if (qfixN == 2) then + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + humtp(1,jc,ic) = 0.0 + end do + end do + qfixN = 1 + end if +end subroutine SetHumidityBCs + +subroutine CreateInitialHumidity + use mpih + integer :: ic, jc, kc + real :: rnum, r, r2, amp + real :: bz(nxm), qz(nxm) + logical :: exists + character(len=30) :: dsetname, filename + + call random_seed() + + filename = trim("drizzle.h5") + + inquire(file=filename, exist=exists) + if (exists) then + if (ismaster) then + dsetname = trim("b") + call HdfSerialReadReal1D(dsetname, filename, bz, nxm) + dsetname = trim("q") + call HdfSerialReadReal1D(dsetname, filename, qz, nxm) + end if + call MPI_BCAST(bz, nxm, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(qz, nxm, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + amp = 1e-3 + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + call random_number(rnum) + temp(kc,jc,ic) = bz(kc) + amp*rnum + call random_number(rnum) + humid(kc,jc,ic) = qz(kc) + amp*rnum + end do + end do + end do + else + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + call random_number(rnum) + temp(kc,jc,ic) = (1.0 - xm(kc))*xm(kc)*rnum*1e-3 + end do + end do + end do + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + call random_number(rnum) + ! r = sqrt((ym(jc) - 0.5*ylen)**2 + (xm(kc) - 0.1)**2 + (zm(ic) - 0.5*zlen)**2) + r2 = (ym(jc) - 0.5*ylen)**2 + (xm(kc) - 0.1)**2 + (zm(ic) - 0.5*zlen)**2 + ! humid(kc,jc,ic) = humbp(1,jc,ic) + (humtp(1,jc,ic) - humbp(1,jc,ic))*xm(kc) + ! humid(kc,jc,ic) = 1.1*qsat(kc,jc,ic)*0.5*(1.0 - tanh(100*(r - 0.1))) + ! humid(kc,jc,ic) = 0.5*(1.0 - tanh(100*(r - 0.1))) + humid(kc,jc,ic) = 5.0*exp(-r2/0.005) + humid(kc,jc,ic) = humid(kc,jc,ic) + 1e-3*rnum + ! humid(kc,jc,ic) = temp(kc,jc,ic) + end do + end do + end do + end if + + call update_halo(temp,lvlhalo) + call UpdateSaturation + call update_halo(humid,lvlhalo) + +end subroutine CreateInitialHumidity + +!> +!! Compute the explicit terms for the specific humidity equation +!! and store the result in rkhumid +subroutine ExplicitHumidity + integer :: ip, jp + integer :: ic, jc, kc + integer :: im, jm + real :: udy, udz, itau + real :: dyyq, dzzq, condensation, hqx, hqy, hqz + + udy = 0.5*dy + udz = 0.5*dz + + itau = 1.0/tau_q + + do ic=xstart(3),xend(3) + im = ic-1 + ip = ic+1 + do jc=xstart(2),xend(2) + jm = jc-1 + jp = jc+1 + do kc=1,nxm + + ! x-advection d/dx(vx * q) + if (kc==1) then + hqx = ( & + vx(kc+1,jc,ic)*(humid(kc+1,jc,ic) + humid(kc,jc,ic)) & + - vx(kc ,jc,ic)*2.0*humbp(1,jc,ic) & + )*0.5*udx3m(kc) + elseif (kc==nxm) then + hqx = ( & + vx(kc+1,jc,ic)*2.0*humtp(1,jc,ic) & + - vx(kc ,jc,ic)*(humid(kc,jc,ic) + humid(kc-1,jc,ic)) & + )*0.5*udx3m(kc) + else + hqx = ( & + vx(kc+1,jc,ic)*(humid(kc+1,jc,ic) + humid(kc ,jc,ic)) & + - vx(kc ,jc,ic)*(humid(kc ,jc,ic) + humid(kc-1,jc,ic)) & + )*0.5*udx3m(kc) + end if + + ! y-advection d/dy(vy * q) + hqy = ( & + vy(kc,jp,ic)*(humid(kc,jp,ic) + humid(kc,jc,ic)) & + - vy(kc,jc,ic)*(humid(kc,jc,ic) + humid(kc,jm,ic)) & + )*udy + + ! z-advection d/dz(vz * q) + hqz = ( & + vz(kc,jc,ip)*(humid(kc,jc,ip) + humid(kc,jc,ic)) & + - vz(kc,jc,ic)*(humid(kc,jc,ic) + humid(kc,jc,im)) & + )*udz + + ! yy second derivative of humidity + dyyq = (humid(kc,jp,ic) - 2.0*humid(kc,jc,ic) + humid(kc,jm,ic))*dyq + + ! zz second derivative of humidity + dzzq = (humid(kc,jc,ip) - 2.0*humid(kc,jc,ic) + humid(kc,jc,im))*dzq + + rkhumid(kc,jc,ic) = kappa_q*(dyyq + dzzq) - (hqx + hqy + hqz) + + ! add -(q-q_s)/tau H(q-q_s) to rhs + if (humid(kc,jc,ic) > qsat(kc,jc,ic)) then + condensation = (qsat(kc,jc,ic) - humid(kc,jc,ic))*itau + rkhumid(kc,jc,ic) = rkhumid(kc,jc,ic) + condensation + end if + ! condensation = (qsat(kc,jc,ic) - humid(kc,jc,ic))*itau & + ! ! tanh approximation to Heaviside function + ! *0.5*(1.0 + tanh(1e5*(humid(kc,jc,ic) - qsat(kc,jc,ic)))) + ! rkhumid(kc,jc,ic) = rkhumid(kc,jc,ic) + condensation + + end do + end do + end do + +end subroutine ExplicitHumidity + +!> Add the effect of condensation to the buoyancy equation +!! gamma/tau*(q - q_s) H(q-q_s) +subroutine AddCondensation + integer :: ic,jc,kc + real :: condensation, gtau + + gtau = gamma_q/tau_q + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + if (humid(kc,jc,ic) > qsat(kc,jc,ic)) then + condensation = gtau*(humid(kc,jc,ic) - qsat(kc,jc,ic)) + hro(kc,jc,ic) = hro(kc,jc,ic) + condensation + end if + ! condensation = (humid(kc,jc,ic) - qsat(kc,jc,ic))*gtau & + ! ! tanh approximation to Heaviside function + ! *0.5*(1.0 + tanh(1e5*(humid(kc,jc,ic) - qsat(kc,jc,ic)))) + ! hro(kc,jc,ic) = hro(kc,jc,ic) + condensation + + end do + end do + end do + +end subroutine AddCondensation + +!> Update the saturation variable qsat using the buoyancy field temp +subroutine UpdateSaturation + integer :: ic, jc, kc + + ! Recall for this model we use the variable temp for the *buoyancy* field, not the temperature + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + qsat(kc,jc,ic) = exp(alpha_q*(temp(kc,jc,ic) - beta_q*xm(kc))) + end do + end do + end do +end subroutine UpdateSaturation + +!> Compute the implicit terms for the humidity evolution +subroutine ImplicitHumidity + integer :: ic, jc, kc + real :: dxxq, alpec + + alpec = al/pecq + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + ! Second xx derivative + ! Apply lower BC + if (kc==1) then + dxxq = humid(kc+1,jc,ic)*ap3ssk(kc) & + + humid(kc ,jc,ic)*ac3ssk(kc) & + - (ap3ssk(kc) + ac3ssk(kc))*humbp(1,jc,ic)*qfixS + ! Apply upper BC + elseif (kc==nxm) then + dxxq = -(am3ssk(kc) + ac3ssk(kc))*humtp(1,jc,ic)*qfixN & + + humid(kc ,jc,ic)*ac3ssk(kc) & + + humid(kc-1,jc,ic)*am3ssk(kc) + ! Compute dxxq in the interior + else + dxxq = humid(kc+1,jc,ic)*ap3ssk(kc) & + + humid(kc ,jc,ic)*ac3ssk(kc) & + + humid(kc-1,jc,ic)*am3ssk(kc) + end if + + rhs(kc,jc,ic) = (ga*rkhumid(kc,jc,ic) + ro*hkhumid(kc,jc,ic) + alpec*dxxq)*dt + + hkhumid(kc,jc,ic) = rkhumid(kc,jc,ic) + + end do + end do + end do + + call SolveImpEqnUpdate_Humidity + +end subroutine ImplicitHumidity + +!> Solve the implicit system for the humidity +!! and update the global variable +subroutine SolveImpEqnUpdate_Humidity + integer :: ic, jc, kc, info, ipkv(nxm), nrhs + real :: amkl(nxm), ackl(nxm), apkl(nxm), ackl_b + real :: amkq(nxm-1), ackq(nxm), apkq(nxm-1), appk(nxm-2), betadx + + betadx = 0.5d0*al*dt/pecq + + ! Construct tridiagonal matrix for LHS + do kc=1,nxm + ackl_b = 1.d0/(1.d0 - ac3ssk(kc)*betadx) + amkl(kc) = -am3ssk(kc)*betadx*ackl_b + ackl(kc) = 1.d0 + apkl(kc) = -ap3ssk(kc)*betadx*ackl_b + end do + + amkq = amkl(2:nxm) + ackq = ackl(1:nxm) + apkq = apkl(1:(nxm-1)) + + call dgttrf(nxm, amkq, ackq, apkq, appk, ipkv, info) + + nrhs = (xend(3)-xstart(3)+1)*(xend(2)-xstart(2)+1) + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + ackl_b = 1.0/(1.0 - ac3ssk(kc)*betadx) + rhs(kc,jc,ic) = rhs(kc,jc,ic)*ackl_b + end do + end do + end do + + call dgttrs('N', nxm, nrhs, amkq, ackq, apkq, appk, ipkv, rhs(1:nxm,:,:), nxm, info) + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + humid(kc,jc,ic) = humid(kc,jc,ic) + rhs(kc,jc,ic) + end do + end do + end do + +end subroutine SolveImpEqnUpdate_Humidity + +!> Calculate vertical profiles of moisture-related statistics and +!! store them in means.h5 +subroutine CalcMoistStats + + real, dimension(nxm) :: qbar !! Horizontally-averaged specific humidity + real, dimension(nxm) :: qrms !! Horizontally-averaged rms humidity + real, dimension(nxm) :: qsbar !! Horizontally-averaged saturation humidity + real, dimension(nxm) :: qrel !! Horizontally-averaged relative humidity (q/qs) + + real, dimension(nxm) :: vxq !! Advective flux of moisture (x) + real, dimension(nxm) :: vyq !! Advective flux of moisture (y) + real, dimension(nxm) :: vzq !! Advective flux of moisture (z) + + real :: inyzm !! 1.0/nym/nzm + + character(30) :: dsetname !! Dataset name for HDF5 file + character(30) :: filename !! HDF5 file name for statistic storage + character( 5) :: nstat !! Character string of statistic index + + integer :: i, j, k, ip, jp + + inyzm = 1.0/nym/nzm + + filename = trim("outputdir/means.h5") + + qbar(:) = 0.0; qsbar(:) = 0.0; qrms(:) = 0.0; qrel(:) = 0.0 + vxq(:) = 0.0; vyq(:) = 0.0; vzq(:) = 0.0 + + do i=xstart(3),xend(3) + ip = i + 1 + do j=xstart(2),xend(2) + jp = j + 1 + do k=1,nxm + qbar(k) = qbar(k) + humid(k,j,i) + qsbar(k) = qsbar(k) + qsat(k,j,i) + qrms(k) = qrms(k) + humid(k,j,i)**2 + qrel(k) = qrel(k) + humid(k,j,i)/qsat(k,j,i) + vxq(k) = vxq(k) + 0.5*(vx(k,j,i) + vx(k+1,j,i))*humid(k,j,i) + vyq(k) = vyq(k) + 0.5*(vy(k,j,i) + vy(k, jp,i))*humid(k,j,i) + vzq(k) = vzq(k) + 0.5*(vz(k,j,i) + vz(k,j, ip))*humid(k,j,i) + end do + end do + end do + + call MpiSumReal1D(qbar,nxm) + call MpiSumReal1D(qsbar,nxm) + call MpiSumReal1D(qrms,nxm) + call MpiSumReal1D(qrel,nxm) + call MpiSumReal1D(vxq,nxm) + call MpiSumReal1D(vyq,nxm) + call MpiSumReal1D(vzq,nxm) + + do k=1,nxm + qbar(k) = qbar(k)*inyzm + qsbar(k) = qsbar(k)*inyzm + qrms(k) = sqrt(qrms(k)*inyzm) + qrel(k) = qrel(k)*inyzm + vxq(k) = vxq(k)*inyzm + vyq(k) = vyq(k)*inyzm + vzq(k) = vzq(k)*inyzm + end do + + ! Store index as character string + write(nstat,"(i5.5)")nint(time/tout) + + + if (ismaster) then + dsetname = trim("qbar/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, qbar, nxm) + dsetname = trim("qsbar/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, qsbar, nxm) + dsetname = trim("qrms/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, qrms, nxm) + dsetname = trim("qrel/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, qrel, nxm) + dsetname = trim("vxq/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, vxq, nxm) + dsetname = trim("vyq/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, vyq, nxm) + dsetname = trim("vzq/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, vzq, nxm) + end if + + call MpiBarrier + +end subroutine CalcMoistStats + +!> Create the groups in the means.h5 file to store the +!! moist-related statistics +subroutine CreateMoistH5Groups(filename) + use HDF5 + + character(30), intent(in) :: filename + integer(HID_T) :: file_id, group_id + integer :: hdf_error + + call h5fopen_f(filename, H5F_ACC_RDWR_F, file_id, hdf_error) + + call h5gcreate_f(file_id, "qbar", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "qsbar", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "qrms", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "qrel", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "vxq", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "vyq", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "vzq", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + + call h5fclose_f(file_id, hdf_error) + +end subroutine CreateMoistH5Groups + +end module afid_moisture \ No newline at end of file diff --git a/src/multires/CreateMgrdGrid.F90 b/src/multires/CreateMgrdGrid.F90 index d90c1f9b..8fbd1395 100644 --- a/src/multires/CreateMgrdGrid.F90 +++ b/src/multires/CreateMgrdGrid.F90 @@ -12,12 +12,11 @@ subroutine CreateMgrdGrid use param use AuxiliaryRoutines use GridModule + use afid_salinity, only: SfixS, SfixN, PraS, ap3sskr, ac3sskr, am3sskr + use afid_phasefield, only: ap3spkr, ac3spkr, am3spkr implicit none - real :: x1,x2,x3 - - integer :: i, j, kc, km, kp - logical :: fexist + integer :: kc do kc=1,nxmr kmvr(kc)=kc-1 @@ -72,6 +71,11 @@ subroutine CreateMgrdGrid if (istr3r==4) call tanh_grid(xcr(1:nxr),xmr(1:nxmr),nxmr,alx3,str3) +! +! OPTION 5: SCALLOP-FOCUSED LOWER WALL CLUSTERING +! + if (istr3==5) call scallop_grid(xcr(1:nxr), xmr(1:nxmr), nxmr, alx3, dPdy, 0.5) + ! ! OPTION 6: CLIPPED CHEBYCHEV-TYPE CLUSTERING ! @@ -154,9 +158,9 @@ subroutine CreateMgrdGrid ! ! Salinity differentiation - call second_derivative_coeff(ap3sskr, ac3sskr, am3sskr, xmr(1:nxmr), alx3, SfixN, SfixS) + if (salinity) call second_derivative_coeff(ap3sskr, ac3sskr, am3sskr, xmr(1:nxmr), alx3, SfixN, SfixS) ! Phase-field differentiation (ensuring zero gradient at boundaries) - call second_derivative_coeff(ap3spkr, ac3spkr, am3spkr, xmr(1:nxmr), alx3, 0, 0) + if (phasefield) call second_derivative_coeff(ap3spkr, ac3spkr, am3spkr, xmr(1:nxmr), alx3, 0, 0) return end subroutine CreateMgrdGrid \ No newline at end of file diff --git a/src/multires/CreateMgrdStencil.F90 b/src/multires/CreateMgrdStencil.F90 index b7a9d856..bffa0e93 100644 --- a/src/multires/CreateMgrdStencil.F90 +++ b/src/multires/CreateMgrdStencil.F90 @@ -78,4 +78,29 @@ subroutine CreateMgrdStencil end if return -end subroutine CreateMgrdStencil \ No newline at end of file +end subroutine CreateMgrdStencil + + +subroutine CreatePFStencil + use param + use mgrd_arrays + use HermiteInterpolations + + implicit none + + integer, dimension(0:nyr) :: jrange !< temporary index counter in y + integer, dimension(0:nzr) :: krange !< temporary index counter in z + + ! Save indices for fast interpolation + call interpolation_indices(yc_to_ymr, yc(1:nym), ymr(1:nymr), ylen) + call interpolation_indices(zc_to_zmr, zc(1:nzm), zmr(1:nzmr), zlen) + + ! Create temporary index counters for stencil creation + call interpolation_indices(jrange, ymr(1:nymr), yc(1:nym), ylen) + call interpolation_indices(krange, zmr(1:nzmr), zc(1:nzm), zlen) + + ! Create stencils + call construct_stencil(cych, ymr(1:nymr), yc(1:nym), ylen, jrange, "y") + call construct_stencil(czch, zmr(1:nzmr), zc(1:nzm), zlen, krange, "z") + +end subroutine CreatePFStencil \ No newline at end of file diff --git a/src/multires/DeallocateMgrdVariables.F90 b/src/multires/DeallocateMgrdVariables.F90 index abe136e1..7ffeb26e 100644 --- a/src/multires/DeallocateMgrdVariables.F90 +++ b/src/multires/DeallocateMgrdVariables.F90 @@ -28,17 +28,17 @@ subroutine DeallocateMgrdVariables call DestroyReal1DArray(udx3cr) call DestroyReal1DArray(udx3mr) - call DestroyReal1DArray(ap3ckr) - call DestroyReal1DArray(ac3ckr) - call DestroyReal1DArray(am3ckr) + ! call DestroyReal1DArray(ap3ckr) + ! call DestroyReal1DArray(ac3ckr) + ! call DestroyReal1DArray(am3ckr) - call DestroyReal1DArray(ap3sskr) - call DestroyReal1DArray(ac3sskr) - call DestroyReal1DArray(am3sskr) + ! call DestroyReal1DArray(ap3sskr) + ! call DestroyReal1DArray(ac3sskr) + ! call DestroyReal1DArray(am3sskr) - call DestroyReal1DArray(ap3spkr) - call DestroyReal1DArray(ac3spkr) - call DestroyReal1DArray(am3spkr) + ! call DestroyReal1DArray(ap3spkr) + ! call DestroyReal1DArray(ac3spkr) + ! call DestroyReal1DArray(am3spkr) call DestroyInt1dArray(kmcr) call DestroyInt1dArray(kpcr) @@ -61,6 +61,9 @@ subroutine DeallocateMgrdVariables call DestroyInt1dArray(jrangb) !CS mgrd call DestroyInt1dArray(krangb) !CS mgrd + call DestroyInt1dArray(yc_to_ymr) + call DestroyInt1dArray(zc_to_zmr) + call DestroyReal2DArray(cxvx) !CS mgrd call DestroyReal2DArray(cxvy) !CS mgrd call DestroyReal2DArray(cxvz) !CS mgrd @@ -85,9 +88,12 @@ subroutine DeallocateMgrdVariables call DestroyReal2DArray(cyphic) call DestroyReal2DArray(czphic) - call DestroyReal3DArray(vxr) - call DestroyReal3DArray(vyr) - call DestroyReal3DArray(vzr) + if (IBM) then + if (phasefield) then + call DestroyReal2DArray(cych) + call DestroyReal2DArray(czch) + end if + end if call DestroyReal3DArray(tpdv) call DestroyReal3DArray(tpdvr) !CS BUG: ERROR WHILE DEALLOCATING diff --git a/src/multires/HermiteInterpolations.F90 b/src/multires/HermiteInterpolations.F90 index d7d11562..6cf3dcf7 100644 --- a/src/multires/HermiteInterpolations.F90 +++ b/src/multires/HermiteInterpolations.F90 @@ -112,16 +112,21 @@ subroutine construct_stencil(cx, x_old, x_new, x_len, idx, axis) dlp = xo(io+2) - xo(io+1) do in=max(idx(io),1),min(idx(io+1)-1,n_new) t = (xn(in) - xo(io))/dlc - h00 = (1.0 + 2.0*t)*(1.0 - t)**2 - h10 = t*(1.0 - t)**2 - h01 = (1.0 + 2.0*(1.0 - t))*t**2 - h11 = -(1.0 - t)*t**2 - cx(1,in) = -h10*dlc**2/dlm/(dlc + dlm) - cx(2,in) = h00 - h11*dlp/(dlp + dlc) & - + h10*(dlc - dlm)/dlm - cx(3,in) = h01 + h10*dlm/(dlm + dlc) & - + h11*(dlp - dlc)/dlp - cx(4,in) = h11*dlc**2/dlp/(dlp + dlc) + if (t<1e-12) then + cx(:,in) = 0.0 + cx(2,in) = 1.0 + else + h00 = (1.0 + 2.0*t)*(1.0 - t)**2 + h10 = t*(1.0 - t)**2 + h01 = (1.0 + 2.0*(1.0 - t))*t**2 + h11 = -(1.0 - t)*t**2 + cx(1,in) = -h10*dlc**2/dlm/(dlc + dlm) + cx(2,in) = h00 - h11*dlp/(dlp + dlc) & + + h10*(dlc - dlm)/dlm + cx(3,in) = h01 + h10*dlm/(dlm + dlc) & + + h11*(dlp - dlc)/dlp + cx(4,in) = h11*dlc**2/dlp/(dlp + dlc) + end if end do end if end do diff --git a/src/multires/IC_interpolation/CreateOldGrid.F90 b/src/multires/IC_interpolation/CreateOldGrid.F90 index 89b42dd3..04e6f0da 100644 --- a/src/multires/IC_interpolation/CreateOldGrid.F90 +++ b/src/multires/IC_interpolation/CreateOldGrid.F90 @@ -15,13 +15,6 @@ subroutine CreateOldGrid use GridModule implicit none - real :: x1, x2, x3 - real :: delet, etain, tstr3, z2dp - - real, allocatable, dimension(:) :: etaz, etazm - - integer :: i, j, k, nxp, nclip - nxmo = nxo - 1 nymo = nyo - 1 nzmo = nzo - 1 @@ -65,6 +58,22 @@ subroutine CreateOldGrid if (istr3ro==1) call centre_focus_grid(xcro(1:nxro),xmro(1:nxmro),nxmro,alx3,str3o) end if + ! Option 2: Natural turb BL clustering + + if (istr3o==2) call natural_BL_grid(xco(1:nxo),xmo(1:nxmo),nxmo,alx3) + + if (multires) then + if (istr3ro==2) call natural_BL_grid(xcro(1:nxro),xmro(1:nxmro),nxmro,alx3) + end if + + ! Option 3: Symmetric natural turb BL clustering + + if (istr3o==3) call sym_natural_BL_grid(xco(1:nxo),xmo(1:nxmo),nxmo,alx3, 1.0) + + if (multires) then + if (istr3ro==3) call sym_natural_BL_grid(xcro(1:nxro),xmro(1:nxmro),nxmro,alx3, 1.0) + end if + ! Option 4: Hyperbolic tangent-type clustering if (istr3o==4) call tanh_grid(xco(1:nxo),xmo(1:nxmo),nxmo,alx3,str3o) @@ -73,6 +82,14 @@ subroutine CreateOldGrid if (istr3ro==4) call tanh_grid(xcro(1:nxro),xmro(1:nxmro),nxmro,alx3,str3o) end if + ! Option 5: Scallop-focused wall clustering + + if (istr3==5) call scallop_grid(xco(1:nxo), xmo(1:nxmo), nxmo, alx3, dPdy, 0.5) + + if (multires) then + if (istr3ro==5) call scallop_grid(xcro(1:nxro),xmro(1:nxmro),nxmro,alx3,dPdy, 0.5) + end if + ! Option 6: Clipped Chebychev-type clustering if (istr3o==6) call cheb_grid(xco(1:nxo),xmo(1:nxmo),nxmo,alx3,str3o) diff --git a/src/multires/IC_interpolation/InterpInputPhi.F90 b/src/multires/IC_interpolation/InterpInputPhi.F90 index e9fc3ba9..ba6f020e 100644 --- a/src/multires/IC_interpolation/InterpInputPhi.F90 +++ b/src/multires/IC_interpolation/InterpInputPhi.F90 @@ -2,18 +2,15 @@ subroutine InterpInputPhi use param use input_grids - use mgrd_arrays + ! use mgrd_arrays + use afid_phasefield, only: phi use mpih use decomp_2d use AuxiliaryRoutines use HermiteInterpolations, only: interpolate_xyz_old_to_new_ref implicit none - integer :: ic,jc,kc, ip,jp,kp, icr,jcr,kcr - - real,dimension(4,4,4) :: qv3 - real,dimension(4,4) :: qv2 - real,dimension(4) :: qv1 + integer :: ic,jc real, allocatable, dimension(:,:,:) :: phio diff --git a/src/multires/IC_interpolation/InterpInputSal.F90 b/src/multires/IC_interpolation/InterpInputSal.F90 index 0655c8fc..b10b0821 100644 --- a/src/multires/IC_interpolation/InterpInputSal.F90 +++ b/src/multires/IC_interpolation/InterpInputSal.F90 @@ -2,18 +2,14 @@ subroutine InterpInputSal use param use input_grids - use mgrd_arrays + use afid_salinity use mpih use decomp_2d use AuxiliaryRoutines use HermiteInterpolations, only: interpolate_xyz_old_to_new_ref implicit none - integer :: ic,jc,kc, ip,jp,kp, icr,jcr,kcr - - real,dimension(4,4,4) :: qv3 - real,dimension(4,4) :: qv2 - real,dimension(4) :: qv1 + integer :: ic,jc real :: Sup, Slo diff --git a/src/multires/IC_interpolation/InterpInputVel.F90 b/src/multires/IC_interpolation/InterpInputVel.F90 index 4cc879d5..419dfba1 100644 --- a/src/multires/IC_interpolation/InterpInputVel.F90 +++ b/src/multires/IC_interpolation/InterpInputVel.F90 @@ -12,9 +12,8 @@ subroutine InterpInputVel integer :: ic,jc,kc, ip,jp,kp, icr,jcr,kcr integer :: jc0,jcr0, ic0,icr0 - integer :: comm_col,comm_row,comm,ierror + integer :: comm_col,comm_row,ierror - real,dimension(4,4,4) :: qv3 real,dimension(4,4) :: qv2 real,dimension(4) :: qv1 @@ -27,7 +26,7 @@ subroutine InterpInputVel !-- Allocate temporary arrays for velocities and gradients real,allocatable,dimension(:,:) :: dvyloc,dvybot, dvzloc,dvzbot - call AllocateReal3DArray(tpdvo,-1,nx+1, & + call AllocateReal3DArray(tpdvo,-1,nxo+1, & xstarto(2)-2,xendo(2)+2,xstarto(3)-2,xendo(3)+2) tpdvo(:,:,:) = 0.d0 ! Temporary gradient array - old @@ -114,8 +113,8 @@ subroutine InterpInputVel !-- Boundary points, enforce zero velocity !CJH e.g. v=0 on x=0 => dv/dy=0 on x=0 - tpdvo(0,jc,ic) = -tpdvo(1,jc,ic) - tpdvo(nxo,jc,ic) = -tpdvo(nxmo,jc,ic) + tpdvo(0,jc,ic) = (1-2*inslwS)*tpdvo(1,jc,ic) + tpdvo(nxo,jc,ic) = (1-2*inslwN)*tpdvo(nxmo,jc,ic) enddo enddo @@ -136,8 +135,8 @@ subroutine InterpInputVel vyxzc(kc,ic)=vyo(kc,jc0,ic) !CS Halo updates can be optimised. Otherwise lvlhalo=2 required enddo ! x boundaries - vyxzc(0,ic) = -vyxzc(1,ic) - vyxzc(nxo,ic) = -vyxzc(nxmo,ic) + vyxzc(0,ic) = (1-2*inslwS)*vyxzc(1,ic) + vyxzc(nxo,ic) = (1-2*inslwN)*vyxzc(nxmo,ic) enddo do ic=xstarto(3)-1,xendo(3) diff --git a/src/multires/InitMgrdVariables.F90 b/src/multires/InitMgrdVariables.F90 index fbdf18bf..275e6e4e 100644 --- a/src/multires/InitMgrdVariables.F90 +++ b/src/multires/InitMgrdVariables.F90 @@ -37,17 +37,17 @@ subroutine InitMgrdVariables call AllocateReal1DArray(udx3cr,1,nxr) call AllocateReal1DArray(udx3mr,1,nxr) - call AllocateReal1DArray(ap3ckr,1,nxr) - call AllocateReal1DArray(ac3ckr,1,nxr) - call AllocateReal1DArray(am3ckr,1,nxr) + ! call AllocateReal1DArray(ap3ckr,1,nxr) + ! call AllocateReal1DArray(ac3ckr,1,nxr) + ! call AllocateReal1DArray(am3ckr,1,nxr) - call AllocateReal1DArray(ap3sskr,1,nxr) - call AllocateReal1DArray(ac3sskr,1,nxr) - call AllocateReal1DArray(am3sskr,1,nxr) + ! call AllocateReal1DArray(ap3sskr,1,nxr) + ! call AllocateReal1DArray(ac3sskr,1,nxr) + ! call AllocateReal1DArray(am3sskr,1,nxr) - call AllocateReal1DArray(ap3spkr,1,nxr) - call AllocateReal1DArray(ac3spkr,1,nxr) - call AllocateReal1DArray(am3spkr,1,nxr) + ! call AllocateReal1DArray(ap3spkr,1,nxr) + ! call AllocateReal1DArray(ac3spkr,1,nxr) + ! call AllocateReal1DArray(am3spkr,1,nxr) call AllocateInt1dArray(kmcr,1,nxr) call AllocateInt1dArray(kpcr,1,nxr) @@ -70,6 +70,9 @@ subroutine InitMgrdVariables call AllocateInt1dArray(jrangr,0,nyr) call AllocateInt1dArray(krangr,0,nzr) + call AllocateInt1dArray(yc_to_ymr,0,ny) + call AllocateInt1dArray(zc_to_zmr,0,nz) + call AllocateReal2DArray(cxvx,1,4,1,nxr) call AllocateReal2DArray(cxvy,1,4,1,nxmr) call AllocateReal2DArray(cxvz,1,4,1,nxmr) @@ -94,6 +97,13 @@ subroutine InitMgrdVariables call AllocateReal2DArray(cysalc,1,4,1,nym) call AllocateReal2DArray(czsalc,1,4,1,nzm) + if (IBM) then + if (phasefield) then + call AllocateReal2DArray(cych,1,4,1,nym) + call AllocateReal2DArray(czch,1,4,1,nzm) + end if + end if + !------------------------------------------------- ! Arrays with ghost cells !------------------------------------------------- @@ -103,6 +113,6 @@ subroutine InitMgrdVariables call AllocateReal3DArray(tpdvr,-1,nxr+1,xstartr(2)-2,xendr(2)+2,xstartr(3)-2,xendr(3)+2) ! RHS array without ghost cells - call AllocateReal3DArray(rhsr,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) + call AllocateReal3DArray(rhsr,1,nxmr,xstartr(2),xendr(2),xstartr(3),xendr(3)) return end subroutine InitMgrdVariables \ No newline at end of file diff --git a/src/multires/InterpVelMgrd.F90 b/src/multires/InterpVelMgrd.F90 index e1cdcac6..fa4f1a7f 100644 --- a/src/multires/InterpVelMgrd.F90 +++ b/src/multires/InterpVelMgrd.F90 @@ -3,6 +3,7 @@ subroutine InterpVelMgrd use param use local_arrays, only: vx,vy,vz use mgrd_arrays + use afid_salinity, only: vxr, vyr, vzr use mpih use decomp_2d use AuxiliaryRoutines @@ -10,11 +11,9 @@ subroutine InterpVelMgrd implicit none integer :: ic,jc,kc, ip,jp,kp, icr,jcr,kcr - integer :: icc,jcc,kcc integer :: jc0,jcr0, ic0,icr0 - integer :: comm_col,comm_row,comm,ierror,chk + integer :: comm_col,comm_row,ierror - real,dimension(4,4,4) :: qv3 real,dimension(4,4) :: qv2 real,dimension(4) :: qv1 diff --git a/src/flow_solver/CalcLocalDivergence.F90 b/src/old/CalcLocalDivergence.F90 similarity index 100% rename from src/flow_solver/CalcLocalDivergence.F90 rename to src/old/CalcLocalDivergence.F90 diff --git a/src/flow_solver/CorrectPressure.F90 b/src/old/CorrectPressure.F90 similarity index 100% rename from src/flow_solver/CorrectPressure.F90 rename to src/old/CorrectPressure.F90 diff --git a/src/flow_solver/InitPressureSolver.F90 b/src/old/InitPressureSolver.F90 similarity index 100% rename from src/flow_solver/InitPressureSolver.F90 rename to src/old/InitPressureSolver.F90 diff --git a/src/flow_solver/SolvePressureCorrection.F90 b/src/old/SolvePressureCorrection.F90 similarity index 100% rename from src/flow_solver/SolvePressureCorrection.F90 rename to src/old/SolvePressureCorrection.F90 diff --git a/src/multires/phase-field/AddLatentHeat.F90 b/src/old/phase-field/AddLatentHeat.F90 similarity index 100% rename from src/multires/phase-field/AddLatentHeat.F90 rename to src/old/phase-field/AddLatentHeat.F90 diff --git a/src/multires/phase-field/CreateICPF.F90 b/src/old/phase-field/CreateICPF.F90 similarity index 93% rename from src/multires/phase-field/CreateICPF.F90 rename to src/old/phase-field/CreateICPF.F90 index bafc4a73..a343dd2d 100644 --- a/src/multires/phase-field/CreateICPF.F90 +++ b/src/old/phase-field/CreateICPF.F90 @@ -18,7 +18,7 @@ subroutine CreateICPF implicit none integer :: i,j,k,kmid - real :: r, x0, lambda, h0, t0, A, B, alpha + real :: r, x0, h0, t0, A, B, alpha real, dimension(11) :: yh, zh if (pf_IC == 1) then ! 1D freezing validation @@ -129,6 +129,17 @@ subroutine CreateICPF end if end if + ! if (IBM) then + ! do i=xstartr(3),xendr(3) + ! do j=xstartr(2),xendr(2) + ! h0 = 0.25 + (ymr(j) - 0.5)**2 + ! do k=1,nxmr + ! phi(k,j,i) = 0.5*(1.0 - tanh((xmr(k) - h0)/2/pf_eps)) + ! end do + ! end do + ! end do + ! end if + return end subroutine CreateICPF diff --git a/src/multires/phase-field/DeallocatePFVariables.F90 b/src/old/phase-field/DeallocatePFVariables.F90 similarity index 100% rename from src/multires/phase-field/DeallocatePFVariables.F90 rename to src/old/phase-field/DeallocatePFVariables.F90 diff --git a/src/multires/phase-field/ExplicitTermsPhi.F90 b/src/old/phase-field/ExplicitTermsPhi.F90 similarity index 100% rename from src/multires/phase-field/ExplicitTermsPhi.F90 rename to src/old/phase-field/ExplicitTermsPhi.F90 diff --git a/src/multires/phase-field/ImmersedBoundary.F90 b/src/old/phase-field/ImmersedBoundary.F90 similarity index 100% rename from src/multires/phase-field/ImmersedBoundary.F90 rename to src/old/phase-field/ImmersedBoundary.F90 diff --git a/src/multires/phase-field/ImplicitAndUpdatePhi.F90 b/src/old/phase-field/ImplicitAndUpdatePhi.F90 similarity index 96% rename from src/multires/phase-field/ImplicitAndUpdatePhi.F90 rename to src/old/phase-field/ImplicitAndUpdatePhi.F90 index f9b4e8f2..6a6ba92d 100644 --- a/src/multires/phase-field/ImplicitAndUpdatePhi.F90 +++ b/src/old/phase-field/ImplicitAndUpdatePhi.F90 @@ -13,13 +13,11 @@ subroutine ImplicitAndUpdatePhi use param use mgrd_arrays, only: phi,hphi,rhsr,ruphi - use decomp_2d, only: xstartr,xendr,nrank + use decomp_2d, only: xstartr,xendr use mpih implicit none integer :: jc,kc,ic - integer :: km,kp real :: alpec,dxxp - real :: app,acc,amm alpec=al*pf_D diff --git a/src/multires/phase-field/InitPFVariables.F90 b/src/old/phase-field/InitPFVariables.F90 similarity index 77% rename from src/multires/phase-field/InitPFVariables.F90 rename to src/old/phase-field/InitPFVariables.F90 index 832acc3f..221b8b63 100644 --- a/src/multires/phase-field/InitPFVariables.F90 +++ b/src/old/phase-field/InitPFVariables.F90 @@ -27,5 +27,12 @@ subroutine InitPFVariables ! Coarse array for phi or d(phi)/dt call AllocateReal3DArray(phic,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + if (IBM) then + call AllocateReal2DArray(solid_height,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + call AllocateReal2DArray(height_vx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal2DArray(height_vy,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal2DArray(height_vz,xstart(2),xend(2),xstart(3),xend(3)) + end if + return end subroutine InitPFVariables \ No newline at end of file diff --git a/src/multires/phase-field/InterpPhiMgrd.F90 b/src/old/phase-field/InterpPhiMgrd.F90 similarity index 87% rename from src/multires/phase-field/InterpPhiMgrd.F90 rename to src/old/phase-field/InterpPhiMgrd.F90 index e2a5c735..69a0284c 100644 --- a/src/multires/phase-field/InterpPhiMgrd.F90 +++ b/src/old/phase-field/InterpPhiMgrd.F90 @@ -11,18 +11,14 @@ subroutine InterpPhiMgrd use param - use mgrd_arrays, only: phi,phic,cxsalc,cysalc,czsalc,irangr,jrangr,krangr,tpdvr + use mgrd_arrays, only: phi,phic,tpdvr use mpih use decomp_2d use AuxiliaryRoutines use HermiteInterpolations, only: interpolate_xyz_to_coarse, interpolate_xyz_to_coarse_fast implicit none - integer :: ic,jc,kc, icr,jcr,kcr - - real,dimension(4,4,4) :: qv3 - real,dimension(4,4) :: qv2 - real,dimension(4) :: qv1 + integer :: ic,jc,kc ! Interpolate to coarse grid here. A better option is to apply ! a box filter. diff --git a/src/multires/phase-field/InterpTempMgrd.F90 b/src/old/phase-field/InterpTempMgrd.F90 similarity index 86% rename from src/multires/phase-field/InterpTempMgrd.F90 rename to src/old/phase-field/InterpTempMgrd.F90 index bce20087..c9d38818 100644 --- a/src/multires/phase-field/InterpTempMgrd.F90 +++ b/src/old/phase-field/InterpTempMgrd.F90 @@ -11,18 +11,14 @@ subroutine InterpTempMgrd use param use local_arrays, only: temp - use mgrd_arrays, only: tempr,cxrs,cyrs,czrs,irangs,jrangs,krangs,tpdv + use mgrd_arrays, only: tempr,tpdv use mpih use decomp_2d use AuxiliaryRoutines use HermiteInterpolations, only: interpolate_xyz_to_refined implicit none - integer :: ic,jc,kc,icr,jcr,kcr - - real,dimension(4,4,4) :: qv3 - real,dimension(4,4) :: qv2 - real,dimension(4) :: qv1 + integer :: ic,jc,kc tempr(:,:,:) = 0.d0 diff --git a/src/multires/phase-field/SolveImpEqnUpdate_Phi.F90 b/src/old/phase-field/SolveImpEqnUpdate_Phi.F90 similarity index 95% rename from src/multires/phase-field/SolveImpEqnUpdate_Phi.F90 rename to src/old/phase-field/SolveImpEqnUpdate_Phi.F90 index 5c6a67c7..7bd5243b 100644 --- a/src/multires/phase-field/SolveImpEqnUpdate_Phi.F90 +++ b/src/old/phase-field/SolveImpEqnUpdate_Phi.F90 @@ -50,7 +50,7 @@ subroutine SolveImpEqnUpdate_Phi end do end do - call dgttrs('N',nxmr,nrhs,amkT,ackT,apkT,appk,ipkv,rhsr(1:nxmr,:,:),nxmr,info) + call dgttrs('N',nxmr,nrhs,amkT,ackT,apkT,appk,ipkv,rhsr,nxmr,info) do ic=xstartr(3),xendr(3) do jc=xstartr(2),xendr(2) diff --git a/src/old/phase-field/UpdateIBMLocation.F90 b/src/old/phase-field/UpdateIBMLocation.F90 new file mode 100644 index 00000000..6ed72d59 --- /dev/null +++ b/src/old/phase-field/UpdateIBMLocation.F90 @@ -0,0 +1,28 @@ +!! Subroutine to update location of the immersed boundary from the phase field variable. +!! Assumes that the solid phase lies beneath the +subroutine UpdateIBMLocation + use mgrd_arrays, only: phi + use ibm_param + use IBMTools + + ! Integrate phi vertically to calculate height profile for solid + call calc_interface_height(phi,solid_height) + ! if (ismaster) write(*,*) 'height at first index:', solid_height(1,1) + + ! Interpolate solid height profile from refined grid to each velocity grid + call interp_height_to_vel_grid(solid_height, height_vx, height_vy, height_vz) + + ! Build `ibmaskX` variables from height profile + ! call mask_below_height(height_vx, ibmaskx, 'x') + ! call mask_below_height(height_vy, ibmasky, 'y') + ! call mask_below_height(height_vz, ibmaskz, 'z') + call mask_above_height(height_vx, ibmaskx, 'x') + call mask_above_height(height_vy, ibmasky, 'y') + call mask_above_height(height_vz, ibmaskz, 'z') + + ! For each velocity grid, store the interpolation values for the boundary points + call calc_IBM_interpolation(height_vx, ibmaskx, distx, 'x') + call calc_IBM_interpolation(height_vy, ibmasky, disty, 'y') + call calc_IBM_interpolation(height_vz, ibmaskz, distz, 'z') + +end subroutine UpdateIBMLocation \ No newline at end of file diff --git a/src/multires/salinity/CreateICSal.F90 b/src/old/salinity/CreateICSal.F90 similarity index 99% rename from src/multires/salinity/CreateICSal.F90 rename to src/old/salinity/CreateICSal.F90 index 9cb3647e..dc3b5ef1 100644 --- a/src/multires/salinity/CreateICSal.F90 +++ b/src/old/salinity/CreateICSal.F90 @@ -15,7 +15,7 @@ subroutine CreateICSal use mpih implicit none integer :: i,k,j, kmid - real :: xxx,yyy,eps,varptb,amp + real :: eps,varptb,amp real :: gamma, t0, x0, h0, A, B, alpha real, dimension(11) :: yh, zh diff --git a/src/multires/salinity/DeallocateSalVariables.F90 b/src/old/salinity/DeallocateSalVariables.F90 similarity index 100% rename from src/multires/salinity/DeallocateSalVariables.F90 rename to src/old/salinity/DeallocateSalVariables.F90 diff --git a/src/multires/salinity/ExplicitTermsSal.F90 b/src/old/salinity/ExplicitTermsSal.F90 similarity index 100% rename from src/multires/salinity/ExplicitTermsSal.F90 rename to src/old/salinity/ExplicitTermsSal.F90 diff --git a/src/multires/salinity/ImplicitAndUpdateSal.F90 b/src/old/salinity/ImplicitAndUpdateSal.F90 similarity index 75% rename from src/multires/salinity/ImplicitAndUpdateSal.F90 rename to src/old/salinity/ImplicitAndUpdateSal.F90 index 21ed4f0a..f5edd5cd 100644 --- a/src/multires/salinity/ImplicitAndUpdateSal.F90 +++ b/src/old/salinity/ImplicitAndUpdateSal.F90 @@ -17,10 +17,7 @@ subroutine ImplicitAndUpdateSal use ibm_param implicit none integer :: jc,kc,ic - integer :: km,kp,ke real :: alpec,dxxs - real :: app,acc,amm - real :: usaldto,sale alpec=al/pecs @@ -40,24 +37,24 @@ subroutine ImplicitAndUpdateSal ! Calculate second derivative of salinty in the x-direction. ! This is the only term calculated implicitly for salinity. if (kc.eq.1) then !CJH Apply lower BC - dxxs = sal(kc+1,jc,ic)*ap3sskr(kc) & - + sal(kc,jc,ic)*ac3sskr(kc) & - - (ap3sskr(kc) + ac3sskr(kc))*salbp(1,jc,ic)*SfixS + dxxs= sal(kc+1,jc,ic)*ap3sskr(kc) & + + sal(kc ,jc,ic)*ac3sskr(kc) & + - (ap3sskr(kc) + ac3sskr(kc))*salbp(1,jc,ic)*SfixS elseif(kc.eq.nxmr) then !CJH Apply upper BC - dxxs = sal(kc,jc,ic)*ac3sskr(kc) & - + sal(kc-1,jc,ic)*am3sskr(kc) & - - (am3sskr(kc) + ac3sskr(kc))*saltp(1,jc,ic)*SfixN + dxxs= sal(kc ,jc,ic)*ac3sskr(kc) & + + sal(kc-1,jc,ic)*am3sskr(kc) & + - (am3sskr(kc) + ac3sskr(kc))*saltp(1,jc,ic)*SfixN else - dxxs = sal(kc+1,jc,ic)*ap3sskr(kc) & - + sal(kc ,jc,ic)*ac3sskr(kc) & - + sal(kc-1,jc,ic)*am3sskr(kc) + dxxs= sal(kc+1,jc,ic)*ap3sskr(kc) & + + sal(kc ,jc,ic)*ac3sskr(kc) & + + sal(kc-1,jc,ic)*am3sskr(kc) end if ! Calculate right hand side of Eq. 5 (VO96) rhsr(kc,jc,ic) = (ga*hsal(kc,jc,ic) + ro*rusal(kc,jc,ic) & - + alpec*dxxs)*dt + + alpec*dxxs)*dt ! Store the non-linear terms for the calculation of ! the next timestep diff --git a/src/multires/salinity/InitSalVariables.F90 b/src/old/salinity/InitSalVariables.F90 similarity index 100% rename from src/multires/salinity/InitSalVariables.F90 rename to src/old/salinity/InitSalVariables.F90 diff --git a/src/multires/salinity/InterpSalMgrd.F90 b/src/old/salinity/InterpSalMgrd.F90 similarity index 86% rename from src/multires/salinity/InterpSalMgrd.F90 rename to src/old/salinity/InterpSalMgrd.F90 index 69a3530d..4e7d9fd5 100644 --- a/src/multires/salinity/InterpSalMgrd.F90 +++ b/src/old/salinity/InterpSalMgrd.F90 @@ -1,18 +1,14 @@ subroutine InterpSalMgrd use param - use mgrd_arrays, only: sal,salc,cxsalc,cysalc,czsalc,irangr,jrangr,krangr,tpdvr + use mgrd_arrays, only: sal,salc,tpdvr use mpih use decomp_2d use AuxiliaryRoutines use HermiteInterpolations, only: interpolate_xyz_to_coarse, interpolate_xyz_to_coarse_fast implicit none - integer :: ic,jc,kc, icr,jcr,kcr - - real,dimension(4,4,4) :: qv3 - real,dimension(4,4) :: qv2 - real,dimension(4) :: qv1 + integer :: icr,jcr,kcr ! Interpolate to coarse grid here. A better option is to apply ! a box filter. diff --git a/src/multires/salinity/SetSalBCs.F90 b/src/old/salinity/SetSalBCs.F90 similarity index 100% rename from src/multires/salinity/SetSalBCs.F90 rename to src/old/salinity/SetSalBCs.F90 diff --git a/src/multires/salinity/SolveImpEqnUpdate_Sal.F90 b/src/old/salinity/SolveImpEqnUpdate_Sal.F90 similarity index 96% rename from src/multires/salinity/SolveImpEqnUpdate_Sal.F90 rename to src/old/salinity/SolveImpEqnUpdate_Sal.F90 index 256775ea..376e8f9c 100644 --- a/src/multires/salinity/SolveImpEqnUpdate_Sal.F90 +++ b/src/old/salinity/SolveImpEqnUpdate_Sal.F90 @@ -56,7 +56,7 @@ subroutine SolveImpEqnUpdate_Sal end do end do - call dgttrs('N',nxmr,nrhs,amkT,ackT,apkT,appk,ipkv,rhsr(1:nxmr,:,:),nxmr,info) + call dgttrs('N',nxmr,nrhs,amkT,ackT,apkT,appk,ipkv,rhsr,nxmr,info) do ic=xstartr(3),xendr(3) do jc=xstartr(2),xendr(2) diff --git a/src/multires/salinity/UpdateScalarBCs.F90 b/src/old/salinity/UpdateScalarBCs.F90 similarity index 100% rename from src/multires/salinity/UpdateScalarBCs.F90 rename to src/old/salinity/UpdateScalarBCs.F90 diff --git a/src/phasefield.F90 b/src/phasefield.F90 new file mode 100644 index 00000000..32f44595 --- /dev/null +++ b/src/phasefield.F90 @@ -0,0 +1,806 @@ +!> Module adding phase-field model to AFiD +!! to simulate the melting of solid objects +module afid_phasefield + use param + use mgrd_arrays + use decomp_2d, only: xstart, xend, xstartr, xendr, update_halo + use AuxiliaryRoutines + use HermiteInterpolations, only: interpolate_xyz_to_coarse, interpolate_xyz_to_coarse_fast, interpolate_xyz_to_refined + use ibm_param, only: solidr + implicit none + + real, allocatable, dimension(:,:,:) :: phi !! Phase-field variable + real, allocatable, dimension(:,:,:) :: ruphi !! RK storage array for phase-field (previous substep) + real, allocatable, dimension(:,:,:) :: hphi !! RK storage array for phase-field + real, allocatable, dimension(:,:,:) :: phic !! Interpolated phase-field on coarse grid (also used to store d(phi)/dt) + real, allocatable, dimension(:,:,:) :: tempr !! Interpolated temperature field on refined grid + + real :: pf_A !! Phase-field Gibbs-Thomson parameter + real :: pf_eps !! Phase-field interface thickness + real :: pf_D !! Dimensionless diffusivity of the phase-field + real :: pf_S !! Stefan number + real :: pf_Tm !! Equilibrium melting temperature + real :: pf_Lambda !! Dimensionless liquidus slope + + real, allocatable, dimension(:) :: ap3spkr !! Upper diagonal derivative coefficient for salinity + real, allocatable, dimension(:) :: ac3spkr !! Diagonal derivative coefficient for salinity + real, allocatable, dimension(:) :: am3spkr !! Lower diagonal derivative coefficient for salinity + + real :: pf_direct_force = 0.9 !! phi contour above which velocity is forced exactly to zero pre-pressure solve + +contains + +!> Initialise and allocate memory for phase-field variables +subroutine InitPFVariables + + ! Main array with ghost cells + call AllocateReal3DArray(phi,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + ! Refined temperature array + call AllocateReal3DArray(tempr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + + ! Arrays without ghost cells + call AllocateReal3DArray(ruphi,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) + call AllocateReal3DArray(hphi,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) + + ! Coarse array for phi or d(phi)/dt + call AllocateReal3DArray(phic,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + + ! Second derivative coefficients + call AllocateReal1DArray(ap3spkr,1,nxr) + call AllocateReal1DArray(ac3spkr,1,nxr) + call AllocateReal1DArray(am3spkr,1,nxr) + +end subroutine InitPFVariables + +!> Deallocate the variables used for evolving the phase-field +subroutine DeallocatePFVariables + + ! Main array + call DestroyReal3DArray(phi) + + ! Array for refined temperature + call DestroyReal3DArray(tempr) + + ! Arrays without ghost cells + call DestroyReal3DArray(ruphi) + call DestroyReal3DArray(hphi) + + ! Coarse array for phi or d(phi)/dt + call DestroyReal3DArray(phic) + + ! Second derivative coefficients + call DestroyReal1DArray(ap3spkr) + call DestroyReal1DArray(ac3spkr) + call DestroyReal1DArray(am3spkr) + +end subroutine DeallocatePFVariables + +!> Set the initial state of the phase field +subroutine CreateInitialPhase + use afid_salinity, only: sal, PraS + real :: A, B, alpha, t0, x0, h0 + integer :: i, j, k + + if (salinity) then + !! Ice above salty water (1D_DDMelting example) + call set_multicomponent_interface(0.8, h0) + call set_flat_interface(h0, .true.) + else + !! 1D freezing/moving example + !! (RayT > 0: melting; RayT < 0: freezing) + if (pf_IC==1) then + h0 = 0.1 + call set_flat_interface(h0, .false.) + call set_temperature_interface(h0, .false.) + + !! 2D axisymmetric melting example (disc radius 0.1) + else if (pf_IC==2) then + call set_ice_disc(r0=0.1) + + !! Favier (2019) appendix A.3 validation case + elseif (pf_IC==3) then + call set_flat_interface(0.5, .true.) + call add_temperature_mode(amp=0.1, ymode=2, zmode=2) + + !! 1D supercooling example + elseif (pf_IC==5) then + h0 = 0.02 + call set_flat_interface(h0, .false.) + call set_temperature_interface(h0, .true.) + end if + end if + +end subroutine CreateInitialPhase + +!> Read a simple input file with three parameters definining +!! the initial state for multicomponent melting +subroutine read_phase_field_params(A, B, alpha) + real, intent(out) :: A, B, alpha + + integer :: io + logical :: exists + + inquire(file="pfparam.in", exist=exists) + if (exists) then + open(newunit=io, file="pfparam.in", status="old", action="read") + read(io, *) A, B, alpha + close(io) + else + A = 1.132 + B = 0.3796 + alpha = 3.987e-2 + end if +end subroutine read_phase_field_params + +!> Impose a flat interface at height x=h0 +!! with ice phase for x>h0 if ice_above=.true. and vice versa +subroutine set_flat_interface(h0, ice_above) + real, intent(in) :: h0 !! height (in x) of interface + logical, intent(in) :: ice_above !! flag to determine whether ice phase is above or below interface + + integer :: i, j, k + + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + if (ice_above) then + phi(k,j,i) = 0.5*(1.0 + tanh((xmr(k) - h0)/2/pf_eps)) + else + phi(k,j,i) = 0.5*(1.0 - tanh((xmr(k) - h0)/2/pf_eps)) + end if + end do + end do + end do +end subroutine set_flat_interface + +!> Impose a diffusive boundary layer profile for temperature +!! on one side of the flat interface at height x=h0 +subroutine set_temperature_interface(h0, diffuse_above) + use local_arrays, only: temp + real, intent(in) :: h0 !! Interface position (at sim time t=0) + logical, intent(in) :: diffuse_above !! Flag: set diffusion above interface (or not) + + real :: t0, Lambda + integer :: i, j, k + + ! Set normalising similarity value (see docs for justification) + if (diffuse_above) then + Lambda = 0.060314 + else + Lambda = 0.620063 + end if + + ! Effective time assuming that initial interface diffused from + ! x=0 at t=0 + t0 = PecT*(h0/2.0/Lambda)**2 + + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + if (diffuse_above) then + !! For the 1D supercooling example + if (xm(k) > h0) then + temp(k,j,i) = erfc(xm(k)*sqrt(pect/t0)/2.0)/erfc(Lambda) + else + temp(k,j,i) = 1.0 + end if + else + !! For the 1D freezing example + if (xm(k) < h0) then + temp(k,j,i) = erf(xm(k)*sqrt(pect/t0)/2)/erf(Lambda) + else + temp(k,j,i) = 1.0 + end if + !! For the 1D melting case + if (RayT > 0) temp(k,j,i) = 1.0 - temp(k,j,i) + end if + end do + end do + end do +end subroutine set_temperature_interface + +!> Set temperature and salinity profiles to a diffusive boundary layer below +!! an ice interface that began as a step profile at height x0 +subroutine set_multicomponent_interface(x0, h0) + use local_arrays, only: temp + use afid_salinity, only: sal, PraS + real, intent(in) :: x0 !! initial interface height (at diffusion start) + real, intent(out) :: h0 !! interface height at simulation time 0 + + real :: t0, A, B, alpha + integer :: i, j, k + + call read_phase_field_params(A, B, alpha) + t0 = 1e-3 + h0 = x0 + 2*alpha*sqrt(t0) + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + sal(k,j,i) = 1.0 - B*erfc((x0 - xmr(k))/sqrt(PraT/PraS*t0)/2.0) + end do + end do + end do + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + if (xm(k) <= h0) then + temp(k,j,i) = 1 - A*erfc((x0 - xm(k))/sqrt(t0)/2.0) + else + temp(k,j,i) = 1 - A*erfc(-alpha) + end if + end do + end do + end do +end subroutine set_multicomponent_interface + +!> Add a modal perturbation to the temperature field in the lower half of the domain +!! Intended for use with phase-field validation case of 2D Melting RBC from Favier et al (2019) +!! Vertical structure is of the form sin^2(2 pi x) +subroutine add_temperature_mode(amp, ymode, zmode) + use local_arrays, only: temp + real, intent(in) :: amp !! amplitude of perturbation + integer, intent(in) :: ymode !! mode number in y + integer, intent(in) :: zmode !! mode number in z (only used if nzm>1) + integer :: i, j, k, kmid + real :: xxx, yyy, zzz + + kmid = nxm/2 + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + if (nzm > 1) then + do k=1,kmid ! If domain 3D, add in z perturbation too + xxx = xm(k) + yyy = ym(j)/ylen + zzz = zm(i)/zlen + + temp(k,j,i) = temp(k,j,i) & + + amp*sin(2.0*pi*ymode*yyy)*cos(2.0*pi*zmode*zzz)*sin(2.0*pi*xxx)**2 + end do + else + do k=1,kmid + xxx = xm(k) + yyy = ym(j)/ylen + + temp(k,j,i) = temp(k,j,i) & + + amp*sin(2.0*pi*ymode*yyy)*sin(2.0*pi*xxx)**2 + end do + end if + end do + end do +end subroutine add_temperature_mode + +!> Add an ice disc at the centre of the domain of radius r0 +!! with a corresponding temperature profile +subroutine set_ice_disc(r0) + use local_arrays, only: temp + real, intent(in) :: r0 !! Initial disc radius + + real :: r + integer :: i, j, k + + ! Temperature field (0 in disc, 1 out of disc, tanh interface width approx 1e-2) + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + r = sqrt((xm(k) - 0.5*alx3)**2 + (ym(j) - 0.5*ylen)**2) + temp(k,j,i) = 0.5*(1.0 + tanh(100.0*(r - r0))) + end do + end do + end do + ! Phase-field (0 out of disc, 1 in disc) + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + r = sqrt((xmr(k) - 0.5*alx3)**2 + (ymr(j) - 0.5*ylen)**2) + phi(k,j,i) = 0.5*(1.0 - tanh(0.5*(r - r0)/pf_eps)) + end do + end do + end do +end subroutine set_ice_disc + +!> Compute the explicit terms for the phase-field evolution +!! and store the result in `hphi` +subroutine ExplicitPhase + integer :: ic, jc, kc + integer :: im, jm + integer :: ip, jp + + real :: udyrq, udzrq + real :: pf_B, nlphi + real :: dyyp, dzzp + + ! Nonlinear term prefactor + pf_B = pf_D/(pf_eps)**2 + + ! Diffusion coefficients + udzrq = pf_D*dzqr + udyrq = pf_D*dyqr + + do ic=xstartr(3),xendr(3) + im = ic - 1 + ip = ic + 1 + do jc=xstartr(2),xendr(2) + jm = jc - 1 + jp = jc + 1 + do kc=1,nxmr + ! yy second derivative of phi + dyyp = (phi(kc,jp,ic) - 2.0*phi(kc,jc,ic) + phi(kc,jm,ic))*udyrq + ! zz second derivative of phi + dzzp = (phi(kc,jc,ip) - 2.0*phi(kc,jc,ic) + phi(kc,jc,im))*udzrq + ! Extra nonlinear terms + nlphi = pf_B*phi(kc,jc,ic)*(1.0 - phi(kc,jc,ic)) & + *(1.0 - 2.0*phi(kc,jc,ic) + pf_A*(tempr(kc,jc,ic) - pf_Tm)) + + hphi(kc,jc,ic) = dyyp + dzzp - nlphi + end do + end do + end do + +end subroutine ExplicitPhase + +!> Force the velocity field to be zero in the ice phase +!! Apply this only above the phase-field level pf_direct_force +subroutine ForceIceVelZero + use local_arrays, only: vx, vy, vz + integer :: i, j, k, im, jm + + do i=xstart(3),xend(3) + im = i - 1 + do j=xstart(2),xend(2) + jm = j - 1 + do k=2,nxm + if (0.5*(phic(k,j,i) + phic(k-1,j,i)) > pf_direct_force) then + vx(k,j,i) = 0.0 + end if + end do + do k=1,nxm + if (0.5*(phic(k,j,i) + phic(k,jm,i)) > pf_direct_force) then + vy(k,j,i) = 0.0 + end if + if (0.5*(phic(k,j,i) + phic(k,j,im)) > pf_direct_force) then + vz(k,j,i) = 0.0 + end if + end do + end do + end do + +end subroutine ForceIceVelZero + +!> Compute the implicit terms for the phase-field evolution +subroutine ImplicitPhase + integer :: jc,kc,ic + real :: alpec,dxxp + + alpec=al*pf_D + + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + + ! Second xx derivative + ! Apply lower BC (d/dx(phi)=0) + if (kc.eq.1) then + dxxp = phi(kc+1,jc, ic)*ap3spkr(kc) & + + phi(kc ,jc, ic)*ac3spkr(kc) + ! Apply upper BC (d/dx(phi)=0) + elseif(kc.eq.nxmr) then + dxxp = phi(kc ,jc,ic)*ac3spkr(kc) & + + phi(kc-1,jc,ic)*am3spkr(kc) + else + dxxp = phi(kc+1,jc,ic)*ap3spkr(kc) & + + phi(kc ,jc,ic)*ac3spkr(kc) & + + phi(kc-1,jc,ic)*am3spkr(kc) + end if + + rhsr(kc,jc,ic) = (ga*hphi(kc,jc,ic) + ro*ruphi(kc,jc,ic) + alpec*dxxp)*dt + + ruphi(kc,jc,ic) = hphi(kc,jc,ic) + + end do + end do + end do + +! Solve equation and update salinity + call SolveImpEqnUpdate_Phi + +end subroutine ImplicitPhase + +!> Solve the implicit system for the phase-field +!! and update the global variable phi +subroutine SolveImpEqnUpdate_Phi + integer :: jc,kc,info,ipkv(nxr),ic,nrhs + real :: betadx,ackl_b + real :: amkT(nxmr-1),ackT(nxmr),apkT(nxmr-1),appk(nxmr-2) + + betadx=0.5d0*al*dt*pf_D + + ! Construct tridiagonal matrix for LHS + ! (normalised to prevent floating point errors) + do kc=1,nxmr + ackl_b=1.0d0/(1.-ac3spkr(kc)*betadx) + if (kc > 1) amkT(kc-1) = -am3spkr(kc)*betadx*ackl_b + ackT(kc)=1.0d0 + if (kc < nxmr) apkT(kc) = -ap3spkr(kc)*betadx*ackl_b + end do + + ! Factor the tridiagonal matrix + call dgttrf(nxmr,amkT,ackT,apkT,appk,ipkv,info) + + ! Rescale RHS to match rescaling of LHS + nrhs=(xendr(3)-xstartr(3)+1)*(xendr(2)-xstartr(2)+1) + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + ackl_b=1.0/(1.0-ac3spkr(kc)*betadx) + rhsr(kc,jc,ic)=rhsr(kc,jc,ic)*ackl_b + end do + end do + end do + + ! Solve tridiagonal system + call dgttrs('N',nxmr,nrhs,amkT,ackT,apkT,appk,ipkv,rhsr,nxmr,info) + + ! Update global variable + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + phi(kc,jc,ic) = phi(kc,jc,ic) + rhsr(kc,jc,ic) + end do + end do + end do + +end subroutine SolveImpEqnUpdate_Phi + +!> Interpolate the phase-field onto the coarse grid to provide +!! volume penalty forcing for the momentum equation +subroutine InterpPhiMultigrid + integer :: ic,jc,kc + + ! Set coarse phase-field array to zero + phic(:,:,:) = 0.d0 + + ! Construct temporary array with extended range for interpolation + ! (using dphi/dx = 0 BC) + tpdvr(:,:,:) = 0.d0 + do ic=xstartr(3)-lvlhalo,xendr(3)+lvlhalo + do jc=xstartr(2)-lvlhalo,xendr(2)+lvlhalo + tpdvr(0,jc,ic) = phi(1,jc,ic) + do kc=1,nxmr + tpdvr(kc,jc,ic) = phi(kc,jc,ic) + end do + tpdvr(nxr,jc,ic) = phi(nxmr,jc,ic) + end do + end do + + ! Interpolate the refined field to the coarse grid, storing in phic + if ((xmr(1) < xm(1)) .and. (xmr(nxmr) > xm(nxm))) then + call interpolate_xyz_to_coarse_fast(tpdvr, phic(1:nxm,:,:), "phi") + else + call interpolate_xyz_to_coarse(tpdvr, phic(1:nxm,:,:)) + end if + +end subroutine InterpPhiMultigrid + +!> Interpolate the temperature field onto the refined grid +!! for calculation of nonlinear terms in phase-field equation +subroutine InterpTempMultigrid + use local_arrays, only: temp, tempbp, temptp + integer :: ic,jc,kc + + tempr(:,:,:) = 0.d0 + tpdv(:,:,:) = 0.d0 + + ! Fill temporary array with temperature field and BCs + do ic=xstart(3)-lvlhalo,xend(3)+lvlhalo + do jc=xstart(2)-lvlhalo,xend(2)+lvlhalo + tpdv(0,jc,ic) = 2.0*tempbp(1,jc,ic) - temp(1,jc,ic) + tpdv(nx,jc,ic) = 2.0*temptp(1,jc,ic) - temp(nxm,jc,ic) + do kc=1,nxm + tpdv(kc,jc,ic) = temp(kc,jc,ic) + end do + end do + end do + + ! Interpolate to the refined grid, storing in tempr + call interpolate_xyz_to_refined(tpdv, tempr(1:nxmr,:,:)) + +end subroutine InterpTempMultigrid + +!> Add volume penalty forcing to each component of the momentum +!! equation by modifying the RK explicit terms +subroutine AddVolumePenalty + use local_arrays, only: qcap, dph, dq, vx, vy, vz + real :: pf_eta, heta, volpen + integer :: ic, jc, kc, jmm, imm + + ! Set optimal damping factor following Hester et al (2020) + pf_eta = ren*(1.51044285*pf_eps)**2 + heta = 0.5d0/pf_eta + + ! x-component + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=2,nxm + volpen = (phic(kc,jc,ic) + phic(kc-1,jc,ic))*vx(kc,jc,ic)*heta + qcap(kc,jc,ic) = qcap(kc,jc,ic) - volpen + end do + end do + end do + + ! y-component + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + jmm = jc - 1 + do kc=1,nxm + volpen = (phic(kc,jc,ic) + phic(kc,jmm,ic))*vy(kc,jc,ic)*heta + dph(kc,jc,ic) = dph(kc,jc,ic) - volpen + end do + end do + end do + + ! z-component + do ic=xstart(3),xend(3) + imm = ic - 1 + do jc=xstart(2),xend(2) + do kc=1,nxm + volpen = (phic(kc,jc,ic) + phic(kc,jc,imm))*vz(kc,jc,ic)*heta + dq(kc,jc,ic) = dq(kc,jc,ic) - volpen + end do + end do + end do + +end subroutine AddVolumePenalty + +!> Add the normal interfacial salt flux term to the +!! explicit terms of the salinity evolution when using +!! the phase-field method +!! -1/Pe_S (grad(phi).grad(S)/(1 - phi + delta)) +subroutine AddSaltFluxInterface + use afid_salinity + integer :: ic, jc, kc + integer :: im, jm, km + integer :: ip, jp, kp + real :: sdx(nxmr), sqdx(nxmr), sqdy, sqdz, kaps + real :: dpdsx, dpdsy, dpdsz, phfrac + real :: pf_delta = 1e-6 + + do kc=1,nxmr + sdx(kc) = 0.5*dxr/g3rmr(kc) + sqdx(kc) = sdx(kc)**2 + end do + sqdy = (0.5d0*dyr)**2 + sqdz = (0.5d0*dzr)**2 + kaps = 1.d0/pecs + + do ic=xstartr(3),xendr(3) + im = ic - 1 + ip = ic + 1 + do jc=xstartr(2),xendr(2) + jm = jc - 1 + jp = jc + 1 + do kc=1,nxmr + km = kc - 1 + kp = kc + 1 + + if (kc==1) then + dpdsx = sqdx(kc)*(phi(kp,jc,ic) - phi(kc,jc,ic))& + *(sal(kp,jc,ic) - sal(kc,jc,ic) & + + 2.0*SfixS*(sal(kc,jc,ic) - salbp(1,jc,ic))) + elseif (kc==nxmr) then + dpdsx = sqdx(kc)*(phi(kc,jc,ic) - phi(km,jc,ic)) & + *(sal(kc,jc,ic) - sal(km,jc,ic) & + + 2.0*SfixN*(saltp(1,jc,ic) - sal(kc,jc,ic))) + else + dpdsx = sqdx(kc)*(phi(kp,jc,ic) - phi(km,jc,ic)) & + *(sal(kp,jc,ic) - sal(km,jc,ic)) + end if + dpdsy = sqdy*(phi(kc,jp,ic) - phi(kc,jm,ic))& + *(sal(kc,jp,ic) - sal(kc,jm,ic)) + dpdsz = sqdz*(phi(kc,jc,ip) - phi(kc,jc,im))& + *(sal(kc,jc,ip) - sal(kc,jc,im)) + + phfrac = 1.0/(1.0 - phi(kc,jc,ic) + pf_delta) + + hsal(kc,jc,ic) = hsal(kc,jc,ic) - phfrac*kaps*(dpdsx + dpdsy + dpdsz) + end do + end do + end do +end subroutine AddSaltFluxInterface + +!> Add forcing terms to rhs of phase-field equation when +!! using salinity to reflect change in local melting temperature +!! with salt concentration +!! -D/eps^2 phi(1 - phi) A\Lambda S +subroutine AdjustMeltPoint + use afid_salinity, only: sal + integer :: ic, jc, kc + real :: bcl + + ! Nonlinear term prefactor + bcl = pf_D/(pf_eps)**2*pf_A*pf_Lambda + + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + hphi(kc,jc,ic) = hphi(kc,jc,ic) & + - bcl*sal(kc,jc,ic)*phi(kc,jc,ic)*(1.0 - phi(kc,jc,ic)) + end do + end do + end do +end subroutine AdjustMeltPoint + +!> Use the result from the implicit solution of phi to compute a latent heat +!! term, interpolate it onto the coarse grid, and add it to the temperature +!! equation. Interpolation is explicitly written out here, rather than using a +!! function since we are adding to the RK array hro +!! When this routine is called, rhsr should contain phi(l+1)-phi(l) = d/dt(phi)*dt +subroutine AddLatentHeat + use local_arrays, only: hro + integer :: ic,jc,kc,icr,jcr,kcr + + real, dimension(4,4,4) :: qv3 + real, dimension(4,4) :: qv2 + real, dimension(4) :: qv1 + + real :: phi_rhs, aldt + + phi_rhs = 0.d0 + aldt = 1.0/al/dt + + tpdvr(:,:,:) = 0.d0 ! Temporary array with extended range for interpolation + + ! Fill temporary array with dphi/dt + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + !CJH Note if d/dx(phi)==0 on boundary, then + ! d/dx(dphi/dt)=0 on boundary + tpdvr(0,jc,ic) = rhsr(1,jc,ic) + do kc=1,nxmr + tpdvr(kc,jc,ic) = rhsr(kc,jc,ic) + end do + tpdvr(nxr,jc,ic) = rhsr(nxmr,jc,ic) + end do + end do + + ! Fill in halo values + call update_halo(tpdvr,lvlhalo) + + ! If the grid is sufficiently fine at the boundary, use a more efficient interpolation + if ((xmr(1) < xm(1)) .and. (xmr(nxmr) > xm(nxm))) then + do ic=xstart(3),xend(3) + icr = krangs(ic) + do jc=xstart(2),xend(2) + jcr = jrangs(jc) + do kc=1,nxm + kcr = irangs(kc) + + qv3 = tpdvr(kcr-2:kcr+1,jcr-2:jcr+1,icr-2:icr+1) + qv2(:,:) = qv3(:,:,1)*czsalc(1,ic) + qv3(:,:,2)*czsalc(2,ic) & + + qv3(:,:,3)*czsalc(3,ic) + qv3(:,:,4)*czsalc(4,ic) + qv1(:) = qv2(:,1)*cysalc(1,jc) + qv2(:,2)*cysalc(2,jc) & + + qv2(:,3)*cysalc(3,jc) + qv2(:,4)*cysalc(4,jc) + + phi_rhs = sum(qv1(1:4)*cxsalc(1:4,kc)) + hro(kc,jc,ic) = hro(kc,jc,ic) + pf_S*phi_rhs*aldt + end do + end do + end do + else + do icr=xstartr(3)-1,xendr(3) + do jcr=xstartr(2)-1,xendr(2) + do kcr=0,nxmr + + qv3 = tpdvr(kcr-1:kcr+2,jcr-1:jcr+2,icr-1:icr+2) + do ic=max(krangr(icr),xstart(3)),min(krangr(icr+1)-1,xend(3)) + qv2(:,:) = qv3(:,:,1)*czsalc(1,ic) + qv3(:,:,2)*czsalc(2,ic)& + +qv3(:,:,3)*czsalc(3,ic) + qv3(:,:,4)*czsalc(4,ic) + do jc=max(jrangr(jcr),xstart(2)),min(jrangr(jcr+1)-1,xend(2)) + qv1(:) = qv2(:,1)*cysalc(1,jc) + qv2(:,2)*cysalc(2,jc) & + +qv2(:,3)*cysalc(3,jc) + qv2(:,4)*cysalc(4,jc) + do kc=max(irangr(kcr),1),min(irangr(kcr+1)-1,nxm) + phi_rhs = sum(qv1(1:4)*cxsalc(1:4,kc)) + + hro(kc,jc,ic) = hro(kc,jc,ic) + pf_S*phi_rhs*aldt + end do + end do + end do + end do + end do + end do + end if +end subroutine AddLatentHeat + +!> Add "latent salt" term to the RK forcing array for salinity (hsal), +!! having calculated d/dt(phi) from the implicit solve and stored it in rhsr +subroutine AddLatentSalt + use afid_salinity, only: sal, hsal + real :: aldt, phfrac + real :: pf_delta = 1e-6 + integer :: ic, jc, kc + + aldt = 1.0/al/dt + + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + phfrac = 1.0/(1.0 - phi(kc,jc,ic) + pf_delta) + hsal(kc,jc,ic) = hsal(kc,jc,ic) & + + phfrac*sal(kc,jc,ic)*rhsr(kc,jc,ic)*aldt + end do + end do + end do +end subroutine AddLatentSalt + +!> Calculate vertical profiles of phase-field-related statistics and +!! store them in means.h5 +subroutine CalcPhiStats + + real, dimension(nxmr) :: phibar !! Horizontally-averaged phase-field variable + real, dimension(nxmr) :: phirms !! Horizontally-averaged rms phase-field + + real :: inyzm !! 1.0/nymr/nzmr + + character(30) :: dsetname !! Dataset name for HDF5 file + character(30) :: filename !! HDF5 file name for statistic storage + character( 5) :: nstat !! Character string of statistic index + + integer :: i, j, k + + inyzm = 1.0/nymr/nzmr + + filename = trim("outputdir/means.h5") + + phibar(:) = 0.0; phirms(:) = 0.0; + + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + phibar(k) = phibar(k) + phi(k,j,i) + phirms(k) = phirms(k) + phi(k,j,i)**2 + end do + end do + end do + + call MpiSumReal1D(phibar,nxmr) + call MpiSumReal1D(phirms,nxmr) + + do k=1,nxmr + phibar(k) = phibar(k)*inyzm + phirms(k) = sqrt(phirms(k)*inyzm) + end do + + ! Store index as character string + write(nstat,"(i5.5)")nint(time/tout) + + if (ismaster) then + dsetname = trim("phibar/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, phibar, nxmr) + dsetname = trim("phirms/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, phirms, nxmr) + end if + + call MpiBarrier + +end subroutine CalcPhiStats + +!> Create the groups in the means.h5 file to store the +!! statistics related to the phase-field +subroutine CreatePhaseH5Groups(filename) + use HDF5 + + character(30), intent(in) :: filename + integer(HID_T) :: file_id, group_id + integer :: hdf_error + + call h5fopen_f(filename, H5F_ACC_RDWR_F, file_id, hdf_error) + + call h5gcreate_f(file_id, "phibar", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "phirms", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + + call h5fclose_f(file_id, hdf_error) + +end subroutine CreatePhaseH5Groups + +end module afid_phasefield \ No newline at end of file diff --git a/src/pressure.F90 b/src/pressure.F90 new file mode 100644 index 00000000..547eb5b9 --- /dev/null +++ b/src/pressure.F90 @@ -0,0 +1,380 @@ +!> Module containing all subroutines and arrays needed to compute the +!! pressure correction needed to ensure a divergence-free velocity field +module afid_pressure + use param + use fftw_params + use local_arrays, only: pr, dph, dphhalo + use decomp_2d + use decomp_2d_fft + use AuxiliaryRoutines + use mpih + implicit none + + logical :: sidewall !! Flag to determine whether to impose sidewalls in y (using a DCT) + + real, allocatable, dimension(:) :: ak1 !! Modified wavenumber in z + real, allocatable, dimension(:) :: ak2 !! Modified wavenumber in y + + real, allocatable, dimension(:) :: amphk !! Lower diagonal Poisson coefficient for pressure solve + real, allocatable, dimension(:) :: acphk !! Centre diagonal Poisson coefficient for pressure solve + real, allocatable, dimension(:) :: apphk !! Upper diagonal Poisson coefficient for pressure solve + +contains + +!> Allocate the arrays used for the solution of the pressure field +subroutine InitPressureVars + + !! Modified wavenumbers + call AllocateReal1DArray(ak1,1,nz) + call AllocateReal1DArray(ak2,1,ny) + + !! Poisson solver coefficients + call AllocateReal1DArray(amphk,1,nx) + call AllocateReal1DArray(acphk,1,nx) + call AllocateReal1DArray(apphk,1,nx) + +end subroutine InitPressureVars + +!> Free up memory from arrays used in pressure correction +subroutine DeallocatePressureVars + + !! Modified wavenumbers + call DestroyReal1DArray(ak1) + call DestroyReal1DArray(ak2) + + !! Poisson solver coefficients + call DestroyReal1DArray(amphk) + call DestroyReal1DArray(acphk) + call DestroyReal1DArray(apphk) + + !! Temporary mid-FFT arrays + if(allocated(dphc)) deallocate(dphc) + if(allocated(rz1)) deallocate(rz1) + if(allocated(cz1)) deallocate(cz1) + if(allocated(ry1)) deallocate(ry1) + if(allocated(cy1)) deallocate(cy1) + +end subroutine DeallocatePressureVars + +!> Compute the metric terms and modified wavenumbers used in the +!! solution of the pressure correction +subroutine InitPressureSolver + integer :: nymh, nymp, i, j, nzmh, nzmp + integer :: kc, km, kp + real :: ugmmm, a33icc, a33icp + real :: ao(1:nz), ap(1:ny) + + !! Index counters for wavenumbers + nymh=nym/2+1 + nymp=nymh+1 + nzmh=nzm/2+1 + nzmp=nzmh+1 + + ! call AllocateReal1DArray(ao,1,nz) + ! call AllocateReal1DArray(ap,1,ny) + + !! Construct modified wavenumbers + !! z + do i=1,nzmh + ao(i)=(i-1)*2.d0*pi + end do + do i=nzmp,nzm + ao(i)=-(nzm-i+1)*2.d0*pi + end do + do i=1,nzm + ak1(i)=2.d0*(1.d0-dcos(ao(i)/nzm))*(float(nzm)/zlen)**2 + end do + + !! y + do j=1,nymh + ap(j)=(j-1)*2.d0*pi + end do + do j=nymp,nym + ap(j)=-(nym-j+1)*2.d0*pi + end do + do j=1,nym + ak2(j)=2.d0*(1.d0-dcos(ap(j)/nym))*(float(nym)/ylen)**2 + end do + + ! call DestroyReal1DArray(ao) + ! call DestroyReal1DArray(ap) + + !! Initialise tridiagonal matrices for Poisson solver + do kc=1,nxm + km = kmv(kc) + kp = kpv(kc) + + a33icc = kmc(kc)*dxq/d3xm(km) + a33icp = kpc(kc)*dxq/d3xm(kc) + ugmmm = 1.0d0/d3xc(kc) + + amphk(kc) = a33icc*ugmmm + apphk(kc) = a33icp*ugmmm + acphk(kc) = -(amphk(kc)+apphk(kc)) + end do + + !! Initialise pencil transposes for the pressure solver + call decomp_2d_fft_init + + !! Allocate temporary mid-FFT arrays + allocate(ry1(ph%yst(1):ph%yen(1), & + ph%yst(2):ph%yen(2), & + ph%yst(3):ph%yen(3))) + allocate(rz1(ph%zst(1):ph%zen(1), & + ph%zst(2):ph%zen(2), & + ph%zst(3):ph%zen(3))) + allocate(cy1(sp%yst(1):sp%yen(1), & + sp%yst(2):sp%yen(2), & + sp%yst(3):sp%yen(3))) + allocate(cz1(sp%zst(1):sp%zen(1), & + sp%zst(2):sp%zen(2), & + sp%zst(3):sp%zen(3))) + allocate(dphc(sp%xst(1):sp%xen(1), & + sp%xst(2):sp%xen(2), & + sp%xst(3):sp%xen(3))) + +end subroutine InitPressureSolver + +!> Compute the local divergence of the intermediate velocity field at +!! every point for the pressure correction step +!! div(u)/al/dt is stored in the array dph after this step +subroutine CalcLocalDivergence + use local_arrays, only: vx, vy, vz + + integer :: ic, ip, jc, jp, kc, kp + real :: usdtal, dqcap + + usdtal = 1.d0/(dt*al) + + do ic=xstart(3),xend(3) + ip = ic+1 + do jc=xstart(2),xend(2) + jp = jc+1 + do kc=1,nxm + kp = kc+1 + dqcap = (vz(kc,jc,ip)-vz(kc,jc,ic))*dz & + +(vy(kc,jp,ic)-vy(kc,jc,ic))*dy & + +(vx(kp,jc,ic)-vx(kc,jc,ic))*udx3m(kc) + dph(kc,jc,ic)=dqcap*usdtal + end do + end do + end do +end subroutine CalcLocalDivergence + +!> Compute the pressure correction by solving a Poisson equation +subroutine SolvePressureCorrection + integer :: nymh + + nymh=nym/2+1 + + !! Compute a 2D Fourier transform in y and z + call transpose_x_to_y(dph,ry1,ph) + ! If calling for the first time, plan the Fourier transform + if (.not. planned) call PlanFourierTransform + call dfftw_execute_dft_r2c(fwd_guruplan_y,ry1,cy1) + call transpose_y_to_z(cy1,cz1,sp) + call dfftw_execute_dft(fwd_guruplan_z,cz1,cz1) + + ! Normalise the transformed array. FFTW does not do this automatically + cz1 = cz1 / (nzm*nym) + + ! Transpose back to an x-pencil before the tridiagonal solve + call transpose_z_to_x(cz1,dphc,sp) + + ! Solve the tridiagonal system with complex coefficients + call SolveTridiagonalPressure + + call transpose_x_to_z(dphc,cz1,sp) + + call dfftw_execute_dft(bwd_guruplan_z,cz1,cz1) + + call transpose_z_to_y(cz1,cy1,sp) + + call dfftw_execute_dft_c2r(bwd_guruplan_y,cy1,ry1) + + call transpose_y_to_x(ry1,dph,ph) + +end subroutine SolvePressureCorrection + +!> Solve the tridiagonal system for the x-derivatives of the Poisson equation +!! for the Fourier transformed pressure +subroutine SolveTridiagonalPressure + integer :: i, j, k, info + complex :: acphT_b + complex :: appph(nxm-2) + complex :: acphT(nxm), apph(nxm-1), amph(nxm-1) + integer :: phpiv(nxm) + + ! Construct tridiagonal system for each Fourier mode + do i=sp%xst(3),sp%xen(3) + do j=sp%xst(2),sp%xen(2) + do k=1,nxm + ! Normalise RHS & LHS to avoid floating point errors + acphT_b = 1.0/(acphk(k) - ak2(j) - ak1(i)) + dphc(k,j,i) = dphc(k,j,i)*acphT_b + if (k < nxm) apph(k ) = apphk(k)*acphT_b + if (k > 1) amph(k-1) = amphk(k)*acphT_b + ! Small perturbation needed to prevent singular matrix + ! when using uniform grid: + acphT(k) = 1.0d0 + 1.0d-15 + end do + + ! Factor the tridiagonal matrix + call zgttrf(nxm, amph, acphT, apph, appph, phpiv, info) + + if (info.gt.0) then + print*,'Singular value found in LAPACK routine zgttrf: info=',info + print*,'Please try to adjust either NX or STR3 in bou.in' + call MPI_Abort(MPI_COMM_WORLD,1,ierr) + endif + + ! Solve the tridiagonal system + call zgttrs('N',nxm,1,amph,acphT,apph,appph,phpiv,dphc(1:nxm,j,i),nxm,info) + + end do + end do +end subroutine SolveTridiagonalPressure + +!> Remove the divergent component of velocity from the velocity field +!! using the solved pressure correction +subroutine RemoveDivergence + use local_arrays, only: vx, vy, vz + real :: usukm, udy, udz, locdph + integer :: ic, im, jc, jm, kc, km + + udy = al*dt*dy + udz = al*dt*dz + + do ic=xstart(3),xend(3) + im = ic-1 + do jc=xstart(2),xend(2) + jm = jc-1 + do kc=1,nxm + km = kmv(kc) + + usukm = al*dt*udx3c(kc) + locdph = dphhalo(kc,jc,ic) + + vx(kc,jc,ic) = vx(kc,jc,ic) & + - (locdph - dphhalo(km,jc,ic))*usukm + vy(kc,jc,ic) = vy(kc,jc,ic) & + - (locdph - dphhalo(kc,jm,ic))*udy + vz(kc,jc,ic) = vz(kc,jc,ic) & + - (locdph - dphhalo(kc,jc,im))*udz + end do + end do + end do +end subroutine RemoveDivergence + +!> Update the total pressure using the solved pressure correction +subroutine CorrectPressure + integer :: ic, jc, kc + integer :: im, jm, km + integer :: ip, jp, kp + real :: be, amm, acc, app + + be = al*beta + do ic=xstart(3),xend(3) + im = ic - 1 + ip = ic + 1 + do jc=xstart(2),xend(2) + jm = jc - 1 + jp = jc + 1 + do kc=1,nxm + kp = kpv(kc) + km = kmv(kc) + + amm = amphk(kc) + acc = acphk(kc) + app = apphk(kc) + + pr(kc,jc,ic) = pr(kc,jc,ic) + dphhalo(kc,jc,ic) - be*( & + ! dzz(p') + (dphhalo(kc,jc,ip) - 2.0*dphhalo(kc,jc,ic) + dphhalo(kc,jc,im))*dzq + & + ! dyy(p') + (dphhalo(kc,jp,ic) - 2.0*dphhalo(kc,jc,ic) + dphhalo(kc,jm,ic))*dyq + & + ! dxx(p') + (dphhalo(kp,jc,ic)*app + dphhalo(kc,jc,ic)*acc + dphhalo(km,jc,ic)*amm)) + end do + end do + end do +end subroutine CorrectPressure + +!> Create FFTW plan +subroutine PlanFourierTransform + + type(fftw_iodim),dimension(1) :: iodim + type(fftw_iodim),dimension(2) :: iodim_howmany + integer :: info + + iodim(1) % n = nzm + iodim(1) % is = (sp%zen(1)-sp%zst(1)+1)*(sp%zen(2)-sp%zst(2)+1) + iodim(1) % os = (sp%zen(1)-sp%zst(1)+1)*(sp%zen(2)-sp%zst(2)+1) + + iodim_howmany(1) % n = (sp%zen(1)-sp%zst(1)+1) + iodim_howmany(1) % is = 1 + iodim_howmany(1) % os = 1 + iodim_howmany(2) % n = (sp%zen(2)-sp%zst(2)+1) + iodim_howmany(2) % is = (sp%zen(1)-sp%zst(1)+1) + iodim_howmany(2) % os = (sp%zen(1)-sp%zst(1)+1) + + ! Construct forward plan for z transform + fwd_guruplan_z = fftw_plan_guru_dft(1, iodim, & + 2, iodim_howmany, cz1, cz1, & + FFTW_FORWARD, FFTW_ESTIMATE) + + iodim(1) % n = nzm + ! Construct backward plan for z transform + bwd_guruplan_z = fftw_plan_guru_dft(1, iodim, & + 2,iodim_howmany,cz1,cz1, & + FFTW_BACKWARD,FFTW_ESTIMATE) + + if (.not.c_associated(bwd_guruplan_z)) then + if (ismaster) print*,'Failed to create guru plan. You should' + if (ismaster) print*,'link with FFTW3 before MKL' + if (ismaster) print*,'Please check linking order.' + call MPI_Abort(MPI_COMM_WORLD,1,info) + end if + + iodim(1) % n = nym + iodim(1) % is = ph%yen(1)-ph%yst(1)+1 + iodim(1) % os = sp%yen(1)-sp%yst(1)+1 + + iodim_howmany(1) % n = (ph%yen(1)-ph%yst(1)+1) + iodim_howmany(1) % is = 1 + iodim_howmany(1) % os = 1 + iodim_howmany(2) % n = (ph%yen(3)-ph%yst(3)+1) + iodim_howmany(2) % is = (ph%yen(1)-ph%yst(1)+1) & + * (ph%yen(2)-ph%yst(2)+1) + iodim_howmany(2) % os = (sp%yen(1)-sp%yst(1)+1) & + * (sp%yen(2)-sp%yst(2)+1) + + ! Construct forward plan for y transform + fwd_guruplan_y = fftw_plan_guru_dft_r2c(1, iodim, & + 2, iodim_howmany, ry1, cy1, & + FFTW_ESTIMATE) + + iodim(1) % n = nym + iodim(1) % is = sp%yen(1)-sp%yst(1)+1 + iodim(1) % os = ph%yen(1)-ph%yst(1)+1 + + iodim_howmany(1) % n=(sp%yen(1)-sp%yst(1)+1) + iodim_howmany(1) % is=1 + iodim_howmany(1) % os=1 + iodim_howmany(2) % n=(sp%yen(3)-sp%yst(3)+1) + iodim_howmany(2) % is=(sp%yen(1)-sp%yst(1)+1) & + * (sp%yen(2)-sp%yst(2)+1) + iodim_howmany(2) % os=(ph%yen(1)-ph%yst(1)+1) & + * (ph%yen(2)-ph%yst(2)+1) + + ! Construct backward plan for y transform + bwd_guruplan_y = fftw_plan_guru_dft_c2r(1,iodim, & + 2,iodim_howmany,cy1,ry1, & + FFTW_ESTIMATE) + + ! Save that we have planned the FFT + planned=.true. + +end subroutine PlanFourierTransform + +end module afid_pressure \ No newline at end of file diff --git a/src/salinity.F90 b/src/salinity.F90 new file mode 100644 index 00000000..6563650e --- /dev/null +++ b/src/salinity.F90 @@ -0,0 +1,712 @@ +!> Module adding evolution of salinity to AFiD +!! Salinity is evolved on a refined grid, and therefore +!! this module depends on the multiple-resolution and +!! interpolation modules +module afid_salinity + use param + use mgrd_arrays + use decomp_2d, only: xstart, xend, xstartr, xendr, update_halo + use AuxiliaryRoutines + use HermiteInterpolations, only: interpolate_xyz_to_coarse, interpolate_xyz_to_coarse_fast + use ibm_param, only: solidr + implicit none + + real, allocatable, dimension(:,:,:) :: sal !! Salinity field + real, allocatable, dimension(:,:,:) :: rusal !! RK storage array for salinity (previous substep) + real, allocatable, dimension(:,:,:) :: hsal !! RK storage array for salinity + real, allocatable, dimension(:,:,:) :: salc !! Interpolated salinity field on coarse grid + real, allocatable, dimension(:,:,:) :: vxr !! Velocity interpolated to refined grid (x component) + real, allocatable, dimension(:,:,:) :: vyr !! Velocity interpolated to refined grid (y component) + real, allocatable, dimension(:,:,:) :: vzr !! Velocity interpolated to refined grid (z component) + + real :: rays !! Solutal Rayleigh number + real :: pras !! Schmidt number (solutal Prandtl number) + real :: pecs !! Solutal Peclet number + real :: bycs !! Buoyancy prefactor for salinity + + integer :: SfixS !! Flag for whether salinity is fixed at lower plate + integer :: SfixN !! Flag for whether salinity is fixed at upper plate + + real, allocatable, dimension(:,:,:) :: salbp !! Salinity boundary value (lower plate) + real, allocatable, dimension(:,:,:) :: saltp !! Salinity boundary value (upper plate) + + real, allocatable, dimension(:) :: ap3sskr !! Upper diagonal derivative coefficient for salinity + real, allocatable, dimension(:) :: ac3sskr !! Diagonal derivative coefficient for salinity + real, allocatable, dimension(:) :: am3sskr !! Lower diagonal derivative coefficient for salinity + +contains + +!> Subroutine to allocate memory for salinity-related variables +subroutine InitSalVariables + + ! Boundary planes + call AllocateReal3DArray(salbp,1,1,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + call AllocateReal3DArray(saltp,1,1,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + + ! Main arrays with ghost cells + call AllocateReal3DArray(sal,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + call AllocateReal3DArray(vxr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + call AllocateReal3DArray(vyr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + call AllocateReal3DArray(vzr,1,nxr,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + + ! Runge-Kutta storage arrays (without ghost cells) + call AllocateReal3DArray(rusal,1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) + call AllocateReal3DArray(hsal, 1,nxr,xstartr(2),xendr(2),xstartr(3),xendr(3)) + + ! Coarse array + call AllocateReal3DArray(salc,1,nx,xstart(2)-lvlhalo,xend(2)+lvlhalo,xstart(3)-lvlhalo,xend(3)+lvlhalo) + + !CJH Needed for melt boundary condition + if (melt) then + call AllocateReal3DArray(Tplaner,1,1,xstartr(2)-lvlhalo,xendr(2)+lvlhalo,xstartr(3)-lvlhalo,xendr(3)+lvlhalo) + end if + + ! Second derivative coefficients + call AllocateReal1DArray(ap3sskr,1,nxr) + call AllocateReal1DArray(ac3sskr,1,nxr) + call AllocateReal1DArray(am3sskr,1,nxr) + +end subroutine InitSalVariables + +!> Deallocate the variables used for evolving salinity +subroutine DeallocateSalVariables + + ! Boundary planes + call DestroyReal3DArray(salbp) + call DestroyReal3DArray(saltp) + + ! Main array + call DestroyReal3DArray(sal) + call DestroyReal3DArray(vxr) + call DestroyReal3DArray(vyr) + call DestroyReal3DArray(vzr) + + ! RK arrays + call DestroyReal3DArray(rhsr) + call DestroyReal3DArray(rusal) + call DestroyReal3DArray(hsal) + + ! Coarse array + call DestroyReal3DArray(salc) + + ! Extra T slice for melt condition + if (melt) then + call DestroyReal3DArray(Tplaner) + end if + + ! Second derivative coefficients + call DestroyReal1DArray(ap3sskr) + call DestroyReal1DArray(ac3sskr) + call DestroyReal1DArray(am3sskr) + +end subroutine DeallocateSalVariables + +!> Set the values for the boundary planes of salinity +subroutine SetSalBCs + integer :: i, j + + if (rays>=0) then ! unstable S gradient + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + saltp(1,j,i) = 0.5d0 + salbp(1,j,i) = -0.5d0 + end do + end do + else ! stable S gradient + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + saltp(1,j,i) = -0.5d0 + salbp(1,j,i) = 0.5d0 + end do + end do + end if + if (phasefield) then + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + saltp(1,j,i) = 0.d0 + salbp(1,j,i) = 1.d0 + end do + end do + end if + + ! Update halo for interpolation routine + call update_halo(saltp,lvlhalo) + call update_halo(salbp,lvlhalo) + +end subroutine SetSalBCs + +!> Set initial conditions for salinity field +!! N.B. This can get overwritten by CreateInitialPhase if also using phase-field +subroutine CreateInitialSalinity + + !! Rayleigh-Taylor setup for pore-scale simulation + if (IBM) then + call SetSaltTwoLayer(h0=0.5*alx3, eps=1e-7, stable=.false.) + call AddSalinityNoise(amp=0.1, localised=.true., h0=0.5*alx3) + + !! Bounded double-diffusive convection (begin with small amplitude noise) + else if ((active_S==1) .and. (active_T==1) .and. (gAxis==1)) then + call SetZeroSalinity + call AddSalinityNoise(amp=5e-3, localised=.false.) + + !! Stratified shear layer setup + else if ((RayS < 0) .and. (RayT < 0)) then + call SetSaltTwoLayer(h0=0.5*alx3, eps=1.0, stable=.true.) + call AddSalinityNoise(amp=1e-2, localised=.true., h0=0.5*alx3) + + !! Default: linear profile + small noise (e.g. RBC, VC) + else + call SetLinearSalinity + call AddSalinityNoise(amp=5e-3, localised=.false.) + end if + +end subroutine CreateInitialSalinity + +!> Set salinity variable to linear profile between boundary values +subroutine SetLinearSalinity + integer :: i, j, k + + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + sal(k,j,i) = salbp(1,j,i) - (salbp(1,j,i) - saltp(1,j,i))*xmr(k)/xcr(nxr) + end do + end do + end do +end subroutine SetLinearSalinity + +!> Set salinity variable to zero everywhere +subroutine SetZeroSalinity + integer :: i, j, k + + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + sal(k,j,i) = 0.0 + end do + end do + end do +end subroutine SetZeroSalinity + +!> Set the salinity field up as a two-layer system with a tanh profile +!! with interface thickness eps +subroutine SetSaltTwoLayer(h0, eps, stable, mode, mode_amp) + real, intent(in) :: h0 !! Mean height of interface + real, intent(in) :: eps !! Width of tanh interface + logical, intent(in) :: stable !! Flag determining gravitational stability of profile + integer, intent(in), optional :: mode !! Optional mode number to perturb interface + real, intent(in), optional :: mode_amp !! Amplitude of optional modal perturbation + + real :: x0 + integer :: i, j, k + + x0 = h0 + + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + ! Use pf_IC input parameter as mode number for initial perturbation + if (present(mode)) x0 = h0 + mode_amp*sin(mode*2.0*pi*ymr(j)/ylen) + sal(k,j,i) = 0.5*tanh((xmr(k) - x0)/eps) + if (stable) sal(k,j,i) = -sal(k,j,i) + end do + end do + end do + +end subroutine SetSaltTwoLayer + +!> Add random noise to the salinity field, either locally at an interface +!! or uniformly. In both cases, noise is limited such that the absolute value +!! of salinity does not exceed 0.5 +subroutine AddSalinityNoise(amp, localised, h0) + real, intent(in) :: amp !! Amplitude of random noise + logical, intent(in) :: localised !! Flag determining whether to add noise around interface or everywhere + real, intent(in), optional :: h0 !! Height of interface if using localised noise + + integer :: i, j, k + real :: a2, varptb + + call random_seed() + + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + call random_number(varptb) + !! Add noise locally + if (localised) then + sal(k,j,i) = sal(k,j,i) + amp/cosh((xmr(k) - h0)/0.01)**2*varptb + ! Restrict initial salinity field to [-0.5,0.5] + sal(k,j,i) = min(0.5, sal(k,j,i)) + sal(k,j,i) = max(-0.5, sal(k,j,i)) + !! Add noise everywhere uniformly + else + ! Prevent values of |S| exceeding 0.5 by restricting noise amplitude locally + if (abs(sal(k,j,i)) + amp > 0.5) then + a2 = 0.5 - abs(sal(k,j,i)) + sal(k,j,i) = sal(k,j,i) + a2*(2.d0*varptb - 1.d0) + else + sal(k,j,i) = sal(k,j,i) + amp*(2.d0*varptb - 1.d0) + end if + end if + end do + end do + end do + +end subroutine AddSalinityNoise + +!> Compute the explicit terms for the salinity evolution +!! and store the result in hsal +subroutine ExplicitSalinity + integer :: ic, jc, kc + integer :: im, jm, km + integer :: ip, jp, kp + + real :: udyr, udzr, udyrq, udzrq + real :: aldt + real, dimension(1:nxmr) :: sdx + real :: hsx, hsy, hsz + real :: dyys, dzzs + + ! Advection coefficients + udyr = 0.5d0*dyr + udzr = 0.5d0*dzr + ! Diffusion coefficients + udyrq = dyqr/pecs + udzrq = dzqr/pecs + + ! x-advection coefficients + do kc=1,nxmr + sdx(kc) = 0.5*dxr/g3rmr(kc) + end do + ! Time advancing pre-factor + aldt = 1.0/al/dt + + do ic=xstartr(3),xendr(3) + im = ic - 1 + ip = ic + 1 + do jc=xstartr(2),xendr(2) + jm = jc - 1 + jp = jc + 1 + do kc=1,nxmr + km = kc - 1 + kp = kc + 1 + + ! x-advection d/dx (vx * S) + if (kc==1) then + hsx = ( & + vxr(kp,jc,ic)*(sal(kp,jc,ic) + sal(kc,jc,ic)) & + - vxr(kc,jc,ic)*2.d0*salbp(1,jc,ic) & + )*udx3mr(kc)*0.5d0 + elseif (kc==nxmr) then + hsx = ( & + vxr(kp,jc,ic)*2.d0*saltp(1,jc,ic) & + - vxr(kc,jc,ic)*(sal(kc,jc,ic) + sal(km,jc,ic)) & + )*udx3mr(kc)*0.5d0 + else + hsx = ( & + vxr(kp,jc,ic)*(sal(kp,jc,ic) + sal(kc,jc,ic)) & + - vxr(kc,jc,ic)*(sal(kc,jc,ic) + sal(km,jc,ic)) & + )*udx3mr(kc)*0.5d0 + end if + + ! y-advection d/dy(vy * S) + hsy = ( & + vyr(kc,jp,ic)*(sal(kc,jp,ic) + sal(kc,jc,ic)) & + - vyr(kc,jc,ic)*(sal(kc,jc,ic) + sal(kc,jm,ic)) & + )*udyr + + ! z-advection d/dz(vz * S) + hsz = ( & + vzr(kc,jc,ip)*(sal(kc,jc,ip) + sal(kc,jc,ic)) & + - vzr(kc,jc,ic)*(sal(kc,jc,ic) + sal(kc,jc,im)) & + )*udzr + + !! If using immersed boundary, enforce zero lateral gradient at interface + if (IBM) then + !!! ADD THIS TO IBM MODULE, CALL AS SUBROUTINE + ! yy second derivative of salinity + if (solidr(kc,jp,ic)) then + dyys = (sal(kc,jm,ic) - sal(kc,jc,ic))*udyrq + elseif (solidr(kc,jm,ic)) then + dyys = (sal(kc,jp,ic) - sal(kc,jc,ic))*udyrq + else + dyys = (sal(kc,jp,ic) - 2.0*sal(kc,jc,ic) + sal(kc,jm,ic))*udyrq + end if + + ! zz second derivative of salinity + if (solidr(kc,jc,ip)) then + dzzs = (sal(kc,jc,im) - sal(kc,jc,ic))*udzrq + elseif (solidr(kc,jc,im)) then + dzzs = (sal(kc,jc,ip) - sal(kc,jc,ic))*udzrq + else + dzzs = (sal(kc,jc,ip) - 2.0*sal(kc,jc,ic) + sal(kc,jc,im))*udzrq + end if + else + ! yy second derivative of salinity + dyys = (sal(kc,jp,ic) - 2.0*sal(kc,jc,ic) + sal(kc,jm,ic))*udyrq + ! zz second derivative of salinity + dzzs = (sal(kc,jc,ip) - 2.0*sal(kc,jc,ic) + sal(kc,jc,im))*udzrq + end if + + ! Sum explicit terms + hsal(kc,jc,ic) = -(hsx + hsy + hsz) + dyys + dzzs + end do + end do + end do + +end subroutine ExplicitSalinity + +!> Compute the implicit terms for the salinity evolution +subroutine ImplicitSalinity + integer :: ic, jc, kc + real :: dxxs, alpec + + alpec = al/pecs + + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + + ! Second xx derivative + ! Apply lower BC + if (kc==1) then + dxxs= sal(kc+1,jc,ic)*ap3sskr(kc) & + + sal(kc ,jc,ic)*ac3sskr(kc) & + - (ap3sskr(kc) + ac3sskr(kc))*salbp(1,jc,ic)*SfixS + ! Apply upper BC + elseif (kc==nxmr) then + dxxs= sal(kc ,jc,ic)*ac3sskr(kc) & + + sal(kc-1,jc,ic)*am3sskr(kc) & + - (am3sskr(kc) + ac3sskr(kc))*saltp(1,jc,ic)*SfixN + else + dxxs= sal(kc+1,jc,ic)*ap3sskr(kc) & + + sal(kc ,jc,ic)*ac3sskr(kc) & + + sal(kc-1,jc,ic)*am3sskr(kc) + end if + + rhsr(kc,jc,ic) = (ga*hsal(kc,jc,ic) + ro*rusal(kc,jc,ic) + alpec*dxxs)*dt + + rusal(kc,jc,ic) = hsal(kc,jc,ic) + end do + end do + end do + + if (IBM) then + call SolveImpEqnUpdate_Sal_ibm + else + call SolveImpEqnUpdate_Sal + end if + +end subroutine ImplicitSalinity + +!> Solve the implicit system for the salinity +!! and update the global variable sal +subroutine SolveImpEqnUpdate_Sal + real :: betadx, ackl_b + integer :: ic, jc, kc, nrhs, ipkv(nxr), info + real :: amkT(nxmr-1), ackT(nxmr), apkT(nxmr-1), appk(nxmr-2) + + betadx = 0.5d0*al*dt/pecs + + ! Construct tridiagonal matrix for LHS + do kc=1,nxmr + ackl_b = 1.0d0/(1. - ac3sskr(kc)*betadx) + if (kc > 1) amkT(kc-1) = -am3sskr(kc)*betadx*ackl_b + ackT(kc) = 1.d0 + if (kc < nxmr) apkT(kc) = -ap3sskr(kc)*betadx*ackl_b + end do + + ! Factor the tridiagonal matrix + call dgttrf(nxmr,amkT,ackT,apkT,appk,ipkv,info) + + ! Rescale RHS to match rescaling of LHS + nrhs=(xendr(3)-xstartr(3)+1)*(xendr(2)-xstartr(2)+1) + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + ackl_b = 1.0/(1.0 - ac3sskr(kc)*betadx) + rhsr(kc,jc,ic) = rhsr(kc,jc,ic)*ackl_b + end do + end do + end do + + ! Solve tridiagonal system + call dgttrs('N',nxmr,nrhs,amkT,ackT,apkT,appk,ipkv,rhsr,nxmr,info) + + ! Update global variable + do ic=xstartr(3),xendr(3) + do jc=xstartr(2),xendr(2) + do kc=1,nxmr + sal(kc,jc,ic) = sal(kc,jc,ic) + rhsr(kc,jc,ic) + end do + end do + end do + +end subroutine SolveImpEqnUpdate_Sal + +!> Interpolate the salinity field onto the coarse grid to +!! provide buoyancy forcing to the momentum equation +subroutine InterpSalMultigrid + integer :: icr, jcr, kcr + + ! Set coarse salinity array to zero + salc(:,:,:) = 0.d0 + + ! Extend refined array in wall-normal direction to give sufficient points + ! for cubic interpolation + do icr=xstartr(3)-lvlhalo,xendr(3)+lvlhalo + do jcr=xstartr(2)-lvlhalo,xendr(2)+lvlhalo + do kcr=1,nxmr + tpdvr(kcr,jcr,icr) = sal(kcr,jcr,icr) + end do + if (SfixS==1) then + tpdvr(0,jcr,icr) = 2.0*salbp(1,jcr,icr) - sal(1,jcr,icr) + else + tpdvr(0,jcr,icr) = sal(1,jcr,icr) + end if + if (SfixN==1) then + tpdvr(nxr,jcr,icr) = 2.0*saltp(1,jcr,icr) - sal(nxmr,jcr,icr) + else + tpdvr(nxr,jcr,icr) = sal(nxmr,jcr,icr) + end if + end do + end do + + ! Interpolate the refined field to the coarse grid, storing in salc + if ((xmr(1) < xm(1)) .and. (xmr(nxmr) > xm(nxm))) then + call interpolate_xyz_to_coarse_fast(tpdvr, salc(1:nxm,:,:), "sal") + else + call interpolate_xyz_to_coarse(tpdvr, salc(1:nxm,:,:)) + end if + +end subroutine InterpSalMultigrid + +!> Add buoyancy contribution from the salinity to one of the +!! momentum forcing arrays +subroutine AddSalBuoyancy(rkv) + real, dimension(:,xstart(2):,xstart(3):), intent(inout) :: rkv + integer :: ic, jc, kc + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + do kc=1,nxm + rkv(kc,jc,ic) = rkv(kc,jc,ic) - bycs*salc(kc,jc,ic) + end do + end do + end do +end subroutine + +!> Calculate and save vertical profiles related to the salinity +!! field, storing the data in means.h5 +subroutine CalcSalStats + real, dimension(nxmr) :: Sbar !! Horizontally-averaged salinity + real, dimension(nxmr) :: Srms !! Horizontally-averaged rms salinity + real, dimension(nxmr) :: chiS !! Horizontally-averaged solutal dissipation rate (nu grad(S)^2) + + real, dimension(nxmr) :: vxS !! Advective flux of salinity (x) + real, dimension(nxmr) :: vyS !! Advective flux of salinity (y) + real, dimension(nxmr) :: vzS !! Advective flux of salinity (z) + + real :: inyzmr !! 1.0/nymr/nzmr + + character(30) :: dsetname !! Dataset name for HDF5 file + character(30) :: filename !! HDF5 file name for statistic storage + character( 5) :: nstat !! Character string of statistic index + + integer :: i, j, k + + inyzmr = 1.0/nymr/nzmr + + filename = trim("outputdir/means.h5") + + Sbar(:) = 0.0; Srms(:) = 0.0; chiS(:) = 0.0 + vxS(:) = 0.0; vyS(:) = 0.0; vzS(:) = 0.0 + + if (IBM) then + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + ! Only record data from fluid phase + if (.not. solidr(k,j,i)) then + Sbar(k) = Sbar(k) + sal(k,j,i) + Srms(k) = Srms(k) + sal(k,j,i)**2 + end if + end do + end do + end do + else + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + Sbar(k) = Sbar(k) + sal(k,j,i) + Srms(k) = Srms(k) + sal(k,j,i)**2 + end do + end do + end do + end if + + ! Since velocities are zero in solid, no need to use if statement for IBM here + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + vxS(k) = vxS(k) + 0.5*(vxr(k,j,i)+vxr(k+1,j,i))*sal(k,j,i) + vyS(k) = vyS(k) + 0.5*(vyr(k,j,i)+vyr(k,j+1,i))*sal(k,j,i) + vzS(k) = vzS(k) + 0.5*(vzr(k,j,i)+vzr(k,j,i+1))*sal(k,j,i) + end do + end do + end do + + call CalcDissipationSal(chiS) + + call MpiSumReal1D(Sbar, nxmr) + call MpiSumReal1D(Srms, nxmr) + call MpiSumReal1D(vxS, nxmr) + call MpiSumReal1D(vyS, nxmr) + call MpiSumReal1D(vzS, nxmr) + call MpiSumReal1D(chiS, nxmr) + + ! Turn sums into averages + ! (and root Srms & scale chi) + do k=1,nxmr + Sbar(k) = Sbar(k)*inyzmr + Srms(k) = sqrt(Srms(k)*inyzmr) + vxS(k) = vxS(k)*inyzmr + vyS(k) = vyS(k)*inyzmr + vzS(k) = vzS(k)*inyzmr + chiS(k) = chiS(k)/pecs*inyzmr + end do + + ! Store index as character string + write(nstat,"(i5.5)")nint(time/tout) + + if (ismaster) then + dsetname = trim("Sbar/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, Sbar, nxmr) + dsetname = trim("Srms/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, Srms, nxmr) + dsetname = trim("vxS/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, vxS, nxmr) + dsetname = trim("vyS/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, vyS, nxmr) + dsetname = trim("vzS/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, vzS, nxmr) + dsetname = trim("chiS/"//nstat) + call HdfSerialWriteReal1D(dsetname, filename, chiS, nxmr) + end if + + call MpiBarrier + +end subroutine CalcSalStats + +!> Calculate dissipation rate for salinity (on local process, no MPI action here) +subroutine CalcDissipationSal(chiS) + real, dimension(:), intent(out) :: chiS + + integer :: i, ip, im + integer :: j, jp, jm + integer :: k + + real, dimension(1:nxmr) :: tdxr + + do k=1,nxmr + tdxr(k) = 0.5*dxr/g3rmr(k) + end do + + ! If we are using the IBM, do not add contributions from the solid phase, + ! and enforce zero gradient at the boundary points + if (IBM) then + do i=xstartr(3),xendr(3) + ip = i + 1 + im = i - 1 + do j=xstartr(2),xendr(2) + jp = j + 1 + jm = j - 1 + do k=1,nxmr + if (.not. solidr(k,j,i)) then + if (solidr(k,j,ip)) then + chiS(k) = chiS(k) + ((sal(k,j,i ) - sal(k,j,im))*0.5*dzr)**2 + elseif (solidr(k,j,im)) then + chiS(k) = chiS(k) + ((sal(k,j,ip) - sal(k,j,i ))*0.5*dzr)**2 + else + chiS(k) = chiS(k) + ((sal(k,j,ip) - sal(k,j,im))*0.5*dzr)**2 + end if + if (solidr(k,jp,i)) then + chiS(k) = chiS(k) + ((sal(k,j ,i) - sal(k,jm,i))*0.5*dyr)**2 + elseif (solidr(k,jm,i)) then + chiS(k) = chiS(k) + ((sal(k,jp,i) - sal(k,j ,i))*0.5*dyr)**2 + else + chiS(k) = chiS(k) + ((sal(k,jp,i) - sal(k,jm,i))*0.5*dyr)**2 + end if + end if + end do + if (.not. solidr(1,j,i)) then + chiS(1) = chiS(1) + (( & + sal(2,j,i) - sal(1,j,i) + 2.0*SfixS*(sal(1,j,i)-salbp(1,j,i))& + )*tdxr(1))**2 + end if + do k=2,nxmr-1 + if (.not. solidr(k,j,i)) then + chiS(k) = chiS(k) + ((sal(k+1,j,i) - sal(k-1,j,i))*tdxr(k))**2 + end if + end do + if (.not. solidr(nxmr,j,i)) then + chiS(nxmr) = chiS(nxmr) + (( & + sal(nxmr,j,i) - sal(nxmr-1,j,i) + 2.0*SfixN*(saltp(1,j,i)-sal(nxmr,j,i)) & + )*tdxr(nxmr))**2 + end if + end do + end do + else + do i=xstartr(3),xendr(3) + ip = i + 1 + im = i - 1 + do j=xstartr(2),xendr(2) + jp = j + 1 + jm = j - 1 + do k=1,nxmr + chiS(k) = chiS(k) + ((sal(k,j,ip)-sal(k,j,im))*0.5*dzr)**2 + chiS(k) = chiS(k) + ((sal(k,jp,i)-sal(k,jm,i))*0.5*dyr)**2 + end do + chiS(1) = chiS(1) + (( & + sal(2,j,i) - sal(1,j,i) + 2.0*SfixS*(sal(1,j,i) - salbp(1,j,i)) & + )*tdxr(1))**2 + do k=2,nxmr-1 + chiS(k) = chiS(k) + ((sal(k+1,j,i) - sal(k-1,j,i))*tdxr(k))**2 + end do + chiS(nxmr) = chiS(nxmr) + (( & + sal(nxmr,j,i) - sal(nxmr-1,j,i) + 2.0*SfixN*(saltp(1,j,i)-sal(nxmr,j,i)) & + )*tdxr(nxmr))**2 + end do + end do + end if +end subroutine CalcDissipationSal + +!> Create the groups in the means.h5 file to store the +!! salinity-related statistics +subroutine CreateSalinityH5Groups(filename) + use HDF5 + + character(30), intent(in) :: filename + integer(HID_T) :: file_id, group_id + integer :: hdf_error + + call h5fopen_f(filename, H5F_ACC_RDWR_F, file_id, hdf_error) + + call h5gcreate_f(file_id, "Sbar", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "Srms", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "vxS", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "vyS", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "vzS", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + call h5gcreate_f(file_id, "chiS", group_id, hdf_error) + call h5gclose_f(group_id, hdf_error) + + call h5fclose_f(file_id, hdf_error) + +end subroutine CreateSalinityH5Groups + +end module afid_salinity \ No newline at end of file diff --git a/tools/afidtools/xdmf.py b/tools/afidtools/xdmf.py index d1b39180..b7ad966a 100644 --- a/tools/afidtools/xdmf.py +++ b/tools/afidtools/xdmf.py @@ -219,33 +219,33 @@ def generate_field_xmf(folder, var): with open(folder+"/outputdir/"+var+"_fields.xmf","w") as f: f.write(formatted_xmf.toprettyxml(indent=" ")) -def interpolate_field_to_uniform(folder, var): +def interpolate_field_to_uniform(folder, var, scale=2): # Create directory to store uniform-gridded snapshots os.makedirs(folder+"/outputdir/viz", exist_ok=True) # Create uniform x-grid to interpolate to grid = Grid(folder) inputs = InputParams(folder) nxm = grid.xm.size - nxu = nxm//2 + nxu = nxm//scale xu = linspace(0, inputs.alx3, nxu+1) xu = 0.5*(xu[:-1] + xu[1:]) # Pick corresponding grid for flow variable if var=="vx": xs = grid.xc - nyu, nzu = grid.ym.size//2, grid.zm.size//2 + nyu, nzu = grid.ym.size//scale, grid.zm.size//scale elif var=="sal" or var=="phi": xs = grid.xmr - nyu, nzu = grid.ymr.size//2, grid.zmr.size//2 + nyu, nzu = grid.ymr.size//scale, grid.zmr.size//scale else: xs = grid.xm - nyu, nzu = grid.ym.size//2, grid.zm.size//2 + nyu, nzu = grid.ym.size//scale, grid.zm.size//scale filelist = sorted(os.listdir(folder+"/outputdir/fields")) fvlist = list(filter(lambda fname: var in fname, filelist)) Funi = zeros((nzu, nyu, nxu), dtype=float32) for fname in fvlist: for k in range(nzu): with h5py.File(folder+"/outputdir/fields/"+fname, 'r') as f: - F = f['var'][2*k,::2,:] + F = f['var'][scale*k,::scale,:] if var=="vx": itp = interp1d(xs, F, kind='cubic', axis=-1) else: @@ -254,7 +254,7 @@ def interpolate_field_to_uniform(folder, var): with h5py.File(folder+"/outputdir/viz/"+fname, 'a') as f: f['var'] = Funi -def generate_uniform_xmf(folder, var): +def generate_uniform_xmf(folder, var, scale=2): """ Generates an xmf file in the Xdmf format to allow reading of the 3D fields in ParaView. This function produces a 3D array @@ -271,12 +271,12 @@ def generate_uniform_xmf(folder, var): nxmr, nymr, nzmr = grid.xmr.size, grid.ymr.size, grid.zmr.size # Store the appropriate grid sizes and names based on the variable - nxu = nxm//2 + nxu = nxm//scale fulldims = (nzm, nym, nxu) if var in "phisal": - nyu, nzu = nymr//2, nzmr//2 + nyu, nzu = nymr//scale, nzmr//scale else: - nyu, nzu = nym//2, nzm//2 + nyu, nzu = nym//scale, nzm//scale dx, dy, dz = grid.xc[-1]/nxu, grid.yc[-1]/nyu, grid.zc[-1]/nzu fulldims = (nzu, nyu, nxu) dims = fulldims diff --git a/tools/drizzle.py b/tools/drizzle.py new file mode 100644 index 00000000..d45398f3 --- /dev/null +++ b/tools/drizzle.py @@ -0,0 +1,64 @@ +import h5py +from scipy.special import lambertw +import os +import afidtools as afid +import numpy as np + +fld = os.environ['WORK']+'/RainyBenard/drizzle_IC' +inputs = afid.InputParams(fld) +nxm = inputs.nxm +print(nxm) + +saturated_top = False + +alpha = 3.0 +beta = 2.0 +if saturated_top: + gamma = beta/(1 - np.exp(-alpha)) +else: + gamma = beta +print(alpha, beta, gamma) + +def W(z): + return np.real(lambertw(z)) + +# Define vertical coordinate grid +zc = np.linspace(0, 1, nxm+1) +z = 0.5*(zc[1:] + zc[:-1]) + +# Set moist static energy +P = gamma +if saturated_top: + Q = beta - 1 - gamma*(1 - np.exp(-alpha)) +else: + Q = beta - 1 - gamma +m = P + Q*z + +# Define temperature profile +C = P + (Q - beta)*z +T = C - W(alpha*gamma*np.exp(alpha*C))/alpha + +# Infer buoyancy profile +b = T + beta*z + +# Compute humidity profile +q = (m - b)/gamma + +if not saturated_top: + # Calculate height above which humidity not at saturation + D = alpha*(beta + 1)*np.exp(1-alpha)/(1 + alpha*beta*np.exp(1-alpha)) + B = beta*D - 1 + zc = 1 + 1/alpha/(B-beta) + # Prescribe linear profiles above this level + b[z > zc] = beta - 1 + B*(z[z > zc] - 1) + q[z > zc] = D*(1 - z[z > zc]) + +import matplotlib.pyplot as plt +fig, ax = plt.subplots(layout='constrained') +ax.plot(b, z) +ax.plot(q, z) +fig.savefig('drizzle.png') + +with h5py.File(fld+'/drizzle.h5','w') as f: + f['b'] = b + f['q'] = q \ No newline at end of file