diff --git a/.gitignore b/.gitignore index 22752d026..4ccdc3b3d 100644 --- a/.gitignore +++ b/.gitignore @@ -65,3 +65,10 @@ source/third_party_open/netcdf/x64 source/sfincs/sfincs.opt.yaml /source/sfincs_lib/*.yaml /source/third_party_open/netcdf/netcdf-fortran-4.6.1/Debug +/source/build_nvfortran_gpu_h7.sh +/source/build_nvfortran_gpu.sh +/source/Singularityfile-gpu.def +/source/Dockerfile.gpu.update01 +/source/Dockerfile.gpu.test +/source/Dockerfile.gpu.25.5.ccall +/source/Dockerfile.gpu diff --git a/source/Dockerfile.gpu b/source/Dockerfile.gpu deleted file mode 100644 index 24da778bc..000000000 --- a/source/Dockerfile.gpu +++ /dev/null @@ -1,31 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:20.11-devel-cuda_multi-ubuntu20.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y dos2unix -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/20.11/REDIST/compilers && tar -czvhf /root/nvidia-20.11.tar.gz * && cd - -COPY . /usr/src/sfincs -WORKDIR /usr/src/sfincs - -RUN find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -RUN find . -name \*.f90|xargs dos2unix -RUN find . -name \*.F90|xargs dos2unix -RUN find . -name \*.am|xargs dos2unix -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=8" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -FROM nvidia/cuda:11.2.2-runtime-ubuntu20.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /root/nvidia-20.11.tar.gz . -# copy software -COPY --from=0 /usr/local /usr/local -RUN mkdir -p /opt/nvidia/compilers -RUN tar -xzf nvidia-20.11.tar.gz -C /opt/nvidia/compilers -# add the REDIST libs to the environment -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] - diff --git a/source/Dockerfile.gpu.25.5.ccall b/source/Dockerfile.gpu.25.5.ccall deleted file mode 100644 index fb8e25234..000000000 --- a/source/Dockerfile.gpu.25.5.ccall +++ /dev/null @@ -1,32 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:25.5-devel-cuda_multi-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y dos2unix -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/2025/REDIST/compilers && tar -czvhf /root/nvidia-25.5.tar.gz * && cd - -COPY . /usr/src/sfincs -WORKDIR /usr/src/sfincs - -RUN find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -RUN find . -name \*.f90|xargs dos2unix -RUN find . -name \*.F90|xargs dos2unix -RUN find . -name \*.am|xargs dos2unix -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=8" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install -#RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -FROM nvidia/cuda:12.5.0-runtime-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /root/nvidia-25.5.tar.gz . -# copy software -COPY --from=0 /usr/local /usr/local -RUN mkdir -p /opt/nvidia/compilers -RUN tar -xzf nvidia-25.5.tar.gz -C /opt/nvidia/compilers -# add the REDIST libs to the environment -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] - diff --git a/source/Dockerfile.gpu.test b/source/Dockerfile.gpu.test deleted file mode 100644 index 87870c2ac..000000000 --- a/source/Dockerfile.gpu.test +++ /dev/null @@ -1,15 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:20.11-devel-cuda_multi-ubuntu20.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/20.11/REDIST/compilers && tar -czvhf /root/nvidia-20.11.tar.gz * && cd - -WORKDIR /usr/src/sfincs -COPY . /usr/src/sfincs -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -# create a debug version -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -g -O0 -gpu=ccall" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] diff --git a/source/Dockerfile.gpu.update01 b/source/Dockerfile.gpu.update01 deleted file mode 100644 index bd54b6412..000000000 --- a/source/Dockerfile.gpu.update01 +++ /dev/null @@ -1,32 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:24.3-devel-cuda_multi-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y dos2unix -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/2024/REDIST/compilers && tar -czvhf /root/nvidia-24.5.tar.gz * && cd - -COPY . /usr/src/sfincs -WORKDIR /usr/src/sfincs - -RUN find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -RUN find . -name \*.f90|xargs dos2unix -RUN find . -name \*.F90|xargs dos2unix -RUN find . -name \*.am|xargs dos2unix -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=8" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install -#RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -FROM nvidia/cuda:12.5.0-runtime-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /root/nvidia-24.5.tar.gz . -# copy software -COPY --from=0 /usr/local /usr/local -RUN mkdir -p /opt/nvidia/compilers -RUN tar -xzf nvidia-24.5.tar.gz -C /opt/nvidia/compilers -# add the REDIST libs to the environment -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] - diff --git a/source/Singularityfile-gpu.def b/source/Singularityfile-gpu.def deleted file mode 100644 index 69ab0a258..000000000 --- a/source/Singularityfile-gpu.def +++ /dev/null @@ -1,58 +0,0 @@ -BootStrap: library -From: library://library/default/ubuntu:jammy -Stage: build - -# -# The source files of sfincs -# - -%files - . /usr/src/sfincs - -# -# Compiler flags -# -fallow-argument-mismatch needed for https://github.com/Unidata/netcdf-fortran/issues/212 -# - -%environment - FCFLAGS="-fopenmp -O3 -fallow-argument-mismatch -w" - FFLAGS="-fopenmp -O3 -fallow-argument-mismatch -w" - -# -# Installing compilers -# Compiling sfincs -# - -%post - apt clean && apt autoclean && apt -y update --fix-missing && apt -y upgrade -y && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config gfortran gdb m4 - -# Install software into /usr/local by default - - cd /usr/src/sfincs - export CONFIG_SHELL=/bin/bash - autoreconf -vif - ./autogen.sh - ./configure FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=999" FC=nvfortran --disable-shared --disable-openmp --program-suffix="_async" - - make clean - make - make install - -# Install binary in a much smaller image - -BootStrap: library -From: library://library/default/ubuntu:jammy -Stage: final - -# install binary from stage one -%files from build - /usr/local/bin/sfincs /usr/local/bin/sfincs - -%files - run-sfincs.sh /usr/local/bin/run-sfincs.sh - -%post - apt clean && apt autoclean && apt update --fix-missing && apt upgrade -y - apt install -y libnetcdf19 libgfortran5 libgomp1 - apt clean && apt autoclean - chmod +x /usr/local/bin/run-sfincs.sh diff --git a/source/build_nvfortran_gpu.sh b/source/build_nvfortran_gpu.sh deleted file mode 100644 index a38ae8a7e..000000000 --- a/source/build_nvfortran_gpu.sh +++ /dev/null @@ -1,35 +0,0 @@ -#!/bin/bash - -LOGFILE=build_$(date +%Y%m%d_%H%M%S).log -exec > >(tee "$LOGFILE") 2>&1 - -echo "Starting build at $(date)" -echo "Saving log to $LOGFILE" - -find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -find . -name \*.f90|xargs dos2unix -find . -name \*.F90|xargs dos2unix -find . -name \*.am|xargs dos2unix -find . -name \*.sh|xargs dos2unix - -MANPATH=$MANPATH:/opt/nvidia/hpc_sdk/Linux_x86_64/24.5/compilers/man; export MANPATH -PATH=/opt/nvidia/hpc_sdk/Linux_x86_64/24.5/compilers/bin:$PATH; export PATH - -LD_LIBRARY_PATH=/usr/lib/wsl/lib:$LD_LIBRARY_PATH; export LD_LIBRARY_PATH - -apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -export CONFIG_SHELL=/bin/bash - -autoreconf -vif - -./autogen.sh - -./configure FCFLAGS="-acc=gpu -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=999" FC=nvfortran --disable-shared --disable-openmp --prefix=/usr/local/bin/sfincs/nvfortran_gpu_ccall - -make clean - -make - -make install - -echo "Build finished at $(date)" diff --git a/source/build_nvfortran_gpu_h7.sh b/source/build_nvfortran_gpu_h7.sh deleted file mode 100644 index c50fb36a6..000000000 --- a/source/build_nvfortran_gpu_h7.sh +++ /dev/null @@ -1,21 +0,0 @@ -#!/bin/sh - -module load nvidia/nvhpc/24.1 -module load netcdf - -export LDFLAGS="-L${NETCDF_C_LIBRARY}" - -find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -find . -name \*.f90|xargs dos2unix -find . -name \*.F90|xargs dos2unix -find . -name \*.am|xargs dos2unix - -./autogen.sh - -./configure FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=999" FC=nvfortran --disable-shared --disable-openmp --prefix /u/${USER}/bin/sfincs_nvfortran_gpu - -make clean - -make - -make install diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 43a3918b9..ed9b4d918 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -589,7 +589,10 @@ module sfincs_data real*4, dimension(:), allocatable :: df real*4, dimension(:), allocatable :: dwig real*4, dimension(:), allocatable :: dfig - real*4, dimension(:), allocatable :: cg + real*4, dimension(:), allocatable :: cg + real*4, dimension(:), allocatable :: cgig + real*4, dimension(:), allocatable :: qb + real*4, dimension(:), allocatable :: gam real*4, dimension(:), allocatable :: betamean real*4, dimension(:), allocatable :: srcig real*4, dimension(:), allocatable :: alphaig diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index f4b8fe1b1..724b028e6 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2581,6 +2581,12 @@ subroutine initialize_hydro() dfig = 0.0 allocate(cg(np)) cg = 0.0 + allocate(cgig(np)) + cgig = 0.0 + allocate(qb(np)) + qb = 0.0 + allocate(gam(np)) + gam = 0.0 allocate(betamean(np)) betamean = 0.0 allocate(srcig(np)) diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index d1881be4d..8f3ffd4bf 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -92,8 +92,8 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.1 mt. Faber" - build_date = "$Date: 2025-12-18" + build_revision = "$Rev: v2.3.0+ mt. Faber:branch-281" + build_date = "$Date: 2026-04-01" ! call write_log('', 1) call write_log('------------ Welcome to SFINCS ------------', 1) diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 3de049327..e82491b99 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -54,7 +54,7 @@ module sfincs_ncoutput integer :: patm_varid, wind_speed_varid, wind_dir_varid integer :: inp_varid, total_runtime_varid, average_dt_varid, status_varid integer :: hm0_varid, hm0ig_varid, zsm_varid, tp_varid, tpig_varid, wavdir_varid, dirspr_varid - integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, beta_varid, srcig_varid, alphaig_varid + integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, cgig_varid, beta_varid, srcig_varid, alphaig_varid, qb_varid, gam_varid integer :: runup_gauge_name_varid, runup_gauge_zs_varid ! end type @@ -1849,6 +1849,27 @@ subroutine ncoutput_his_init() NF90(nf90_put_att(his_file%ncid, his_file%cg_varid, 'long_name', 'wave group velocity')) NF90(nf90_put_att(his_file%ncid, his_file%cg_varid, 'coordinates', 'station_id station_name point_x point_y')) ! + NF90(nf90_def_var(his_file%ncid, 'cgig', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%cgig_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'units', 'm/s')) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'standard_name', 'infragravity_wave_velocity')) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'long_name', 'infragravity wave velocity')) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'coordinates', 'station_id station_name point_x point_y')) + ! + NF90(nf90_def_var(his_file%ncid, 'qb', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%qb_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'units', '-')) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'standard_name', 'fraction_breaking_waves')) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'long_name', 'fraction breaking incident waves')) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'coordinates', 'station_id station_name point_x point_y')) + ! + NF90(nf90_def_var(his_file%ncid, 'gam', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%gam_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'units', '-')) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'standard_name', 'local_wave_height_water_depth_ratio')) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'long_name', 'local wave height water depth ratio')) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'coordinates', 'station_id station_name point_x point_y')) + ! NF90(nf90_def_var(his_file%ncid, 'beta', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%beta_varid)) ! time-varying water level point NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, 'units', '-')) @@ -2915,6 +2936,9 @@ subroutine ncoutput_update_his(t,nthisout) real*4, dimension(nobs) :: dwigobs real*4, dimension(nobs) :: dfigobs real*4, dimension(nobs) :: cgobs + real*4, dimension(nobs) :: cgigobs + real*4, dimension(nobs) :: qbobs + real*4, dimension(nobs) :: gamobs real*4, dimension(nobs) :: betaobs real*4, dimension(nobs) :: srcigobs real*4, dimension(nobs) :: alphaigobs @@ -2938,6 +2962,9 @@ subroutine ncoutput_update_his(t,nthisout) dwobs = FILL_VALUE dfobs = FILL_VALUE cgobs = FILL_VALUE + cgigobs = FILL_VALUE + qbobs = FILL_VALUE + gamobs = FILL_VALUE betaobs = FILL_VALUE srcigobs = FILL_VALUE alphaigobs = FILL_VALUE @@ -3038,7 +3065,10 @@ subroutine ncoutput_update_his(t,nthisout) dfobs(iobs) = df(nm) dwigobs(iobs) = dwig(nm) dfigobs(iobs) = dfig(nm) - cgobs(iobs) = cg(nm) + cgobs(iobs) = cg(nm) + cgigobs(iobs) = cgig(nm) + qbobs(iobs) = qb(nm) + gamobs(iobs) = gam(nm) betaobs(iobs) = betamean(nm) srcigobs(iobs) = srcig(nm) alphaigobs(iobs) = alphaig(nm) @@ -3101,6 +3131,9 @@ subroutine ncoutput_update_his(t,nthisout) NF90(nf90_put_var(his_file%ncid, his_file%dfig_varid, dfigobs, (/1, nthisout/))) ! NF90(nf90_put_var(his_file%ncid, his_file%cg_varid, cgobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%cgig_varid, cgigobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%qb_varid, qbobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%gam_varid, gamobs, (/1, nthisout/))) ! NF90(nf90_put_var(his_file%ncid, his_file%beta_varid, betaobs, (/1, nthisout/))) NF90(nf90_put_var(his_file%ncid, his_file%srcig_varid, srcigobs, (/1, nthisout/))) @@ -3690,7 +3723,7 @@ subroutine ncoutput_add_params(ncid, varid) ! ! Because of overlapping names, only important specific values from snapwave_data use snapwave_data, only: gamma, gammax, alpha, hmin, fw0, fw0_ig, dt, tol, dtheta, crit, nr_sweeps, baldock_opt, baldock_ratio, & - igwaves_opt, alpha_ig, gamma_ig, shinc2ig, alphaigfac, baldock_ratio_ig, ig_opt, herbers_opt, tpig_opt, eeinc2ig, tinc2ig, & + igwaves_opt, alpha_ig, gamma_ig, gamma_fac_br, shinc2ig, alphaigfac, baldock_ratio_ig, ig_opt, herbers_opt, tpig_opt, eeinc2ig, tinc2ig, & snapwave_jonswapfile, snapwave_encfile, snapwave_bndfile, snapwave_bhsfile, snapwave_btpfile, snapwave_bwdfile, snapwave_bdsfile, upwfile, gridfile ! @@ -3883,7 +3916,8 @@ subroutine ncoutput_add_params(ncid, varid) ! NF90(nf90_put_att(ncid, varid, 'snapwave_igwaves',igwaves_opt)) NF90(nf90_put_att(ncid, varid, 'snapwave_alpha_ig',alpha_ig)) - NF90(nf90_put_att(ncid, varid, 'snapwave_gammaig',gamma_ig)) + NF90(nf90_put_att(ncid, varid, 'snapwave_gammaig',gamma_ig)) + NF90(nf90_put_att(ncid, varid, 'snapwave_gamma_fac_br',gamma_fac_br)) NF90(nf90_put_att(ncid, varid, 'snapwave_shinc2ig',shinc2ig)) NF90(nf90_put_att(ncid, varid, 'snapwave_alphaigfac',alphaigfac)) NF90(nf90_put_att(ncid, varid, 'snapwave_baldock_ratio_ig',baldock_ratio_ig)) diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 05281e4d2..2a9a3edae 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -26,6 +26,9 @@ module sfincs_snapwave real*4, dimension(:), allocatable :: snapwave_Dwig real*4, dimension(:), allocatable :: snapwave_Dfig real*4, dimension(:), allocatable :: snapwave_cg + real*4, dimension(:), allocatable :: snapwave_cgig + real*4, dimension(:), allocatable :: snapwave_qb + real*4, dimension(:), allocatable :: snapwave_gam real*4, dimension(:), allocatable :: snapwave_beta real*4, dimension(:), allocatable :: snapwave_srcig real*4, dimension(:), allocatable :: snapwave_alphaig @@ -297,7 +300,10 @@ subroutine update_wave_field(t, tloop) real*4, dimension(:), allocatable :: df0 real*4, dimension(:), allocatable :: dwig0 real*4, dimension(:), allocatable :: dfig0 - real*4, dimension(:), allocatable :: cg0 + real*4, dimension(:), allocatable :: cg0 + real*4, dimension(:), allocatable :: cgig0 + real*4, dimension(:), allocatable :: qb0 + real*4, dimension(:), allocatable :: gam0 real*4, dimension(:), allocatable :: beta0 real*4, dimension(:), allocatable :: srcig0 real*4, dimension(:), allocatable :: alphaig0 @@ -312,7 +318,10 @@ subroutine update_wave_field(t, tloop) allocate(df0(np)) allocate(dwig0(np)) allocate(dfig0(np)) - allocate(cg0(np)) + allocate(cg0(np)) + allocate(cgig0(np)) + allocate(qb0(np)) + allocate(gam0(np)) allocate(beta0(np)) allocate(srcig0(np)) allocate(alphaig0(np)) @@ -324,6 +333,9 @@ subroutine update_wave_field(t, tloop) dwig0 = 0.0 dfig0 = 0.0 cg0 = 0.0 + cgig0 = 0.0 + qb0 = 0.0 + gam0 = 0.0 beta0 = 0.0 srcig0 = 0.0 alphaig0 = 0.0 @@ -418,6 +430,9 @@ subroutine update_wave_field(t, tloop) dwig0(nm) = snapwave_Dwig(ip) dfig0(nm) = snapwave_Dfig(ip) cg0(nm) = snapwave_cg(ip) + cgig0(nm) = snapwave_cgig(ip) + qb0(nm) = snapwave_qb(ip) + gam0(nm) = snapwave_gam(ip) beta0(nm) = snapwave_beta(ip) srcig0(nm) = snapwave_srcig(ip) alphaig0(nm) = snapwave_alphaig(ip) @@ -441,6 +456,9 @@ subroutine update_wave_field(t, tloop) dwig0(nm) = 0.0 dfig0(nm) = 0.0 cg0(nm) = 0.0 + cgig0(nm) = 0.0 + qb0(nm) = 0.0 + gam0(nm) = 0.0 beta0(nm) = 0.0 srcig0(nm) = 0.0 alphaig0(nm) = 0.0 @@ -459,7 +477,10 @@ subroutine update_wave_field(t, tloop) df(nm) = df0(nm) dwig(nm) = dwig0(nm) dfig(nm) = dfig0(nm) - cg(nm) = cg0(nm) + cg(nm) = cg0(nm) + cgig(nm) = cgig0(nm) + qb(nm) = qb0(nm) + gam(nm) = gam0(nm) betamean(nm) = beta0(nm) srcig(nm) = srcig0(nm) alphaig(nm) = alphaig0(nm) @@ -536,6 +557,9 @@ subroutine compute_snapwave(t) snapwave_Dwig = Dw_ig snapwave_Dfig = Df_ig snapwave_cg = cg + snapwave_cgig = cg_ig + snapwave_qb = qb + snapwave_gam = gam snapwave_beta = beta snapwave_srcig = srcig snapwave_alphaig = alphaig @@ -552,10 +576,13 @@ subroutine compute_snapwave(t) snapwave_mean_direction(k) = 0.0 snapwave_directional_spreading(k) = 0.0 snapwave_cg(k) = 0.0 + snapwave_qb(k) = 0.0 + snapwave_gam(k) = 0.0 endif ! if (snapwave_H_ig(k) <= 0.0) then - snapwave_Tp_ig(k) = 0.0 + snapwave_Tp_ig(k) = 0.0 + snapwave_cgig(k) = 0.0 endif enddo ! @@ -621,14 +648,23 @@ subroutine read_snapwave_input() ! ! Settings related to IG waves: call read_int_input(500,'snapwave_igwaves',igwaves_opt,1) - call read_real_input(500,'snapwave_alpha_ig',alpha_ig,1.0) !TODO choose whether snapwave_alphaig or snapwave_gamma_ig - call read_real_input(500,'snapwave_gammaig',gamma_ig,0.2) + call read_real_input(500,'snapwave_alpha_ig',alpha_ig,1.0) !FIXME choose whether snapwave_alphaig or snapwave_gamma_ig + call read_real_input(500,'snapwave_gammaig',gamma_ig,0.7) !FIXME choose whether snapwave_alphaig or snapwave_gamma_ig + !call read_real_input(500,'snapwave_gamma_fac_br',gamma_fac_br,2.0/3.0) ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point + call read_real_input(500,'snapwave_gamma_fac_br',gamma_fac_br,0.45) ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point call read_real_input(500,'snapwave_shinc2ig',shinc2ig,1.0) ! Ratio of how much of the calculated IG wave source term, is subtracted from the incident wave energy (0-1, 1=default=all energy as sink) call read_real_input(500,'snapwave_alphaigfac',alphaigfac,1.0) ! Multiplication factor for IG shoaling source/sink term call read_real_input(500,'snapwave_baldock_ratio_ig',baldock_ratio_ig,0.2) call read_int_input(500,'snapwave_ig_opt',ig_opt,1) call read_int_input(500,'snapwave_iterative_srcig',iterative_srcig_opt,0) ! Option whether to calculate IG source/sink term in iterative lower (better, but potentially slower, 1=default), or effectively based on previous timestep (faster, potential mismatch, =0) ! + ! IG steep slope related: + call read_real_input(500,'snapwave_steep_fac1',steep_fac1,0.0) ! Cut-off gamma, below this alphaig_steep = 0 + call read_real_input(500,'snapwave_steep_fac2',steep_fac2,0.1) ! Multiplication factor + call read_real_input(500,'snapwave_steep_fac3',steep_fac3,0.07) ! Cut-off beta, below this alphaig_steep = 0, and above this it increases with beta + call read_real_input(500,'snapwave_steep_fac4',steep_fac4,0.6) ! Exponent + call read_real_input(500,'snapwave_steep_fac5',steep_fac5,1.0) ! Cut-off gamma, above this alphaig_steep = 0 + ! ! IG boundary conditions options: call read_int_input(500,'snapwave_use_herbers',herbers_opt,1) ! Choice whether you want IG Hm0&Tp be calculated by herbers (=1, default), or want to specify user defined values (0> then snapwave_eeinc2ig & snapwave_Tinc2ig are used) call read_int_input(500,'snapwave_tpig_opt',tpig_opt,1) ! IG wave period option based on Herbers calculated spectrum, only used if snapwave_use_herbers = 1. Options are: 1=Tm01 (default), 2=Tpsmooth, 3=Tp, 4=Tm-1,0 diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index d90cc331d..7f5d267c5 100644 --- a/source/src/snapwave/snapwave_data.f90 +++ b/source/src/snapwave/snapwave_data.f90 @@ -73,6 +73,8 @@ module snapwave_data real*4, dimension(:), allocatable :: beta real*4, dimension(:), allocatable :: srcig real*4, dimension(:), allocatable :: alphaig + real*4, dimension(:), allocatable :: qb + real*4, dimension(:), allocatable :: gam ! integer*4, dimension(:), allocatable :: index_snapwave_in_quadtree integer*4, dimension(:), allocatable :: index_quadtree_in_snapwave @@ -258,6 +260,7 @@ module snapwave_data ! integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) real*4 :: alpha_ig,gamma_ig ! coefficients in Baldock wave breaking dissipation model for IG waves + real*4 :: gamma_fac_br ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point real*4 :: shinc2ig ! Ratio of how much of the calculated IG wave source term, is subtracted from the incident wave energy (0-1, 0=default) real*4 :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0 real*4 :: eeinc2ig ! ratio of incident wave energy as first estimate of IG wave energy at boundary @@ -269,6 +272,7 @@ module snapwave_data ! ... or just a priori based on effectively incident wave energy from previous timestep only integer :: herbers_opt ! Choice whether you want IG Hm0&Tp be calculated by herbers (=1, default), or want to specify user defined values (0> then snapwave_eeinc2ig & snapwave_Tinc2ig are used) integer :: tpig_opt ! IG wave period option based on Herbers calculated spectrum, only used if herbers_opt = 1. Options are: 1=Tm01 (default), 2=Tpsmooth, 3=Tp, 4=Tm-1,0 + real*4 :: steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5 ! ! Switches logical :: igwaves ! switch whether include IG or not diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index bea58fa33..988cc72c3 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -136,6 +136,8 @@ subroutine initialize_snapwave_domain() allocate(beta(no_nodes)) allocate(srcig(no_nodes)) allocate(alphaig(no_nodes)) + allocate(qb(no_nodes)) + allocate(gam(no_nodes)) ! allocate(uorb(no_nodes)) allocate(ctheta(ntheta,no_nodes)) allocate(ctheta_ig(ntheta,no_nodes)) @@ -200,6 +202,7 @@ subroutine initialize_snapwave_domain() prev360 = 0 H = 0.0 H_ig = 0.0 + Dw = 0.0 aa = 0.0 sig = 0.0 WsorE = 0.0 diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 409709f55..234604403 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -86,7 +86,8 @@ subroutine compute_wave_field() enddo if (igwaves) then do k = 1, no_nodes - Hmx_ig(k) = 0.88/kwav_ig(k)*tanh(gamma_ig*kwav_ig(k)*depth(k)/0.88) ! Note - uses gamma_ig + !Hmx_ig(k) = 0.88/kwav_ig(k)*tanh(gamma_ig*kwav_ig(k)*depth(k)/0.88) ! Note - uses gamma_ig + Hmx_ig(k) = gamma_ig*depth(k) enddo else Hmx_ig = 0.0 @@ -133,9 +134,10 @@ subroutine compute_wave_field() aa, sig, jadcgdx, sigmin, sigmax,& c_dispT, WsorE, WsorA, SwE, SwA, Tpini, & igwaves,kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & - beta, srcig, alphaig, Dw_ig, Df_ig, & + beta, srcig, alphaig, Dw_ig, Df_ig, qb, gam, & + steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5, & vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, & - zb, nwav, ig_opt, alpha_ig, gamma_ig, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) + zb, nwav, ig_opt, alpha_ig, gamma_ig, gamma_fac_br, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) ! call timer(t3) ! @@ -158,9 +160,10 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & aa, sig, jadcgdx, sigmin, sigmax,& c_dispT, WsorE, WsorA, SwE, SwA, Tpini, & igwaves,kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & - betamean, srcig, alphaig, Dw_ig, Df_ig, & + betamean, srcig, alphaig, Dw_ig, Df_ig, qb, gam, & + steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5, & vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, & - zb, nwav, ig_opt, alfa_ig, gamma_ig, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) + zb, nwav, ig_opt, alfa_ig, gamma_ig, gamma_fac_br, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) ! use snapwave_windsource !use snapwave_ncoutput ! TL: removed, we don't use this in SF+SW @@ -188,13 +191,15 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4, dimension(no_nodes), intent(in) :: kwav_ig ! wave number real*4, dimension(no_nodes), intent(inout) :: cg ! group velocity real*4, dimension(ntheta,no_nodes), intent(inout):: ctheta ! refractioon speed - real*4, dimension(no_nodes), intent(in) :: cg_ig ! group velocity + real*4, dimension(no_nodes), intent(inout) :: cg_ig ! group velocity real*4, dimension(no_nodes), intent(in) :: nwav ! wave number n real*4, dimension(ntheta,no_nodes), intent(inout):: ee ! real*4, dimension(ntheta,no_nodes), intent(inout):: ee_ig ! real*4, dimension(ntheta,no_nodes), intent(in) :: ctheta_ig ! refractioon speed real*4, dimension(no_nodes), intent(in) :: fw ! wave friction factor real*4, dimension(no_nodes), intent(in) :: fw_ig ! wave friction factor + real*4, dimension(no_nodes), intent(out) :: qb ! Fraction of breaking waves according to Baldock's formulation + real*4, dimension(no_nodes), intent(out) :: gam ! Local incident wave height water depth ratio real*4, dimension(no_nodes), intent(out) :: betamean ! Mean local bed slope parameter real*4, dimension(no_nodes), intent(out) :: srcig ! Directionally averaged incident wave sink/infragravity source term real*4, dimension(no_nodes), intent(out) :: alphaig ! Mean IG shoaling parameter alpha @@ -205,7 +210,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4, intent(in) :: baldock_ratio_ig ! option controlling from what depth wave breaking should take place for IG waves: (Hk_ig>baldock_ratio_ig*Hmx_ig(k)), default baldock_ratio_ig=0.2 real*4, dimension(no_nodes), intent(inout) :: H ! wave height - TODO - TL - CHECK > inout needed to have updated 'H' for determining srcig real*4, dimension(no_nodes), intent(out) :: H_ig ! wave height - real*4, dimension(no_nodes), intent(out) :: Dw ! wave breaking dissipation + real*4, dimension(no_nodes), intent(inout) :: Dw ! wave breaking dissipation real*4, dimension(no_nodes), intent(out) :: Dw_ig ! wave breaking dissipation IG real*4, dimension(no_nodes), intent(out) :: F ! wave force Dw/C/rho/h real*4, dimension(no_nodes), intent(out) :: Df ! wave friction dissipation @@ -219,7 +224,8 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4, dimension(no_nodes), intent(in) :: u10 ! wind speed and direction integer, intent(in) :: niter ! max number of iterations real*4, intent(in) :: crit ! relative accuracy for stopping criterion - integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) + integer, intent(in) :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) + real*4, intent(in) :: steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5 ! ! wind source vars ! @@ -257,7 +263,9 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4, dimension(:,:), allocatable :: srcig_local ! Energy source/sink term because of IG wave energy transfer from incident waves real*4, dimension(:,:), allocatable :: beta_local ! Local bed slope based on bed level per direction real*4, dimension(:,:), allocatable :: alphaig_local ! Local infragravity wave shoaling parameter alpha - real*4, dimension(:,:), allocatable :: depthprev ! water depth at upwind intersection point per direction + real*4, dimension(:,:), allocatable :: depthprev ! water depth at upwind intersection point per direction + real*4, dimension(:,:), allocatable :: qb_local ! local percentage of breaking waves Qb + real*4, dimension(:,:), allocatable :: gam_local ! local incident wave height over water depth ratio real*4, dimension(:), allocatable :: dee ! difference with energy previous iteration real*4, dimension(:), allocatable :: eeprev, cgprev ! energy density and group velocity at upwind intersection point real*4, dimension(:), allocatable :: eeprev_ig, cgprev_ig ! energy density and group velocity at upwind intersection point @@ -292,6 +300,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4 :: Ek_ig real*4 :: Hk_ig real*4 :: alfa_ig,gamma_ig ! coefficients in Baldock wave breaking dissipation model for IG waves + real*4 :: gamma_fac_br ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point real*4 :: eeinc2ig ! ratio of incident wave energy as first estimate of IG wave energy at boundary real*4 :: Tinc2ig ! ratio compared to period Tinc to estimate Tig real*4 :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0 @@ -340,6 +349,9 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & allocate(depthprev(ntheta,no_nodes)); depthprev=0.0 allocate(beta_local(ntheta,no_nodes)); beta_local=0.0 allocate(alphaig_local(ntheta,no_nodes)); alphaig_local=0.0 + allocate(qb_local(ntheta,no_nodes)); qb_local=0.0 + allocate(gam_local(ntheta,no_nodes)); gam_local=0.0 + endif ! if (wind) then @@ -359,7 +371,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & enddo ! df = 0.0 - dw = 0.0 + !dw = 0.0 ! TODO - TL: CHeck > needed for restart for IG > set to 0 now in snapwave_domain.f90 F = 0.0 ! ok = 0 @@ -471,10 +483,10 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! As defined in Leijnse, van Ormondt, van Dongeren, Aerts & Muis et al. 2024 ! ! Actual determining of source term: - ! - call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + ! + call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, Dw, Hmx, qb_local, gam_local, gamma, gamma_fac_br, steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) ! - ! inout: alphaig_local, srcig_local - eeprev, eeprev_ig, cgprev, beta_local + ! inout: alphaig_local, srcig_local, cgprev, beta_local, qb_local, gam_local ! in: the rest ! ! NOTE - This is based on the energy in the precious SnapWave timestep 'ee' and 'ee_ig', and waveheight 'H', which should therefore be made available. @@ -555,7 +567,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! ! Actual determining of source term - every first sweep of iteration ! - call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, Dw, Hmx, qb_local, gam_local, gamma, gamma_fac_br, steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) ! endif ! @@ -595,7 +607,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! enddo ! - Ek = sum(eeprev)*dtheta ! to check + Ek = sum(eeprev)*dtheta ! to check - Maarten: should not be 'eeprev'? ! depthlimfac = max(1.0, (sqrt(Ek/rhog8)/(gammax*depth(k)))**2.0) Hk = min(sqrt(Ek/rhog8), gamma*depth(k)) @@ -661,26 +673,29 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! do itheta = 1, ntheta ! - R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcig_local(itheta, k) * shinc2ig + R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcig_local(itheta, k) * shinc2ig ! enddo ! do itheta = 2, ntheta - 1 ! A(itheta) = -ctheta(itheta - 1, k)*oneover2dtheta - B(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverE(k) + ! + B(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverE(k) + ! C(itheta) = ctheta(itheta + 1, k)*oneover2dtheta ! enddo ! A(1) = -ctheta(ntheta, k)*oneover2dtheta - B(1) = oneoverdt + cg(k)/ds(1,k) + DoverE(k) C(1) = ctheta(2, k)*oneover2dtheta ! A(ntheta) = -ctheta(ntheta - 1, k)*oneover2dtheta - B(ntheta) = oneoverdt + cg(k)/ds(ntheta,k) + DoverE(k) C(ntheta) = ctheta(1, k)*oneover2dtheta ! + B(1) = oneoverdt + cg(k)/ds(1,k) + DoverE(k) + B(ntheta) = oneoverdt + cg(k)/ds(ntheta,k) + DoverE(k) + ! ! Solve tridiagonal system per point ! if (wind) then @@ -742,7 +757,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! IG ! if (igwaves) then - Ek_ig = sum(eeprev_ig)*dtheta + Ek_ig = sum(eeprev_ig)*dtheta ! to check - Maarten: should not be 'eeprev_ig'? !Hk_ig = sqrt(Ek_ig/rhog8) !org trunk Hk_ig = min(sqrt(Ek_ig/rhog8), gamma_ig*depth(k)) !TL: Question - why not this one? Ek_ig = rhog8*Hk_ig**2 @@ -764,14 +779,16 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! do itheta = 1, ntheta ! - R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta)*eeprev_ig(itheta)/ds(itheta, k) + srcig_local(itheta, k) !TL: new version + R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta)*eeprev_ig(itheta)/ds(itheta, k) + srcig_local(itheta, k) ! enddo ! do itheta = 2, ntheta - 1 ! A_ig(itheta) = -ctheta_ig(itheta - 1, k)*oneover2dtheta + ! B_ig(itheta) = oneoverdt + cg_ig(k)/ds(itheta,k) + DoverE_ig(k) + ! C_ig(itheta) = ctheta_ig(itheta + 1, k)*oneover2dtheta ! enddo @@ -779,20 +796,27 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & if (ctheta_ig(1,k)<0) then A_ig(1) = 0.0 B_ig(1) = oneoverdt - ctheta_ig(1, k)/dtheta + cg_ig(k)/ds(1, k) + DoverE_ig(k) + ! C_ig(1) = ctheta_ig(2, k)/dtheta else A_ig(1)=0.0 + ! B_ig(1)=1.0/dt + cg_ig(k)/ds(1, k) + DoverE_ig(k) + ! C_ig(1)=0.0 endif ! if (ctheta_ig(ntheta, k)>0) then A_ig(ntheta) = -ctheta_ig(ntheta - 1, k)/dtheta + ! B_ig(ntheta) = oneoverdt + ctheta_ig(ntheta, k)/dtheta + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k) + ! C_ig(ntheta) = 0.0 else A_ig(ntheta) = 0.0 - B_ig(ntheta) = oneoverdt + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k) + ! + B_ig(ntheta) = oneoverdt + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k) + ! C_ig(ntheta) = 0.0 endif ! @@ -938,6 +962,16 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! srcig(k) = sum(srcig_local(:,k)) /ntheta ! real mean ! + !qb(k) = sum(qb_local(:,k)) /ntheta ! real mean + qb(k) = maxval(qb_local(:,k)) ! max + ! + !gam(k) = sum(gam_local(:,k)) /ntheta ! real mean + gam(k) = maxval(gam_local(:,k)) ! max + ! + !if (qb(k) > 0.01) then + ! write(*,*)'k qb gam',k, qb_local(:,k), gam_local(:,k) + !endif + ! endif ! if (wind) then @@ -1017,8 +1051,8 @@ subroutine baldock (rho,g,alfa,gamma,depth,H,T,opt,Dw,Hmax) endif ! end subroutine baldock - - subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + + subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, Dw, Hmx, qb_local, gam_local, gamma, gamma_fac_br, steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) ! implicit none ! @@ -1028,22 +1062,29 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d real*4, dimension(2,ntheta,no_nodes),intent(in) :: w ! weights of upwind grid points, 2 per grid point and per wave direction real*4, dimension(ntheta,no_nodes), intent(in) :: ds ! distance to interpolated upwind point, per grid point and direction integer, dimension(2,ntheta,no_nodes),intent(in) :: prev ! two upwind grid points per grid point and wave direction - real*4, dimension(no_nodes), intent(in) :: cg_ig ! group velocity + real*4, dimension(no_nodes), intent(inout) :: cg_ig ! group velocity real*4, dimension(no_nodes), intent(in) :: nwav ! wave number n real*4, dimension(no_nodes), intent(in) :: depth ! water depth real*4, dimension(no_nodes), intent(in) :: zb ! actual bed level - real*4, dimension(no_nodes), intent(in) :: H ! wave height + real*4, dimension(no_nodes), intent(in) :: H ! wave height real*4, dimension(ntheta,no_nodes), intent(in) :: ee ! energy density - real*4, dimension(ntheta,no_nodes), intent(in) :: ee_ig ! energy density infragravity waves + real*4, dimension(ntheta,no_nodes), intent(in) :: ee_ig ! energy density infragravity waves integer, intent(in) :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) real*4, intent(in) :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0 + real*4, intent(in) :: dtheta ! directional resolution + real*4, intent(in) :: gamma ! coefficients in Baldock wave breaking dissipation + real*4, intent(in) :: gamma_fac_br ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point + real*4, dimension(no_nodes), intent(in) :: Dw ! wave breaking dissipation + real*4, dimension(no_nodes), intent(in) :: Hmx ! Hmax + real*4, intent(in) :: steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5 ! ! Inout variables real*4, dimension(:,:), intent(inout) :: alphaig_local ! Local infragravity wave shoaling parameter alpha real*4, dimension(:,:), intent(inout) :: srcig_local ! Energy source/sink term because of IG wave shoaling - real*4, dimension(:), intent(inout) :: eeprev, cgprev ! energy density and group velocity at upwind intersection point - real*4, dimension(:), intent(inout) :: eeprev_ig ! energy density at upwind intersection point - real*4, dimension(ntheta,no_nodes), intent(inout):: beta_local ! Local bed slope based on bed level per direction + real*4, dimension(:), intent(inout) :: cgprev ! group velocity at upwind intersection point + real*4, dimension(ntheta,no_nodes), intent(inout):: beta_local ! Local bed slope based on bed level per direction + real*4, dimension(ntheta,no_nodes), intent(inout):: qb_local ! local percentage of breaking waves Qb + real*4, dimension(ntheta,no_nodes), intent(inout):: gam_local ! local incident wave height over water depth ratio ! ! Internal variables integer :: itheta ! directional counter @@ -1051,22 +1092,73 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d integer :: k1,k2 ! upwind counters (k is grid index) real*4 :: gam ! local gamma (Hinc / depth ratio) real*4, dimension(ntheta,no_nodes) :: depthprev ! water depth at upwind intersection point - real*4, dimension(ntheta,no_nodes) :: Sxx ! Radiation Stress + real*4, dimension(no_nodes) :: Sxx ! Radiation Stress real*4, dimension(:), allocatable :: Sxxprev ! radiation stress at upwind intersection point real*4, dimension(:), allocatable :: Hprev ! Incident wave height at upwind intersection point + real*4, dimension(:), allocatable :: Eprev ! Mean incident wave energy at upwind intersection point + real*4, dimension(:), allocatable :: Eprev_ig ! Mean infragravity wave energy at upwind intersection point + real*4, dimension(no_nodes) :: E_local ! mean wave energy waves - just local + real*4, dimension(no_nodes) :: E_ig_local ! mean wave energy infragravity waves - just local real*4 :: dSxx ! difference in Radiation stress real*4 :: Sxx_cons ! conservative estimate of radiation stress using conservative shoaling + real*4 :: delta_Dw ! difference of Dw compared to upwind point, to get sign for max breaking point + real*4 :: Qb ! Percentage of breaking incident waves + real*4 :: transition_factor ! Transition factor for letting srcig go to zero smoothly, around gamma*gamma_fac_br + real*4 :: transition_factor_width_1 ! Width factor of generalized (Fermi–Dirac style) transfer function with adjustable midpoint and width + real*4 :: transition_factor_width_2 ! Width factor of generalized (Fermi–Dirac style) transfer function with adjustable midpoint and width + real*4 :: gamma_fac_br_transition ! Transitioned version of gamma_fac_br, so that for steep slopes it remains 1.0 + real*4 :: beta_limit_1 ! Cut-off beta_local for end of validity alphaig formulation of Leijnse et al. 2024 + real*4 :: beta_limit_2 ! Beta_local limit for transition function ! ! Allocate internal variables allocate(Sxxprev(ntheta)) allocate(Hprev(ntheta)) + allocate(Eprev(ntheta)) + allocate(Eprev_ig(ntheta)) ! Sxx = 0.0 + Hprev = 0.0 + Eprev = 0.0 + Eprev_ig = 0.0 + ! + E_local = 0.0 + E_ig_local = 0.0 + ! + ! Used is generalized (Fermi–Dirac style) transfer function with adjustable midpoint and width + ! + transition_factor_width_1 = 0.005 + transition_factor_width_2 = 0.002 + beta_limit_1 = 0.07 + !beta_limit_2 = beta_limit_1 - 0.01 + beta_limit_2 = beta_limit_1 - 0.02 + ! + ! Precompute all Sxx - FIXME - add parallellisation ! do k = 1, no_nodes ! - if (inner(k)) then + if (inner(k)) then + ! + ! Update E (not saved from previous timestep) ! + E_local(k) = sum(ee(:,k))*dtheta + ! + ! Update E_ig (not saved from previous timestep) + E_ig_local(k) = sum(ee_ig(:, k))*dtheta + ! + endif + ! + ! Calculate Sxx for all points - FIXME double check correct also for not inner + ! + Sxx(k) = (((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) * E_local(k)) + ! + enddo + ! + ! Actual computation of srcig + ! + do k = 1, no_nodes + ! + if (inner(k)) then + ! ! Compute exchange source term inc to ig waves - per direction ! do itheta = 1, ntheta @@ -1081,6 +1173,9 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! beta_local(itheta,k) = max((w(1, itheta, k)*(zb(k) - zb(k1)) + w(2, itheta, k)*(zb(k) - zb(k2)))/ds(itheta, k), 0.0) ! + ! FIXME - shorter, but also same result? + !beta_local(itheta,k) = max(zb(k) - (w(1, itheta, k) * zb(k1) + w(2, itheta, k) * zb(k2))/ds(itheta, k), 0.0) + ! ! Notes: ! - use actual bed level now for slope, because depth changes because of wave setup/tide/surge ! - in zb, depth is negative > therefore zb(k) minus zb(k1) @@ -1088,30 +1183,66 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! !betan_local(itheta,k) = (beta/sigm_ig)*sqrt(9.81/max(depth(k), hmin)) ! TL: in case in the future we would need the normalised bed slope again ! - ! TL - Note: cg_ig = cg + ! Fraction of breaking waves, based on H(k) + !Qb = min(max(exp(-(Hmx(k)/H(k))**2), 0.0), 1.0) ! Qb percentage of breaking waves according to Baldock's formulation, between 0 and 1 + ! Base on upwind point: + Qb = min(max(exp(-((w(1, itheta, k)*Hmx(k1) + w(2, itheta, k)*Hmx(k2)) / Hprev(itheta))**2), 0.0), 1.0) ! Qb percentage of breaking waves according to Baldock's formulation, between 0 and 1 + ! + qb_local(itheta, k) = Qb + ! cgprev(itheta) = w(1, itheta, k)*cg_ig(k1) + w(2, itheta, k)*cg_ig(k2) - ! - Sxx(itheta,k1) = ((2.0 * max(0.0,min(1.0,nwav(k1)))) - 0.5) * ee(itheta, k1) ! limit so value of nwav is between 0 and 1 - Sxx(itheta,k2) = ((2.0 * max(0.0,min(1.0,nwav(k2)))) - 0.5) * ee(itheta, k2) ! limit so value of nwav is between 0 and 1 ! - Sxxprev(itheta) = w(1, itheta, k)*Sxx(itheta,k1) + w(2, itheta, k)*Sxx(itheta,k2) + Sxxprev(itheta) = w(1, itheta, k)*Sxx(k1) + w(2, itheta, k)*Sxx(k2) ! - eeprev(itheta) = w(1, itheta, k)*ee(itheta, k1) + w(2, itheta, k)*ee(itheta, k2) - eeprev_ig(itheta) = w(1, itheta, k)*ee_ig(itheta, k1) + w(2, itheta, k)*ee_ig(itheta, k2) + Eprev(itheta) = w(1, itheta, k)*E_local(k1) + w(2, itheta, k)*E_local(k2) + Eprev_ig(itheta) = w(1, itheta, k)*E_ig_local(k1) + w(2, itheta, k)*E_ig_local(k2) ! - Hprev(itheta) = w(1, itheta, k)*H(k1) + w(2, itheta, k)*H(k2) + Hprev(itheta) = w(1, itheta, k)*H(k1) + w(2, itheta, k)*H(k2) ! ! Determine relative waterdepth 'gam' ! gam = max(0.5*(Hprev(itheta)/depthprev(itheta,k) + H(k)/depth(k)), 0.0) ! mean gamma over current and upwind point ! + gam_local(itheta, k) = gam + ! + ! Free waves and no IG source/sink term if incident waves start breaking + ! + ! Adjust cg_ig for free infragravity waves release in surfzone + ! TL - Note: cg_ig = cg + !if (ig_opt == X) then + ! ! + ! if (gam > (gamma_fac_br * gamma)) then + ! ! + ! cg_ig(k) = sqrt(9.81 * depth(k)) + ! ! + ! endif + ! ! + !endif + ! ! Determine dSxx and IG source/sink term 'srcig' ! - if (ig_opt == 1 .or. ig_opt == 2) then + if (ig_opt == 1 .or. ig_opt == 2 .or. ig_opt == 11 .or. ig_opt == 12 .or. ig_opt == 13 .or. ig_opt == 14 .or. ig_opt == 15) then ! ! Calculate shoaling parameter alpha_ig following Leijnse et al. (2024) ! + if (ig_opt == 11 .or. ig_opt == 12 .or. ig_opt == 13) then + ! + ! Limit beta to max 0.07 (=beta_limit_1) before going into alphaig parametrisation + ! + beta_local(itheta,k) = min(beta_local(itheta,k), beta_limit_1) + ! + endif + ! call estimate_shoaling_parameter_alphaig(beta_local(itheta,k), gam, alphaig_local(itheta,k)) ! [input, input, output] + ! + ! Steep slope addition + ! + if (ig_opt == 14 .or. ig_opt == 15) then + ! + call estimate_shoaling_parameter_alphaig_steep_slopes(beta_local(itheta,k), gam, alphaig_local(itheta,k), steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) + ! [input, input, inout, input, input, input, input, input] + ! + endif ! ! Now calculate source term component ! @@ -1121,26 +1252,69 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d srcig_local(itheta, k) = 0.0 !Avoid big jumps in dSxx that can happen if a upwind point is a boundary point with Hinc=0 ! else - ! - if (ig_opt == 1) then ! Option using conservative shoaling for dSxx/dx + ! + if (ig_opt == 1 .or. ig_opt == 11 .or. ig_opt == 12 .or. ig_opt == 13 .or. ig_opt == 14 .or. ig_opt == 15) then ! Option using conservative shoaling for dSxx/dx + ! + ! Calculate Sxx based on conservative shoaling of upwind point's energy: + ! Sxx_cons = E(i-1) * Cg(i-1) / Cg * (2 * n(i) - 0.5) + Sxx_cons = Eprev(itheta) * cgprev(itheta) / cg_ig(k) * ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) + ! Note - limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite + ! + dSxx = Sxx_cons - Sxxprev(itheta) + ! + elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx + ! + dSxx = Sxx(itheta) - Sxxprev(itheta) + endif ! - ! Calculate Sxx based on conservative shoaling of upwind point's energy: - ! Sxx_cons = E(i-1) * Cg(i-1) / Cg * (2 * n(i) - 0.5) - Sxx_cons = eeprev(itheta) * cgprev(itheta) / cg_ig(k) * ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) - ! Note - limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite + dSxx = max(dSxx, 0.0) ! - dSxx = Sxx_cons - Sxxprev(itheta) + !if (ig_opt == 1 .or. ig_opt == 2.or. ig_opt == 12 .or. ig_opt == 13) then ! - elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx + ! Base on E_prev_ig instead of eeprev_ig(itheta) > no bins but total energy + ! NOTE - already here multiplied with ee(itheta,k), for direct inclusion in 'R'-term + srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(Eprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) /max(E_local(k), 1.0e-6) * ee(itheta,k) ! - dSxx = Sxx(itheta,k) - Sxxprev(itheta) + !elseif (ig_opt == 20) then + ! + ! NOTE - in main script this is multiplied with ee(itheta,k) to get directional energy, for direct inclusion in 'B'-term + ! + !srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(Eprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) /max(E_local(k), 1.0e-6) !* ee(itheta,k) + !endif + ! + ! Limit srcig to 0 after waves start (significantly) breaking, as defined here as gam=Hrms,inc / h > (gamma_fac_br * gamma) + ! + ! Ergo, it is assumed that after this point IG waves are free, and no bound wave forcing is happening anymore, so srcig should be 0 from here on + ! + if (ig_opt == 12) then + ! + ! Let srcig transition to 0 more smoothly using fac_transition that reduced from 1 to 0 around gamma_fac_br * snapwave_gamma + ! Similar as before, but then smooth: + ! ! Note - gam is in Hrms + ! if (gam > (gamma_fac_br * gamma)) then + ! ! + ! srcig_local(itheta, k) = 0.0 + ! ! + ! endif + ! + transition_factor = 1.0 - (1.0 / (1.0 + exp(- (gam - (gamma_fac_br * gamma)) / transition_factor_width_1))) + ! + srcig_local(itheta, k) = transition_factor * srcig_local(itheta, k) + ! + elseif (ig_opt == 13 .or. ig_opt == 15) then + ! + ! Let srcig transition to 0 more smoothly using fac_transition that reduced from 1 to 0 around gamma_fac_br * snapwave_gamma + ! + ! But, only for beta_local < 0.07, so adjust based on beta_local so that transition_factor = 1.0 for Beta_local = 0.07 + ! + gamma_fac_br_transition = gamma_fac_br + ((1-gamma_fac_br) / (1 + exp(- (beta_local(itheta,k) - beta_limit_2) / transition_factor_width_2))) + ! + transition_factor = 1.0 - (1.0 / (1.0 + exp(- (gam - (gamma_fac_br_transition * gamma)) / transition_factor_width_1))) + ! + srcig_local(itheta, k) = transition_factor * srcig_local(itheta, k) + ! + endif ! - endif - ! - dSxx = max(dSxx, 0.0) - ! - srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(eeprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) - ! endif ! else ! TL: option to add future parameterisations here for e.g. coral reef type coasts @@ -1176,6 +1350,7 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) beta3 = 17.7104 beta4 = 1 beta5 = 0.7 + !beta5 = 0.5 beta6 = 0.11841 beta7 = 0.34037 ! @@ -1190,8 +1365,10 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) alphaig = exp(-beta3 * beta ** beta4) * ((beta5 - gam) * beta6 + (beta7 - gam) * (beta1 / beta ** beta2)) ! elseif (gam >= beta7) then ! shallow water - for gam>0.7 the fit automatically goes to 0 + !elseif (gam >= beta7 .and. gam < 0.5) then ! shallow water - for gam>0.7 the fit automatically goes to 0 ! - alphaig = exp(-beta3 * beta ** beta4) * (max(beta5 - gam, 0.0)) * beta6 + alphaig = exp(-beta3 * beta ** beta4) * (max(beta5 - gam, 0.0)) * beta6 + !alphaig = exp(-beta3 * beta ** beta4) * (max(0.5 - gam, 0.0)) * beta6 ! else ! for safety, but negative gamma should not occur ! @@ -1208,6 +1385,72 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) ! end subroutine estimate_shoaling_parameter_alphaig + + subroutine estimate_shoaling_parameter_alphaig_steep_slopes(beta, gam, alphaig, fac1, fac2, fac3, fac4, fac5) + ! [input, input, inout, input, input, input, input, input] + real*4, intent(in) :: beta + real*4, intent(in) :: gam + real*4, intent(inout) :: alphaig + ! + real*4, intent(in) :: fac1, fac2, fac3, fac4, fac5 + real*4 :: alphaig_steep + ! + ! Estimate shoaling parameter alphaig - for steep slopes with beta > 0.07 + ! These were not covered in training dataset of Leijnse et al. 2024 + ! + !alphaig_total = alphaig Eq11 + alphaig_steep + !alphaig_steep = 0.1 * max(beta-0.07, 0)^0.6 * max(1.0-gamma, 0) + !alphaig_steep = fac2 * max(beta-fac3, 0)**fac4 * max(fac5-gam, 0) + ! + ! Determine constants + ! + !fac1 = 0.3 !Cut-off gamma, below this alphaig_steep = 0 + !fac1 = 0.0 !Cut-off gamma, below this alphaig_steep = 0 + !fac2 = 0.1 ! Multiplication factor + !fac3 = 0.07 ! Cut-off beta, below this alphaig_steep = 0, and above this it increases with beta + !fac4 = 0.6 ! Exponent + !fac5 = 1.0 ! Cut-off gamma, above this alphaig_steep = 0 + ! + ! For deep water or negative slope, alphaig_steep = 0 + ! + alphaig_steep = 0.0 + ! + !if (beta > 0.0) then + ! write(*,*)'Starting alphaig steep slope addition, beta, gam, alphaig before', beta, gam, alphaig + !endif + ! + ! If positively increasing local bed slope beta + ! + if (beta > 0.0) then + ! + if (gam >= fac1) then ! shallow(er) water - for gam>1.0 (=fac5) the fit automatically goes to 0 + ! + alphaig_steep = fac2 * (max(beta - fac3, 0.0) ** fac4) * (max(fac5 - gam, 0.0)) + ! + else ! for gam < fac1 + ! + alphaig_steep = 0.0 + ! + endif + ! + endif + ! + !if (beta > 0.0) then + ! write(*,*)'Calculated alphaig_steep addition', alphaig_steep + !endif + ! + ! Combine + ! + alphaig = alphaig + alphaig_steep + ! + ! Limit total alphaig between [0, 1] to prevent large overshoots in case of low gamma and very small beta + ! + alphaig = max(alphaig, 0.0) + alphaig = min(alphaig, 1.0) + ! + end subroutine estimate_shoaling_parameter_alphaig_steep_slopes + + subroutine hpsort_eps_epw (n, ra, ind, eps) !--------------------------------------------------------------------- ! sort an array ra(1:n) into ascending order using heapsort algorithm,