diff --git a/src/Makefile.openmpi b/src/Makefile.openmpi new file mode 100644 index 0000000..b65fce0 --- /dev/null +++ b/src/Makefile.openmpi @@ -0,0 +1,265 @@ +# FC: Specify Fortran90 compiler command. +FC=mpif90 + +# Added for displacement. +# PRJ4_DIR: Specify PROJ.4 install directory. +# CC&CFLAGS: Specify C compiler command and flags. +PROJ4_DIR=/home/baba/proj4 +CC=mpicc +CFLAGS=-g -I$(PROJ4_DIR)/include + +# BASE: Specify basic compiler options. +BASE=-cpp + +# FFTW3_INCLUDE_DIR: Specify FFTW3 include directory. +FFTW3_INCLUDE_DIR=/home/baba/fftw3/include + +# FFTW3_LIB: Specify linker options to link FFTW3. +FFTW3_LIB=-L/home/baba/fftw3/lib -lfftw3 + +# OPT: Specify compiler options about optimization level. +#OPT=-O2 -fopenmp -I$(FFTW3_INCLUDE_DIR) -march=native +OPT=-O2 -fopenmp -I$(FFTW3_INCLUDE_DIR) +#OPT=-O2 -I$(FFTW3_INCLUDE_DIR) + +# NETCDF: Specify the path to NetCDF library. +NETCDF=/home/baba/netcdf-fortran + +# LIBS: Specify linker options. +LIBS=-L$(NETCDF)/lib -lnetcdff -lnetcdf -L/home/baba/hdf5/lib -lhdf5_hl -lhdf5 -lcurl -lsz -L$(PROJ4_DIR)/lib -lproj $(FFTW3_LIB) # for NetCDF4! + +# EXEC: Specify the name of executable. +EXEC=jagurs + +################################################################################ +### Set parameters here! ### +################################################################################ +# PREC=REAL_DBLE: All calc. (except file I/O) is performed with double precision. +# DBLE_MATH: Only intrinsic math functions are performed with double precision. +# Else, all calc. is performed with single precision. +#PREC=DBLE_MATH +PREC=REAL_DBLE + +# MPI=ON: MPI version is made. +# Else, serial version is made. +#MPI=ON + +# USE_ALLTOALLV=ON: Use MPI_Alltoallv for inter-nest communications. +# Else, MPI_Allreduce is used. +#USE_ALLTOALLV=ON + +# A2A3D=ON: Use A2A3D for inter-nest communications. +# This is valid only if USE_ALLTOALLV=ON. +#A2A3D=ON + +# SINGLE_A2A=ON: Use SINGLE_A2A for inter-nest communications. +# This is valid only if USE_ALLTOALLV=ON. +#SINGLE_A2A=ON + +# TIMER=ON: Built-in timer is available. +# DETAIL: More detail. +# Else, built-in timer is disabled. +TIMER=ON +#TIMER=DETAIL + +# CONV_CHECK=ON: Convergence check is available on dispersive calc. +# Else, convergence check is disabled. +#CONV_CHECK=ON + +# OUTPUT=NCDIO: Output is put into single NetCDF file for each domain/MPI processes. +# DIROUT: GMT snapshot files are written into directories. +# Else, original GMT output. +#OUTPUT=DIROUT +#OUTPUT=NCDIO + +# MULTI=ON: Multi-members on a job is supported. +#MULTI=ON + +# SINGLE_TGS=ON: Station (tgs) output is written into single file per process. +# This file can be split into files for each station by "splittgs.sh". +#SINGLE_TGS=ON + +# CARTESIAN=ON: Cartesian axis is utilized. +#CARTESIAN=ON + +# UPWIND3=ON: Difference scheme for advective term is changed into 3rd-order upwind +# from original 1st-order upwind. +#UPWIND3=ON + +# SKIP_MAX_VEL=ON: Do not compute max velocity (saves time, most noticeably for +# linear computation) +#SKIP_MAX_VEL=ON + +# LESS_CC=ON: Convergence check on dispersive calc. is done only every 10 steps. +LESS_CC=ON + +# ONEFILE=ON: MPI parallelized version can read/write non-splitted GMT file. +ONEFILE=ON + +# BANKFILE=ON: Bank line data becomes available. +#BANKFILE=ON + +# HZMINOUT=ON: Min. height is output. +#HZMINOUT=ON + +#DUMP1D=ON + +################################################################################ +################################################################################ +### DON'T CHANGE THE FOLLOWINGS!!!! ### +################################################################################ +################################################################################ + +################################################################################ +### Parameters are processed here. ### +################################################################################ +FFLAGS=$(BASE) $(OPT) -I$(NETCDF)/include + +# Precision +ifeq ($(PREC),REAL_DBLE) + FFLAGS+=-DREAL_DBLE +else +ifeq ($(PREC),DBLE_MATH) + FFLAGS+=-DDBLE_MATH +endif +endif + +# MPI +ifeq ($(MPI),ON) + FFLAGS+=-DMPI -lmpi +ifeq ($(CARTESIAN),ON) + OBJS=mapproject.o okada.sub.o \ + mod_multi.o \ + mod_a2a3d.o mod_truncation.o mod_timer.o mod_mpi_fatal.o mod_grid.o mod_params.o mod_mpi.o mod_onefile.o mod_mygmt_gridio.o \ + mod_fxy.o mod_fxy_cartesian.o mod_nest.o mod_interpolation.o mod_fxy_disp.o mod_fxy_disp_cartesian.o mod_hxy.o mod_hxy_cartesian.o mod_bank.o \ + mod_displacement.o mod_rwg.o \ + mod_tgs.o mod_newsubs.o mod_ncdio.o mod_init_disp_gaussian.o mod_init_disp_sinwave.o mod_restart.o mod_dump1d.o JAGURS.o +else + OBJS=mapproject.o okada.sub.o \ + mod_multi.o \ + mod_a2a3d.o mod_truncation.o mod_timer.o mod_mpi_fatal.o mod_grid.o mod_params.o mod_mpi.o mod_onefile.o mod_mygmt_gridio.o \ + mod_fxy.o mod_fxy_cartesian.o mod_nest.o mod_interpolation.o mod_fxy_disp.o mod_fxy_disp_cartesian.o mod_hxy.o mod_hxy_cartesian.o mod_bank.o \ + mod_fxy_linear_Coriolis.o mod_fxy_linear_Coriolis_disp.o \ + mod_density.o mod_displacement.o mod_loading.o mod_rwg.o \ + mod_tgs.o mod_newsubs.o mod_ncdio.o mod_init_disp_gaussian.o mod_init_disp_sinwave.o mod_restart.o mod_dump1d.o JAGURS.o +endif +ifeq ($(USE_ALLTOALLV),ON) + FFLAGS+=-DUSE_ALLTOALLV +ifeq ($(A2A3D),ON) + FFLAGS+=-DA2A3D +endif +ifeq ($(SINGLE_A2A),ON) + FFLAGS+=-DSINGLE_A2A +endif +endif +else +ifeq ($(CARTESIAN),ON) + OBJS=mapproject.o okada.sub.o \ + mod_multi.o \ + mod_truncation.o mod_timer.o mod_grid.o mod_params.o mod_mygmt_gridio.o \ + mod_fxy.o mod_fxy_cartesian.o mod_nest.o mod_interpolation.o mod_fxy_disp.o mod_fxy_disp_cartesian.o mod_hxy.o mod_hxy_cartesian.o mod_bank.o \ + mod_displacement.o mod_rwg.o \ + mod_tgs.o mod_newsubs.o mod_ncdio.o mod_init_disp_gaussian.o mod_init_disp_sinwave.o mod_restart.o mod_dump1d.o JAGURS.o +else + OBJS=mapproject.o okada.sub.o \ + mod_multi.o \ + mod_truncation.o mod_timer.o mod_grid.o mod_params.o mod_mygmt_gridio.o \ + mod_fxy.o mod_fxy_cartesian.o mod_nest.o mod_interpolation.o mod_fxy_disp.o mod_fxy_disp_cartesian.o mod_hxy.o mod_hxy_cartesian.o mod_bank.o \ + mod_fxy_linear_Coriolis.o mod_fxy_linear_Coriolis_disp.o \ + mod_density.o mod_displacement.o mod_loading.o mod_rwg.o \ + mod_tgs.o mod_newsubs.o mod_ncdio.o mod_init_disp_gaussian.o mod_init_disp_sinwave.o mod_restart.o mod_dump1d.o JAGURS.o +endif +endif + +# Timer +ifeq ($(TIMER),DETAIL) + FFLAGS+=-DTIMER -DTIMER_DETAIL +else +ifeq ($(TIMER),ON) + FFLAGS+=-DTIMER +endif +endif + +# Convergence check +ifeq ($(CONV_CHECK),ON) + FFLAGS+=-DCONV_CHECK +endif + +# File output +ifeq ($(OUTPUT),NCDIO) + FFLAGS+=-DNCDIO +else +ifeq ($(OUTPUT),DIROUT) + FFLAGS+=-DDIROUT +endif +endif + +# Multi-members +ifeq ($(MULTI),ON) +ifeq ($(MPI),ON) + FFLAGS+=-DMULTI +else + FFLAGS+=-DMULTI -lmpi +endif +endif + +# Single tgs output +ifeq ($(SINGLE_TGS),ON) + FFLAGS+=-DSINGLE_TGS +endif + +# Cartesian axis +ifeq ($(CARTESIAN),ON) + FFLAGS+=-DCARTESIAN +endif + +# 3rd-order upwind +ifeq ($(UPWIND3),ON) + FFLAGS+=-DUPWIND3 +endif + +# Skip max velocity computation +ifeq ($(SKIP_MAX_VEL),ON) + FFLAGS+=-DSKIP_MAX_VEL +endif + +# Convergence check on dispersive calc. is done only every 10 steps. +ifeq ($(LESS_CC),ON) + FFLAGS+=-DLESS_CC +endif + +ifeq ($(ONEFILE),ON) + FFLAGS+=-DONEFILE +endif + +# Bank line data becomes available. +ifeq ($(BANKFILE),ON) + FFLAGS+=-DBANKFILE +endif + +# Min. height is output. +ifeq ($(HZMINOUT),ON) + FFLAGS+=-DHZMINOUT +endif + +ifeq ($(DUMP1D),ON) + FFLAGS+=-DDUMP1D +endif + +################################################################################ +### Make rules ### +################################################################################ +$(EXEC): $(OBJS) + $(FC) -o $@ $(FFLAGS) $^ $(LIBS) + +%.o: %.f90 + $(FC) -c $(FFLAGS) $^ + +%.o: %.f + $(FC) -c $(FFLAGS) $^ + +%.o: %.c + $(CC) -c $(CFLAGS) $^ + +clean: + rm -f *.o *.mod *__genmod.f90 $(EXEC) diff --git a/src/mod_bank.f90 b/src/mod_bank.f90 index 96daad2..65dec93 100644 --- a/src/mod_bank.f90 +++ b/src/mod_bank.f90 @@ -1,4 +1,5 @@ !#define OLDFORMAT +#define BANKFILEREAL #ifdef DBLE_MATH #include "dble_math.h" #endif @@ -44,7 +45,13 @@ subroutine read_bank_file(dg) integer(kind=4), allocatable, dimension(:) :: x, y, val #ifdef CARTESIAN +#ifndef BANKFILEREAL integer(kind=4) :: n, i, j, num_lines, xtmp, ytmp, valtmp +#else + integer(kind=4) :: n, i, j, num_lines, valtmp + real(kind=REAL_BYTE) :: xtmp, ytmp + real(kind=REAL_BYTE), allocatable, dimension(:) :: xin, yin +#endif #else integer(kind=4) :: n, i, j, num_lines, valtmp #endif @@ -53,7 +60,13 @@ subroutine read_bank_file(dg) real(kind=REAL_BYTE), allocatable, dimension(:) :: height #ifdef CARTESIAN +#ifndef BANKFILEREAL integer(kind=4) :: n, i, j, num_lines, xtmp, ytmp, irtmp +#else + integer(kind=4) :: n, i, j, num_lines, irtmp + real(kind=REAL_BYTE) :: xtmp, ytmp + real(kind=REAL_BYTE), allocatable, dimension(:) :: xin, yin +#endif #else integer(kind=4) :: n, i, j, num_lines, irtmp #endif @@ -91,6 +104,10 @@ subroutine read_bank_file(dg) #endif allocate(x(num_lines)) allocate(y(num_lines)) +#ifdef BANKFILEREAL + allocate(xin(num_lines)) + allocate(yin(num_lines)) +#endif #ifdef OLDFORMAT allocate(val(num_lines)) #else @@ -101,6 +118,7 @@ subroutine read_bank_file(dg) rewind(1) do n = 1, num_lines #ifdef CARTESIAN +#ifndef BANKFILEREAL #ifdef OLDFORMAT read(1,*) x(n), y(n), val(n) #else @@ -112,6 +130,13 @@ subroutine read_bank_file(dg) #else read(1,*) yin(n), xin(n), irread(n), height(n) #endif +#endif +#else +#ifdef OLDFORMAT + read(1,*) xin(n), yin(n), val(n) +#else + read(1,*) yin(n), xin(n), irread(n), height(n) +#endif #endif end do @@ -164,8 +189,13 @@ subroutine read_bank_file(dg) ! x(n) = int((x(n) - mlon0 + 0.5d0)/dh) + 1 ! y(n) = int((y(n) - mlat0 + 0.5d0)/dh) + 1 #ifdef CARTESIAN +#ifndef BANKFILEREAL x(n) = floor((x(n) - mlon0 + 0.5d0)/dh) + 1 y(n) = floor((y(n) - mlat0 + 0.5d0)/dh) + 1 +#else + x(n) = floor((xin(n) - mlon0 + 0.5d0*dh)/dh) + 1 + y(n) = floor((yin(n) - mlat0 + 0.5d0*dh)/dh) + 1 +#endif y(n) = ny - y(n) + 1 #else x(n) = floor((xin(n) - mlon0)/dh) + 1 @@ -228,10 +258,10 @@ subroutine read_bank_file(dg) do i = -1, nx+1 if(ir(i,j) == 0) then if(abs(btx(i,j)) > tiny(bh)) then - write(0,*) 'Somthing is wrong on btx!!!', btx(i,j) + write(0,*) 'Something is wrong on btx!!!', btx(i,j) end if if(abs(bty(i,j-1)) > tiny(bh)) then - write(0,*) 'Somthing is wrong on bty!!!', bty(i,j-1) + write(0,*) 'Something is wrong on bty!!!', bty(i,j-1) end if else write(400000,'(3i6,2f15.3)') i, j, ir(i,j), btx(i,j), bty(i,j-1) @@ -246,6 +276,10 @@ subroutine read_bank_file(dg) #endif deallocate(x) deallocate(y) +#ifdef BANKFILEREAL + deallocate(xin) + deallocate(yin) +#endif #ifdef OLDFORMAT deallocate(val) #else diff --git a/src/mod_displacement.f90 b/src/mod_displacement.f90 index b9235d8..8a44d30 100644 --- a/src/mod_displacement.f90 +++ b/src/mod_displacement.f90 @@ -572,6 +572,7 @@ subroutine displacement_calc_displacement(dg, ig) u1, u2, u3) #endif uh1(i,j) = u1*sin(strike(n)*DEG2RAD) - u2*cos(strike(n)*DEG2RAD) + ! note that downward is positive along Y direction in JAGURS uh2(i,j) = -u1*cos(strike(n)*DEG2RAD) - u2*sin(strike(n)*DEG2RAD) uz(i,j) = u3 #ifndef CARTESIAN @@ -601,7 +602,9 @@ subroutine displacement_calc_displacement(dg, ig) do j = 1, nlat do i = 1, nlon uh1(i,j) = uh1(i,j)*0.01d0 - uh2(i,j) = -uh2(i,j)*0.01d0 + !a bug fix on 14/04/2019, Baba + !uh2(i,j) = -uh2(i,j)*0.01d0 + uh2(i,j) = uh2(i,j)*0.01d0 uz(i,j) = uz(i,j) *0.01d0 end do end do diff --git a/src/mod_fxy.f90 b/src/mod_fxy.f90 index 0910674..9e9a7eb 100644 --- a/src/mod_fxy.f90 +++ b/src/mod_fxy.f90 @@ -815,18 +815,17 @@ subroutine fxynl_rwg(wfld,dfld,ffld,ifz,cfs,cfl,cflag,dt,th0,dth,joff,nlon,nlat, else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddx_tmp**(-1.0d0/3.0d0) end if + ! explicit ! fric = dt*bcf*fx_old(i,j)*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 - fric = dt*bcf*(half*(fx(i,j)+fx_old(i,j)))*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fx_old(i,j)*fx_old(i,j)+fybar*fybar)/ddx_tmp/ddx_tmp + fx(i,j) = (fx(i,j)-fric*fx_old(i,j))/(1.0d0+fric) else - fric = zap + fx(i,j) = zap +! fric = zap end if - fx(i,j) = fx(i,j) - fric +! fx(i,j) = fx(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fx(i,j) > 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! else if(fx(i,j) < -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) d = max(0.0d0, half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) @@ -863,18 +862,17 @@ subroutine fxynl_rwg(wfld,dfld,ffld,ifz,cfs,cfl,cflag,dt,th0,dth,joff,nlon,nlat, else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddy_tmp**(-1.0d0/3.0d0) end if + ! explicit ! fric = dt*bcf*fy_old(i,j)*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 - fric = dt*bcf*(half*(fy(i,j)+fy_old(i,j)))*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fy_old(i,j)*fy_old(i,j)+fxbar*fxbar)/ddy_tmp/ddy_tmp + fy(i,j) = (fy(i,j)-fric*fy_old(i,j))/(1.0d0+fric) else - fric = zap + fy(i,j) = zap +! fric = zap end if - fy(i,j) = fy(i,j) - fric +! fy(i,j) = fy(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fy(i,j) > 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! else if(fy(i,j) < -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) d = max(0.0d0, half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) diff --git a/src/mod_fxy_cartesian.f90 b/src/mod_fxy_cartesian.f90 index 6347c2b..a9e6b5e 100644 --- a/src/mod_fxy_cartesian.f90 +++ b/src/mod_fxy_cartesian.f90 @@ -749,18 +749,17 @@ subroutine fxynl_rwg(wfld,dfld,ffld,ifz,cfs,cfl,dt,dxdy,nlon,nlat,gflag,smallh,b else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddx_tmp**(-1.0d0/3.0d0) end if + ! explicit ! fric = dt*bcf*fx_old(i,j)*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 - fric = dt*bcf*(half*(fx(i,j)+fx_old(i,j)))*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fx_old(i,j)*fx_old(i,j)+fybar*fybar)/ddx_tmp/ddx_tmp + fx(i,j) = (fx(i,j)-fric*fx_old(i,j))/(1.0d0+fric) else - fric = zap + fx(i,j) = zap +! fric = zap end if - fx(i,j) = fx(i,j) - fric +! fx(i,j) = fx(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fx(i,j) > 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! else if(fx(i,j) < -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) d = max(0.0d0, half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) @@ -797,18 +796,17 @@ subroutine fxynl_rwg(wfld,dfld,ffld,ifz,cfs,cfl,dt,dxdy,nlon,nlat,gflag,smallh,b else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddy_tmp**(-1.0d0/3.0d0) end if + ! explicit ! fric = dt*bcf*fy_old(i,j)*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 - fric = dt*bcf*(half*(fy(i,j)+fy_old(i,j)))*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fy_old(i,j)*fy_old(i,j)+fxbar*fxbar)/ddy_tmp/ddy_tmp + fy(i,j) = (fy(i,j)-fric*fy_old(i,j))/(1.0d0+fric) else - fric = zap + fy(i,j) = zap +! fric = zap end if - fy(i,j) = fy(i,j) - fric +! fy(i,j) = fy(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fy(i,j) > 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! else if(fy(i,j) < -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) d = max(0.0d0, half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) diff --git a/src/mod_fxy_disp.f90 b/src/mod_fxy_disp.f90 index feb4798..4a569f6 100644 --- a/src/mod_fxy_disp.f90 +++ b/src/mod_fxy_disp.f90 @@ -1272,18 +1272,17 @@ subroutine fxynl_rwg_disp(wfld,dfld,ffld,ifz,cfs,cfl,cflag,dt,th0,dth,joff,nlon, else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddx_tmp**(-1.0d0/3.0d0) end if + ! explicit ! fric = dt*bcf*fx_old(i,j)*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 - fric = dt*bcf*(half*(fx(i,j)+fx_old(i,j)))*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fx_old(i,j)*fx_old(i,j)+fybar*fybar)/ddx_tmp/ddx_tmp + fx(i,j) = (fx(i,j)-fric*fx_old(i,j))/(1.0d0+fric) else - fric = zap +! fric = zap + fx(i,j) = zap end if - fx(i,j) = fx(i,j) - fric +! fx(i,j) = fx(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fx(i,j) > 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! else if(fx(i,j) < -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) d = max(0.0d0, half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) @@ -1384,18 +1383,17 @@ subroutine fxynl_rwg_disp(wfld,dfld,ffld,ifz,cfs,cfl,cflag,dt,th0,dth,joff,nlon, else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddy_tmp**(-1.0d0/3.0d0) end if + ! explicit ! fric = dt*bcf*fy_old(i,j)*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 - fric = dt*bcf*(half*(fy(i,j)+fy_old(i,j)))*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fy_old(i,j)*fy_old(i,j)+fxbar*fxbar)/ddy_tmp/ddy_tmp + fy(i,j) = (fy(i,j)-fric*fy_old(i,j))/(1.0d0+fric) else - fric = zap + fy(i,j) = zap +! fric = zap end if - fy(i,j) = fy(i,j) - fric +! fy(i,j) = fy(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fy(i,j) > 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! else if(fy(i,j) < -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) d = max(0.0d0, half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) diff --git a/src/mod_fxy_disp_cartesian.f90 b/src/mod_fxy_disp_cartesian.f90 index ce90ac2..cdbc052 100644 --- a/src/mod_fxy_disp_cartesian.f90 +++ b/src/mod_fxy_disp_cartesian.f90 @@ -1170,18 +1170,17 @@ subroutine fxynl_rwg_disp(wfld,dfld,ffld,ifz,cfs,cfl,dt,dxdy,nlon,nlat, & else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddx_tmp**(-1.0d0/3.0d0) end if + ! explict ! fric = dt*bcf*fx_old(i,j)*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 - fric = dt*bcf*(half*(fx(i,j)+fx_old(i,j)))*sqrt(fx_old(i,j)*fx_old(i,j) + fybar*fybar)/ddx_tmp**2 + ! semi-implicit added on 14/04/2019, Baba + fric = dt*bcf*half*sqrt(fx_old(i,j)*fx_old(i,j)+fybar*fybar)/ddx_tmp/ddx_tmp + fx(i,j) = (fx(i,j)-fric*fx_old(i,j))/(1.0d0+fric) else - fric = zap + fx(i,j) = zap +! fric = zap end if - fx(i,j) = fx(i,j) - fric +! fx(i,j) = fx(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fx(i,j) > 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = 15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! else if(fx(i,j) < -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) then -! fx(i,j) = -15.0d0*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j)) d = max(0.0d0, half*(dz(i+1,j)+dz(i,j)+hz_old(i+1,j)+hz_old(i,j))) @@ -1262,18 +1261,16 @@ subroutine fxynl_rwg_disp(wfld,dfld,ffld,ifz,cfs,cfl,dt,dxdy,nlon,nlat, & else ! Manning's roughness coefficent bcf = cf*cf*9.8d0*ddy_tmp**(-1.0d0/3.0d0) end if + ! explict ! fric = dt*bcf*fy_old(i,j)*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 - fric = dt*bcf*(half*(fy(i,j)+fy_old(i,j)))*sqrt(fy_old(i,j)*fy_old(i,j) + fxbar*fxbar)/ddy_tmp**2 + fric = dt*bcf*half*sqrt(fy_old(i,j)*fy_old(i,j)+fxbar*fxbar)/ddy_tmp/ddy_tmp + fy(i,j) = (fy(i,j)-fric*fy_old(i,j))/(1.0d0+fric) else - fric = zap + fy(i,j) = zap +! fric = zap end if - fy(i,j) = fy(i,j) - fric +! fy(i,j) = fy(i,j) - fric ! === Limiter with max Froude number. ========================================== -! if(fy(i,j) > 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = 15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! else if(fy(i,j) < -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j))) then -! fy(i,j) = -15.0d0*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) -! end if ! === To prevent sqrt of negative numbers. ===================================== ! d = half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)) d = max(0.0d0, half*(dz(i,j+1)+dz(i,j)+hz_old(i,j+1)+hz_old(i,j)))