diff --git a/Makefile b/Makefile deleted file mode 100644 index b68b84676..000000000 --- a/Makefile +++ /dev/null @@ -1,38 +0,0 @@ -include Makefile.in - -all: - make metis - make parms - make cleanfesom - make fesom_ini - make cleanfesom - make fesom - -metis: - make -C ./lib/metis-5.1.0 - -parms: - make -C ./lib/parms - -fesom_ini: - make -C ./src run_ini - -fesom: - make -C ./src run - -cleanparms: - make -C ./lib/parms cleanall - -cleanmetis: - make -C ./lib/metis-5.1.0 clean - -cleanfesom: - make -C ./src clean - -cleanall: - make cleanparms - make cleanmetis - make cleanfesom - -clean: - make cleanfesom diff --git a/Makefile.in_gnu_impi b/Makefile.in_gnu_impi deleted file mode 100644 index cdf9a6792..000000000 --- a/Makefile.in_gnu_impi +++ /dev/null @@ -1,45 +0,0 @@ - -# Compiler and Linker -CC = gcc -MPICC = mpicc -MPIFC = mpif90 -LD = $(MPIFC) - -# Optimization specs for compiler - used for FESOM an parms. -# Metis falls back to conservative and non-compiler specific -O2, as -# it is only used in the less runtime-critical setup phase. -OPT = -O3 -march=native -funroll-loops -ffree-line-length-none -COPT = -O3 -march=native -funroll-loops -DLINUX -DUNDER_ -DMPI2 -Iinclude -## Options for debugging: -# OPT = -g -fbounds-check -fbacktrace - -# Preprocessor flags for FESOM2 -CPP_DEFS = -DDISABLE_MULTITHREADING - -# NetCDF library and include definitions -NCFDIR = /global/AWIsoft/netcdf/4.6.1_gnu -NC_LIB = -L$(NCFDIR)/lib -Wl,-rpath,$(NCFDIR)/lib -lnetcdff -NC_INC = -I$(NCFDIR)/include - -# Definitions for MPI, if not set by compiler wrappers -MPIROOT = -MPI_LIB = -MPI_INC = - -## For partitioning, the FESOM initialization offers an interface to Metis 5. -## The option -DPART_WEIGHTED distributes 2D _and_ 3D-nodes equally, but the -## partitioning might be quite irregular. -## Without DPART_WEIGHTED, only the number of 2D nodes is considered for load balancing. -## You can start Metis 5 with different random seeds -> different partitions. - -METIS_DEFS = -DMETIS_VERSION=5 -DPART_WEIGHTED -DMETISRANDOMSEED=35243 -# METISRANDDOMSEED: any integer - -# Used to compile the fesom wrapper to parms -PARMS_DEFS = -DUSE_MPI -DREAL=double -DDBL -DVOID_POINTER_SIZE_8 - -####################################################### -# for pARMS and metis Library, archive and options - -AR = ar qv -RANLIB = ranlib diff --git a/Makefile.in_intel b/Makefile.in_intel deleted file mode 100644 index 1bb5c4800..000000000 --- a/Makefile.in_intel +++ /dev/null @@ -1,45 +0,0 @@ - -# Compiler and Linker -CC = icc -MPICC = mpiicc -MPIFC = mpiifort -LD = $(MPIFC) - -# Optimization specs for compiler - used for FESOM an parms. -# Metis falls back to conservative and non-compiler specific -O2, as -# it is only used in the less runtime-critical setup phase. -OPT = -O3 -r8 -fp-model precise -xHost -ip -implicitnone -static-intel -COPT = -O3 -xHost -DLINUX -DUNDER_ -DMPI2 -Iinclude -## Options for debugging: -# OPT = -g -xHost -traceback -fpe:0 -check all -static-intel - -# Preprocessor flags for FESOM2 -CPP_DEFS = -DDISABLE_MULTITHREADING - -# NetCDF library and include definitions -NCFDIR = /global/AWIsoft/netcdf/4.6.1_intel_18.0.3 -NC_LIB = -L$(NCFDIR)/lib -Wl,-rpath,$(NCFDIR)/lib -lnetcdff -NC_INC = -I$(NCFDIR)/include - -# Definitions for MPI, if not set by compiler wrappers -MPIROOT = -MPI_LIB = -MPI_INC = - -## For partitioning, the FESOM initialization offers an interface to Metis 5. -## The option -DPART_WEIGHTED distributes 2D _and_ 3D-nodes equally, but the -## partitioning might be quite irregular. -## Without DPART_WEIGHTED, only the number of 2D nodes is considered for load balancing. -## You can start Metis 5 with different random seeds -> different partitions. - -METIS_DEFS = -DMETIS_VERSION=5 -DPART_WEIGHTED -DMETISRANDOMSEED=35243 -# METISRANDDOMSEED: any integer - -# Used to compile the fesom wrapper to parms -PARMS_DEFS = -DUSE_MPI -DREAL=double -DDBL -DVOID_POINTER_SIZE_8 - -####################################################### -# for pARMS and metis Library, archive and options - -AR = ar qv -RANLIB = ranlib diff --git a/src/Makefile b/src/Makefile deleted file mode 100755 index 598dc8d39..000000000 --- a/src/Makefile +++ /dev/null @@ -1,165 +0,0 @@ -###################################################### -# Makefile -###################################################### - -###################################################### -# Include architecture-specific definitions - -include ../Makefile.in - -# Definition of metis include and library -METIS_DIR = ../lib/metis-5.1.0 -LIB_METIS = -L$(METIS_DIR)/lib -lmetis -METIS_INC = -I$(METIS_DIR)/include - -# Definition of pARMS include and library -PARMS_DIR = ../lib/parms -PARMS_INC = -I$(PARMS_DIR)/include -I$(PARMS_DIR)/src/include -LIB_PARMS = -L$(PARMS_DIR)/lib -lparms -CPP_SOL = -DPARMS - -###### Objects for Mesh Partitioning ################################################ -OBJ_INI = fort_part.o oce_modules.o MOD_READ_BINARY_ARRAYS.o MOD_WRITE_BINARY_ARRAYS.o MOD_MESH.o MOD_PARTIT.o \ - gen_modules_partitioning.o gen_modules_config.o \ - gen_modules_rotate_grid.o oce_local.o gen_comm.o fvom_init.o - -# objects - -OBJECTS = fortran_utils.o \ - oce_modules.o \ - info_module.o \ - command_line_options.o \ - MOD_READ_BINARY_ARRAYS.o \ - MOD_WRITE_BINARY_ARRAYS.o \ - MOD_MESH.o \ - MOD_PARTIT.o \ - MOD_DYN.o \ - gen_modules_partitioning.o \ - MOD_ICE.o \ - MOD_TRACER.o \ - ice_modules.o \ - gen_modules_config.o \ - gen_modules_clock.o \ - gen_modules_rotate_grid.o \ - gen_modules_read_NetCDF.o \ - gen_modules_forcing.o \ - gen_halo_exchange.o \ - gen_support.o \ - oce_ale_mixing_kpp.o \ - oce_adv_tra_hor.o \ - oce_adv_tra_ver.o \ - oce_adv_tra_fct.o \ - oce_adv_tra_driver.o \ - gen_modules_diag.o \ - psolve.o \ - oce_ale_mixing_pp.o \ - oce_tracer_mod.o \ - cvmix_kinds_and_types.o \ - cvmix_utils.o \ - cvmix_put_get.o \ - cvmix_tke.o \ - cvmix_idemix.o \ - gen_modules_cvmix_idemix.o \ - gen_modules_cvmix_tke.o \ - cvmix_math.o \ - cvmix_kpp.o \ - gen_modules_cvmix_kpp.o \ - cvmix_tidal.o \ - gen_modules_cvmix_tidal.o \ - cvmix_shear.o \ - gen_modules_cvmix_pp.o \ - async_threads_module.o \ - forcing_provider_netcdf_module.o \ - forcing_lookahead_reader_module.o \ - forcing_provider_async_module.o \ - io_gather.o \ - mpi_topology_module.o \ - io_netcdf_workaround_module.o \ - io_data_strategy.o \ - fesom_version_info.o \ - io_netcdf_attribute_module.o \ - io_netcdf_file_module.o \ - io_scatter.o \ - io_meandata.o \ - io_restart_derivedtype.o \ - io_fesom_file.o \ - io_restart_file_group.o \ - io_restart.o \ - io_blowup.o \ - io_mesh_info.o \ - oce_ale_pressure_bv.o \ - gen_ic3d.o \ - gen_surface_forcing.o \ - gen_modules_gpot.o \ - toy_channel_soufflet.o \ - gen_modules_backscatter.o \ - solver.o \ - oce_ale_vel_rhs.o \ - write_step_info.o \ - oce_fer_gm.o \ - oce_ale_tracer.o \ - oce_ale.o \ - oce_setup_step.o \ - gen_comm.o \ - oce_mesh.o \ - oce_dyn.o \ - oce_muscl_adv.o \ - oce_shortwave_pene.o \ - cavity_param.o \ - ice_EVP.o \ - ice_maEVP.o \ - ice_fct.o \ - ice_thermo_oce.o \ - ice_setup_step.o \ - ice_oce_coupling.o \ - gen_model_setup.o \ - gen_forcing_init.o \ - gen_bulk_formulae.o \ - gen_forcing_couple.o \ - gen_interpolation.o \ - gen_events.o \ - oce_mo_conv.o \ - oce_spp.o \ - fesom_module.o \ - fesom_main.o - - -# oce_pressure_bv.o \ -# oce_redi_gm.o \ -# Name of executables - -EXE = fesom.x -EXE_INI = fesom_ini.x - -# TARGETS - -default: run - -run: $(MODULES) $(OBJECTS) - @echo "======= Building FESOM ==========" - $(LD) $(OPT) -o $(EXE) $(OBJECTS) \ - $(MPI_LIB) $(LIB_LAP) $(NC_LIB) $(LIB_PARMS) -# cp -pf $(EXE) ../bin/. - -run_ini: CPP_DEFS+=-DFVOM_INIT -run_ini: cleanomod $(OBJ_INI) - @echo "======= Building FESOM partioning program ==========" - $(LD) $(OPT) -o $(EXE_INI) $(OBJ_INI) \ - $(MPI_LIB) $(LIB_METIS) $(NC_LIB) - cp -pf $(EXE_INI) ../bin/. - -.SUFFIXES: .c .F90 .o - -.c.o : - $(MPICC) $(COPT) $(CPP_SOL) $(METIS_DEFS) $(PARMS_DEFS) $(CPP_DEFS) $(METIS_INC) $(PARMS_INC) $(MPI_INC) -c $*.c - -.F90.o : - @echo $@ - $(MPIFC) $(CPP_DEFS) $(CPP_SOL) $(NC_INC) $(OPT) $(OASIS3_INC) $(PARMS_INC) $(MPI_INC) -c $*.F90 - -clean : - rm -f *.o *.mod *~ *.f90 $(EXE) $(EXE_INI) - -cleanomod: - rm -f *.o *.mod - CPP_DEFS=; export CPP_DEFS diff --git a/src/Makefile_hlrn b/src/Makefile_hlrn deleted file mode 100755 index 1f893b580..000000000 --- a/src/Makefile_hlrn +++ /dev/null @@ -1,104 +0,0 @@ -###################################################### -# Makefile -###################################################### - -###################################################### -# Include architecture-specific definitions - -include Makefile.in_hlrn - -###### Objects for Mesh Partitioning ################################################ -# modules -MOD_INI = fort_part.o oce_modules.o gen_modules_config.o gen_modules_rotate_grid.o - -OBJ_INI = fvom_init.o \ - oce_local.o \ - gen_comm.o - -# objects - -MODULES = oce_modules.o \ - ice_modules.o \ - gen_modules_config.o \ - gen_modules_clock.o \ - gen_modules_rotate_grid.o \ - gen_modules_read_NetCDF.o \ - gen_modules_forcing.o \ - gen_halo_exchange.o \ - psolve.o \ - oce_tracer_mod.o \ - gen_str2unstr.o \ - fort_part.o \ - gen_input.o - -OBJECTS= fvom_main.o \ - oce_setup_step.o \ - oce_mesh.o \ - oce_dyn.o \ - oce_vel_rhs.o \ - oce_vel_rhs_vinv.o \ - oce_pressure_bv.o \ - oce_fer_gm.o \ - oce_muscl_adv.o \ - oce_fct_muscl_adv.o \ - oce_mixing.o \ - oce_ice_init_state.o \ - oce_shortwave_pene.o \ - ice_setup_step.o \ - ice_evp.o \ - ice_fct.o \ - ice_oce_coupling.o \ - ice_thermodynamics.o \ - gen_comm.o \ - gen_model_setup.o \ - gen_forcing_init.o \ - gen_bulk_formulae.o \ - gen_forcing_couple.o \ - gen_forcing_ocean.o \ - gen_interpolation.o \ - gen_output_mean.o \ - gen_output_restart.o - - - -# oce_pressure_bv.o \ -# oce_redi_gm.o \ -# Name of executables - -EXE = fvom.x -EXE_INI = fvom_ini.x - -# TARGETS - -default: run - -run: $(MODULES) $(OBJECTS) - @echo "======= Building FESOM ==========" - $(LD) $(OPT) -o $(EXE) $(FOS_INC) $(MODULES) $(OBJECTS) \ - $(MPI_LIB) $(LIB_LAP) $(LIB_PARMS) $(LIB_METIS) $(OASIS3_LIB) $(NCLIB) - cp -pf $(EXE) ../bin/. - -run_ini: CPP_DEFS+=-DFVOM_INIT -run_ini: cleanomod $(MOD_INI) $(OBJ_INI) - @echo "======= Building FESOM paritioning program ==========" - $(LD) $(OPT) -o $(EXE_INI) $(MOD_INI) $(OBJ_INI) \ - $(MPI_LIB) $(LIB_METIS) $(OASIS3_LIB) $(NCLIB) - cp -pf $(EXE_INI) ../bin/. - -.SUFFIXES: .c .F90 .o - -.c.o : - $(CC) $(COPT) $(METIS_DEFS) $(CPP_DEFS) $(PARMS_DEFS) $(METIS_INC) $(PARMS_INC) $(MPI_INCLUDE) -c $*.c - -.F90.o : - @echo $@ - $(CPP) $(CPP_DEFS) $(CPP_SOL) $(PETSC_INC) $(PARMS_INC) $(PETSCCONF_INC) $(FOS_INC) $(MPI_INC) $(NCINC) $(LIBS_SLV) < $*.F90 > $*.f90 - $(FC) $(OPT) $(CPP_SOL) $(CPP_DEFS) $(OASIS3_INC) $(FOS_INC)\ - $(PARMS_INC) $(MPI_INCLUDE) $(NCINC) -c $*.f90 - -clean : - rm -f *.o *.mod *~ *.f90 fvom.x fvom_ini.x - -cleanomod: - rm -f *.o *.mod - CPP_DEFS=; export CPP_DEFS diff --git a/src/Makefile_ollie b/src/Makefile_ollie deleted file mode 100755 index bd8eb9c16..000000000 --- a/src/Makefile_ollie +++ /dev/null @@ -1,117 +0,0 @@ -###################################################### -# Makefile -###################################################### - -###################################################### -# Include architecture-specific definitions - -include Makefile.in - -###### Objects for Mesh Partitioning ################################################ -# modules -MOD_INI = fort_part.o oce_modules.o gen_modules_config.o gen_modules_partitioning.o gen_modules_rotate_grid.o - -OBJ_INI = fvom_init.o \ - oce_local.o \ - gen_comm.o - -# objects - -MODULES = oce_modules.o \ - ice_modules.o \ - gen_modules_config.o \ - gen_modules_partitioning.o \ - gen_modules_clock.o \ - gen_modules_rotate_grid.o \ - gen_modules_read_NetCDF.o \ - gen_modules_forcing.o \ - gen_halo_exchange.o \ - gen_support.o \ - oce_ale_mixing_kpp.o \ - gen_modules_diag.o \ - psolve.o \ - oce_ale_mixing_pp.o \ - oce_tracer_mod.o \ - fort_part.o \ - gen_input.o \ - io_meandata.o \ - io_restart.o \ - io_blowup.o \ - io_mesh_info.o - -OBJECTS= fvom_main.o \ - gen_comm.o \ - oce_setup_step.o \ - oce_mesh.o \ - oce_dyn.o \ - oce_ale_vel_rhs.o \ - oce_vel_rhs_vinv.o \ - oce_ale_pressure_bv.o \ - oce_fer_gm.o \ - oce_muscl_adv.o \ - oce_ale_fct_3d_adv.o \ - oce_ice_init_state.o \ - oce_shortwave_pene.o \ - oce_ale.o \ - oce_ale_tracer.o \ - ice_evp.o \ - ice_maEVP.o \ - ice_setup_step.o \ - ice_fct.o \ - ice_oce_coupling.o \ - ice_thermo_oce.o \ - gen_model_setup.o \ - gen_forcing_init.o \ - gen_bulk_formulae.o \ - gen_forcing_couple.o \ - gen_interpolation.o \ - gen_events.o \ - write_step_info.o - - - -# oce_pressure_bv.o \ -# oce_redi_gm.o \ -# Name of executables - -EXE = fesom.x -EXE_INI = fvom_ini.x - -# TARGETS - -default: run - -run: $(MODULES) $(OBJECTS) - @echo "======= Building FESOM ==========" - $(LD) $(OPT) -o $(EXE) $(FOS_INC) $(MODULES) $(OBJECTS) \ - $(MPI_LIB) $(LIB_LAP) $(LIB_PARMS) $(LIB_METIS) $(NCLIB) - cp -pf $(EXE) ../bin/. - -run_ini: CPP_DEFS+=-DFVOM_INIT -run_ini: cleanomod $(MOD_INI) $(OBJ_INI) - @echo "======= Building FESOM paritioning program ==========" - $(LD) $(OPT) -o $(EXE_INI) $(MOD_INI) $(OBJ_INI) \ - $(MPI_LIB) $(LIB_METIS) $(NCLIB) - cp -pf $(EXE_INI) ../bin/. - -.SUFFIXES: .c .F90 .o - -.c.o : - $(CC) $(COPT) $(METIS_DEFS) $(CPP_DEFS) $(PARMS_DEFS) $(METIS_INC) $(PARMS_INC) $(MPI_INCLUDE) -c $*.c - -#.F90.o : -# @echo $@ -# $(CPP) $(CPP_DEFS) $(CPP_SOL) $(PETSC_INC) $(PARMS_INC) $(PETSCCONF_INC) $(FOS_INC) $(MPI_INC) $(NCINC) $(LIBS_SLV) < $*.F90 > $*.f90 -# $(FC) $(OPT) $(CPP_SOL) $(CPP_DEFS) $(OASIS3_INC) $(FOS_INC)\ -# $(PARMS_INC) $(MPI_INCLUDE) $(NCINC) -c $*.f90 -.F90.o : - @echo $@ - $(FC) $(CPP_DEFS) $(CPP_SOL) $(PETSC_INC) $(PARMS_INC) $(PETSCCONF_INC) $(FOS_INC) $(MPI_INC) $(NCINC) $(LIBS_SLV) $(OPT) $(CPP_SOL) $(CPP_DEFS) $(OASIS3_INC) $(FOS_INC)\ - $(PARMS_INC) $(MPI_INCLUDE) $(NCINC) -c $*.F90 - -clean : - rm -f *.o *.mod *~ *.f90 fvom.x fvom_ini.x - -cleanomod: - rm -f *.o *.mod - CPP_DEFS=; export CPP_DEFS diff --git a/src/temp/MOD_MESH.F90 b/src/temp/MOD_MESH.F90 deleted file mode 100644 index 4eb0c23e1..000000000 --- a/src/temp/MOD_MESH.F90 +++ /dev/null @@ -1,329 +0,0 @@ -!========================================================== -MODULE MOD_MESH -USE O_PARAM -USE MOD_WRITE_BINARY_ARRAYS -USE MOD_READ_BINARY_ARRAYS -USE, intrinsic :: ISO_FORTRAN_ENV -IMPLICIT NONE -SAVE -integer, parameter :: MAX_ADJACENT=32 ! Max allowed number of adjacent nodes - -TYPE SPARSE_MATRIX - integer :: nza - integer :: dim - real(kind=WP), allocatable, dimension(:) :: values - integer(int32), allocatable, dimension(:) :: colind - integer(int32), allocatable, dimension(:) :: rowptr - integer(int32), allocatable, dimension(:) :: colind_loc - integer(int32), allocatable, dimension(:) :: rowptr_loc -END TYPE SPARSE_MATRIX - -TYPE T_MESH -integer :: nod2D ! the number of 2D nodes -real(kind=WP) :: ocean_area, ocean_areawithcav -real(kind=WP), allocatable, dimension(:,:) :: coord_nod2D, geo_coord_nod2D -integer :: edge2D ! the number of 2D edges -integer :: edge2D_in ! the number of internal 2D edges -integer :: elem2D ! the number of 2D elements -integer, allocatable, dimension(:,:) :: elem2D_nodes ! elem2D_nodes(:,n) lists; 3 nodes of element n -integer, allocatable, dimension(:,:) :: edges ! edge(:,n) lists 2 nodes; edge n -integer, allocatable, dimension(:,:) :: edge_tri ! edge_tri(:,n) lists 2 - ! elements containing edge n: the first one is to left - ! of the line directed to the second node -integer, allocatable, dimension(:,:) :: elem_edges ! elem_edges(:,n) are edges of element n. -real(kind=WP), allocatable, dimension(:) :: elem_area -real(kind=WP), allocatable, dimension(:,:) :: edge_dxdy, edge_cross_dxdy -real(kind=WP), allocatable, dimension(:) :: elem_cos, metric_factor -integer, allocatable, dimension(:,:) :: elem_neighbors -integer, allocatable, dimension(:,:) :: nod_in_elem2D -real(kind=WP), allocatable, dimension(:,:) :: x_corners, y_corners ! cornes for the scalar points -integer, allocatable, dimension(:) :: nod_in_elem2D_num -real(kind=WP), allocatable, dimension(:) :: depth ! depth(n) is the depths at node n -real(kind=WP), allocatable, dimension(:,:) :: gradient_vec - ! coefficients of linear reconstruction - ! of velocities on elements -real(kind=WP), allocatable, dimension(:,:) :: gradient_sca ! Coefficients to compute gradient of scalars - ! on elements -INTEGER, ALLOCATABLE, DIMENSION(:) :: bc_index_nod2D(:) - ! vertical structure -! -! -!___vertical mesh info__________________________________________________________ -! total number of layers -integer :: nl - -! initial layer, mid-depth layer and element depth -real(kind=WP), allocatable, dimension(:) :: zbar, Z,elem_depth - -! upper boudnary index of all vertical vertice/element loops, default==1 but when -! cavity is used becomes index of cavity-ocean boundary at vertices and elements -integer, allocatable, dimension(:) :: ulevels, ulevels_nod2D, ulevels_nod2D_max - -! number of levels at elem and vertices considering bottom topography -integer, allocatable, dimension(:) :: nlevels, nlevels_nod2D, nlevels_nod2D_min - -! -! -!___horizontal mesh info________________________________________________________ -real(kind=WP), allocatable, dimension(:,:) :: area, area_inv, areasvol, areasvol_inv -real(kind=WP), allocatable, dimension(:) :: mesh_resolution - -! -! -!___cavity mesh info____________________________________________________________ -! level index of cavity-ocean boundary at vertices and elements -! --> see: ulevels, ulevels_nod2D (fvom_main) - -! vertice/element yes=1/no=0 flag if cavity exists -integer, allocatable, dimension(:) :: cavity_flag_n, cavity_flag_e - -! depth of cavity-ocean interface -real(kind=WP), allocatable, dimension(:) :: cavity_depth - - -real(kind=WP), allocatable, dimension(:,:) :: cavity_nrst_cavlpnt_xyz - -! -! -!___Elevation stiffness matrix__________________________________________________ -type(sparse_matrix) :: ssh_stiff - -!#if defined (__oasis) -real(kind=WP), allocatable, dimension(:) :: lump2d_south, lump2d_north -integer, allocatable, dimension(:) :: ind_south, ind_north -!#endif - -integer :: nn_size -integer, allocatable, dimension(:) :: nn_num -integer, allocatable, dimension(:,:) :: nn_pos - -!_______________________________________________________________________________ -! Arrays added for ALE implementation: -! --> layer thinkness at node and depthlayer for t=n and t=n+1 -real(kind=WP), allocatable,dimension(:,:) :: hnode, hnode_new, zbar_3d_n, Z_3d_n - -! --> layer thinkness at elements, interpolated from hnode -real(kind=WP), allocatable,dimension(:,:) :: helem - -! --> thinkness of bottom elem (important for partial cells) -real(kind=WP), allocatable,dimension(:) :: bottom_elem_thickness -real(kind=WP), allocatable,dimension(:) :: bottom_node_thickness - -! --> The increment of total fluid depth on elements. It is used to update the matrix -real(kind=WP), allocatable,dimension(:) :: dhe - -! --> hbar, hbar_old: correspond to the elevation, but on semi-integer time steps. -real(kind=WP), allocatable,dimension(:) :: hbar, hbar_old - -! --> auxiliary array to store depth of layers and depth of mid level due to changing -! layer thinkness at every node -real(kind=WP), allocatable,dimension(:) :: zbar_n, Z_n - -! new bottom depth at node and element due to partial cells -real(kind=WP), allocatable,dimension(:) :: zbar_n_bot -real(kind=WP), allocatable,dimension(:) :: zbar_e_bot - -! new depth of cavity-ocean interface at node and element due to partial cells -real(kind=WP), allocatable,dimension(:) :: zbar_n_srf -real(kind=WP), allocatable,dimension(:) :: zbar_e_srf - -character(:), allocatable :: representative_checksum - -contains - procedure write_t_mesh - procedure read_t_mesh - generic :: write(unformatted) => write_t_mesh - generic :: read(unformatted) => read_t_mesh -END TYPE T_MESH - -contains - -! Unformatted writing for t_mesh -subroutine write_t_mesh(mesh, unit, iostat, iomsg) - IMPLICIT NONE - class(t_mesh), intent(in) :: mesh - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: i, j, k - integer :: s1, s2, s3 - ! write records (giving sizes for the allocation for arrays) - write(unit, iostat=iostat, iomsg=iomsg) mesh%nod2D - write(unit, iostat=iostat, iomsg=iomsg) mesh%ocean_area - write(unit, iostat=iostat, iomsg=iomsg) mesh%ocean_areawithcav - write(unit, iostat=iostat, iomsg=iomsg) mesh%edge2D - write(unit, iostat=iostat, iomsg=iomsg) mesh%edge2D_in - write(unit, iostat=iostat, iomsg=iomsg) mesh%elem2D - call write_bin_array(mesh%elem2D_nodes, unit, iostat, iomsg) - call write_bin_array(mesh%edges, unit, iostat, iomsg) - call write_bin_array(mesh%edge_tri, unit, iostat, iomsg) - call write_bin_array(mesh%elem_edges, unit, iostat, iomsg) - call write_bin_array(mesh%elem_area, unit, iostat, iomsg) - call write_bin_array(mesh%edge_dxdy, unit, iostat, iomsg) - - call write_bin_array(mesh%edge_cross_dxdy, unit, iostat, iomsg) - call write_bin_array(mesh%elem_cos, unit, iostat, iomsg) - call write_bin_array(mesh%metric_factor, unit, iostat, iomsg) - call write_bin_array(mesh%elem_neighbors, unit, iostat, iomsg) - call write_bin_array(mesh%nod_in_elem2D, unit, iostat, iomsg) - call write_bin_array(mesh%x_corners, unit, iostat, iomsg) - call write_bin_array(mesh%y_corners, unit, iostat, iomsg) - call write_bin_array(mesh%nod_in_elem2D_num, unit, iostat, iomsg) - call write_bin_array(mesh%depth, unit, iostat, iomsg) - call write_bin_array(mesh%gradient_vec, unit, iostat, iomsg) - call write_bin_array(mesh%gradient_sca, unit, iostat, iomsg) - call write_bin_array(mesh%bc_index_nod2D, unit, iostat, iomsg) - - write(unit, iostat=iostat, iomsg=iomsg) mesh%nl - - call write_bin_array(mesh%zbar, unit, iostat, iomsg) - call write_bin_array(mesh%Z, unit, iostat, iomsg) - call write_bin_array(mesh%elem_depth, unit, iostat, iomsg) - call write_bin_array(mesh%ulevels, unit, iostat, iomsg) - call write_bin_array(mesh%ulevels_nod2D, unit, iostat, iomsg) - call write_bin_array(mesh%ulevels_nod2D_max, unit, iostat, iomsg) - call write_bin_array(mesh%nlevels, unit, iostat, iomsg) - call write_bin_array(mesh%nlevels_nod2D, unit, iostat, iomsg) - call write_bin_array(mesh%nlevels_nod2D_min, unit, iostat, iomsg) - call write_bin_array(mesh%area, unit, iostat, iomsg) - call write_bin_array(mesh%area_inv, unit, iostat, iomsg) - call write_bin_array(mesh%areasvol, unit, iostat, iomsg) - call write_bin_array(mesh%areasvol_inv, unit, iostat, iomsg) - call write_bin_array(mesh%mesh_resolution, unit, iostat, iomsg) - - call write_bin_array(mesh%cavity_flag_n, unit, iostat, iomsg) - call write_bin_array(mesh%cavity_flag_e, unit, iostat, iomsg) - call write_bin_array(mesh%cavity_depth, unit, iostat, iomsg) - call write_bin_array(mesh%cavity_nrst_cavlpnt_xyz, unit, iostat, iomsg) - - write(unit, iostat=iostat, iomsg=iomsg) mesh%ssh_stiff%dim - write(unit, iostat=iostat, iomsg=iomsg) mesh%ssh_stiff%nza - - call write_bin_array(mesh%ssh_stiff%rowptr, unit, iostat, iomsg) - call write_bin_array(mesh%ssh_stiff%colind, unit, iostat, iomsg) - call write_bin_array(mesh%ssh_stiff%values, unit, iostat, iomsg) - call write_bin_array(mesh%ssh_stiff%colind_loc, unit, iostat, iomsg) - call write_bin_array(mesh%ssh_stiff%rowptr_loc, unit, iostat, iomsg) - - call write_bin_array(mesh%lump2d_south, unit, iostat, iomsg) - call write_bin_array(mesh%lump2d_north, unit, iostat, iomsg) - call write_bin_array(mesh%ind_south, unit, iostat, iomsg) - call write_bin_array(mesh%ind_north, unit, iostat, iomsg) - write(unit, iostat=iostat, iomsg=iomsg) mesh%nn_size - call write_bin_array(mesh%nn_num, unit, iostat, iomsg) - call write_bin_array(mesh%nn_pos, unit, iostat, iomsg) - call write_bin_array(mesh%hnode, unit, iostat, iomsg) - call write_bin_array(mesh%hnode_new, unit, iostat, iomsg) - call write_bin_array(mesh%zbar_3d_n, unit, iostat, iomsg) - call write_bin_array(mesh%Z_3d_n, unit, iostat, iomsg) - call write_bin_array(mesh%helem, unit, iostat, iomsg) - call write_bin_array(mesh%bottom_elem_thickness, unit, iostat, iomsg) - call write_bin_array(mesh%bottom_node_thickness, unit, iostat, iomsg) - call write_bin_array(mesh%dhe, unit, iostat, iomsg) - call write_bin_array(mesh%hbar, unit, iostat, iomsg) - call write_bin_array(mesh%hbar_old, unit, iostat, iomsg) - call write_bin_array(mesh%zbar_n, unit, iostat, iomsg) - call write_bin_array(mesh%Z_n, unit, iostat, iomsg) - call write_bin_array(mesh%zbar_n_bot, unit, iostat, iomsg) - call write_bin_array(mesh%zbar_e_bot, unit, iostat, iomsg) - call write_bin_array(mesh%zbar_n_srf, unit, iostat, iomsg) - call write_bin_array(mesh%zbar_e_srf, unit, iostat, iomsg) -! call write_bin_array(mesh%representative_checksum, unit, iostat, iomsg) -end subroutine write_t_mesh - -! Unformatted reading for t_mesh -subroutine read_t_mesh(mesh, unit, iostat, iomsg) - IMPLICIT NONE - class(t_mesh), intent(inout) :: mesh - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: i, j, k - integer :: s1, s2, s3 - ! write records (giving sizes for the allocation for arrays) - read(unit, iostat=iostat, iomsg=iomsg) mesh%nod2D - read(unit, iostat=iostat, iomsg=iomsg) mesh%ocean_area - read(unit, iostat=iostat, iomsg=iomsg) mesh%ocean_areawithcav - read(unit, iostat=iostat, iomsg=iomsg) mesh%edge2D - read(unit, iostat=iostat, iomsg=iomsg) mesh%edge2D_in - read(unit, iostat=iostat, iomsg=iomsg) mesh%elem2D - - call read_bin_array(mesh%elem2D_nodes, unit, iostat, iomsg) - call read_bin_array(mesh%edges, unit, iostat, iomsg) - call read_bin_array(mesh%edge_tri, unit, iostat, iomsg) - call read_bin_array(mesh%elem_edges, unit, iostat, iomsg) - call read_bin_array(mesh%elem_area, unit, iostat, iomsg) - call read_bin_array(mesh%edge_dxdy, unit, iostat, iomsg) - - call read_bin_array(mesh%edge_cross_dxdy, unit, iostat, iomsg) - call read_bin_array(mesh%elem_cos, unit, iostat, iomsg) - call read_bin_array(mesh%metric_factor, unit, iostat, iomsg) - call read_bin_array(mesh%elem_neighbors, unit, iostat, iomsg) - call read_bin_array(mesh%nod_in_elem2D, unit, iostat, iomsg) - call read_bin_array(mesh%x_corners, unit, iostat, iomsg) - call read_bin_array(mesh%y_corners, unit, iostat, iomsg) - call read_bin_array(mesh%nod_in_elem2D_num, unit, iostat, iomsg) - call read_bin_array(mesh%depth, unit, iostat, iomsg) - call read_bin_array(mesh%gradient_vec, unit, iostat, iomsg) - call read_bin_array(mesh%gradient_sca, unit, iostat, iomsg) - call read_bin_array(mesh%bc_index_nod2D, unit, iostat, iomsg) - - read(unit, iostat=iostat, iomsg=iomsg) mesh%nl - - call read_bin_array(mesh%zbar, unit, iostat, iomsg) - call read_bin_array(mesh%Z, unit, iostat, iomsg) - call read_bin_array(mesh%elem_depth, unit, iostat, iomsg) - call read_bin_array(mesh%ulevels, unit, iostat, iomsg) - call read_bin_array(mesh%ulevels_nod2D, unit, iostat, iomsg) - call read_bin_array(mesh%ulevels_nod2D_max, unit, iostat, iomsg) - call read_bin_array(mesh%nlevels, unit, iostat, iomsg) - call read_bin_array(mesh%nlevels_nod2D, unit, iostat, iomsg) - call read_bin_array(mesh%nlevels_nod2D_min, unit, iostat, iomsg) - call read_bin_array(mesh%area, unit, iostat, iomsg) - call read_bin_array(mesh%area_inv, unit, iostat, iomsg) - call read_bin_array(mesh%areasvol, unit, iostat, iomsg) - call read_bin_array(mesh%areasvol_inv, unit, iostat, iomsg) - call read_bin_array(mesh%mesh_resolution, unit, iostat, iomsg) - - call read_bin_array(mesh%cavity_flag_n, unit, iostat, iomsg) - call read_bin_array(mesh%cavity_flag_e, unit, iostat, iomsg) - call read_bin_array(mesh%cavity_depth, unit, iostat, iomsg) - call read_bin_array(mesh%cavity_nrst_cavlpnt_xyz, unit, iostat, iomsg) - - read(unit, iostat=iostat, iomsg=iomsg) mesh%ssh_stiff%dim - read(unit, iostat=iostat, iomsg=iomsg) mesh%ssh_stiff%nza - - call read_bin_array(mesh%ssh_stiff%rowptr, unit, iostat, iomsg) - call read_bin_array(mesh%ssh_stiff%colind, unit, iostat, iomsg) - call read_bin_array(mesh%ssh_stiff%values, unit, iostat, iomsg) - call read_bin_array(mesh%ssh_stiff%colind_loc, unit, iostat, iomsg) - call read_bin_array(mesh%ssh_stiff%rowptr_loc, unit, iostat, iomsg) - - call read_bin_array(mesh%lump2d_south, unit, iostat, iomsg) - call read_bin_array(mesh%lump2d_north, unit, iostat, iomsg) - call read_bin_array(mesh%ind_south, unit, iostat, iomsg) - call read_bin_array(mesh%ind_north, unit, iostat, iomsg) - read(unit, iostat=iostat, iomsg=iomsg) mesh%nn_size - call read_bin_array(mesh%nn_num, unit, iostat, iomsg) - call read_bin_array(mesh%nn_pos, unit, iostat, iomsg) - call read_bin_array(mesh%hnode, unit, iostat, iomsg) - call read_bin_array(mesh%hnode_new, unit, iostat, iomsg) - call read_bin_array(mesh%zbar_3d_n, unit, iostat, iomsg) - call read_bin_array(mesh%Z_3d_n, unit, iostat, iomsg) - call read_bin_array(mesh%helem, unit, iostat, iomsg) - call read_bin_array(mesh%bottom_elem_thickness, unit, iostat, iomsg) - call read_bin_array(mesh%bottom_node_thickness, unit, iostat, iomsg) - call read_bin_array(mesh%dhe, unit, iostat, iomsg) - call read_bin_array(mesh%hbar, unit, iostat, iomsg) - call read_bin_array(mesh%hbar_old, unit, iostat, iomsg) - call read_bin_array(mesh%zbar_n, unit, iostat, iomsg) - call read_bin_array(mesh%Z_n, unit, iostat, iomsg) - call read_bin_array(mesh%zbar_n_bot, unit, iostat, iomsg) - call read_bin_array(mesh%zbar_e_bot, unit, iostat, iomsg) - call read_bin_array(mesh%zbar_n_srf, unit, iostat, iomsg) - call read_bin_array(mesh%zbar_e_srf, unit, iostat, iomsg) -! call read_bin_array(mesh%representative_checksum, unit, iostat, iomsg) -end subroutine read_t_mesh -end module MOD_MESH -!========================================================== - diff --git a/src/temp/MOD_PARTIT.F90 b/src/temp/MOD_PARTIT.F90 deleted file mode 100644 index 818b46541..000000000 --- a/src/temp/MOD_PARTIT.F90 +++ /dev/null @@ -1,189 +0,0 @@ -!========================================================== -! Variables to organize parallel work -module MOD_PARTIT -USE O_PARAM -USE, intrinsic :: ISO_FORTRAN_ENV -USE MOD_WRITE_BINARY_ARRAYS -USE MOD_READ_BINARY_ARRAYS -USE mpi -IMPLICIT NONE -SAVE -integer, parameter :: MAX_LAENDERECK=16 -integer, parameter :: MAX_NEIGHBOR_PARTITIONS=32 - - -type com_struct - integer :: rPEnum ! the number of PE I receive info from - integer, dimension(MAX_NEIGHBOR_PARTITIONS) :: rPE ! their list - integer, dimension(MAX_NEIGHBOR_PARTITIONS+1) :: rptr ! allocatables to the list of nodes - integer, dimension(:), allocatable :: rlist ! the list of nodes - integer :: sPEnum ! send part - integer, dimension(MAX_NEIGHBOR_PARTITIONS) :: sPE - integer, dimension(MAX_NEIGHBOR_PARTITIONS) :: sptr - integer, dimension(:), allocatable :: slist - integer, dimension(:), allocatable :: req ! request for MPI_Wait - integer :: nreq ! number of requests for MPI_Wait - ! (to combine halo exchange of several fields) - contains - procedure WRITE_T_COM_STRUCT - procedure READ_T_COM_STRUCT - generic :: write(unformatted) => WRITE_T_COM_STRUCT - generic :: read(unformatted) => READ_T_COM_STRUCT -end type com_struct - -TYPE T_PARTIT - integer :: MPI_COMM_FESOM ! FESOM communicator (for ocean only runs if often a copy of MPI_COMM_WORLD) - - type(com_struct) :: com_nod2D - type(com_struct) :: com_elem2D - type(com_struct) :: com_elem2D_full - - ! MPI Datatypes for interface exchange - ! Element fields (2D; 2D integer; 3D with nl-1 or nl levels, 1 - 4 values) - ! small halo and / or full halo - !!! s(r)_mpitype_* are constructed during the runtime ans shall not be dumped!!! - integer, allocatable :: s_mpitype_elem2D(:,:), r_mpitype_elem2D(:,:) - integer, allocatable :: s_mpitype_elem2D_full_i(:), r_mpitype_elem2D_full_i(:) - integer, allocatable :: s_mpitype_elem2D_full(:,:), r_mpitype_elem2D_full(:,:) - integer, allocatable :: s_mpitype_elem3D(:,:,:), r_mpitype_elem3D(:,:,:) - integer, allocatable :: s_mpitype_elem3D_full(:,:,:),r_mpitype_elem3D_full(:,:,:) - - ! Nodal fields (2D; 2D integer; 3D with nl-1 or nl levels, one, two, or three values) - integer, allocatable :: s_mpitype_nod2D(:), r_mpitype_nod2D(:) - integer, allocatable :: s_mpitype_nod2D_i(:), r_mpitype_nod2D_i(:) - integer, allocatable :: s_mpitype_nod3D(:,:,:), r_mpitype_nod3D(:,:,:) - - ! general MPI part - integer :: MPIERR - integer :: npes - integer :: mype - integer :: maxPEnum=100 - integer, allocatable, dimension(:) :: part - - ! Mesh partition - integer :: myDim_nod2D, eDim_nod2D - integer, allocatable, dimension(:) :: myList_nod2D - integer :: myDim_elem2D, eDim_elem2D, eXDim_elem2D - integer, allocatable, dimension(:) :: myList_elem2D - integer :: myDim_edge2D, eDim_edge2D - integer, allocatable, dimension(:) :: myList_edge2D - - integer :: pe_status = 0 ! if /=0 then something is wrong - !!! remPtr_* are constructed during the runtime ans shall not be dumped!!! - integer, allocatable :: remPtr_nod2D(:), remList_nod2D(:) - integer, allocatable :: remPtr_elem2D(:), remList_elem2D(:) - - logical :: elem_full_flag - contains - procedure WRITE_T_PARTIT - procedure READ_T_PARTIT - generic :: write(unformatted) => WRITE_T_PARTIT - generic :: read(unformatted) => READ_T_PARTIT -END TYPE T_PARTIT -contains - -! Unformatted writing for COM_STRUCT TYPE -subroutine WRITE_T_COM_STRUCT(tstruct, unit, iostat, iomsg) - IMPLICIT NONE - class(COM_STRUCT), intent(in) :: tstruct - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - write(unit, iostat=iostat, iomsg=iomsg) tstruct%rPEnum - call write1d_int_static(tstruct%rPE, unit, iostat, iomsg) - call write1d_int_static(tstruct%rptr, unit, iostat, iomsg) - call write_bin_array(tstruct%rlist, unit, iostat, iomsg) - write(unit, iostat=iostat, iomsg=iomsg) tstruct%sPEnum - call write1d_int_static(tstruct%sPE, unit, iostat, iomsg) - call write1d_int_static(tstruct%sptr, unit, iostat, iomsg) - call write_bin_array(tstruct%slist, unit, iostat, iomsg) - ! req is constructed during the runtime - ! call write_bin_array(tstruct%req, unit, iostat, iomsg) - write(unit, iostat=iostat, iomsg=iomsg) tstruct%nreq -end subroutine WRITE_T_COM_STRUCT - -subroutine READ_T_COM_STRUCT(tstruct, unit, iostat, iomsg) - IMPLICIT NONE - class(COM_STRUCT), intent(inout) :: tstruct - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - read(unit, iostat=iostat, iomsg=iomsg) tstruct%rPEnum - call read1d_int_static(tstruct%rPE, unit, iostat, iomsg) - call read1d_int_static(tstruct%rptr, unit, iostat, iomsg) - call read_bin_array(tstruct%rlist, unit, iostat, iomsg) - read(unit, iostat=iostat, iomsg=iomsg) tstruct%sPEnum - call read1d_int_static(tstruct%sPE, unit, iostat, iomsg) - call read1d_int_static(tstruct%sptr, unit, iostat, iomsg) - call read_bin_array(tstruct%slist, unit, iostat, iomsg) -! req is constructed during the runtime -! call read_bin_array(tstruct%req, unit, iostat, iomsg) - read(unit, iostat=iostat, iomsg=iomsg) tstruct%nreq -end subroutine READ_T_COM_STRUCT - -! Unformatted writing for T_PARTIT -subroutine WRITE_T_PARTIT(partit, unit, iostat, iomsg) - IMPLICIT NONE - class(T_PARTIT), intent(in) :: partit - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - write(unit, iostat=iostat, iomsg=iomsg) partit%com_nod2D - write(unit, iostat=iostat, iomsg=iomsg) partit%com_elem2D - write(unit, iostat=iostat, iomsg=iomsg) partit%com_elem2D_full - - write(unit, iostat=iostat, iomsg=iomsg) partit%npes - write(unit, iostat=iostat, iomsg=iomsg) partit%mype - write(unit, iostat=iostat, iomsg=iomsg) partit%maxPEnum - call write_bin_array(partit%part, unit, iostat, iomsg) - - write(unit, iostat=iostat, iomsg=iomsg) partit%myDim_nod2D - write(unit, iostat=iostat, iomsg=iomsg) partit%eDim_nod2D - call write_bin_array(partit%myList_nod2D, unit, iostat, iomsg) - - write(unit, iostat=iostat, iomsg=iomsg) partit%myDim_elem2D - write(unit, iostat=iostat, iomsg=iomsg) partit%eDim_elem2D - write(unit, iostat=iostat, iomsg=iomsg) partit%eXDim_elem2D - call write_bin_array(partit%myList_elem2D, unit, iostat, iomsg) - - write(unit, iostat=iostat, iomsg=iomsg) partit%myDim_edge2D - write(unit, iostat=iostat, iomsg=iomsg) partit%eDim_edge2D - call write_bin_array(partit%myList_edge2D, unit, iostat, iomsg) - write(unit, iostat=iostat, iomsg=iomsg) partit%pe_status -end subroutine WRITE_T_PARTIT -! Unformatted reading for T_PARTIT -subroutine READ_T_PARTIT(partit, unit, iostat, iomsg) - IMPLICIT NONE - class(T_PARTIT), intent(inout) :: partit - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - read(unit, iostat=iostat, iomsg=iomsg) partit%com_nod2D - read(unit, iostat=iostat, iomsg=iomsg) partit%com_elem2D - read(unit, iostat=iostat, iomsg=iomsg) partit%com_elem2D_full - - read(unit, iostat=iostat, iomsg=iomsg) partit%npes - read(unit, iostat=iostat, iomsg=iomsg) partit%mype - read(unit, iostat=iostat, iomsg=iomsg) partit%maxPEnum - call read_bin_array(partit%part, unit, iostat, iomsg) - - read(unit, iostat=iostat, iomsg=iomsg) partit%myDim_nod2D - read(unit, iostat=iostat, iomsg=iomsg) partit%eDim_nod2D - call read_bin_array(partit%myList_nod2D, unit, iostat, iomsg) - - read(unit, iostat=iostat, iomsg=iomsg) partit%myDim_elem2D - read(unit, iostat=iostat, iomsg=iomsg) partit%eDim_elem2D - read(unit, iostat=iostat, iomsg=iomsg) partit%eXDim_elem2D - call read_bin_array(partit%myList_elem2D, unit, iostat, iomsg) - - read(unit, iostat=iostat, iomsg=iomsg) partit%myDim_edge2D - read(unit, iostat=iostat, iomsg=iomsg) partit%eDim_edge2D - call read_bin_array(partit%myList_edge2D, unit, iostat, iomsg) - read(unit, iostat=iostat, iomsg=iomsg) partit%pe_status -end subroutine READ_T_PARTIT - -end module MOD_PARTIT \ No newline at end of file diff --git a/src/temp/MOD_READ_BINARY_ARRAYS.F90 b/src/temp/MOD_READ_BINARY_ARRAYS.F90 deleted file mode 100644 index 87f0b2389..000000000 --- a/src/temp/MOD_READ_BINARY_ARRAYS.F90 +++ /dev/null @@ -1,118 +0,0 @@ -!========================================================== -! -!------------------------------------------------------------------------------------------ -! useful interface (read_bin_array) for reading arbitary binary arrays into an opened file -MODULE MOD_READ_BINARY_ARRAYS -use o_PARAM -private -public :: read_bin_array, read1d_int_static -INTERFACE read_bin_array - MODULE PROCEDURE read1d_real, read1d_int, read1d_char, read2d_real, read2d_int, read3d_real, read3d_int -END INTERFACE -contains -subroutine read1d_real(arr, unit, iostat, iomsg) - real(kind=WP), intent(inout), allocatable :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - read(unit, iostat=iostat, iomsg=iomsg) s1 - if (s1==0) return - allocate(arr(s1)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) -end subroutine read1d_real - -subroutine read1d_int(arr, unit, iostat, iomsg) - integer, intent(inout), allocatable :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - read(unit, iostat=iostat, iomsg=iomsg) s1 - if (s1==0) return - allocate(arr(s1)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) -end subroutine read1d_int - -subroutine read1d_char(arr, unit, iostat, iomsg) - character, intent(inout), allocatable :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - read(unit, iostat=iostat, iomsg=iomsg) s1 - if (s1==0) return - allocate(arr(s1)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) -end subroutine read1d_char - -subroutine read1d_int_static(arr, unit, iostat, iomsg) - IMPLICIT NONE - integer, intent(inout) :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - read(unit, iostat=iostat, iomsg=iomsg) s1 - if (s1==0) return - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) -end subroutine read1d_int_static - -subroutine read2d_real(arr, unit, iostat, iomsg) - real(kind=WP), intent(inout), allocatable :: arr(:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2 - - read(unit, iostat=iostat, iomsg=iomsg) s1, s2 - if ((s1==0) .or. (s2==0)) return - allocate(arr(s1, s2)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2) -end subroutine read2d_real - -subroutine read2d_int(arr, unit, iostat, iomsg) - integer, intent(inout), allocatable :: arr(:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2 - - read(unit, iostat=iostat, iomsg=iomsg) s1, s2 - if ((s1==0) .or. (s2==0)) return - allocate(arr(s1, s2)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2) -end subroutine read2d_int - -subroutine read3d_real(arr, unit, iostat, iomsg) - real(kind=WP), intent(inout), allocatable :: arr(:,:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2, s3 - - read(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3 - if ((s1==0) .or. (s2==0) .or. (s3==0)) return - allocate(arr(s1,s2,s3)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3) -end subroutine read3d_real - -subroutine read3d_int(arr, unit, iostat, iomsg) - integer, intent(inout), allocatable :: arr(:,:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2, s3 - - read(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3 - if ((s1==0) .or. (s2==0) .or. (s3==0)) return - allocate(arr(s1,s2,s3)) - read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3) -end subroutine read3d_int -end module MOD_READ_BINARY_ARRAYS -!========================================================== - diff --git a/src/temp/MOD_TRACER.F90 b/src/temp/MOD_TRACER.F90 deleted file mode 100644 index 8e8247830..000000000 --- a/src/temp/MOD_TRACER.F90 +++ /dev/null @@ -1,228 +0,0 @@ -!========================================================== -MODULE MOD_TRACER -USE O_PARAM -USE, intrinsic :: ISO_FORTRAN_ENV -USE MOD_WRITE_BINARY_ARRAYS -USE MOD_READ_BINARY_ARRAYS -IMPLICIT NONE -SAVE - -TYPE T_TRACER_DATA -real(kind=WP), allocatable, dimension(:,:) :: values, valuesAB ! instant values & Adams-Bashfort interpolation -logical :: smooth_bh_tra=.false. -real(kind=WP) :: gamma0_tra, gamma1_tra, gamma2_tra -logical :: i_vert_diff =.false. -character(20) :: tra_adv_hor, tra_adv_ver, tra_adv_lim ! type of the advection scheme for this tracer -real(kind=WP) :: tra_adv_ph = 1. ! a parameter to be used in horizontal advection (for MUSCL it is the fraction of fourth-order contribution in the solution) -real(kind=WP) :: tra_adv_pv = 1. ! a parameter to be used in horizontal advection (for QR4C it is the fraction of fourth-order contribution in the solution) -integer :: ID - -contains - procedure WRITE_T_TRACER_DATA - procedure READ_T_TRACER_DATA - generic :: write(unformatted) => WRITE_T_TRACER_DATA - generic :: read(unformatted) => READ_T_TRACER_DATA -END TYPE T_TRACER_DATA - - -TYPE T_TRACER_WORK -!auxuary arrays to work with tracers: -real(kind=WP), allocatable :: del_ttf(:,:) -real(kind=WP), allocatable :: del_ttf_advhoriz(:,:),del_ttf_advvert(:,:) -!_______________________________________________________________________________ -! in case ldiag_DVD=.true. --> calculate discrete variance decay (DVD) -real(kind=WP), allocatable :: tr_dvd_horiz(:,:,:), tr_dvd_vert(:,:,:) -! The fct part -real(kind=WP),allocatable,dimension(:,:) :: fct_LO ! Low-order solution -real(kind=WP),allocatable,dimension(:,:) :: adv_flux_hor ! Antidif. horiz. contrib. from edges / backup for iterafive fct scheme -real(kind=WP),allocatable,dimension(:,:) :: adv_flux_ver ! Antidif. vert. fluxes from nodes / backup for iterafive fct scheme - -real(kind=WP),allocatable,dimension(:,:) :: fct_ttf_max,fct_ttf_min -real(kind=WP),allocatable,dimension(:,:) :: fct_plus,fct_minus -! MUSCL type reconstruction -integer,allocatable,dimension(:) :: nboundary_lay -integer,allocatable,dimension(:,:) :: edge_up_dn_tri -real(kind=WP),allocatable,dimension(:,:,:) :: edge_up_dn_grad - -contains - procedure WRITE_T_TRACER_WORK - procedure READ_T_TRACER_WORK - generic :: write(unformatted) => WRITE_T_TRACER_WORK - generic :: read(unformatted) => READ_T_TRACER_WORK -END TYPE T_TRACER_WORK - -! auxury type for reading namelist.tra -TYPE NML_TRACER_LIST_TYPE - INTEGER :: ID =-1 - CHARACTER(len=4) :: adv_hor ='NONE' - CHARACTER(len=4) :: adv_ver ='NONE' - CHARACTER(len=4) :: adv_lim ='NONE' - REAL(kind=WP) :: adv_ph =1. - REAL(kind=WP) :: adv_pv =1. -END TYPE NML_TRACER_LIST_TYPE - -TYPE T_TRACER -! total number of tracers: -integer :: num_tracers=2 -type(t_tracer_data), allocatable :: data(:) -type(t_tracer_work) :: work -! general options for all tracers (can be moved to T_TRACER is needed) -! bharmonic diffusion for tracers. We recommend to use this option in very high resolution runs (Redi is generally off there). -logical :: smooth_bh_tra = .false. -real(kind=WP) :: gamma0_tra = 0.0005 -real(kind=WP) :: gamma1_tra = 0.0125 -real(kind=WP) :: gamma2_tra = 0. -logical :: i_vert_diff = .true. - -contains -procedure WRITE_T_TRACER -procedure READ_T_TRACER -generic :: write(unformatted) => WRITE_T_TRACER -generic :: read(unformatted) => READ_T_TRACER -END TYPE T_TRACER - -contains - -! Unformatted writing for T_TRACER_DATA -subroutine WRITE_T_TRACER_DATA(tdata, unit, iostat, iomsg) - IMPLICIT NONE - class(T_TRACER_DATA), intent(in) :: tdata - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - call write_bin_array(tdata%values, unit, iostat, iomsg) - call write_bin_array(tdata%valuesAB, unit, iostat, iomsg) - write(unit, iostat=iostat, iomsg=iomsg) tdata%smooth_bh_tra - write(unit, iostat=iostat, iomsg=iomsg) tdata%gamma0_tra - write(unit, iostat=iostat, iomsg=iomsg) tdata%gamma1_tra - write(unit, iostat=iostat, iomsg=iomsg) tdata%gamma2_tra - write(unit, iostat=iostat, iomsg=iomsg) tdata%i_vert_diff - write(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_hor - write(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_ver - write(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_lim - write(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_ph - write(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_pv - write(unit, iostat=iostat, iomsg=iomsg) tdata%ID -end subroutine WRITE_T_TRACER_DATA - -! Unformatted reading for T_TRACER_DATA -subroutine READ_T_TRACER_DATA(tdata, unit, iostat, iomsg) - IMPLICIT NONE - class(T_TRACER_DATA), intent(inout) :: tdata - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - call read_bin_array(tdata%values, unit, iostat, iomsg) - call read_bin_array(tdata%valuesAB, unit, iostat, iomsg) - read(unit, iostat=iostat, iomsg=iomsg) tdata%smooth_bh_tra - read(unit, iostat=iostat, iomsg=iomsg) tdata%gamma0_tra - read(unit, iostat=iostat, iomsg=iomsg) tdata%gamma1_tra - read(unit, iostat=iostat, iomsg=iomsg) tdata%gamma2_tra - read(unit, iostat=iostat, iomsg=iomsg) tdata%i_vert_diff - read(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_hor - read(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_ver - read(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_lim - read(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_ph - read(unit, iostat=iostat, iomsg=iomsg) tdata%tra_adv_pv - read(unit, iostat=iostat, iomsg=iomsg) tdata%ID -end subroutine READ_T_TRACER_DATA - -! Unformatted writing for T_TRACER_WORK -subroutine WRITE_T_TRACER_WORK(twork, unit, iostat, iomsg) - IMPLICIT NONE - class(T_TRACER_WORK), intent(in) :: twork - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - call write_bin_array(twork%del_ttf, unit, iostat, iomsg) - call write_bin_array(twork%del_ttf_advhoriz, unit, iostat, iomsg) - call write_bin_array(twork%del_ttf_advvert, unit, iostat, iomsg) - call write_bin_array(twork%tr_dvd_horiz, unit, iostat, iomsg) - call write_bin_array(twork%tr_dvd_vert, unit, iostat, iomsg) - call write_bin_array(twork%fct_LO, unit, iostat, iomsg) - call write_bin_array(twork%adv_flux_hor, unit, iostat, iomsg) - call write_bin_array(twork%adv_flux_ver, unit, iostat, iomsg) - call write_bin_array(twork%fct_ttf_max, unit, iostat, iomsg) - call write_bin_array(twork%fct_ttf_min, unit, iostat, iomsg) - call write_bin_array(twork%fct_plus, unit, iostat, iomsg) - call write_bin_array(twork%fct_minus, unit, iostat, iomsg) - call write_bin_array(twork%nboundary_lay, unit, iostat, iomsg) - call write_bin_array(twork%edge_up_dn_tri, unit, iostat, iomsg) - call write_bin_array(twork%edge_up_dn_grad, unit, iostat, iomsg) -end subroutine WRITE_T_TRACER_WORK - -! Unformatted reading for T_TRACER_WORK -subroutine READ_T_TRACER_WORK(twork, unit, iostat, iomsg) - IMPLICIT NONE - class(T_TRACER_WORK), intent(inout) :: twork - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - - call read_bin_array(twork%del_ttf, unit, iostat, iomsg) - call read_bin_array(twork%del_ttf_advhoriz, unit, iostat, iomsg) - call read_bin_array(twork%del_ttf_advvert, unit, iostat, iomsg) - call read_bin_array(twork%tr_dvd_horiz, unit, iostat, iomsg) - call read_bin_array(twork%tr_dvd_vert, unit, iostat, iomsg) - call read_bin_array(twork%fct_LO, unit, iostat, iomsg) - call read_bin_array(twork%adv_flux_hor, unit, iostat, iomsg) - call read_bin_array(twork%adv_flux_ver, unit, iostat, iomsg) - call read_bin_array(twork%fct_ttf_max, unit, iostat, iomsg) - call read_bin_array(twork%fct_ttf_min, unit, iostat, iomsg) - call read_bin_array(twork%fct_plus, unit, iostat, iomsg) - call read_bin_array(twork%fct_minus, unit, iostat, iomsg) - call read_bin_array(twork%nboundary_lay, unit, iostat, iomsg) - call read_bin_array(twork%edge_up_dn_tri, unit, iostat, iomsg) - call read_bin_array(twork%edge_up_dn_grad, unit, iostat, iomsg) -end subroutine READ_T_TRACER_WORK - -! Unformatted writing for T_TRACER -subroutine WRITE_T_TRACER(tracer, unit, iostat, iomsg) - IMPLICIT NONE - class(T_TRACER), intent(in) :: tracer - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: i - - write(unit, iostat=iostat, iomsg=iomsg) tracer%num_tracers - do i=1, tracer%num_tracers - write(unit, iostat=iostat, iomsg=iomsg) tracer%data(i) - end do - write(unit, iostat=iostat, iomsg=iomsg) tracer%work - write(unit, iostat=iostat, iomsg=iomsg) tracer%smooth_bh_tra - write(unit, iostat=iostat, iomsg=iomsg) tracer%gamma0_tra - write(unit, iostat=iostat, iomsg=iomsg) tracer%gamma1_tra - write(unit, iostat=iostat, iomsg=iomsg) tracer%gamma2_tra - write(unit, iostat=iostat, iomsg=iomsg) tracer%i_vert_diff -end subroutine WRITE_T_TRACER - -! Unformatted reading for T_TRACER -subroutine READ_T_TRACER(tracer, unit, iostat, iomsg) - IMPLICIT NONE - class(T_TRACER), intent(inout) :: tracer - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: i - - read(unit, iostat=iostat, iomsg=iomsg) tracer%num_tracers -! write(*,*) 'number of tracers to read: ', tracer%num_tracers - allocate(tracer%data(tracer%num_tracers)) - do i=1, tracer%num_tracers - read(unit, iostat=iostat, iomsg=iomsg) tracer%data(i) -! write(*,*) 'tracer info:', tracer%data(i)%ID, TRIM(tracer%data(i)%tra_adv_hor), TRIM(tracer%data(i)%tra_adv_ver), TRIM(tracer%data(i)%tra_adv_lim) - end do - read(unit, iostat=iostat, iomsg=iomsg) tracer%work - read(unit, iostat=iostat, iomsg=iomsg) tracer%smooth_bh_tra - read(unit, iostat=iostat, iomsg=iomsg) tracer%gamma0_tra - read(unit, iostat=iostat, iomsg=iomsg) tracer%gamma1_tra - read(unit, iostat=iostat, iomsg=iomsg) tracer%gamma2_tra - read(unit, iostat=iostat, iomsg=iomsg) tracer%i_vert_diff -end subroutine READ_T_TRACER -end module MOD_TRACER -!========================================================== - diff --git a/src/temp/MOD_WRITE_BINARY_ARRAYS.F90 b/src/temp/MOD_WRITE_BINARY_ARRAYS.F90 deleted file mode 100644 index 4f03b5cea..000000000 --- a/src/temp/MOD_WRITE_BINARY_ARRAYS.F90 +++ /dev/null @@ -1,160 +0,0 @@ -!========================================================== -! -!------------------------------------------------------------------------------------------ -! useful interface (write_bin_array) for writing arbitary binary arrays into an opened file -MODULE MOD_WRITE_BINARY_ARRAYS -use o_PARAM -private -public :: write_bin_array, write1d_int_static -INTERFACE write_bin_array - MODULE PROCEDURE write1d_real, write1d_int, write1d_char, write2d_real, write2d_int, write3d_real, write3d_int -END INTERFACE -contains - -subroutine write1d_real(arr, unit, iostat, iomsg) - real(kind=WP), intent(in), allocatable :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - if (allocated(arr)) then - s1=size(arr, 1) - write(unit, iostat=iostat, iomsg=iomsg) s1 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) - else - s1=0 - write(unit, iostat=iostat, iomsg=iomsg) s1 - end if -end subroutine write1d_real - -subroutine write1d_int(arr, unit, iostat, iomsg) - integer, intent(in), allocatable :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - if (allocated(arr)) then - s1=size(arr, 1) - write(unit, iostat=iostat, iomsg=iomsg) s1 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) - else - s1=0 - write(unit, iostat=iostat, iomsg=iomsg) s1 - end if -end subroutine write1d_int - -subroutine write1d_char(arr, unit, iostat, iomsg) - character, intent(in), allocatable :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - if (allocated(arr)) then - s1=size(arr, 1) - write(unit, iostat=iostat, iomsg=iomsg) s1 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) - else - s1=0 - write(unit, iostat=iostat, iomsg=iomsg) s1 - end if -end subroutine write1d_char - -subroutine write1d_int_static(arr, unit, iostat, iomsg) - IMPLICIT NONE - integer, intent(in) :: arr(:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1 - - s1=size(arr, 1) - write(unit, iostat=iostat, iomsg=iomsg) s1 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1) -end subroutine write1d_int_static - -subroutine write2d_real(arr, unit, iostat, iomsg) - real(kind=WP), intent(in), allocatable :: arr(:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2 - - if (allocated(arr)) then - s1=size(arr, 1) - s2=size(arr, 2) - write(unit, iostat=iostat, iomsg=iomsg) s1, s2 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2) - else - s1=0 - s2=0 - write(unit, iostat=iostat, iomsg=iomsg) s1, s2 - end if -end subroutine write2d_real - -subroutine write2d_int(arr, unit, iostat, iomsg) - integer, intent(in), allocatable :: arr(:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2 - - if (allocated(arr)) then - s1=size(arr, 1) - s2=size(arr, 2) - write(unit, iostat=iostat, iomsg=iomsg) s1, s2 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2) - else - s1=0 - s2=0 - write(unit, iostat=iostat, iomsg=iomsg) s1, s2 - end if -end subroutine write2d_int - - -subroutine write3d_real(arr, unit, iostat, iomsg) - real(kind=WP), intent(in), allocatable :: arr(:,:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2, s3 - - if (allocated(arr)) then - s1=size(arr, 1) - s2=size(arr, 2) - s3=size(arr, 3) - write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3) - else - s1=0 - s2=0 - s3=0 - write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3 - end if -end subroutine write3d_real - -subroutine write3d_int(arr, unit, iostat, iomsg) - integer, intent(in), allocatable :: arr(:,:,:) - integer, intent(in) :: unit - integer, intent(out) :: iostat - character(*), intent(inout) :: iomsg - integer :: s1, s2, s3 - - if (allocated(arr)) then - s1=size(arr, 1) - s2=size(arr, 2) - s3=size(arr, 3) - write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3 - write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3) - else - s1=0 - s2=0 - s3=0 - write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3 - end if -end subroutine write3d_int -end module MOD_WRITE_BINARY_ARRAYS -!========================================================== - diff --git a/src/temp/gen_halo_exchange.F90 b/src/temp/gen_halo_exchange.F90 deleted file mode 100755 index 7b9f66e6b..000000000 --- a/src/temp/gen_halo_exchange.F90 +++ /dev/null @@ -1,2381 +0,0 @@ -! ======================================================================== -! Halo exchange routines + broadcast routines that collect information -! on the entire field (needed for output) -! The routines here are very similar, difference is the data type and -! exchange pattern. -! exchange_nod2D_i(arr(myDim_nod2D+eDim_nod2D)) INTEGER -! exchange_nod2D(arr(myDim_nod2D+eDim_nod2D)) WP -! exchange_nod3D(arr(nl-1,myDim_nod2D+eDim_nod2D)) WP -! exchange_nod3D_full(arr(nl,myDim_nod2D+eDim_nod2D)) WP -! exchange_edge2D(edge_array2D) WP not used currently !!! no buffer!!! -! exchange_edge3D(edge_array3D) WP not used currently !!! no buffer!!! -! exchange_elem3D(elem_array3D) WP -! exchange_elem2d_full -! exchange_elem2d_full_i -! ======================================================================== - -module g_comm - - use, intrinsic :: ISO_FORTRAN_ENV - - implicit none - -contains - -#ifdef DEBUG -! General version of the communication routine for 2D nodal fields -! Only needed in debug mode -subroutine check_mpi_comm(rn, sn, r_mpitype, s_mpitype, rPE, sPE, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer, intent(in) :: sn, rn, r_mpitype(:), s_mpitype(:), rPE(:), sPE(:) -integer :: n, sdebug, rdebug, status(MPI_STATUS_SIZE), request -#include "associate_part_def.h" -#include "associate_part_ass.h" -DO n=1,rn - CALL MPI_TYPE_SIZE(r_mpitype(n), rdebug, MPIerr) - CALL MPI_ISEND(rdebug, 1, MPI_INTEGER, rPE(n), 10, MPI_COMM_FESOM, request, MPIerr) -END DO -DO n=1, sn - call MPI_RECV(sdebug, 1, MPI_INTEGER, sPE(n), 10, MPI_COMM_FESOM, & - status, MPIerr) - call MPI_TYPE_SIZE(s_mpitype(n), rdebug, MPIerr) - if (sdebug /= rdebug) then - print *, "Mismatching MPI send/recieve message lengths." - print *,"Send/receive process numbers: ", mype, '/', sPE(n) - print *,"Number of send/receive bytes: ", sdebug, '/', rdebug - call MPI_ABORT( MPI_COMM_FESOM, 1 ) - end if -END DO -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -END SUBROUTINE check_mpi_comm -#endif - - -subroutine exchange_nod2D_i(nod_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer, intent(inout) :: nod_array2D(:) -#include "associate_part_def.h" -#include "associate_part_ass.h" -if (npes > 1) then - call exchange_nod2D_i_begin(nod_array2D, partit) - call exchange_nod_end(partit) -endif -END SUBROUTINE exchange_nod2D_i - -!============================================================================= -! General version of the communication routine for 2D nodal fields -subroutine exchange_nod2D_i_begin(nod_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer, intent(inout) :: nod_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod2D_i, s_mpitype_nod2D_i, & - com_nod2D%rPE, com_nod2D%sPE) -#endif - - DO n=1,rn - - call MPI_IRECV(nod_array2D, 1, r_mpitype_nod2D_i(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(n), MPIerr) - END DO - - DO n=1, sn - - call MPI_ISEND(nod_array2D, 1, s_mpitype_nod2D_i(n), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(rn+n), MPIerr) - END DO - - com_nod2D%nreq = rn+sn - - endif -END SUBROUTINE exchange_nod2D_i_begin - -! ======================================================================== -! General version of the communication routine for 2D nodal fields -subroutine exchange_nod2D(nod_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array2D(:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - call exchange_nod2D_begin(nod_array2D, partit) - call exchange_nod_end(partit) - end if - -END SUBROUTINE exchange_nod2D - -! ======================================================================== -! General version of the communication routine for 2D nodal fields -subroutine exchange_nod2D_begin(nod_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod2D, s_mpitype_nod2D, & - com_nod2D%rPE, com_nod2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(nod_array2D, 1, r_mpitype_nod2D(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(nod_array2D, 1, s_mpitype_nod2D(n), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(rn+n), MPIerr) - END DO - - com_nod2D%nreq = rn+sn - - end if - -END SUBROUTINE exchange_nod2D_begin -!=============================================== -! General version of the communication routine for 2D nodal fields -subroutine exchange_nod2D_2fields(nod1_array2D, nod2_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - - - if (npes > 1) then - call exchange_nod2D_2fields_begin(nod1_array2D, nod2_array2D, partit) - call exchange_nod_end(partit) - end if - -END SUBROUTINE exchange_nod2D_2fields - -! ======================================================================== -! General version of the communication routine for 2D nodal fields -subroutine exchange_nod2D_2fields_begin(nod1_array2D, nod2_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes > 1) then - - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod2D, s_mpitype_nod2D, & - com_nod2D%rPE, com_nod2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(nod1_array2D, 1, r_mpitype_nod2D(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(2*n-1), MPIerr) - - call MPI_IRECV(nod2_array2D, 1, r_mpitype_nod2D(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n)+npes, MPI_COMM_FESOM, com_nod2D%req(2*n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(nod1_array2D, 1, s_mpitype_nod2D(n), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(2*rn+2*n-1), MPIerr) - - call MPI_ISEND(nod2_array2D, 1, s_mpitype_nod2D(n), com_nod2D%sPE(n), & - mype+npes, MPI_COMM_FESOM, com_nod2D%req(2*rn+2*n), MPIerr) - END DO - - com_nod2D%nreq = 2*(rn+sn) - -end if - -END SUBROUTINE exchange_nod2D_2fields_begin - -!=============================================== -subroutine exchange_nod2D_3fields(nod1_array2D, nod2_array2D, nod3_array2D, partit) -! General version of the communication routine for 2D nodal fields -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) -real(real64), intent(inout) :: nod3_array2D(:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - - - if (npes > 1) then - call exchange_nod2D_3fields_begin(nod1_array2D, nod2_array2D, nod3_array2D, partit) - call exchange_nod_end(partit) - end if - -END SUBROUTINE exchange_nod2D_3fields - -! ======================================================================== -subroutine exchange_nod2D_3fields_begin(nod1_array2D, nod2_array2D, nod3_array2D, partit) -! General version of the communication routine for 2D nodal fields -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) -real(real64), intent(inout) :: nod3_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod2D, s_mpitype_nod2D, & - com_nod2D%rPE, com_nod2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(nod1_array2D, 1, r_mpitype_nod2D(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(3*n-2), MPIerr) - - call MPI_IRECV(nod2_array2D, 1, r_mpitype_nod2D(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n)+npes, MPI_COMM_FESOM, com_nod2D%req(3*n-1), MPIerr) - - call MPI_IRECV(nod3_array2D, 1, r_mpitype_nod2D(n), com_nod2D%rPE(n), & - com_nod2D%rPE(n)+2*npes, MPI_COMM_FESOM, com_nod2D%req(3*n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(nod1_array2D, 1, s_mpitype_nod2D(n), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(3*rn+3*n-2), MPIerr) - - call MPI_ISEND(nod2_array2D, 1, s_mpitype_nod2D(n), com_nod2D%sPE(n), & - mype+npes, MPI_COMM_FESOM, com_nod2D%req(3*rn+3*n-1), MPIerr) - - call MPI_ISEND(nod3_array2D, 1, s_mpitype_nod2D(n), com_nod2D%sPE(n), & - mype+2*npes, MPI_COMM_FESOM, com_nod2D%req(3*rn+3*n), MPIerr) - END DO - - com_nod2D%nreq = 3*(rn+sn) - -end if - -END SUBROUTINE exchange_nod2D_3fields_begin - -! ======================================================================== -! General version of the communication routine for 3D nodal fields -! stored in (vertical, horizontal) format -subroutine exchange_nod3D(nod_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:) - -if (partit%npes > 1) then - call exchange_nod3D_begin(nod_array3D, partit) - call exchange_nod_end(partit) -endif - -END SUBROUTINE exchange_nod3D - -! ======================================================================== -! General version of the communication routine for 3D nodal fields -! stored in (vertical, horizontal) format -subroutine exchange_nod3D_begin(nod_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:) -integer :: n, sn, rn -integer :: nz, nl1 -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - nl1=ubound(nod_array3D,1) - - if ((nl1ubound(r_mpitype_nod3D, 2))) then - if (mype==0) then - print *,'Subroutine exchange_nod3D not implemented for',nl1,'layers.' - print *,'Adding the MPI datatypes is easy, see oce_modules.F90.' - endif - call par_ex(partit, 1) - endif - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod3D(:,nl1,1), s_mpitype_nod3D(:,nl1,1), & - com_nod2D%rPE, com_nod2D%sPE) -#endif - DO n=1,rn - call MPI_IRECV(nod_array3D, 1, r_mpitype_nod3D(n,nl1,1), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(nod_array3D, 1, s_mpitype_nod3D(n,nl1,1), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(rn+n), MPIerr) - END DO - com_nod2D%nreq = rn+sn - - endif -END SUBROUTINE exchange_nod3D_begin - -! ======================================================================== -! General version of the communication routine for 3D nodal fields -! stored in (vertical, horizontal) format -subroutine exchange_nod3D_2fields(nod1_array3D,nod2_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array3D(:,:) -real(real64), intent(inout) :: nod2_array3D(:,:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes > 1) then - call exchange_nod3D_2fields_begin(nod1_array3D,nod2_array3D, partit) - call exchange_nod_end(partit) -endif -END SUBROUTINE exchange_nod3D_2fields - -! ======================================================================== -subroutine exchange_nod3D_2fields_begin(nod1_array3D,nod2_array3D, partit) -! General version of the communication routine for 3D nodal fields -! stored in (vertical, horizontal) format -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array3D(:,:) -real(real64), intent(inout) :: nod2_array3D(:,:) -integer :: n, sn, rn -integer :: nz, nl1, nl2 -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - nl1 = ubound(nod1_array3D,1) - - if ((nl1ubound(r_mpitype_nod3D, 2))) then - if (mype==0) then - print *,'Subroutine exchange_nod3D not implemented for',nl1,'layers.' - print *,'Adding the MPI datatypes is easy, see oce_modules.F90.' - endif - call par_ex(1) - endif - - nl2 = ubound(nod2_array3D,1) - if ((nl2ubound(r_mpitype_nod3D, 2))) then - if (mype==0) then - print *,'Subroutine exchange_nod3D not implemented for',nl2,'layers.' - print *,'Adding the MPI datatypes is easy, see oce_modules.F90.' - endif - call par_ex(1) - endif - -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod3D(:,nl1,1), s_mpitype_nod3D(:,nl1,1), & - com_nod2D%rPE, com_nod2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(nod1_array3D, 1, r_mpitype_nod3D(n,nl1,1), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(2*n-1), MPIerr) - - call MPI_IRECV(nod2_array3D, 1, r_mpitype_nod3D(n,nl2,1), com_nod2D%rPE(n), & - com_nod2D%rPE(n)+npes, MPI_COMM_FESOM, com_nod2D%req(2*n ), MPIerr) - END DO - - DO n=1, sn - call MPI_ISEND(nod1_array3D, 1, s_mpitype_nod3D(n,nl1,1), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(2*rn+2*n-1), MPIerr) - - call MPI_ISEND(nod2_array3D, 1, s_mpitype_nod3D(n,nl2,1), com_nod2D%sPE(n), & - mype+npes, MPI_COMM_FESOM, com_nod2D%req(2*rn+2*n), MPIerr) - END DO - - com_nod2D%nreq = 2*(rn+sn) - - endif -END SUBROUTINE exchange_nod3D_2fields_begin -! ======================================================================== -subroutine exchange_nod3D_n(nod_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:,:) -if (partit%npes>1) then - call exchange_nod3D_n_begin(nod_array3D, partit) - call exchange_nod_end(partit) -endif - -END SUBROUTINE exchange_nod3D_n - -!================================================= -! General version of the communication routine for 3D nodal fields -! stored in (vertical, horizontal) format -subroutine exchange_nod3D_n_begin(nod_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:,:) -integer :: n, sn, rn -integer :: nz, nl1, n_val -#include "associate_part_def.h" -#include "associate_part_ass.h" -if (npes>1) then - ! nod_array3D(n_val,nl1,nod2D_size) - nl1 = ubound(nod_array3D,2) - n_val = ubound(nod_array3D,1) - if ((nl1ubound(r_mpitype_nod3D, 2)) .or. (n_val > 3)) then - ! This routine also works for swapped dimensions nod_array3D(nl1,n_val, nod2D_size) - nl1 = ubound(nod_array3D,1) - n_val = ubound(nod_array3D,2) - - if ((nl1ubound(r_mpitype_nod3D, 2)) .or. (n_val > 3)) then - if (mype==0) then - print *,'Subroutine exchange_nod3D_n not implemented for' - print *,nl1,'layers and / or ',n_val,'values per element.' - print *,'Adding the MPI datatypes is easy, see oce_modules.F90.' - endif - call par_ex(1) - endif - endif - sn=com_nod2D%sPEnum - rn=com_nod2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_nod3D(:,nl1,n_val), & - s_mpitype_nod3D(:,nl1,n_val), com_nod2D%rPE, com_nod2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(nod_array3D, 1, r_mpitype_nod3D(n,nl1,n_val), com_nod2D%rPE(n), & - com_nod2D%rPE(n), MPI_COMM_FESOM, com_nod2D%req(n), MPIerr) - END DO - - DO n=1, sn - call MPI_ISEND(nod_array3D, 1, s_mpitype_nod3D(n,nl1,n_val), com_nod2D%sPE(n), & - mype, MPI_COMM_FESOM, com_nod2D%req(rn+n), MPIerr) - END DO - - com_nod2D%nreq = rn+sn - - endif - - -END SUBROUTINE exchange_nod3D_n_begin - -!======================================= -! AND WAITING -!======================================= - -SUBROUTINE exchange_nod_end(partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit - -if (partit%npes > 1) & - call MPI_WAITALL(partit%com_nod2D%nreq, partit%com_nod2D%req, MPI_STATUSES_IGNORE, partit%MPIerr) - -END SUBROUTINE exchange_nod_end - -SUBROUTINE exchange_elem_end(partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes > 1) then - if (elem_full_flag) then - call MPI_WAITALL(com_elem2D_full%nreq, & - com_elem2D_full%req, MPI_STATUSES_IGNORE, MPIerr) - else - call MPI_WAITALL(com_elem2D%nreq, & - com_elem2D%req, MPI_STATUSES_IGNORE, MPIerr) - endif - end if -END SUBROUTINE exchange_elem_end -!============================================================================= -subroutine exchange_elem3D(elem_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - -call exchange_elem3D_begin(elem_array3D, partit) -call exchange_elem_end(partit) - -END SUBROUTINE exchange_elem3D -!=========================================== -! General version of the communication routine for 3D elemental fields -! stored in (vertical, horizontal) format -subroutine exchange_elem3D_begin(elem_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:) -integer :: n, sn, rn, nl1 -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - - nl1=ubound(elem_array3D,1) - - if (ubound(elem_array3D,2)<=myDim_elem2D+eDim_elem2D) then - - elem_full_flag = .false. - - sn=com_elem2D%sPEnum - rn=com_elem2D%rPEnum - - if (nl1==ubound(r_mpitype_elem3D, 2) .or. nl1==ubound(r_mpitype_elem3D, 2)-1) then - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem3D(:,nl1,1), s_mpitype_elem3D(:,nl1,1), & - com_elem2D%rPE, com_elem2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array3D, 1, r_mpitype_elem3D(n,nl1,1), com_elem2D%rPE(n), & - com_elem2D%rPE(n), MPI_COMM_FESOM, & - com_elem2D%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array3D, 1, s_mpitype_elem3D(n,nl1,1), com_elem2D%sPE(n), & - mype, MPI_COMM_FESOM, & - com_elem2D%req(rn+n), MPIerr) - END DO - - elseif (nl1 <= 4) then - ! In fact, this is a 2D-array with up to 4 values, e.g. derivatives - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem2D(:,nl1), s_mpitype_elem2D(:,nl1), & - com_elem2D%rPE, com_elem2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array3D, 1, r_mpitype_elem2D(n,nl1), com_elem2D%rPE(n), & - com_elem2D%rPE(n), MPI_COMM_FESOM, & - com_elem2D%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array3D, 1, s_mpitype_elem2D(n,nl1), com_elem2D%sPE(n), & - mype, MPI_COMM_FESOM, & - com_elem2D%req(rn+n), MPIerr) - END DO - else - if (mype==0) print *,'Sorry, no MPI datatype prepared for',nl1,'values per element (exchange_elem3D)' - call par_ex(1) - endif - - com_elem2D%nreq = rn+sn - - else - - elem_full_flag = .true. - - sn=com_elem2D_full%sPEnum - rn=com_elem2D_full%rPEnum - - if (nl1==ubound(r_mpitype_elem3D_full, 2) .or. nl1==ubound(r_mpitype_elem3D_full, 2)-1) then - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem3D_full(:,nl1,1), & - s_mpitype_elem3D_full(:,nl1,1), com_elem2D_full%rPE, com_elem2D_full%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array3D, 1, r_mpitype_elem3D_full(n,nl1,1), & - com_elem2D_full%rPE(n), & - com_elem2D_full%rPE(n), MPI_COMM_FESOM, & - com_elem2D_full%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array3D, 1, s_mpitype_elem3D_full(n,nl1,1), & - com_elem2D_full%sPE(n), & - mype, MPI_COMM_FESOM, & - com_elem2D_full%req(rn+n), MPIerr) - END DO - elseif (nl1 <= 4) then - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem2D_full(:,nl1), & - s_mpitype_elem2D_full(:,nl1), com_elem2D_full%rPE, com_elem2D_full%sPE) -#endif - - ! In fact, this is a 2D-array with up to 4 values, e.g. derivatives - DO n=1,rn - call MPI_IRECV(elem_array3D, 1, r_mpitype_elem2D_full(n,nl1), & - com_elem2D_full%rPE(n), & - com_elem2D_full%rPE(n), MPI_COMM_FESOM, & - com_elem2D_full%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array3D, 1, s_mpitype_elem2D_full(n,nl1), & - com_elem2D_full%sPE(n), & - mype, MPI_COMM_FESOM, & - com_elem2D_full%req(rn+n), MPIerr) - END DO - else - if (mype==0) print *,'Sorry, no MPI datatype prepared for',nl1,'values per element (exchange_elem3D)' - call par_ex(1) - endif - - com_elem2D_full%nreq = rn+sn - - endif - -endif - -END SUBROUTINE exchange_elem3D_begin - -!============================================================================= -! General version of the communication routine for 3D elemental fields -! stored in (vertical, horizontal) format -subroutine exchange_elem3D_n(elem_array3D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:,:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes> 1) then - call exchange_elem3D_n_begin(elem_array3D, partit) - call exchange_elem_end(partit) - endif -END SUBROUTINE exchange_elem3D_n -!============================================================================= -subroutine exchange_elem3D_n_begin(elem_array3D, partit) -! General version of the communication routine for 3D elemental fields -! stored in (vertical, horizontal) format -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:,:) -integer :: n, sn, rn, n_val, nl1 -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes> 1) then - nl1 = ubound(elem_array3D,2) - n_val = ubound(elem_array3D,1) - - if ((nl1ubound(r_mpitype_elem3D, 2)) .or. (n_val > 4)) then - - ! This routine also works for swapped dimensions elem_array3D(nl1,n_val, elem2D_size) - nl1= ubound(elem_array3D,1) - n_val = ubound(elem_array3D,2) - - if ((nl1ubound(r_mpitype_elem3D, 2)) .or. (n_val > 4)) then - if (mype==0) then - print *,'Subroutine exchange_elem3D_n not implemented for' - print *,nl1,'layers and / or ',n_val,'values per element.' - print *,'Adding the MPI datatypes is easy, see oce_modules.F90.' - endif - call par_ex(1) - endif - endif - - if (ubound(elem_array3D,3)<=myDim_elem2D+eDim_elem2D) then - - elem_full_flag = .false. - - sn=com_elem2D%sPEnum - rn=com_elem2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem3D(:,nl1,n_val), & - s_mpitype_elem3D(:,nl1,n_val), com_elem2D%rPE, com_elem2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array3D, 1, r_mpitype_elem3D(n,nl1,n_val), com_elem2D%rPE(n), & - com_elem2D%rPE(n), MPI_COMM_FESOM, com_elem2D%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array3D, 1, s_mpitype_elem3D(n,nl1,n_val), com_elem2D%sPE(n), & - mype, MPI_COMM_FESOM, com_elem2D%req(rn+n), MPIerr) - END DO - - com_elem2D%nreq = rn+sn - - else - - elem_full_flag = .true. - - sn=com_elem2D_full%sPEnum - rn=com_elem2D_full%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem3D_full(:,nl1,n_val), & - s_mpitype_elem3D_full(:,nl1,n_val), com_elem2D_full%rPE, com_elem2D_full%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array3D, 1, r_mpitype_elem3D_full(n,nl1,n_val), com_elem2D_full%rPE(n), & - com_elem2D_full%rPE(n), MPI_COMM_FESOM, com_elem2D_full%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array3D, 1, s_mpitype_elem3D_full(n,nl1,n_val), com_elem2D_full%sPE(n), & - mype, MPI_COMM_FESOM, com_elem2D_full%req(rn+n), MPIerr) - END DO - - com_elem2D_full%nreq = rn+sn - - end if - - -endif -END SUBROUTINE exchange_elem3D_n_begin -!======================================================================== -! General version of the communication routine for 3D elemental fields -! stored in (vertical, horizontal) format -subroutine exchange_elem2D(elem_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array2D(:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes> 1) then - call exchange_elem2D_begin(elem_array2D, partit) - call exchange_elem_end(partit) - end if - -END SUBROUTINE exchange_elem2D -!======================================================================== -! General version of the communication routine for 3D elemental fields -! stored in (vertical, horizontal) format -subroutine exchange_elem2D_begin(elem_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - - if (ubound(elem_array2D,1)<=myDim_elem2D+eDim_elem2D) then - - elem_full_flag = .false. - - sn=com_elem2D%sPEnum - rn=com_elem2D%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem2D(:,1), s_mpitype_elem2D(:,1), & - com_elem2D%rPE, com_elem2D%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array2D, 1, r_mpitype_elem2D(n,1), com_elem2D%rPE(n), & - com_elem2D%rPE(n), MPI_COMM_FESOM, com_elem2D%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array2D, 1, s_mpitype_elem2D(n,1), com_elem2D%sPE(n), & - mype, MPI_COMM_FESOM, com_elem2D%req(rn+n), MPIerr) - END DO - - com_elem2D%nreq = rn+sn - - else - elem_full_flag = .true. - - sn=com_elem2D_full%sPEnum - rn=com_elem2D_full%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem2D_full(:,1), s_mpitype_elem2D_full(:,1), & - com_elem2D_full%rPE, com_elem2D_full%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array2D, 1, r_mpitype_elem2D_full(n,1), com_elem2D_full%rPE(n), & - com_elem2D_full%rPE(n), MPI_COMM_FESOM, com_elem2D_full%req(n), MPIerr) - END DO - DO n=1, sn - call MPI_ISEND(elem_array2D, 1, s_mpitype_elem2D_full(n,1), com_elem2D_full%sPE(n), & - mype, MPI_COMM_FESOM, com_elem2D_full%req(rn+n), MPIerr) - END DO - - com_elem2D_full%nreq = rn+sn - - end if - -end if - -END SUBROUTINE exchange_elem2D_begin -! ======================================================================== -!Exchange with ALL(!) the neighbours -subroutine exchange_elem2D_i(elem_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer, intent(inout) :: elem_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes> 1) then - call exchange_elem2D_i_begin(elem_array2D, partit) - call exchange_elem_end(partit) -end if - -END SUBROUTINE exchange_elem2D_i -!============================================================================= -!Exchange with ALL(!) the neighbours -subroutine exchange_elem2D_i_begin(elem_array2D, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer, intent(inout) :: elem_array2D(:) -integer :: n, sn, rn -#include "associate_part_def.h" -#include "associate_part_ass.h" - - if (npes> 1) then - - elem_full_flag = .true. - - sn=com_elem2D_full%sPEnum - rn=com_elem2D_full%rPEnum - - ! Check MPI point-to-point communication for consistency -#ifdef DEBUG - call check_mpi_comm(rn, sn, r_mpitype_elem2D_full_i, s_mpitype_elem2D_full_i, & - com_elem2D_full%rPE, com_elem2D_full%sPE) -#endif - - DO n=1,rn - call MPI_IRECV(elem_array2D, 1, r_mpitype_elem2D_full_i(n), com_elem2D_full%rPE(n), & - com_elem2D_full%rPE(n), MPI_COMM_FESOM, com_elem2D_full%req(n), MPIerr) - END DO - - DO n=1, sn - - call MPI_ISEND(elem_array2D, 1, s_mpitype_elem2D_full_i(n), com_elem2D_full%sPE(n), & - mype, MPI_COMM_FESOM, com_elem2D_full%req(rn+n), MPIerr) - END DO - - com_elem2D_full%nreq = rn+sn - -end if - -END SUBROUTINE exchange_elem2D_i_begin -! ======================================================================== -! Broadcast routines -! Many because of different sizes. -! ======================================================================== -subroutine broadcast_nod3D(arr3D, arr3Dglobal, partit) -! Distribute the nodal information available on 0 PE to other PEs -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nz, counter,nl1 -integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -real(real64) :: arr3D(:,:) -real(real64) :: arr3Dglobal(:,:) -integer :: node_size -INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf, recvbuf -#include "associate_part_def.h" -#include "associate_part_ass.h" - -node_size=myDim_nod2D+eDim_nod2D -nl1=ubound(arr3D,1) -IF ( mype == 0 ) THEN - if (npes>1) then - arr3D(:,1:node_size)=arr3Dglobal(:,myList_nod2D(1:node_size)) - end if - DO n = 1, npes-1 - CALL MPI_RECV( nTS, 1, MPI_INTEGER, MPI_ANY_SOURCE, & - 0, MPI_COMM_FESOM, status, MPIerr ) - sender = status(MPI_SOURCE) - ALLOCATE(sendbuf(nTS*nl1), irecvbuf(nTS)) - - CALL MPI_RECV(irecvbuf(1), nTS, MPI_INTEGER, sender, & - 1, MPI_COMM_FESOM, status, MPIerr ) - counter=0 - DO i = 1, nTS - DO nz=1, nl1 - counter=counter+1 - sendbuf(counter) = arr3Dglobal(nz,irecvbuf(i)) - ENDDO - ENDDO - - CALL MPI_SEND(sendbuf(1), nTS*nl1, MPI_DOUBLE_PRECISION, & - sender, 2, MPI_COMM_FESOM, MPIerr ) - - DEALLOCATE(irecvbuf, sendbuf) - ENDDO -ELSE - CALL MPI_SEND( node_size, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( myList_nod2D(1), node_size, MPI_INTEGER, 0, 1, & - MPI_COMM_FESOM, MPIerr ) - - ALLOCATE(recvbuf(node_size*nl1)) - CALL MPI_RECV( recvbuf(1), node_size*nl1, MPI_DOUBLE_PRECISION, 0, & - 2, MPI_COMM_FESOM, status, MPIerr ) - counter=0 - DO n = 1, node_size - DO nz=1, nl1 - counter=counter+1 - arr3D(nz,n)=recvbuf(counter) - ENDDO - ENDDO - - DEALLOCATE(recvbuf) -ENDIF -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -end subroutine broadcast_nod3D -! -!============================================================================ -! -subroutine broadcast_nod2D(arr2D, arr2Dglobal, partit) -! A 2D version of the previous routine -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(in), target :: partit -real(real64) :: arr2D(:) -real(real64) :: arr2Dglobal(:) -integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf -integer :: node_size -#include "associate_part_def.h" -#include "associate_part_ass.h" - -node_size=myDim_nod2D+eDim_nod2D - -IF ( mype == 0 ) THEN - if (npes>1) then - arr2D(1:node_size)=arr2Dglobal(myList_nod2D(1:node_size)) - end if - DO n = 1, npes-1 - CALL MPI_RECV( nTS, 1, MPI_INTEGER, MPI_ANY_SOURCE, & - 0, MPI_COMM_FESOM, status, MPIerr ) - sender = status(MPI_SOURCE) - ALLOCATE(sendbuf(nTS), irecvbuf(nTS)) - - CALL MPI_RECV(irecvbuf(1), nTS, MPI_INTEGER, sender, & - 1, MPI_COMM_FESOM, status, MPIerr ) - DO i = 1, nTS - sendbuf(i) = arr2Dglobal(irecvbuf(i)) - ENDDO - - CALL MPI_SEND(sendbuf(1), nTS, MPI_DOUBLE_PRECISION, & - sender, 2, MPI_COMM_FESOM, MPIerr ) - - DEALLOCATE(irecvbuf, sendbuf) - ENDDO -ELSE - CALL MPI_SEND( node_size, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( myList_nod2D(1), node_size, MPI_INTEGER, 0, 1, & - MPI_COMM_FESOM, MPIerr ) - CALL MPI_RECV( arr2D(1), node_size, MPI_DOUBLE_PRECISION, 0, & - 2, MPI_COMM_FESOM, status, MPIerr ) -ENDIF -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -end subroutine broadcast_nod2D -! -!============================================================================ -! -subroutine broadcast_elem3D(arr3D, arr3Dglobal, partit) -! Distribute the elemental information available on 0 PE to other PEs -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(in), target :: partit -INTEGER :: nz, counter,nl1 -integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -real(real64) :: arr3D(:,:) -real(real64) :: arr3Dglobal(:,:) -integer :: elem_size - -INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf, recvbuf -#include "associate_part_def.h" -#include "associate_part_ass.h" - -elem_size=myDim_elem2D+eDim_elem2D - -nl1=ubound(arr3D,1) -IF ( mype == 0 ) THEN - if (npes>1) then - arr3D(:,1:elem_size)=arr3Dglobal(:,myList_elem2D(1:elem_size)) - end if - DO n = 1, npes-1 - CALL MPI_RECV( nTS, 1, MPI_INTEGER, MPI_ANY_SOURCE, & - 0, MPI_COMM_FESOM, status, MPIerr ) - sender = status(MPI_SOURCE) - ALLOCATE(sendbuf(nTS*nl1), irecvbuf(nTS)) - - CALL MPI_RECV(irecvbuf(1), nTS, MPI_INTEGER, sender, & - 1, MPI_COMM_FESOM, status, MPIerr ) - counter=0 - DO i = 1, nTS - DO nz=1, nl1 - counter=counter+1 - sendbuf(counter) = arr3Dglobal(nz,irecvbuf(i)) - ENDDO - ENDDO - - CALL MPI_SEND(sendbuf(1), nTS*nl1, MPI_DOUBLE_PRECISION, & - sender, 2, MPI_COMM_FESOM, MPIerr ) - - DEALLOCATE(irecvbuf, sendbuf) - ENDDO -ELSE - CALL MPI_SEND( elem_size, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( myList_elem2D(1), elem_size, MPI_INTEGER, 0, 1, & - MPI_COMM_FESOM, MPIerr ) - - ALLOCATE(recvbuf(elem_size*nl1)) - CALL MPI_RECV( recvbuf(1), elem_size*nl1, MPI_DOUBLE_PRECISION, 0, & - 2, MPI_COMM_FESOM, status, MPIerr ) - counter=0 - DO n = 1, elem_size - DO nz=1, nl1 - counter=counter+1 - arr3D(nz,n)=recvbuf(counter) - ENDDO - ENDDO - - DEALLOCATE(recvbuf) -ENDIF -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -end subroutine broadcast_elem3D -! -!============================================================================ -! -subroutine broadcast_elem2D(arr2D, arr2Dglobal, partit) -! A 2D version of the previous routine -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(in), target :: partit -integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -real(real64) :: arr2D(:) -real(real64) :: arr2Dglobal(:) -integer :: elem_size -INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf -#include "associate_part_def.h" -#include "associate_part_ass.h" - -elem_size=myDim_elem2D+eDim_elem2D - -IF ( mype == 0 ) THEN - if (npes>1) then - arr2D(1:elem_size)=arr2Dglobal(myList_elem2D(1:elem_size)) - end if - DO n = 1, npes-1 - CALL MPI_RECV( nTS, 1, MPI_INTEGER, MPI_ANY_SOURCE, & - 0, MPI_COMM_FESOM, status, MPIerr ) - sender = status(MPI_SOURCE) - ALLOCATE(sendbuf(1:nTS), irecvbuf(nTS)) - - CALL MPI_RECV(irecvbuf(1), nTS, MPI_INTEGER, sender, & - 1, MPI_COMM_FESOM, status, MPIerr ) - DO i = 1, nTS - sendbuf(i) = arr2Dglobal(irecvbuf(i)) - ENDDO - - CALL MPI_SEND(sendbuf(1), nTS, MPI_DOUBLE_PRECISION, & - sender, 2, MPI_COMM_FESOM, MPIerr ) - - DEALLOCATE(irecvbuf, sendbuf) - ENDDO -ELSE - CALL MPI_SEND( elem_size, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( myList_elem2D(1), elem_size, MPI_INTEGER, 0, 1, & - MPI_COMM_FESOM, MPIerr ) - CALL MPI_RECV( arr2D(1), elem_size, MPI_DOUBLE_PRECISION, 0, & - 2, MPI_COMM_FESOM, status, MPIerr ) -ENDIF -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -end subroutine broadcast_elem2D -! -!============================================================================ -! Make nodal information available to master PE -! Use only with 3D arrays stored in (vertical, horizontal) way -subroutine gather_nod3D(arr3D, arr3D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -real(real64) :: arr3D(:,:) -real(real64) :: arr3D_global(:,:) -real(real64), allocatable :: recvbuf(:,:) -integer :: req(partit%npes-1) -integer :: start, n3D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(nl1,ubound(arr3D_global,2))) - - do n = 1, npes-1 - n3D = (remPtr_nod2D(n+1) - remPtr_nod2D(n))*nl1 - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(1,start), n3D, MPI_DOUBLE_PRECISION, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_nod2D(1:myDim_nod2D)) = arr3D(1:nl1,1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_nod2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - call MPI_SEND( arr3D, myDim_nod2D*nl1, MPI_DOUBLE_PRECISION, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -end if -end subroutine gather_nod3D -! -!============================================================================ -! -subroutine gather_real4_nod3D(arr3D, arr3D_global, partit) - -! Make nodal information available to master PE -! -! Use only with 3D arrays stored in (vertical, horizontal) way -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -real(real32) :: arr3D(:,:) -real(real32) :: arr3D_global(:,:) -real(real32), allocatable :: recvbuf(:,:) -integer :: req(partit%npes-1) -integer :: start, n3D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(nl1,ubound(arr3D_global,2))) - - do n = 1, npes-1 - n3D = (remPtr_nod2D(n+1) - remPtr_nod2D(n))*nl1 - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(1,start), n3D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_nod2D(1:myDim_nod2D)) = arr3D(1:nl1,1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_nod2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - call MPI_SEND( arr3D, myDim_nod2D*nl1, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -end if -end subroutine gather_real4_nod3D -!======================================================= - -subroutine gather_int2_nod3D(arr3D, arr3D_global, partit) - -! Make nodal information available to master PE -! -! Use only with 3D arrays stored in (vertical, horizontal) way -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -integer(int16) :: arr3D(:,:) -integer(int16) :: arr3D_global(:,:) -integer(int16), allocatable :: recvbuf(:,:) -integer :: req(partit%npes-1) -integer :: start, n3D -#include "associate_part_def.h" -#include "associate_part_ass.h" - - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(nl1,ubound(arr3D_global,2))) - - do n = 1, npes-1 - n3D = (remPtr_nod2D(n+1) - remPtr_nod2D(n))*nl1 - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(1,start), n3D, MPI_SHORT, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_nod2D(1:myDim_nod2D)) = arr3D(1:nl1,1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_nod2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - call MPI_SEND( arr3D, myDim_nod2D*nl1, MPI_SHORT, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -end if -end subroutine gather_int2_nod3D -!============================================== -subroutine gather_nod2D(arr2D, arr2D_global, partit) -! Make nodal information available to master PE -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -real(real64) :: arr2D(:) -real(real64) :: arr2D_global(:) -real(real64), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, n2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(ubound(arr2D_global,1))) - do n = 1, npes-1 - n2D = remPtr_nod2D(n+1) - remPtr_nod2D(n) - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(start), n2D, MPI_DOUBLE_PRECISION, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_nod2D(1:myDim_nod2D)) = arr2D(1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1 : remPtr_nod2D(npes)-1) - deallocate(recvbuf) - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_nod2D, MPI_DOUBLE_PRECISION, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_nod2D -!============================================== -subroutine gather_real4_nod2D(arr2D, arr2D_global, partit) -! Make nodal information available to master PE -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -real(real32) :: arr2D(:) -real(real32) :: arr2D_global(:) -real(real32), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, n2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(ubound(arr2D_global,1))) - do n = 1, npes-1 - n2D = remPtr_nod2D(n+1) - remPtr_nod2D(n) - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(start), n2D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_nod2D(1:myDim_nod2D)) = arr2D(1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1 : remPtr_nod2D(npes)-1) - deallocate(recvbuf) - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_nod2D, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_real4_nod2D - -!============================================== -! Make nodal information available to master PE -subroutine gather_int2_nod2D(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -integer(int16) :: arr2D(:) -integer(int16) :: arr2D_global(:) -integer(int16), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, n2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(ubound(arr2D_global,1))) - do n = 1, npes-1 - n2D = remPtr_nod2D(n+1) - remPtr_nod2D(n) - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(start), n2D, MPI_SHORT, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_nod2D(1:myDim_nod2D)) = arr2D(1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1 : remPtr_nod2D(npes)-1) - deallocate(recvbuf) - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_nod2D, MPI_SHORT, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_int2_nod2D - -!============================================================================ -subroutine gather_elem3D(arr3D, arr3D_global, partit) -! Make element information available to master PE -! -! Use only with 3D arrays stored in (vertical, horizontal) way -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -real(real64) :: arr3D(:,:) -real(real64) :: arr3D_global(:,:) -real(real64), allocatable :: recvbuf(:,:) -integer :: req(partit%npes-1) -integer :: start, e3D, ende, err_alloc -integer :: max_loc_Dim, i, status(MPI_STATUS_SIZE) -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global -! (Carefull with duplicate interface elements, coming from two -! PEs at once!) - -IF ( mype == 0 ) THEN - - if (npes>1) then -! - allocate(recvbuf(nl1,remPtr_elem2D(npes))) - - do n = 1, npes-1 - e3D = (remPtr_elem2D(n+1) - remPtr_elem2D(n))*nl1 - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(1,start), e3D, MPI_DOUBLE_PRECISION, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_elem2D(1:myDim_elem2D)) = arr3D(1:nl1,1:myDim_elem2D) - - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - call MPI_SEND( arr3D, myDim_elem2D*nl1, MPI_DOUBLE_PRECISION, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_elem3D - -!=================================================================== -! Make element information available to master PE -! Use only with 3D arrays stored in (vertical, horizontal) way -subroutine gather_real4_elem3D(arr3D, arr3D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -real(real32) :: arr3D(:,:) -real(real32) :: arr3D_global(:,:) -real(real32), allocatable :: recvbuf(:,:) -integer :: req(partit%npes-1) -integer :: start, e3D, ende, err_alloc -integer :: max_loc_Dim, i, status(MPI_STATUS_SIZE) -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global -! (Carefull with duplicate interface elements, coming from two -! PEs at once!) - -IF ( mype == 0 ) THEN - - if (npes>1) then -! - allocate(recvbuf(nl1,remPtr_elem2D(npes))) - - do n = 1, npes-1 - e3D = (remPtr_elem2D(n+1) - remPtr_elem2D(n))*nl1 - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(1,start), e3D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_elem2D(1:myDim_elem2D)) = arr3D(1:nl1,1:myDim_elem2D) - - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - call MPI_SEND( arr3D, myDim_elem2D*nl1, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_real4_elem3D - - -!=================================================================== -! Make element information available to master PE -! Use only with 3D arrays stored in (vertical, horizontal) way -subroutine gather_int2_elem3D(arr3D, arr3D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -integer(int16) :: arr3D(:,:) -integer(int16) :: arr3D_global(:,:) -integer(int16), allocatable :: recvbuf(:,:) -integer :: req(partit%npes-1) -integer :: start, e3D, ende, err_alloc -integer :: max_loc_Dim, i, status(MPI_STATUS_SIZE) -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global -! (Carefull with duplicate interface elements, coming from two -! PEs at once!) - -IF ( mype == 0 ) THEN - - if (npes>1) then -! - allocate(recvbuf(nl1,remPtr_elem2D(npes))) - - do n = 1, npes-1 - e3D = (remPtr_elem2D(n+1) - remPtr_elem2D(n))*nl1 - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(1,start), e3D, MPI_SHORT, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_elem2D(1:myDim_elem2D)) = arr3D(1:nl1,1:myDim_elem2D) - - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - call MPI_SEND( arr3D, myDim_elem2D*nl1, MPI_SHORT, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_int2_elem3D - - -!============================================== -! Make element information available to master PE -subroutine gather_elem2D(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -real(real64) :: arr2D(:) -real(real64) :: arr2D_global(:) -real(real64), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, e2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - - allocate(recvbuf(remPtr_elem2D(npes))) - - do n = 1, npes-1 - e2D = remPtr_elem2D(n+1) - remPtr_elem2D(n) - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(start), e2D, MPI_DOUBLE_PRECISION, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_elem2D(1:myDim_elem2D)) = arr2D(1:myDim_elem2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_elem2D, MPI_DOUBLE_PRECISION, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF -end if - -end subroutine gather_elem2D - -!================================================ -! Make element information available to master PE -subroutine gather_real4_elem2D(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -real(real32) :: arr2D(:) -real(real32) :: arr2D_global(:) -real(real32), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, e2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - - allocate(recvbuf(remPtr_elem2D(npes))) - - do n = 1, npes-1 - e2D = remPtr_elem2D(n+1) - remPtr_elem2D(n) - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(start), e2D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_elem2D(1:myDim_elem2D)) = arr2D(1:myDim_elem2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_elem2D, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF -end if - -end subroutine gather_real4_elem2D - -!================================================ -! Make element information available to master PE -subroutine gather_int2_elem2D(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -integer(int16) :: arr2D(:) -integer(int16) :: arr2D_global(:) -integer(int16), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, e2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - - allocate(recvbuf(remPtr_elem2D(npes))) - - do n = 1, npes-1 - e2D = remPtr_elem2D(n+1) - remPtr_elem2D(n) - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(start), e2D, MPI_SHORT, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_elem2D(1:myDim_elem2D)) = arr2D(1:myDim_elem2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_elem2D, MPI_SHORT, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF -end if - -end subroutine gather_int2_elem2D - - -!============================================================================ -! Make nodal information available to master PE -! Use only with 3D arrays stored in (vertical, horizontal) way -subroutine gather_real8to4_nod3D(arr3D, arr3D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -real(real64) :: arr3D(:,:) -real(real32) :: arr3D_global(:,:) -integer :: req(partit%npes-1) -integer :: start, n3D, ierr -real(real32), allocatable :: recvbuf(:,:) -real(real32), allocatable :: sendbuf(:,:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(nl1, ubound(arr3D_global,2))) - - do n = 1, npes-1 - n3D = (remPtr_nod2D(n+1) - remPtr_nod2D(n))*nl1 - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(1,start), n3D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_nod2D(1:myDim_nod2D)) = arr3D(1:nl1,1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_nod2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - - allocate(sendbuf(nl1,myDim_nod2D)) - sendbuf(1:nl1,1:myDim_nod2D) = arr3D(1:nl1,1:myDim_nod2D) - - call MPI_SEND(sendbuf, myDim_nod2D*nl1, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - deallocate(sendbuf) - -ENDIF - -end if - -end subroutine gather_real8to4_nod3D -!============================================== -! Make nodal information available to master PE -subroutine gather_real8to4_nod2D(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -real(real64) :: arr2D(:) -real(real32) :: arr2D_global(:) -real(real32) :: sendbuf(partit%myDim_nod2D) -real(real64), allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, n2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -! Consider MPI-datatypes to recv directly into arr2D_global! - - if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(ubound(arr2D_global,1))) - do n = 1, npes-1 - n2D = remPtr_nod2D(n+1) - remPtr_nod2D(n) - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(start), n2D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_nod2D(1:myDim_nod2D)) = arr2D(1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1 : remPtr_nod2D(npes)-1) - deallocate(recvbuf) - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - sendbuf(1:myDim_nod2D) = real(arr2D(1:myDim_nod2D),real32) - - call MPI_SEND(sendbuf, myDim_nod2D, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -end if -end subroutine gather_real8to4_nod2D -!============================================================================ -subroutine gather_real8to4_elem3D(arr3D, arr3D_global, partit) -! Make element information available to master PE -! Use only with 3D arrays stored in (vertical, horizontal) way -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -INTEGER :: nl1 -integer :: n -real(real64) :: arr3D(:,:) -real(real32) :: arr3D_global(:,:) -integer :: req(partit%npes-1) -integer :: start, e3D -real(real32), allocatable :: recvbuf(:,:) -real(real32), allocatable :: sendbuf(:,:) -#include "associate_part_def.h" -#include "associate_part_ass.h" - - -if (npes> 1) then -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -nl1=ubound(arr3D,1) - -! Consider MPI-datatypes to recv directly into arr3D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(nl1,remPtr_elem2D(npes))) - - do n = 1, npes-1 - e3D = (remPtr_elem2D(n+1) - remPtr_elem2D(n))*nl1 - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(1,start), e3D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr3D_global(1:nl1,myList_elem2D(1:myDim_elem2D)) = arr3D(1:nl1,1:myDim_elem2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr3D_global(1:nl1, remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1:nl1, 1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - arr3D_global(:,:) = arr3D(:,:) - endif - -ELSE - allocate(sendbuf(nl1,myDim_elem2D)) - sendbuf(1:nl1,1:myDim_elem2D) = arr3D(1:nl1,1:myDim_elem2D) - - call MPI_SEND(sendbuf, myDim_elem2D*nl1, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - deallocate(sendbuf) -ENDIF - -end if -end subroutine gather_real8to4_elem3D -!================================================ -! Make element information available to master PE -subroutine gather_real8to4_elem2D(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -real(real64) :: arr2D(:) -real(real32) :: arr2D_global(:) -real(real32), allocatable :: recvbuf(:) -real(real32) :: sendbuf(partit%myDim_elem2D) -integer :: req(partit%npes-1) -integer :: start, e2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(remPtr_elem2D(npes))) - - do n = 1, npes-1 - e2D = remPtr_elem2D(n+1) - remPtr_elem2D(n) - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(start), e2D, MPI_REAL, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_elem2D(1:myDim_elem2D)) = arr2D(1:myDim_elem2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1 : remPtr_elem2D(npes)-1) - - deallocate(recvbuf) - - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - sendbuf(1:myDim_elem2D) = real(arr2D(1:myDim_elem2D),real32) - call MPI_SEND(sendbuf, myDim_elem2D, MPI_REAL, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -end if -end subroutine gather_real8to4_elem2D -!============================================== -subroutine gather_elem2D_i(arr2D, arr2D_global, partit) -! Make element information available to master PE -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -integer :: arr2D(:) -integer :: arr2D_global(:) -integer, allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, e2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - ! Consider MPI-datatypes to recv directly into arr2D_global! - IF ( mype == 0 ) THEN - if (npes > 1) then - allocate(recvbuf(remPtr_elem2D(npes))) - do n = 1, npes-1 - e2D = remPtr_elem2D(n+1) - remPtr_elem2D(n) - start = remPtr_elem2D(n) - call MPI_IRECV(recvbuf(start), e2D, MPI_INTEGER, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - arr2D_global(myList_elem2D(1:myDim_elem2D)) = arr2D(1:myDim_elem2D) - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - arr2D_global(remList_elem2D(1 : remPtr_elem2D(npes)-1)) & - = recvbuf(1 : remPtr_elem2D(npes)-1) - deallocate(recvbuf) - else - arr2D_global(:) = arr2D(:) - endif - ELSE - call MPI_SEND(arr2D, myDim_elem2D, MPI_INTEGER, 0, 2, MPI_COMM_FESOM, MPIerr ) - ENDIF -end subroutine gather_elem2D_i -!============================================== -! Make nodal information available to master PE -subroutine gather_nod2D_i(arr2D, arr2D_global, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: n -integer :: arr2D(:) -integer :: arr2D_global(:) -integer, allocatable :: recvbuf(:) -integer :: req(partit%npes-1) -integer :: start, n2D -#include "associate_part_def.h" -#include "associate_part_ass.h" - -if (npes> 1) then - -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) - -! Consider MPI-datatypes to recv directly into arr2D_global! - -IF ( mype == 0 ) THEN - - if (npes>1) then - allocate(recvbuf(ubound(arr2D_global, 1))) - do n = 1, npes-1 - n2D = remPtr_nod2D(n+1) - remPtr_nod2D(n) - start = remPtr_nod2D(n) - call MPI_IRECV(recvbuf(start), n2D, MPI_INTEGER, n, 2, MPI_COMM_FESOM, req(n), MPIerr) - enddo - - arr2D_global(myList_nod2D(1:myDim_nod2D)) = arr2D(1:myDim_nod2D) - - call MPI_WAITALL(npes-1, req, MPI_STATUSES_IGNORE, MPIerr) - - arr2D_global(remList_nod2D(1 : remPtr_nod2D(npes)-1)) & - = recvbuf(1 : remPtr_nod2D(npes)-1) - deallocate(recvbuf) - else - - arr2D_global(:) = arr2D(:) - - endif - -ELSE - - call MPI_SEND( arr2D, myDim_nod2D, MPI_INTEGER, 0, 2, MPI_COMM_FESOM, MPIerr ) - -ENDIF - -endif -end subroutine gather_nod2D_i -!============================================================================ -! A 2D version of the previous routine -subroutine gather_edg2D(arr2D, arr2Dglobal, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(in), target :: partit -real(real64) :: arr2D(:) -real(real64) :: arr2Dglobal(:) -integer :: i, n, buf_size, sender, status(MPI_STATUS_SIZE) -INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf -REAL(real64), ALLOCATABLE, DIMENSION(:) :: rbuf -#include "associate_part_def.h" -#include "associate_part_ass.h" - -IF ( mype == 0 ) THEN - arr2Dglobal(myList_edge2D(1:myDim_edge2D))=arr2D(1:myDim_edge2D) - DO n = 1, npes-1 - CALL MPI_RECV( buf_size, 1, MPI_INTEGER, MPI_ANY_SOURCE, & - 0, MPI_COMM_FESOM, status, MPIerr ) - sender = status(MPI_SOURCE) - ALLOCATE(rbuf(buf_size), ibuf(buf_size)) - - CALL MPI_RECV(ibuf(1), buf_size, MPI_INTEGER, sender, & - 1, MPI_COMM_FESOM, status, MPIerr ) - - CALL MPI_RECV(rbuf(1), buf_size, MPI_DOUBLE_PRECISION, sender, & - 2, MPI_COMM_FESOM, status, MPIerr ) - arr2Dglobal(ibuf)=rbuf - DEALLOCATE(ibuf, rbuf) - ENDDO -ELSE - CALL MPI_SEND( myDim_edge2D, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( myList_edge2D(1), myDim_edge2D, MPI_INTEGER, 0, 1, & - MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( arr2D(1), myDim_edge2D, MPI_DOUBLE_PRECISION, 0, 2,& - MPI_COMM_FESOM, MPIerr ) -ENDIF -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -end subroutine gather_edg2D -! -!============================================================================ -! A 2D version of the previous routine -subroutine gather_edg2D_i(arr2D, arr2Dglobal, partit) -use MOD_MESH -use MOD_PARTIT -IMPLICIT NONE -type(t_partit), intent(inout), target :: partit -integer :: arr2D(:) -integer :: arr2Dglobal(:) -integer :: i, n, buf_size, sender, status(MPI_STATUS_SIZE) -INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf, vbuf -#include "associate_part_def.h" -#include "associate_part_ass.h" - -IF ( mype == 0 ) THEN - arr2Dglobal(myList_edge2D(1:myDim_edge2D))=arr2D(1:myDim_edge2D) - DO n = 1, npes-1 - CALL MPI_RECV( buf_size, 1, MPI_INTEGER, MPI_ANY_SOURCE, & - 0, MPI_COMM_FESOM, status, MPIerr ) - sender = status(MPI_SOURCE) - ALLOCATE(ibuf(buf_size), vbuf(buf_size)) - - CALL MPI_RECV(ibuf(1), buf_size, MPI_INTEGER, sender, & - 1, MPI_COMM_FESOM, status, MPIerr ) - - CALL MPI_RECV(vbuf(1), buf_size, MPI_INTEGER, sender, & - 2, MPI_COMM_FESOM, status, MPIerr ) - arr2Dglobal(ibuf)=vbuf - DEALLOCATE(ibuf, vbuf) - ENDDO -ELSE - CALL MPI_SEND( myDim_edge2D, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( myList_edge2D(1), myDim_edge2D, MPI_INTEGER, 0, 1, & - MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( arr2D(1), myDim_edge2D, MPI_INTEGER, 0, 2,& - MPI_COMM_FESOM, MPIerr ) -ENDIF -CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) -end subroutine gather_edg2D_i -!============================================== - -end module g_comm - - - -module g_comm_auto -use g_comm -implicit none -interface exchange_nod - module procedure exchange_nod2D - module procedure exchange_nod2D_i - module procedure exchange_nod2D_2fields - module procedure exchange_nod2D_3fields - module procedure exchange_nod3D - module procedure exchange_nod3D_2fields - module procedure exchange_nod3D_n -end interface exchange_nod - -interface exchange_nod_begin - module procedure exchange_nod2D_begin - module procedure exchange_nod2D_i_begin - module procedure exchange_nod2D_2fields_begin - module procedure exchange_nod2D_3fields_begin - module procedure exchange_nod3D_begin - module procedure exchange_nod3D_2fields_begin - module procedure exchange_nod3D_n_begin -end interface exchange_nod_begin - -!!$interface exchange_edge -!!$ module procedure exchange_edge2D -!!$! module procedure exchange_edge3D ! not available, not used -!!$end interface exchange_edge - -interface exchange_elem - module procedure exchange_elem3D - module procedure exchange_elem3D_n - module procedure exchange_elem2d - module procedure exchange_elem2d_i -end interface exchange_elem - -interface exchange_elem_begin - module procedure exchange_elem3D_begin - module procedure exchange_elem3D_n_begin - module procedure exchange_elem2d_begin - module procedure exchange_elem2d_i_begin -end interface exchange_elem_begin - - -interface broadcast_nod - module procedure broadcast_nod3D - module procedure broadcast_nod2D -end interface broadcast_nod - -interface broadcast_elem - module procedure broadcast_elem3D - module procedure broadcast_elem2D -end interface broadcast_elem - -interface gather_nod - module procedure gather_nod3D - module procedure gather_nod2D - module procedure gather_real4_nod3D - module procedure gather_real4_nod2D - module procedure gather_int2_nod3D - module procedure gather_int2_nod2D - module procedure gather_real8to4_nod3D - module procedure gather_real8to4_nod2D - module procedure gather_nod2D_i -end interface gather_nod - -interface gather_elem - module procedure gather_elem3D - module procedure gather_elem2D - module procedure gather_real4_elem3D - module procedure gather_real4_elem2D - module procedure gather_int2_elem3D - module procedure gather_int2_elem2D - module procedure gather_real8to4_elem3D - module procedure gather_real8to4_elem2D - module procedure gather_elem2D_i -end interface gather_elem - -interface gather_edge - module procedure gather_edg2D - module procedure gather_edg2D_i -end interface gather_edge - - -private ! hides items not listed on public statement -public :: exchange_nod,exchange_elem,broadcast_nod,broadcast_elem, & - gather_nod, gather_elem, exchange_nod_begin, exchange_nod_end, exchange_elem_begin, & - exchange_elem_end, gather_edge -end module g_comm_auto diff --git a/src/temp/gen_modules_partitioning.F90 b/src/temp/gen_modules_partitioning.F90 deleted file mode 100644 index cc7d3c080..000000000 --- a/src/temp/gen_modules_partitioning.F90 +++ /dev/null @@ -1,508 +0,0 @@ -module par_support_interfaces - interface - subroutine par_init(partit) - USE o_PARAM - USE MOD_PARTIT - implicit none - type(t_partit), intent(inout), target :: partit - end subroutine - - subroutine par_ex(partit, abort) - USE MOD_PARTIT - implicit none - type(t_partit), intent(inout), target :: partit - integer,optional :: abort - end subroutine - - subroutine set_par_support(partit, mesh) - use MOD_MESH - use MOD_PARTIT - implicit none - type(t_partit), intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - end subroutine - - subroutine init_gatherLists(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - implicit none - type(t_mesh), intent(in), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine - end interface -end module - -subroutine par_init(partit) ! initializes MPI - USE o_PARAM - USE MOD_PARTIT - implicit none - type(t_partit), intent(inout), target :: partit - integer :: i - integer :: provided_mpi_thread_support_level - character(:), allocatable :: provided_mpi_thread_support_level_name - -#ifndef __oasis - call MPI_Comm_Size(MPI_COMM_WORLD,partit%npes,i) - call MPI_Comm_Rank(MPI_COMM_WORLD,partit%mype,i) - partit%MPI_COMM_FESOM=MPI_COMM_WORLD -#else - call MPI_Comm_Size(MPI_COMM_FESOM,partit%npes,i) - call MPI_Comm_Rank(MPI_COMM_FESOM,partit%mype,i) -#endif - - if(partit%mype==0) then - call MPI_Query_thread(provided_mpi_thread_support_level, i) - if(provided_mpi_thread_support_level == MPI_THREAD_SINGLE) then - provided_mpi_thread_support_level_name = "MPI_THREAD_SINGLE" - else if(provided_mpi_thread_support_level == MPI_THREAD_FUNNELED) then - provided_mpi_thread_support_level_name = "MPI_THREAD_FUNNELED" - else if(provided_mpi_thread_support_level == MPI_THREAD_SERIALIZED) then - provided_mpi_thread_support_level_name = "MPI_THREAD_SERIALIZED" - else if(provided_mpi_thread_support_level == MPI_THREAD_MULTIPLE) then - provided_mpi_thread_support_level_name = "MPI_THREAD_MULTIPLE" - else - provided_mpi_thread_support_level_name = "unknown" - end if - write(*,*) 'MPI has been initialized, provided MPI thread support level: ', & - provided_mpi_thread_support_level_name,provided_mpi_thread_support_level - write(*, *) 'Running on ', partit%npes, ' PEs' - end if -end subroutine par_init -!================================================================= -subroutine par_ex(partit, abort) ! finalizes MPI -USE MOD_PARTIT -#ifndef __oifs -!For standalone and coupled ECHAM runs -#if defined (__oasis) - use mod_prism -#endif - implicit none - type(t_partit), intent(inout), target :: partit - integer,optional :: abort - -#ifndef __oasis - if (present(abort)) then - if (partit%mype==0) write(*,*) 'Run finished unexpectedly!' - call MPI_ABORT(partit%MPI_COMM_FESOM, 1 ) - else - call MPI_Barrier(partit%MPI_COMM_FESOM,partit%MPIerr) - call MPI_Finalize(partit%MPIerr) - endif -#else - if (.not. present(abort)) then - if (partit%mype==0) print *, 'FESOM calls MPI_Barrier before calling prism_terminate' - call MPI_Barrier(MPI_COMM_WORLD, partit%MPIerr) - end if - call prism_terminate_proto(MPIerr) - if (partit%mype==0) print *, 'FESOM calls MPI_Barrier before calling MPI_Finalize' - call MPI_Barrier(MPI_COMM_WORLD, partit%MPIerr) - - if (partit%mype==0) print *, 'FESOM calls MPI_Finalize' - call MPI_Finalize(MPIerr) -#endif - if (partit%mype==0) print *, 'fesom should stop with exit status = 0' -#endif -#if defined (__oifs) -!OIFS coupling doesnt call prism_terminate_proto and uses MPI_COMM_FESOM - implicit none - integer,optional :: abort - if (present(abort)) then - if (partit%mype==0) write(*,*) 'Run finished unexpectedly!' - call MPI_ABORT( partit%MPI_COMM_FESOM, 1 ) - else - call MPI_Barrier(partit%MPI_COMM_FESOM,partit%MPIerr) - call MPI_Finalize(partit%MPIerr) - endif -#endif - -end subroutine par_ex -!======================================================================= -subroutine set_par_support(partit, mesh) - use MOD_MESH - use MOD_PARTIT - implicit none - - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - integer :: n, offset - integer :: i, max_nb, nb, nini, nend, nl1, n_val - integer, allocatable :: blocklen(:), displace(:) - integer, allocatable :: blocklen_tmp(:), displace_tmp(:) - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - ! - ! In the distributed memory version, most of the job is already done - ! at the initialization phase and is taken into account in read_mesh - ! routine. Here, MPI datatypes are built and buffers for MPI wait requests - ! are allocated. - - if (npes > 1) then - -!================================================ -! MPI REQUEST BUFFERS -!================================================ - if (.not. allocated(com_nod2D%req)) allocate(com_nod2D%req( 3*com_nod2D%rPEnum + 3*com_nod2D%sPEnum)) - if (.not. allocated(com_elem2D%req)) allocate(com_elem2D%req( 3*com_elem2D%rPEnum + 3*com_elem2D%sPEnum)) - if (.not. allocated(com_elem2D_full%req)) allocate(com_elem2D_full%req(3*com_elem2D_full%rPEnum + 3*com_elem2D_full%sPEnum)) -!================================================ -! MPI DATATYPES -!================================================ - ! Build MPI Data types for halo exchange: Elements - allocate(partit%r_mpitype_elem2D(com_elem2D%rPEnum,4)) ! 2D, small halo - allocate(partit%s_mpitype_elem2D(com_elem2D%sPEnum,4)) - allocate(partit%r_mpitype_elem2D_full_i(com_elem2D_full%rPEnum)) ! 2D, wide halo, integer - allocate(partit%s_mpitype_elem2D_full_i(com_elem2D_full%sPEnum)) - allocate(partit%r_mpitype_elem2D_full(com_elem2D_full%rPEnum,4)) ! 2D, wide halo - allocate(partit%s_mpitype_elem2D_full(com_elem2D_full%sPEnum,4)) - allocate(partit%r_mpitype_elem3D(com_elem2D%rPEnum, nl-1:nl,4)) ! 3D, small halo - allocate(partit%s_mpitype_elem3D(com_elem2D%sPEnum, nl-1:nl,4)) - allocate(partit%r_mpitype_elem3D_full(com_elem2D_full%rPEnum, nl-1:nl,4)) ! 3D, wide halo - allocate(partit%s_mpitype_elem3D_full(com_elem2D_full%sPEnum, nl-1:nl,4)) -!after the allocation we just reassotiate ALL pointers again here -#include "associate_part_ass.h" - ! Upper limit for the length of the local interface between the neighbor PEs - max_nb = max( & - maxval(com_elem2D%rptr(2:com_elem2D%rPEnum+1) - com_elem2D%rptr(1:com_elem2D%rPEnum)), & - maxval(com_elem2D%sptr(2:com_elem2D%sPEnum+1) - com_elem2D%sptr(1:com_elem2D%sPEnum)), & - maxval(com_elem2D_full%rptr(2:com_elem2D_full%rPEnum+1) - com_elem2D_full%rptr(1:com_elem2D_full%rPEnum)), & - maxval(com_elem2D_full%sptr(2:com_elem2D_full%sPEnum+1) - com_elem2D_full%sptr(1:com_elem2D_full%sPEnum))) - - allocate(displace(max_nb), blocklen(max_nb)) - allocate(displace_tmp(max_nb), blocklen_tmp(max_nb)) - - - do n=1,com_elem2D%rPEnum - nb = 1 - nini = com_elem2D%rptr(n) - nend = com_elem2D%rptr(n+1) - 1 - displace(:) = 0 - displace(1) = com_elem2D%rlist(nini) -1 ! C counting, start at 0 - blocklen(:) = 1 - do i=nini+1, nend - if (com_elem2D%rlist(i) /= com_elem2D%rlist(i-1) + 1) then - ! New block - nb = nb+1 - displace(nb) = com_elem2D%rlist(i) -1 - else - blocklen(nb) = blocklen(nb)+1 - endif - enddo - - DO n_val=1,4 - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val - displace_tmp(1:nb) = displace(1:nb)*n_val - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - r_mpitype_elem2D(n,n_val), MPIerr) - - call MPI_TYPE_COMMIT(r_mpitype_elem2D(n,n_val), MPIerr) - - DO nl1=nl-1, nl - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 - displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - r_mpitype_elem3D(n,nl1,n_val), MPIerr) - - call MPI_TYPE_COMMIT(r_mpitype_elem3D(n,nl1,n_val), MPIerr) - ENDDO - ENDDO - enddo - - do n=1,com_elem2D%sPEnum - nb = 1 - nini = com_elem2D%sptr(n) - nend = com_elem2D%sptr(n+1) - 1 - displace(:) = 0 - displace(1) = com_elem2D%slist(nini) -1 ! C counting, start at 0 - blocklen(:) = 1 - do i=nini+1, nend - if (com_elem2D%slist(i) /= com_elem2D%slist(i-1) + 1) then - ! New block - nb = nb+1 - displace(nb) = com_elem2D%slist(i) -1 - else - blocklen(nb) = blocklen(nb)+1 - endif - enddo - - DO n_val=1,4 - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val - displace_tmp(1:nb) = displace(1:nb)*n_val - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - s_mpitype_elem2D(n, n_val), MPIerr) - - call MPI_TYPE_COMMIT(s_mpitype_elem2D(n, n_val), MPIerr) - - DO nl1=nl-1, nl - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 - displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - s_mpitype_elem3D(n,nl1,n_val), MPIerr) - - call MPI_TYPE_COMMIT(s_mpitype_elem3D(n,nl1,n_val), MPIerr) - ENDDO - ENDDO - enddo - - do n=1,com_elem2D_full%rPEnum - nb = 1 - nini = com_elem2D_full%rptr(n) - nend = com_elem2D_full%rptr(n+1) - 1 - displace(:) = 0 - displace(1) = com_elem2D_full%rlist(nini) -1 ! C counting, start at 0 - blocklen(:) = 1 - do i=nini+1, nend - if (com_elem2D_full%rlist(i) /= com_elem2D_full%rlist(i-1) + 1) then - ! New block - nb = nb+1 - displace(nb) = com_elem2D_full%rlist(i) -1 - else - blocklen(nb) = blocklen(nb)+1 - endif - enddo - - call MPI_TYPE_INDEXED(nb, blocklen,displace,MPI_INTEGER, r_mpitype_elem2D_full_i(n),MPIerr) - - call MPI_TYPE_COMMIT(r_mpitype_elem2D_full_i(n), MPIerr) - - DO n_val=1,4 - - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & - r_mpitype_elem2D_full(n,n_val), MPIerr) - call MPI_TYPE_COMMIT(r_mpitype_elem2D_full(n, n_val), MPIerr) - - DO nl1=nl-1, nl - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 - displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - r_mpitype_elem3D_full(n,nl1,n_val), MPIerr) - - call MPI_TYPE_COMMIT(r_mpitype_elem3D_full(n,nl1,n_val), MPIerr) - ENDDO - ENDDO - enddo - - do n=1,com_elem2D_full%sPEnum - nb = 1 - nini = com_elem2D_full%sptr(n) - nend = com_elem2D_full%sptr(n+1) - 1 - displace(:) = 0 - displace(1) = com_elem2D_full%slist(nini) -1 ! C counting, start at 0 - blocklen(:) = 1 - do i=nini+1, nend - if (com_elem2D_full%slist(i) /= com_elem2D_full%slist(i-1) + 1) then - ! New block - nb = nb+1 - displace(nb) = com_elem2D_full%slist(i) -1 - else - blocklen(nb) = blocklen(nb)+1 - endif - enddo - - call MPI_TYPE_INDEXED(nb, blocklen,displace,MPI_INTEGER, s_mpitype_elem2D_full_i(n), MPIerr) - - call MPI_TYPE_COMMIT(s_mpitype_elem2D_full_i(n), MPIerr) - - DO n_val=1,4 - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & - s_mpitype_elem2D_full(n,n_val),MPIerr) - call MPI_TYPE_COMMIT(s_mpitype_elem2D_full(n,n_val), MPIerr) - - DO nl1=nl-1, nl - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 - displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - s_mpitype_elem3D_full(n,nl1,n_val), MPIerr) - - call MPI_TYPE_COMMIT(s_mpitype_elem3D_full(n,nl1,n_val), MPIerr) - ENDDO - ENDDO - enddo - - deallocate(displace, blocklen) - deallocate(displace_tmp, blocklen_tmp) - - - ! Build MPI Data types for halo exchange: Nodes - - allocate(partit%r_mpitype_nod2D(com_nod2D%rPEnum)) ! 2D - allocate(partit%s_mpitype_nod2D(com_nod2D%sPEnum)) - allocate(partit%r_mpitype_nod2D_i(com_nod2D%rPEnum)) ! 2D integer - allocate(partit%s_mpitype_nod2D_i(com_nod2D%sPEnum)) - - allocate(partit%r_mpitype_nod3D(com_nod2D%rPEnum,nl-1:nl,3)) ! 3D with nl-1 or nl layers, 1-3 values - allocate(partit%s_mpitype_nod3D(com_nod2D%sPEnum,nl-1:nl,3)) -!after the allocation we just reassotiate ALL pointers again here -#include "associate_part_ass.h" - - ! Upper limit for the length of the local interface between the neighbor PEs - max_nb = max(maxval(com_nod2D%rptr(2:com_nod2D%rPEnum+1) - com_nod2D%rptr(1:com_nod2D%rPEnum)), & - maxval(com_nod2D%sptr(2:com_nod2D%sPEnum+1) - com_nod2D%sptr(1:com_nod2D%sPEnum))) - - allocate(displace(max_nb), blocklen(max_nb)) - allocate(displace_tmp(max_nb), blocklen_tmp(max_nb)) - - do n=1,com_nod2D%rPEnum - nb = 1 - nini = com_nod2D%rptr(n) - nend = com_nod2D%rptr(n+1) - 1 - displace(:) = 0 - displace(1) = com_nod2D%rlist(nini) -1 ! C counting, start at 0 - blocklen(:) = 1 - do i=nini+1, nend - if (com_nod2D%rlist(i) /= com_nod2D%rlist(i-1) + 1) then - ! New block - nb = nb+1 - displace(nb) = com_nod2D%rlist(i) -1 - else - blocklen(nb) = blocklen(nb)+1 - endif - enddo - - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & - r_mpitype_nod2D(n), MPIerr) - - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_INTEGER, & - r_mpitype_nod2D_i(n), MPIerr) - - call MPI_TYPE_COMMIT(r_mpitype_nod2D(n), MPIerr) - call MPI_TYPE_COMMIT(r_mpitype_nod2D_i(n), MPIerr) - - DO nl1=nl-1, nl - DO n_val=1,3 - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 - displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - r_mpitype_nod3D(n,nl1,n_val), MPIerr) - - call MPI_TYPE_COMMIT(r_mpitype_nod3D(n,nl1,n_val), MPIerr) - ENDDO - ENDDO - enddo - - do n=1,com_nod2D%sPEnum - nb = 1 - nini = com_nod2D%sptr(n) - nend = com_nod2D%sptr(n+1) - 1 - displace(:) = 0 - displace(1) = com_nod2D%slist(nini) -1 ! C counting, start at 0 - blocklen(:) = 1 - do i=nini+1, nend - if (com_nod2D%slist(i) /= com_nod2D%slist(i-1) + 1) then - ! New block - nb = nb+1 - displace(nb) = com_nod2D%slist(i) -1 - else - blocklen(nb) = blocklen(nb)+1 - endif - enddo - - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & - s_mpitype_nod2D(n), MPIerr) - - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_INTEGER, & - s_mpitype_nod2D_i(n), MPIerr) - - call MPI_TYPE_COMMIT(s_mpitype_nod2D(n), MPIerr) - call MPI_TYPE_COMMIT(s_mpitype_nod2D_i(n), MPIerr) - - DO nl1=nl-1, nl - DO n_val=1,3 - - blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 - displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & - s_mpitype_nod3D(n,nl1,n_val), MPIerr) - - call MPI_TYPE_COMMIT(s_mpitype_nod3D(n,nl1,n_val), MPIerr) - ENDDO - ENDDO - enddo - - deallocate(blocklen, displace) - deallocate(blocklen_tmp, displace_tmp) - - endif - - call init_gatherLists(partit, mesh) - if(mype==0) write(*,*) 'Communication arrays are set' -end subroutine set_par_support - - -!=================================================================== -subroutine init_gatherLists(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - implicit none - type(t_mesh), intent(in), target :: mesh - type(t_partit), intent(inout), target :: partit - integer :: n2D, e2D, sum_loc_elem2D - integer :: n, estart, nstart -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - if (mype==0) then - - if (npes > 1) then - - allocate(partit%remPtr_nod2D(npes)) - allocate(partit%remPtr_elem2D(npes)) -!reassociate the pointers to the just allocated arrays -#include "associate_part_ass.h" - remPtr_nod2D(1) = 1 - remPtr_elem2D(1) = 1 - - do n=1, npes-1 - call MPI_RECV(n2D, 1, MPI_INTEGER, n, 0, MPI_COMM_FESOM, MPI_STATUS_IGNORE, MPIerr ) - call MPI_RECV(e2D, 1, MPI_INTEGER, n, 1, MPI_COMM_FESOM, MPI_STATUS_IGNORE, MPIerr ) - - remPtr_nod2D(n+1) = remPtr_nod2D(n) + n2D - remPtr_elem2D(n+1) = remPtr_elem2D(n) + e2D - enddo - - allocate(partit%remList_nod2D(remPtr_nod2D(npes))) ! this should be nod2D - myDim_nod2D - allocate(partit%remList_elem2D(remPtr_elem2D(npes))) ! this is > elem2D, because the elements overlap. - ! Consider optimization: avoid multiple communication - ! of the same elem from different PEs. -!reassociate the pointers to the just allocated arrays -#include "associate_part_ass.h" - - do n=1, npes-1 - nstart = remPtr_nod2D(n) - n2D = remPtr_nod2D(n+1) - remPtr_nod2D(n) - call MPI_RECV(remList_nod2D(nstart), n2D, MPI_INTEGER, n, 2, MPI_COMM_FESOM, & - MPI_STATUS_IGNORE, MPIerr ) - estart = remPtr_elem2D(n) - e2D = remPtr_elem2D(n+1) - remPtr_elem2D(n) - call MPI_RECV(remList_elem2D(estart),e2D, MPI_INTEGER, n, 3, MPI_COMM_FESOM, & - MPI_STATUS_IGNORE, MPIerr ) - - enddo - end if - else - - call MPI_SEND(myDim_nod2D, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) - call MPI_SEND(myDim_elem2D, 1, MPI_INTEGER, 0, 1, MPI_COMM_FESOM, MPIerr ) - call MPI_SEND(myList_nod2D, myDim_nod2D, MPI_INTEGER, 0, 2, MPI_COMM_FESOM, MPIerr ) - call MPI_SEND(myList_elem2D, myDim_elem2D, MPI_INTEGER, 0, 3, MPI_COMM_FESOM, MPIerr ) - - endif -end subroutine init_gatherLists diff --git a/src/temp/oce_adv_tra_driver.F90 b/src/temp/oce_adv_tra_driver.F90 deleted file mode 100644 index 511e903b7..000000000 --- a/src/temp/oce_adv_tra_driver.F90 +++ /dev/null @@ -1,278 +0,0 @@ -module oce_adv_tra_driver_interfaces - interface - subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, tracers, partit, mesh) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - real(kind=WP), intent(in), target :: dt - integer, intent(in) :: tr_num - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - type(t_tracer), intent(inout), target :: tracers - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(in), target :: W(mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in), target :: WI(mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in), target :: WE(mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - end subroutine - end interface -end module - -module oce_tra_adv_flux2dtracer_interface - interface - subroutine oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, flux_h, flux_v, partit, mesh, use_lo, ttf, lo) - !update the solution for vertical and horizontal flux contributions - use MOD_MESH - use MOD_PARTIT - real(kind=WP), intent(in), target :: dt - type(t_partit),intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(inout) :: dttf_h(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: dttf_v(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux_h(mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(inout) :: flux_v(mesh%nl, partit%myDim_nod2D) - logical, optional :: use_lo - real(kind=WP), optional :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), optional :: lo (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - end subroutine - end interface -end module -! -! -!=============================================================================== -subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, tracers, partit, mesh) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - use oce_adv_tra_hor_interfaces - use oce_adv_tra_ver_interfaces - use oce_adv_tra_fct_interfaces - use oce_tra_adv_flux2dtracer_interface - implicit none - real(kind=WP), intent(in), target :: dt - integer, intent(in) :: tr_num - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - type(t_tracer), intent(inout), target :: tracers - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(in), target :: W(mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in), target :: WI(mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in), target :: WE(mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - - real(kind=WP), pointer, dimension (:,:) :: pwvel - real(kind=WP), pointer, dimension (:,:) :: ttf, ttfAB, fct_LO - real(kind=WP), pointer, dimension (:,:) :: adv_flux_hor, adv_flux_ver, dttf_h, dttf_v - real(kind=WP), pointer, dimension (:,:) :: fct_ttf_min, fct_ttf_max - real(kind=WP), pointer, dimension (:,:) :: fct_plus, fct_minus - - integer, pointer, dimension (:) :: nboundary_lay - real(kind=WP), pointer, dimension (:,:,:) :: edge_up_dn_grad - - integer :: el(2), enodes(2), nz, n, e - integer :: nl12, nu12, nl1, nl2, nu1, nu2 - real(kind=WP) :: cLO, cHO, deltaX1, deltaY1, deltaX2, deltaY2 - real(kind=WP) :: qc, qu, qd - real(kind=WP) :: tvert(mesh%nl), tvert_e(mesh%nl), a, b, c, d, da, db, dg, vflux, Tupw1 - real(kind=WP) :: Tmean, Tmean1, Tmean2, num_ord - real(kind=WP) :: opth, optv - logical :: do_zero_flux - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - ttf => tracers%data(tr_num)%values - ttfAB => tracers%data(tr_num)%valuesAB - opth = tracers%data(tr_num)%tra_adv_ph - optv = tracers%data(tr_num)%tra_adv_pv - fct_LO => tracers%work%fct_LO - adv_flux_ver => tracers%work%adv_flux_ver - adv_flux_hor => tracers%work%adv_flux_hor - edge_up_dn_grad => tracers%work%edge_up_dn_grad - nboundary_lay => tracers%work%nboundary_lay - fct_ttf_min => tracers%work%fct_ttf_min - fct_ttf_max => tracers%work%fct_ttf_max - fct_plus => tracers%work%fct_plus - fct_minus => tracers%work%fct_minus - dttf_h => tracers%work%del_ttf_advhoriz - dttf_v => tracers%work%del_ttf_advvert - !___________________________________________________________________________ - ! compute FCT horzontal and vertical low order solution as well as lw order - ! part of antidiffusive flux - if (trim(tracers%data(tr_num)%tra_adv_lim)=='FCT') then - ! compute the low order upwind horizontal flux - ! init_zero=.true. : zero the horizontal flux before computation - ! init_zero=.false. : input flux will be substracted - call adv_tra_hor_upw1(vel, ttf, partit, mesh, adv_flux_hor, init_zero=.true.) - ! update the LO solution for horizontal contribution - fct_LO=0.0_WP - do e=1, myDim_edge2D - enodes=edges(:,e) - el=edge_tri(:,e) - nl1=nlevels(el(1))-1 - nu1=ulevels(el(1)) - nl2=0 - nu2=0 - if(el(2)>0) then - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - end if - - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) - - !!PS do nz=1, max(nl1, nl2) - do nz=nu12, nl12 - fct_LO(nz, enodes(1))=fct_LO(nz, enodes(1))+adv_flux_hor(nz, e) - fct_LO(nz, enodes(2))=fct_LO(nz, enodes(2))-adv_flux_hor(nz, e) - end do - end do - ! compute the low order upwind vertical flux (explicit part only) - ! zero the input/output flux before computation - call adv_tra_ver_upw1(we, ttf, partit, mesh, adv_flux_ver, init_zero=.true.) - ! update the LO solution for vertical contribution - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !!PS do nz=1, nlevels_nod2D(n)-1 - do nz= nu1, nl1-1 - fct_LO(nz,n)=(ttf(nz,n)*hnode(nz,n)+(fct_LO(nz,n)+(adv_flux_ver(nz, n)-adv_flux_ver(nz+1, n)))*dt/areasvol(nz,n))/hnode_new(nz,n) - end do - end do - if (w_split) then !wvel/=wvel_e - ! update for implicit contribution (w_split option) - call adv_tra_vert_impl(dt, wi, fct_LO, partit, mesh) - ! compute the low order upwind vertical flux (full vertical velocity) - ! zero the input/output flux before computation - ! --> compute here low order part of vertical anti diffusive fluxes, - ! has to be done on the full vertical velocity w - call adv_tra_ver_upw1(w, ttf, partit, mesh, adv_flux_ver, init_zero=.true.) - end if - call exchange_nod(fct_LO, partit) - end if - - do_zero_flux=.true. - if (trim(tracers%data(tr_num)%tra_adv_lim)=='FCT') do_zero_flux=.false. - !___________________________________________________________________________ - ! do horizontal tracer advection, in case of FCT high order solution - SELECT CASE(trim(tracers%data(tr_num)%tra_adv_hor)) - CASE('MUSCL') - ! compute the untidiffusive horizontal flux (init_zero=.false.: input is the LO horizontal flux computed above) - call adv_tra_hor_muscl(vel, ttfAB, partit, mesh, opth, adv_flux_hor, edge_up_dn_grad, nboundary_lay, init_zero=do_zero_flux) - CASE('MFCT') - call adv_tra_hor_mfct(vel, ttfAB, partit, mesh, opth, adv_flux_hor, edge_up_dn_grad, init_zero=do_zero_flux) - CASE('UPW1') - call adv_tra_hor_upw1(vel, ttfAB, partit, mesh, adv_flux_hor, init_zero=do_zero_flux) - CASE DEFAULT !unknown - if (mype==0) write(*,*) 'Unknown horizontal advection type ', trim(tracers%data(tr_num)%tra_adv_hor), '! Check your namelists!' - call par_ex(partit, 1) - END SELECT - if (trim(tracers%data(tr_num)%tra_adv_lim)=='FCT') then - pwvel=>w - else - pwvel=>we - end if - !___________________________________________________________________________ - ! do vertical tracer advection, in case of FCT high order solution - SELECT CASE(trim(tracers%data(tr_num)%tra_adv_ver)) - CASE('QR4C') - ! compute the untidiffusive vertical flux (init_zero=.false.:input is the LO vertical flux computed above) - call adv_tra_ver_qr4c ( pwvel, ttfAB, partit, mesh, optv, adv_flux_ver, init_zero=do_zero_flux) - CASE('CDIFF') - call adv_tra_ver_cdiff( pwvel, ttfAB, partit, mesh, adv_flux_ver, init_zero=do_zero_flux) - CASE('PPM') - call adv_tra_vert_ppm(dt, pwvel, ttfAB, partit, mesh, adv_flux_ver, init_zero=do_zero_flux) - CASE('UPW1') - call adv_tra_ver_upw1 ( pwvel, ttfAB, partit, mesh, adv_flux_ver, init_zero=do_zero_flux) - CASE DEFAULT !unknown - if (mype==0) write(*,*) 'Unknown vertical advection type ', trim(tracers%data(tr_num)%tra_adv_ver), '! Check your namelists!' - call par_ex(1) - ! --> be aware the vertical implicite part in case without FCT is done in - ! oce_ale_tracer.F90 --> subroutine diff_ver_part_impl_ale(tr_num, partit, mesh) - ! for do_wimpl=.true. - END SELECT - !___________________________________________________________________________ - ! - if (trim(tracers%data(tr_num)%tra_adv_lim)=='FCT') then - !edge_up_dn_grad will be used as an auxuary array here - call oce_tra_adv_fct(dt, ttf, fct_LO, adv_flux_hor, adv_flux_ver, fct_ttf_min, fct_ttf_max, fct_plus, fct_minus, edge_up_dn_grad, partit, mesh) - call oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, adv_flux_hor, adv_flux_ver, partit, mesh, use_lo=.TRUE., ttf=ttf, lo=fct_LO) - else - call oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, adv_flux_hor, adv_flux_ver, partit, mesh) - end if -end subroutine do_oce_adv_tra -! -! -!=============================================================================== -subroutine oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, flux_h, flux_v, partit, mesh, use_lo, ttf, lo) - use MOD_MESH - use o_ARRAYS - use MOD_PARTIT - use g_comm_auto - implicit none - real(kind=WP), intent(in), target :: dt - type(t_partit),intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(inout) :: dttf_h(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: dttf_v(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux_h(mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(inout) :: flux_v(mesh%nl, partit%myDim_nod2D) - logical, optional :: use_lo - real(kind=WP), optional :: lo (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), optional :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - integer :: n, nz, k, elem, enodes(3), num, el(2), nu12, nl12, nu1, nu2, nl1, nl2, edge -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - !___________________________________________________________________________ - ! c. Update the solution - ! Vertical - if (present(use_lo)) then - if (use_lo) then - do n=1, myDim_nod2d - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !!PS do nz=1,nlevels_nod2D(n)-1 - do nz=nu1, nl1-1 - dttf_v(nz,n)=dttf_v(nz,n)-ttf(nz,n)*hnode(nz,n)+LO(nz,n)*hnode_new(nz,n) - end do - end do - end if - end if - - do n=1, myDim_nod2d - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - do nz=nu1,nl1-1 - dttf_v(nz,n)=dttf_v(nz,n) + (flux_v(nz,n)-flux_v(nz+1,n))*dt/areasvol(nz,n) - end do - end do - - - ! Horizontal - do edge=1, myDim_edge2D - enodes(1:2)=edges(:,edge) - el=edge_tri(:,edge) - nl1=nlevels(el(1))-1 - nu1=ulevels(el(1)) - - nl2=0 - nu2=0 - if(el(2)>0) then - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - end if - - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) - - !!PS do nz=1, max(nl1, nl2) - do nz=nu12, nl12 - dttf_h(nz,enodes(1))=dttf_h(nz,enodes(1))+flux_h(nz,edge)*dt/areasvol(nz,enodes(1)) - dttf_h(nz,enodes(2))=dttf_h(nz,enodes(2))-flux_h(nz,edge)*dt/areasvol(nz,enodes(2)) - end do - end do -end subroutine oce_tra_adv_flux2dtracer diff --git a/src/temp/oce_adv_tra_fct.F90 b/src/temp/oce_adv_tra_fct.F90 deleted file mode 100644 index 5eb7993a9..000000000 --- a/src/temp/oce_adv_tra_fct.F90 +++ /dev/null @@ -1,365 +0,0 @@ -module oce_adv_tra_fct_interfaces - interface - subroutine oce_adv_tra_fct_init(twork, partit, mesh) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - type(t_mesh), intent(in), target :: mesh - type(t_partit),intent(inout), target :: partit - type(t_tracer_work), intent(inout), target :: twork - end subroutine - - subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, fct_plus, fct_minus, AUX, partit, mesh) - use MOD_MESH - use MOD_PARTIT - real(kind=WP), intent(in), target :: dt - type(t_partit),intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(inout) :: fct_ttf_min(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: fct_ttf_max(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: lo (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: adf_h(mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(inout) :: adf_v(mesh%nl, partit%myDim_nod2D) - real(kind=WP), intent(inout) :: fct_plus(mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(inout) :: fct_minus(mesh%nl, partit%myDim_nod2D) - real(kind=WP), intent(inout) :: AUX(:,:,:) !a large auxuary array - end subroutine - end interface -end module -! -! -!=============================================================================== -subroutine oce_adv_tra_fct_init(twork, partit, mesh) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - implicit none - integer :: my_size - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer_work), intent(inout), target :: twork -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - my_size=myDim_nod2D+eDim_nod2D - allocate(twork%fct_LO(nl-1, my_size)) ! Low-order solution - allocate(twork%adv_flux_hor(nl-1,partit%myDim_edge2D)) ! antidiffusive hor. contributions / from edges - allocate(twork%adv_flux_ver(nl, partit%myDim_nod2D)) ! antidiffusive ver. fluxes / from nodes - - allocate(twork%fct_ttf_max(nl-1, my_size),twork%fct_ttf_min(nl-1, my_size)) - allocate(twork%fct_plus(nl-1, my_size), twork%fct_minus(nl-1, my_size)) - ! Initialize with zeros: - twork%fct_LO=0.0_WP - twork%adv_flux_hor=0.0_WP - twork%adv_flux_ver=0.0_WP - twork%fct_ttf_max=0.0_WP - twork%fct_ttf_min=0.0_WP - twork%fct_plus=0.0_WP - twork%fct_minus=0.0_WP - - if (mype==0) write(*,*) 'FCT is initialized' -end subroutine oce_adv_tra_fct_init - -! -! -!=============================================================================== -subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, fct_plus, fct_minus, AUX, partit, mesh) - ! - ! 3D Flux Corrected Transport scheme - ! Limits antidiffusive fluxes==the difference in flux HO-LO - ! LO ==Low-order (first-order upwind) - ! HO ==High-order (3rd/4th order gradient reconstruction method) - ! Adds limited fluxes to the LO solution - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - implicit none - real(kind=WP), intent(in), target :: dt - type(t_mesh), intent(in), target :: mesh - type(t_partit),intent(inout), target :: partit - real(kind=WP), intent(inout) :: fct_ttf_min(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: fct_ttf_max(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: lo (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: adf_h(mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(inout) :: adf_v(mesh%nl, partit%myDim_nod2D) - real(kind=WP), intent(inout) :: fct_plus (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: fct_minus(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: AUX(:,:,:) !a large auxuary array, let us use twork%edge_up_dn_grad(1:4, 1:NL-2, 1:partit%myDim_edge2D) to save space - integer :: n, nz, k, elem, enodes(3), num, el(2), nl1, nl2, nu1, nu2, nl12, nu12, edge - real(kind=WP) :: flux, ae,tvert_max(mesh%nl-1),tvert_min(mesh%nl-1) - real(kind=WP) :: flux_eps=1e-16 - real(kind=WP) :: bignumber=1e3 - integer :: vlimit=1 - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - ! -------------------------------------------------------------------------- - ! ttf is the tracer field on step n - ! del_ttf is the increment - ! vlimit sets the version of limiting, see below - ! -------------------------------------------------------------------------- - !___________________________________________________________________________ - ! a1. max, min between old solution and updated low-order solution per node - do n=1,myDim_nod2D + edim_nod2d - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - do nz=nu1, nl1-1 - fct_ttf_max(nz,n)=max(LO(nz,n), ttf(nz,n)) - fct_ttf_min(nz,n)=min(LO(nz,n), ttf(nz,n)) - end do - end do - - !___________________________________________________________________________ - ! a2. Admissible increments on elements - ! (only layers below the first and above the last layer) - ! look for max, min bounds for each element --> AUX here auxilary array - do elem=1, myDim_elem2D - enodes=elem2D_nodes(:,elem) - nu1 = ulevels(elem) - nl1 = nlevels(elem) - do nz=nu1, nl1-1 - AUX(1,nz,elem)=maxval(fct_ttf_max(nz,enodes)) - AUX(2,nz,elem)=minval(fct_ttf_min(nz,enodes)) - end do - if (nl1<=nl-1) then - do nz=nl1,nl-1 - AUX(1,nz,elem)=-bignumber - AUX(2,nz,elem)= bignumber - end do - endif - end do ! --> do elem=1, myDim_elem2D - - !___________________________________________________________________________ - ! a3. Bounds on clusters and admissible increments - ! Vertical1: In this version we look at the bounds on the clusters - ! above and below, which leaves wide bounds because typically - ! vertical gradients are larger. - if(vlimit==1) then - !Horizontal - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - - !___________________________________________________________________ - do nz=nu1,nl1-1 - ! max,min horizontal bound in cluster around node n in every - ! vertical layer - ! nod_in_elem2D --> elem indices of which node n is surrounded - ! nod_in_elem2D_num --> max number of surrounded elem - tvert_max(nz)= maxval(AUX(1,nz,nod_in_elem2D(1:nod_in_elem2D_num(n),n))) - tvert_min(nz)= minval(AUX(2,nz,nod_in_elem2D(1:nod_in_elem2D_num(n),n))) - end do - - !___________________________________________________________________ - ! calc max,min increment of surface layer with respect to low order - ! solution - fct_ttf_max(nu1,n)=tvert_max(nu1)-LO(nu1,n) - fct_ttf_min(nu1,n)=tvert_min(nu1)-LO(nu1,n) - - ! calc max,min increment from nz-1:nz+1 with respect to low order - ! solution at layer nz - do nz=nu1+1,nl1-2 - fct_ttf_max(nz,n)=maxval(tvert_max(nz-1:nz+1))-LO(nz,n) - fct_ttf_min(nz,n)=minval(tvert_min(nz-1:nz+1))-LO(nz,n) - end do - ! calc max,min increment of bottom layer -1 with respect to low order - ! solution - nz=nl1-1 - fct_ttf_max(nz,n)=tvert_max(nz)-LO(nz,n) - fct_ttf_min(nz,n)=tvert_min(nz)-LO(nz,n) - end do - end if - - !___________________________________________________________________________ - ! Vertical2: Similar to the version above, but the vertical bounds are more - ! local - if(vlimit==2) then - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - do nz=nu1,nl1-1 - tvert_max(nz)= maxval(AUX(1,nz,nod_in_elem2D(1:nod_in_elem2D_num(n),n))) - tvert_min(nz)= minval(AUX(2,nz,nod_in_elem2D(1:nod_in_elem2D_num(n),n))) - end do - do nz=nu1+1, nl1-2 - tvert_max(nz)=max(tvert_max(nz),maxval(fct_ttf_max(nz-1:nz+1,n))) - tvert_min(nz)=min(tvert_min(nz),minval(fct_ttf_max(nz-1:nz+1,n))) - end do - do nz=nu1,nl1-1 - fct_ttf_max(nz,n)=tvert_max(nz)-LO(nz,n) - fct_ttf_min(nz,n)=tvert_min(nz)-LO(nz,n) - end do - end do - end if - - !___________________________________________________________________________ - ! Vertical3: Vertical bounds are taken into account only if they are narrower than the - ! horizontal ones - if(vlimit==3) then - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - do nz=nu1, nl1-1 - tvert_max(nz)= maxval(AUX(1,nz,nod_in_elem2D(1:nod_in_elem2D_num(n),n))) - tvert_min(nz)= minval(AUX(2,nz,nod_in_elem2D(1:nod_in_elem2D_num(n),n))) - end do - do nz=nu1+1, nl1-2 - tvert_max(nz)=min(tvert_max(nz),maxval(fct_ttf_max(nz-1:nz+1,n))) - tvert_min(nz)=max(tvert_min(nz),minval(fct_ttf_max(nz-1:nz+1,n))) - end do - do nz=nu1, nl1-1 - fct_ttf_max(nz,n)=tvert_max(nz)-LO(nz,n) - fct_ttf_min(nz,n)=tvert_min(nz)-LO(nz,n) - end do - end do - end if - - !___________________________________________________________________________ - ! b1. Split positive and negative antidiffusive contributions - ! --> sum all positive (fct_plus), negative (fct_minus) antidiffusive - ! horizontal element and vertical node contribution to node n and layer nz - ! see. R. Löhner et al. "finite element flux corrected transport (FEM-FCT) - ! for the euler and navier stoke equation - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - do nz=nu1,nl1-1 - fct_plus(nz,n)=0._WP - fct_minus(nz,n)=0._WP - end do - end do - - !Vertical - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - do nz=nu1,nl1-1 -! fct_plus(nz,n)=fct_plus(nz,n)+ & -! (max(0.0_WP,adf_v(nz,n))+max(0.0_WP,-adf_v(nz+1,n))) & -! /hnode(nz,n) -! fct_minus(nz,n)=fct_minus(nz,n)+ & -! (min(0.0_WP,adf_v(nz,n))+min(0.0_WP,-adf_v(nz+1,n))) & -! /hnode(nz,n) - fct_plus(nz,n) =fct_plus(nz,n) +(max(0.0_WP,adf_v(nz,n))+max(0.0_WP,-adf_v(nz+1,n))) - fct_minus(nz,n)=fct_minus(nz,n)+(min(0.0_WP,adf_v(nz,n))+min(0.0_WP,-adf_v(nz+1,n))) - end do - end do - - !Horizontal - do edge=1, myDim_edge2D - enodes(1:2)=edges(:,edge) - el=edge_tri(:,edge) - nl1=nlevels(el(1))-1 - nu1=ulevels(el(1)) - nl2=0 - nu2=0 - if(el(2)>0) then - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - end if - - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) - - do nz=nu12, nl12 - fct_plus (nz,enodes(1))=fct_plus (nz,enodes(1)) + max(0.0_WP, adf_h(nz,edge)) - fct_minus(nz,enodes(1))=fct_minus(nz,enodes(1)) + min(0.0_WP, adf_h(nz,edge)) - fct_plus (nz,enodes(2))=fct_plus (nz,enodes(2)) + max(0.0_WP,-adf_h(nz,edge)) - fct_minus(nz,enodes(2))=fct_minus(nz,enodes(2)) + min(0.0_WP,-adf_h(nz,edge)) - end do - end do - - !___________________________________________________________________________ - ! b2. Limiting factors - do n=1,myDim_nod2D - nu1=ulevels_nod2D(n) - nl1=nlevels_nod2D(n) - do nz=nu1,nl1-1 - flux=fct_plus(nz,n)*dt/areasvol(nz,n)+flux_eps - fct_plus(nz,n)=min(1.0_WP,fct_ttf_max(nz,n)/flux) - flux=fct_minus(nz,n)*dt/areasvol(nz,n)-flux_eps - fct_minus(nz,n)=min(1.0_WP,fct_ttf_min(nz,n)/flux) - end do - end do - - ! fct_minus and fct_plus must be known to neighbouring PE - call exchange_nod(fct_plus, fct_minus, partit) - - !___________________________________________________________________________ - ! b3. Limiting - !Vertical - do n=1, myDim_nod2D - nu1=ulevels_nod2D(n) - nl1=nlevels_nod2D(n) - - !_______________________________________________________________________ - nz=nu1 - ae=1.0_WP - flux=adf_v(nz,n) - if(flux>=0.0_WP) then - ae=min(ae,fct_plus(nz,n)) - else - ae=min(ae,fct_minus(nz,n)) - end if - adf_v(nz,n)=ae*adf_v(nz,n) - - !_______________________________________________________________________ - do nz=nu1+1,nl1-1 - ae=1.0_WP - flux=adf_v(nz,n) - if(flux>=0._WP) then - ae=min(ae,fct_minus(nz-1,n)) - ae=min(ae,fct_plus(nz,n)) - else - ae=min(ae,fct_plus(nz-1,n)) - ae=min(ae,fct_minus(nz,n)) - end if - adf_v(nz,n)=ae*adf_v(nz,n) - end do - ! the bottom flux is always zero - end do - - call exchange_nod_end(partit) ! fct_plus, fct_minus - - !Horizontal - do edge=1, myDim_edge2D - enodes(1:2)=edges(:,edge) - el=edge_tri(:,edge) - nu1=ulevels(el(1)) - nl1=nlevels(el(1))-1 - nl2=0 - nu2=0 - if(el(2)>0) then - nu2=ulevels(el(2)) - nl2=nlevels(el(2))-1 - end if - - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) - - do nz=nu12, nl12 - ae=1.0_WP - flux=adf_h(nz,edge) - - if(flux>=0._WP) then - ae=min(ae,fct_plus(nz,enodes(1))) - ae=min(ae,fct_minus(nz,enodes(2))) - else - ae=min(ae,fct_minus(nz,enodes(1))) - ae=min(ae,fct_plus(nz,enodes(2))) - endif - - adf_h(nz,edge)=ae*adf_h(nz,edge) - end do - end do -end subroutine oce_tra_adv_fct diff --git a/src/temp/oce_adv_tra_hor.F90 b/src/temp/oce_adv_tra_hor.F90 deleted file mode 100644 index 714eccf68..000000000 --- a/src/temp/oce_adv_tra_hor.F90 +++ /dev/null @@ -1,739 +0,0 @@ -!=============================================================================================================================== -!**************** routines for horizontal tracer advection *********************** -module oce_adv_tra_hor_interfaces - interface -! (low order upwind) -! returns flux given at edges which contributes with -! plus sign into 1st. node and with the minus sign into the 2nd node -! IF init_zero=.TRUE. : flux will be set to zero before computation -! IF init_zero=.FALSE. : flux=flux-input flux -! flux is not multiplied with dt - subroutine adv_tra_hor_upw1(vel, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: ttf( mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(inout) :: flux( mesh%nl-1, partit%myDim_edge2D) - logical, optional :: init_zero - end subroutine -!=============================================================================== -! MUSCL -! returns flux given at edges which contributes with -! plus sign into 1st. node and with the minus sign into the 2nd node -! IF init_zero=.TRUE. : flux will be set to zero before computation -! IF init_zero=.FALSE. : flux=flux-input flux -! flux is not multiplied with dt - subroutine adv_tra_hor_muscl(vel, ttf, partit, mesh, num_ord, flux, edge_up_dn_grad, nboundary_lay, init_zero) - use MOD_MESH - use MOD_PARTIT - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: num_ord ! num_ord is the fraction of fourth-order contribution in the solution - real(kind=WP), intent(in) :: ttf( mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(inout) :: flux( mesh%nl-1, partit%myDim_edge2D) - integer, intent(in) :: nboundary_lay(partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: edge_up_dn_grad(4, mesh%nl-1, partit%myDim_edge2D) - logical, optional :: init_zero - end subroutine -! a not stable version of MUSCL (reconstruction in the vicinity of bottom topography is not upwind) -! it runs with FCT option only - subroutine adv_tra_hor_mfct(vel, ttf, partit, mesh, num_ord, flux, edge_up_dn_grad, init_zero) - use MOD_MESH - use MOD_PARTIT - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: num_ord ! num_ord is the fraction of fourth-order contribution in the solution - real(kind=WP), intent(in) :: ttf( mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(inout) :: flux( mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(in) :: edge_up_dn_grad(4, mesh%nl-1, partit%myDim_edge2D) - logical, optional :: init_zero - end subroutine - end interface -end module -! -! -!=============================================================================== -subroutine adv_tra_hor_upw1(vel, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_PARTIT - use g_comm_auto - implicit none - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: ttf( mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(inout) :: flux( mesh%nl-1, partit%myDim_edge2D) - logical, optional :: init_zero - real(kind=WP) :: deltaX1, deltaY1, deltaX2, deltaY2 - real(kind=WP) :: a, vflux - integer :: el(2), enodes(2), nz, edge - integer :: nu12, nl12, nl1, nl2, nu1, nu2 - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - ! The result is the low-order solution horizontal fluxes - ! They are put into flux - !___________________________________________________________________________ - do edge=1, myDim_edge2D - ! local indice of nodes that span up edge ed - enodes=edges(:,edge) - - ! local index of element that contribute to edge - el=edge_tri(:,edge) - - ! number of layers -1 at elem el(1) - nl1=nlevels(el(1))-1 - - ! index off surface layer in case of cavity !=1 - nu1=ulevels(el(1)) - - ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to - ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) - deltaX1=edge_cross_dxdy(1,edge) - deltaY1=edge_cross_dxdy(2,edge) - a=r_earth*elem_cos(el(1)) - - !_______________________________________________________________________ - ! same parameter but for other element el(2) that contributes to edge ed - ! if el(2)==0 than edge is boundary edge - nl2=0 - nu2=0 - if(el(2)>0) then - deltaX2=edge_cross_dxdy(3,edge) - deltaY2=edge_cross_dxdy(4,edge) - ! number of layers -1 at elem el(2) - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - a=0.5_WP*(a+r_earth*elem_cos(el(2))) - end if - - !_______________________________________________________________________ - ! nl12 ... minimum number of layers -1 between element el(1) & el(2) that - ! contribute to edge ed - ! nu12 ... upper index of layers between element el(1) & el(2) that - ! contribute to edge ed - ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 - ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! - nl12=min(nl1,nl2) - nu12=max(nu1,nu2) - - !_______________________________________________________________________ - ! (A) goes only into this loop when the edge has only facing element - ! el(1) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - do nz=nu1, nu12-1 - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - - !____________________________________________________________________ - ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)) & - )-flux(nz, edge) - end do - - !_______________________________________________________________________ - ! (B) goes only into this loop when the edge has only facing elemenmt - ! el(2) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - if (nu2 > 0) then - do nz=nu2, nu12-1 - !___________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - - !___________________________________________________________ - ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) - end do - end if - - !_______________________________________________________________________ - ! (C) Both segments - ! loop over depth layers from top (nu12) to nl12 - ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than nl12=0 so - ! you wont enter in this loop - do nz=nu12, nl12 - !___________________________________________________________________ - ! 1st. low order upwind solution - ! here already assumed that ed is NOT! a boundary edge so el(2) should exist - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & - +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) - end do - - !_______________________________________________________________________ - ! (D) remaining segments on the left or on the right - do nz=nl12+1, nl1 - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - !____________________________________________________________________ - ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)) & - )-flux(nz, edge) - end do - - !_______________________________________________________________________ - ! (E) remaining segments on the left or on the right - do nz=nl12+1, nl2 - !_______________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - !_______________________________________________________________ - ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) - end do - end do -end subroutine adv_tra_hor_upw1 -! -! -!=============================================================================== -subroutine adv_tra_hor_muscl(vel, ttf, partit, mesh, num_ord, flux, edge_up_dn_grad, nboundary_lay, init_zero) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - implicit none - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: num_ord ! num_ord is the fraction of fourth-order contribution in the solution - real(kind=WP), intent(in) :: ttf( mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(inout) :: flux( mesh%nl-1, partit%myDim_edge2D) - integer, intent(in) :: nboundary_lay(partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: edge_up_dn_grad(4, mesh%nl-1, partit%myDim_edge2D) - logical, optional :: init_zero - real(kind=WP) :: deltaX1, deltaY1, deltaX2, deltaY2 - real(kind=WP) :: Tmean1, Tmean2, cHO - real(kind=WP) :: c_lo(2) - real(kind=WP) :: a, vflux - integer :: el(2), enodes(2), nz, edge - integer :: nu12, nl12, nl1, nl2, nu1, nu2 - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - ! The result is the low-order solution horizontal fluxes - ! They are put into flux - !___________________________________________________________________________ - do edge=1, myDim_edge2D - ! local indice of nodes that span up edge ed - enodes=edges(:,edge) - - ! local index of element that contribute to edge - el=edge_tri(:,edge) - - ! number of layers -1 at elem el(1) - nl1=nlevels(el(1))-1 - - ! index off surface layer in case of cavity !=1 - nu1=ulevels(el(1)) - - ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to - ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) - deltaX1=edge_cross_dxdy(1,edge) - deltaY1=edge_cross_dxdy(2,edge) - a=r_earth*elem_cos(el(1)) - - !_______________________________________________________________________ - ! same parameter but for other element el(2) that contributes to edge ed - ! if el(2)==0 than edge is boundary edge - nl2=0 - nu2=0 - if(el(2)>0) then - deltaX2=edge_cross_dxdy(3,edge) - deltaY2=edge_cross_dxdy(4,edge) - ! number of layers -1 at elem el(2) - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - a=0.5_WP*(a+r_earth*elem_cos(el(2))) - end if - - !_______________________________________________________________________ - ! n2 ... minimum number of layers -1 between element el(1) & el(2) that - ! contribute to edge ed - ! nu12 ... upper index of layers between element el(1) & el(2) that - ! contribute to edge ed - ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 - ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! - nl12=min(nl1,nl2) - nu12=max(nu1,nu2) - - !_______________________________________________________________________ - ! (A) goes only into this loop when the edge has only facing element - ! el(1) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - do nz=nu1, nu12-1 - c_lo(1)=real(max(sign(1, nboundary_lay(enodes(1))-nz), 0),WP) - c_lo(2)=real(max(sign(1, nboundary_lay(enodes(2))-nz), 0),WP) - - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP*c_lo(2) - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP*c_lo(1) - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - - !_______________________________________________________________________ - ! (B) goes only into this loop when the edge has only facing elemenmt - ! el(2) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - if (nu2 > 0) then - do nz=nu2, nu12-1 - c_lo(1)=real(max(sign(1, nboundary_lay(enodes(1))-nz), 0),WP) - c_lo(2)=real(max(sign(1, nboundary_lay(enodes(2))-nz), 0),WP) - - !_______________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP*c_lo(2) - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP*c_lo(1) - - !_______________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - end if - - !_______________________________________________________________________ - ! (C) Both segments - ! loop over depth layers from top to n2 - ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than n2=0 so - ! you wont enter in this loop - do nz=nu12, nl12 - c_lo(1)=real(max(sign(1, nboundary_lay(enodes(1))-nz), 0),WP) - c_lo(2)=real(max(sign(1, nboundary_lay(enodes(2))-nz), 0),WP) - - !___________________________________________________________________ - ! MUSCL-type reconstruction - ! check if upwind or downwind triagle is necessary - ! - ! cross product between velocity vector and cross vector edge-elem-center - ! cross product > 0 --> angle vec_v and (dx,dy) --> [0 180] --> upwind triangle - ! cross product < 0 --> angle vec_v and (dx,dy) --> [180 360] --> downwind triangle - ! - ! o o ! o o - ! / \ / \ ! / \ / \ - ! / \ \ vec_v / \ ! / \ / / \ - ! / up \ \ / dn \ ! / up \ / / dn \ - ! o-------o----+---->o-------o ! o-------o----+---->o-------o - ! 1 / 2 ! 1 \vec_v - ! /vec_v ! \ - ! --> downwind triangle ! --> upwind triangle - ! - ! edge_up_dn_grad(1,nz,edge) ... gradTR_x upwind - ! edge_up_dn_grad(2,nz,edge) ... gradTR_x downwind - ! edge_up_dn_grad(3,nz,edge) ... gradTR_y upwind - ! edge_up_dn_grad(4,nz,edge) ... gradTR_y downwind - - !___________________________________________________________________ - ! use downwind triangle to interpolate Tracer to edge center with - ! fancy scheme --> Linear upwind reconstruction - ! T_n+0.5 = T_n+1 - 1/2*deltax*GRADIENT - ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri - ! T_n+0.5 = T_n+1 - 2/6*(T_n+1-T_n) + 1/6*gradT_down - ! --> edge_up_dn_grad ... contains already elemental tracer gradient - ! of up and dn wind triangle - ! --> Tmean2 ... edge center interpolated Tracer using tracer - ! gradient info from upwind triangle - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP*c_lo(2) - - ! use upwind triangle to interpolate Tracer to edge center with - ! fancy scheme --> Linear upwind reconstruction - ! T_n+0.5 = T_n + 1/2*deltax*GRADIENT - ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri - ! T_n+0.5 = T_n + 2/6*(T_n+1-T_n) + 1/6*gradT_down - ! --> Tmean1 ... edge center interpolated Tracer using tracer - ! gradient info from downwind triangle - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP*c_lo(1) - - !___________________________________________________________________ - ! volume flux along the edge segment ed - ! netto volume flux along segment that comes from edge node 1 and 2 - ! - ! - ! C1 (centroid el(1)) --> (u1,v1) - ! x - ! ^ - ! (dx1,dy1) | - ! |---> vec_n1 (dy1,-dx1)--> project vec_u1 onto vec_n1 --> -v1*dx1+u1*dy1 --> - ! | | - ! enodes(1) o----------O---------o enodes(2) |-> calculate volume flux out of/in - ! vflux_________/| | the volume of enode1(enode2) through - ! |---> vec_n2 (dy2,-dx2)--> project vec_u2 onto vec_n2 --> -v2*dx2+u2*dy2 --> sections of dx1,dy1 and dx2,dy2 - ! (dx2,dy2) | --> vflux - ! v - ! x - ! C2 (centroid el(2)) --> (u2,v2) - - ! here already assumed that ed is NOT! a boundary edge so el(2) should exist - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & - +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - - !___________________________________________________________________ - ! (1-num_ord) is done with 3rd order upwind - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - - !_______________________________________________________________________ - ! (D) remaining segments on the left or on the right - do nz=nl12+1, nl1 - c_lo(1)=real(max(sign(1, nboundary_lay(enodes(1))-nz), 0),WP) - c_lo(2)=real(max(sign(1, nboundary_lay(enodes(2))-nz), 0),WP) - - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP*c_lo(2) - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP*c_lo(1) - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - - !_______________________________________________________________________ - ! (E) remaining segments on the left or on the right - do nz=nl12+1, nl2 - c_lo(1)=real(max(sign(1, nboundary_lay(enodes(1))-nz), 0),WP) - c_lo(2)=real(max(sign(1, nboundary_lay(enodes(2))-nz), 0),WP) - - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP*c_lo(2) - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP*c_lo(1) - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - end do -end subroutine adv_tra_hor_muscl -! -! -!=============================================================================== - subroutine adv_tra_hor_mfct(vel, ttf, partit, mesh, num_ord, flux, edge_up_dn_grad, init_zero) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - implicit none - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: num_ord ! num_ord is the fraction of fourth-order contribution in the solution - real(kind=WP), intent(in) :: ttf( mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: vel(2, mesh%nl-1, partit%myDim_elem2D+partit%eDim_elem2D) - real(kind=WP), intent(inout) :: flux( mesh%nl-1, partit%myDim_edge2D) - real(kind=WP), intent(in) :: edge_up_dn_grad(4, mesh%nl-1, partit%myDim_edge2D) - logical, optional :: init_zero - real(kind=WP) :: deltaX1, deltaY1, deltaX2, deltaY2 - real(kind=WP) :: Tmean1, Tmean2, cHO - real(kind=WP) :: a, vflux - integer :: el(2), enodes(2), nz, edge - integer :: nu12, nl12, nl1, nl2, nu1, nu2 - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - ! The result is the low-order solution horizontal fluxes - ! They are put into flux - !___________________________________________________________________________ - do edge=1, myDim_edge2D - ! local indice of nodes that span up edge ed - enodes=edges(:,edge) - - ! local index of element that contribute to edge - el=edge_tri(:,edge) - - ! number of layers -1 at elem el(1) - nl1=nlevels(el(1))-1 - - ! index off surface layer in case of cavity !=1 - nu1=ulevels(el(1)) - - ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to - ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) - deltaX1=edge_cross_dxdy(1,edge) - deltaY1=edge_cross_dxdy(2,edge) - a=r_earth*elem_cos(el(1)) - - !_______________________________________________________________________ - ! same parameter but for other element el(2) that contributes to edge ed - ! if el(2)==0 than edge is boundary edge - nl2=0 - nu2=0 - if(el(2)>0) then - deltaX2=edge_cross_dxdy(3,edge) - deltaY2=edge_cross_dxdy(4,edge) - ! number of layers -1 at elem el(2) - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - a=0.5_WP*(a+r_earth*elem_cos(el(2))) - end if - - !_______________________________________________________________________ - ! n2 ... minimum number of layers -1 between element el(1) & el(2) that - ! contribute to edge ed - ! nu12 ... upper index of layers between element el(1) & el(2) that - ! contribute to edge ed - ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 - ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! - nl12=min(nl1,nl2) - nu12=max(nu1,nu2) - - !_______________________________________________________________________ - ! (A) goes only into this loop when the edge has only facing element - ! el(1) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - do nz=nu1, nu12-1 - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - - !_______________________________________________________________________ - ! (B) goes only into this loop when the edge has only facing elemenmt - ! el(2) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - if (nu2 > 0) then - do nz=nu2,nu12-1 - !___________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - !___________________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - end if - - !_______________________________________________________________________ - ! (C) Both segments - ! loop over depth layers from top to n2 - ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than n2=0 so - ! you wont enter in this loop - do nz=nu12, nl12 - !___________________________________________________________________ - ! MUSCL-type reconstruction - ! check if upwind or downwind triagle is necessary - ! - ! cross product between velocity vector and cross vector edge-elem-center - ! cross product > 0 --> angle vec_v and (dx,dy) --> [0 180] --> upwind triangle - ! cross product < 0 --> angle vec_v and (dx,dy) --> [180 360] --> downwind triangle - ! - ! o o ! o o - ! / \ / \ ! / \ / \ - ! / \ \ vec_v / \ ! / \ / / \ - ! / up \ \ / dn \ ! / up \ / / dn \ - ! o-------o----+---->o-------o ! o-------o----+---->o-------o - ! 1 / 2 ! 1 \vec_v - ! /vec_v ! \ - ! --> downwind triangle ! --> upwind triangle - ! - ! edge_up_dn_grad(1,nz,edge) ... gradTR_x upwind - ! edge_up_dn_grad(2,nz,edge) ... gradTR_x downwind - ! edge_up_dn_grad(3,nz,edge) ... gradTR_y upwind - ! edge_up_dn_grad(4,nz,edge) ... gradTR_y downwind - - !___________________________________________________________________ - ! use downwind triangle to interpolate Tracer to edge center with - ! fancy scheme --> Linear upwind reconstruction - ! T_n+0.5 = T_n+1 - 1/2*deltax*GRADIENT - ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri - ! T_n+0.5 = T_n+1 - 2/6*(T_n+1-T_n) + 1/6*gradT_down - ! --> edge_up_dn_grad ... contains already elemental tracer gradient - ! of up and dn wind triangle - ! --> Tmean2 ... edge center interpolated Tracer using tracer - ! gradient info from upwind triangle - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - ! use upwind triangle to interpolate Tracer to edge center with - ! fancy scheme --> Linear upwind reconstruction - ! T_n+0.5 = T_n + 1/2*deltax*GRADIENT - ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri - ! T_n+0.5 = T_n + 2/6*(T_n+1-T_n) + 1/6*gradT_down - ! --> Tmean1 ... edge center interpolated Tracer using tracer - ! gradient info from downwind triangle - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - !___________________________________________________________________ - ! volume flux along the edge segment ed - ! netto volume flux along segment that comes from edge node 1 and 2 - ! - ! - ! C1 (centroid el(1)) --> (u1,v1) - ! x - ! ^ - ! (dx1,dy1) | - ! |---> vec_n1 (dy1,-dx1)--> project vec_u1 onto vec_n1 --> -v1*dx1+u1*dy1 --> - ! | | - ! enodes(1) o----------O---------o enodes(2) |-> calculate volume flux out of/in - ! vflux_________/| | the volume of enode1(enode2) through - ! |---> vec_n2 (dy2,-dx2)--> project vec_u2 onto vec_n2 --> -v2*dx2+u2*dy2 --> sections of dx1,dy1 and dx2,dy2 - ! (dx2,dy2) | --> vflux - ! v - ! x - ! C2 (centroid el(2)) --> (u2,v2) - - ! here already assumed that ed is NOT! a boundary edge so el(2) should exist - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & - +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - - !___________________________________________________________________ - ! (1-num_ord) is done with 3rd order upwind - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - - !_______________________________________________________________________ - ! (D) remaining segments on the left or on the right - do nz=nl12+1, nl1 - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - - !_______________________________________________________________________ - ! (E) remaining segments on the left or on the right - do nz=nl12+1, nl2 - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - end do -end subroutine adv_tra_hor_mfct - diff --git a/src/temp/oce_adv_tra_ver.F90 b/src/temp/oce_adv_tra_ver.F90 deleted file mode 100644 index eab9847a8..000000000 --- a/src/temp/oce_adv_tra_ver.F90 +++ /dev/null @@ -1,598 +0,0 @@ -module oce_adv_tra_ver_interfaces - interface -! implicit 1st order upwind vertical advection with to solve for fct_LO -! updates the input tracer ttf - subroutine adv_tra_vert_impl(dt, w, ttf, partit, mesh) - use mod_mesh - use MOD_PARTIT - real(kind=WP), intent(in), target :: dt - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(inout) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - end subroutine -!=============================================================================== -! 1st order upwind (explicit) -! returns flux given at vertical interfaces of scalar volumes -! IF init_zero=.TRUE. : flux will be set to zero before computation -! IF init_zero=.FALSE. : flux=flux-input flux -! flux is not multiplied with dt - subroutine adv_tra_ver_upw1(w, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_PARTIT - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - end subroutine -!=============================================================================== -! QR (4th order centerd) -! returns flux given at vertical interfaces of scalar volumes -! IF init_zero=.TRUE. : flux will be set to zero before computation -! IF init_zero=.FALSE. : flux=flux-input flux -! flux is not multiplied with dt - subroutine adv_tra_ver_qr4c(w, ttf, partit, mesh, num_ord, flux, init_zero) - use MOD_MESH - use MOD_PARTIT - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: num_ord ! num_ord is the fraction of fourth-order contribution in the solution - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - end subroutine -!=============================================================================== -! Vertical advection with PPM reconstruction (5th order) -! returns flux given at vertical interfaces of scalar volumes -! IF init_zero=.TRUE. : flux will be set to zero before computation -! IF init_zero=.FALSE. : flux=flux-input flux -! flux is not multiplied with dt - subroutine adv_tra_vert_ppm(dt, w, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_PARTIT - real(kind=WP), intent(in), target :: dt - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - integer :: n, nz, nl1 - real(kind=WP) :: tvert(mesh%nl), tv - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - end subroutine -! central difference reconstruction (2nd order, use only with FCT) -! returns flux given at vertical interfaces of scalar volumes -! IF init_zero=.TRUE. : flux will be set to zero before computation -! IF init_zero=.FALSE. : flux=flux-input flux -! flux is not multiplied with dt - subroutine adv_tra_ver_cdiff(w, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_PARTIT - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - integer :: n, nz, nl1 - real(kind=WP) :: tvert(mesh%nl), tv - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - end subroutine - end interface -end module -!=============================================================================== -subroutine adv_tra_vert_impl(dt, w, ttf, partit, mesh) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - - implicit none - real(kind=WP), intent(in) , target :: dt - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in) , target :: mesh - real(kind=WP), intent(inout) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP) :: a(mesh%nl), b(mesh%nl), c(mesh%nl), tr(mesh%nl) - real(kind=WP) :: cp(mesh%nl), tp(mesh%nl) - integer :: nz, n, nzmax, nzmin, tr_num - real(kind=WP) :: m, zinv, dt_inv, dz - real(kind=WP) :: c1, v_adv - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - dt_inv=1.0_WP/dt - - !___________________________________________________________________________ - ! loop over local nodes - do n=1,myDim_nod2D - - ! initialise - a = 0.0_WP - b = 0.0_WP - c = 0.0_WP - tr = 0.0_WP - tp = 0.0_WP - cp = 0.0_WP - - ! max. number of levels at node n - nzmax=nlevels_nod2D(n) - - ! upper surface index, in case of cavity !=1 - nzmin=ulevels_nod2D(n) - - !___________________________________________________________________________ - ! Here can not exchange zbar_n & Z_n with zbar_3d_n & Z_3d_n because - ! they be calculate from the actualized mesh with hnode_new - ! calculate new zbar (depth of layers) and Z (mid depths of layers) - ! depending on layer thinkness over depth at node n - ! Be carefull here vertical operation have to be done on NEW vertical mesh !!! - zbar_n=0.0_WP - Z_n=0.0_WP - zbar_n(nzmax)=zbar_n_bot(n) - Z_n(nzmax-1) =zbar_n(nzmax) + hnode_new(nzmax-1,n)/2.0_WP - do nz=nzmax-1,nzmin+1,-1 - zbar_n(nz) = zbar_n(nz+1) + hnode_new(nz,n) - Z_n(nz-1) = zbar_n(nz) + hnode_new(nz-1,n)/2.0_WP - end do - zbar_n(nzmin) = zbar_n(nzmin+1) + hnode_new(nzmin,n) - - !_______________________________________________________________________ - ! Regular part of coefficients: --> surface layer - nz=nzmin - - ! 1/dz(nz) - zinv=1.0_WP*dt ! no .../(zbar(1)-zbar(2)) because of ALE - - !!PS a(nz)=0.0_WP - !!PS v_adv=zinv*areasvol(nz+1,n)/areasvol(nz,n) - !!PS b(nz)= hnode_new(nz,n)+W(nz, n)*zinv-min(0._WP, W(nz+1, n))*v_adv - !!PS c(nz)=-max(0._WP, W(nz+1, n))*v_adv - - a(nz)=0.0_WP - v_adv=zinv*area(nz ,n)/areasvol(nz,n) - b(nz)= hnode_new(nz,n)+W(nz, n)*v_adv - - v_adv=zinv*area(nz+1,n)/areasvol(nz,n) - b(nz)= b(nz)-min(0._WP, W(nz+1, n))*v_adv - c(nz)=-max(0._WP, W(nz+1, n))*v_adv - - !_______________________________________________________________________ - ! Regular part of coefficients: --> 2nd...nl-2 layer - do nz=nzmin+1, nzmax-2 - ! update from the vertical advection - v_adv=zinv*area(nz ,n)/areasvol(nz,n) - a(nz)=min(0._WP, W(nz, n))*v_adv - b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*v_adv - - v_adv=zinv*area(nz+1,n)/areasvol(nz,n) - b(nz)=b(nz)-min(0._WP, W(nz+1, n))*v_adv - c(nz)= -max(0._WP, W(nz+1, n))*v_adv - end do ! --> do nz=2, nzmax-2 - - !_______________________________________________________________________ - ! Regular part of coefficients: --> nl-1 layer - nz=nzmax-1 - ! update from the vertical advection - !!PS a(nz)= min(0._WP, W(nz, n))*zinv - !!PS b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*zinv - !!PS c(nz)=0.0_WP - v_adv=zinv*area(nz ,n)/areasvol(nz,n) - a(nz)= min(0._WP, W(nz, n))*v_adv - b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*v_adv - c(nz)=0.0_WP - - !_______________________________________________________________________ - nz=nzmin - dz=hnode_new(nz,n) ! It would be (zbar(nz)-zbar(nz+1)) if not ALE - tr(nz)=-(b(nz)-dz)*ttf(nz,n)-c(nz)*ttf(nz+1,n) - - do nz=nzmin+1,nzmax-2 - dz=hnode_new(nz,n) - tr(nz)=-a(nz)*ttf(nz-1,n)-(b(nz)-dz)*ttf(nz,n)-c(nz)*ttf(nz+1,n) - end do - nz=nzmax-1 - dz=hnode_new(nz,n) - tr(nz)=-a(nz)*ttf(nz-1,n)-(b(nz)-dz)*ttf(nz,n) - - !_______________________________________________________________________ - nz = nzmin - cp(nz) = c(nz)/b(nz) - tp(nz) = tr(nz)/b(nz) - - ! solve for vectors c-prime and t, s-prime - do nz = nzmin+1,nzmax-1 - m = b(nz)-cp(nz-1)*a(nz) - cp(nz) = c(nz)/m - tp(nz) = (tr(nz)-tp(nz-1)*a(nz))/m - end do - - !_______________________________________________________________________ - ! start with back substitution - tr(nzmax-1) = tp(nzmax-1) - - ! solve for x from the vectors c-prime and d-prime - do nz = nzmax-2, nzmin, -1 - tr(nz) = tp(nz)-cp(nz)*tr(nz+1) - end do - - !_______________________________________________________________________ - ! update tracer - do nz=nzmin,nzmax-1 - ttf(nz,n)=ttf(nz,n)+tr(nz) - end do - end do ! --> do n=1,myDim_nod2D -end subroutine adv_tra_vert_impl -! -! -!=============================================================================== -subroutine adv_tra_ver_upw1(w, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - - implicit none - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP) :: tvert(mesh%nl) - integer :: n, nz, nzmax, nzmin - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - do n=1, myDim_nod2D - !_______________________________________________________________________ - nzmax=nlevels_nod2D(n) - nzmin=ulevels_nod2D(n) - - !_______________________________________________________________________ - ! vert. flux at surface layer - nz=nzmin - flux(nz,n)=-W(nz,n)*ttf(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux at bottom layer --> zero bottom flux - nz=nzmax - flux(nz,n)= 0.0_WP-flux(nz,n) - - !_______________________________________________________________________ - ! Be carefull have to do vertical tracer advection here on old vertical grid - ! also horizontal advection is done on old mesh (see helem contains old - ! mesh information) - !_______________________________________________________________________ - ! vert. flux at remaining levels - do nz=nzmin+1,nzmax-1 - flux(nz,n)=-0.5*( & - ttf(nz ,n)*(W(nz,n)+abs(W(nz,n)))+ & - ttf(nz-1,n)*(W(nz,n)-abs(W(nz,n))))*area(nz,n)-flux(nz,n) - end do - end do -end subroutine adv_tra_ver_upw1 -! -! -!=============================================================================== -subroutine adv_tra_ver_qr4c(w, ttf, partit, mesh, num_ord, flux, init_zero) - use MOD_MESH - use o_ARRAYS - use o_PARAM - use MOD_PARTIT - implicit none - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: num_ord ! num_ord is the fraction of fourth-order contribution in the solution - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - real(kind=WP) :: tvert(mesh%nl) - integer :: n, nz, nzmax, nzmin - real(kind=WP) :: Tmean, Tmean1, Tmean2 - real(kind=WP) :: qc, qu, qd - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - do n=1, myDim_nod2D - !_______________________________________________________________________ - nzmax=nlevels_nod2D(n) - nzmin=ulevels_nod2D(n) - !_______________________________________________________________________ - ! vert. flux at surface layer - nz=nzmin - flux(nz,n)=-ttf(nz,n)*W(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux 2nd layer --> centered differences - nz=nzmin+1 - flux(nz,n)=-0.5_WP*(ttf(nz-1,n)+ttf(nz,n))*W(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux at bottom - 1 layer --> centered differences - nz=nzmax-1 - flux(nz,n)=-0.5_WP*(ttf(nz-1,n)+ttf(nz,n))*W(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux at bottom layer --> zero bottom flux - nz=nzmax - flux(nz,n)= 0.0_WP-flux(nz,n) - - !_______________________________________________________________________ - ! Be carefull have to do vertical tracer advection here on old vertical grid - ! also horizontal advection is done on old mesh (see helem contains old - ! mesh information) - !_______________________________________________________________________ - ! vert. flux at remaining levels - do nz=nzmin+2,nzmax-2 - !centered (4th order) - qc=(ttf(nz-1,n)-ttf(nz ,n))/(Z_3d_n(nz-1,n)-Z_3d_n(nz ,n)) - qu=(ttf(nz ,n)-ttf(nz+1,n))/(Z_3d_n(nz ,n)-Z_3d_n(nz+1,n)) - qd=(ttf(nz-2,n)-ttf(nz-1,n))/(Z_3d_n(nz-2,n)-Z_3d_n(nz-1,n)) - - Tmean1=ttf(nz ,n)+(2*qc+qu)*(zbar_3d_n(nz,n)-Z_3d_n(nz ,n))/3.0_WP - Tmean2=ttf(nz-1,n)+(2*qc+qd)*(zbar_3d_n(nz,n)-Z_3d_n(nz-1,n))/3.0_WP - Tmean =(W(nz,n)+abs(W(nz,n)))*Tmean1+(W(nz,n)-abs(W(nz,n)))*Tmean2 - ! flux(nz,n)=-0.5_WP*(num_ord*(Tmean1+Tmean2)*W(nz,n)+(1.0_WP-num_ord)*Tmean)*area(nz,n)-flux(nz,n) - flux(nz,n)=(-0.5_WP*(1.0_WP-num_ord)*Tmean - num_ord*(0.5_WP*(Tmean1+Tmean2))*W(nz,n))*area(nz,n)-flux(nz,n) - end do - end do -end subroutine adv_tra_ver_qr4c -! -! -!=============================================================================== -subroutine adv_tra_vert_ppm(dt, w, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - implicit none - real(kind=WP), intent(in), target :: dt - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in) , target :: mesh - real(kind=WP), intent(in) :: ttf (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - real(kind=WP) :: tvert(mesh%nl), tv(mesh%nl), aL, aR, aj, x - real(kind=WP) :: dzjm1, dzj, dzjp1, dzjp2, deltaj, deltajp1 - integer :: n, nz, nzmax, nzmin - integer :: overshoot_counter, counter - -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - ! -------------------------------------------------------------------------- - ! Vertical advection - ! -------------------------------------------------------------------------- - ! A piecewise parabolic scheme for uniformly-spaced layers. - ! See Colella and Woodward, JCP, 1984, 174-201. It can be coded so as to to take - ! non-uniformity into account, but this is more cumbersome. This is the version for AB - ! time stepping - ! -------------------------------------------------------------------------- - overshoot_counter=0 - counter =0 - do n=1, myDim_nod2D - !_______________________________________________________________________ - !Interpolate to zbar...depth levels --> all quantities (tracer ...) are - ! calculated on mid depth levels - ! nzmax ... number of depth levels at node n - nzmax=nlevels_nod2D(n) - nzmin=ulevels_nod2D(n) - - ! tracer at surface level - tv(nzmin)=ttf(nzmin,n) - ! tracer at surface+1 level -! tv(2)=-ttf(1,n)*min(sign(1.0, W(2,n)), 0._WP)+ttf(2,n)*max(sign(1.0, W(2,n)), 0._WP) -! tv(3)=-ttf(2,n)*min(sign(1.0, W(3,n)), 0._WP)+ttf(3,n)*max(sign(1.0, W(3,n)), 0._WP) - tv(nzmin+1)=0.5*(ttf(nzmin, n)+ttf(nzmin+1,n)) - ! tacer at bottom-1 level - tv(nzmax-1)=-ttf(nzmax-2,n)*min(sign(1.0_wp, W(nzmax-1,n)), 0._WP)+ttf(nzmax-1,n)*max(sign(1.0_wp, W(nzmax-1,n)), 0._WP) -! tv(nzmax-1)=0.5_WP*(ttf(nzmax-2,n)+ttf(nzmax-1,n)) - ! tracer at bottom level - tv(nzmax)=ttf(nzmax-1,n) - - !_______________________________________________________________________ - ! calc tracer for surface+2 until depth-2 layer - ! see Colella and Woodward, JCP, 1984, 174-201 --> equation (1.9) - ! loop over layers (segments) - !!PS do nz=3, nzmax-3 - do nz=nzmin+1, nzmax-3 - !___________________________________________________________________ - ! for uniform spaced vertical grids --> piecewise parabolic method (ppm) - ! equation (1.9) - ! tv(nz)=(7.0_WP*(ttf(nz-1,n)+ttf(nz,n))-(ttf(nz-2,n)+ttf(nz+1,n)))/12.0_WP - - !___________________________________________________________________ - ! for non-uniformity spaced vertical grids --> piecewise parabolic - ! method (ppm) see see Colella and Woodward, JCP, 1984, 174-201 - ! --> full equation (1.6), (1.7) and (1.8) - dzjm1 = hnode_new(nz-1,n) - dzj = hnode_new(nz ,n) - dzjp1 = hnode_new(nz+1,n) - dzjp2 = hnode_new(nz+2,n) - ! Be carefull here vertical operation have to be done on NEW vertical mesh !!! - - !___________________________________________________________________ - ! equation (1.7) - ! --> Here deltaj is the average slope in the jth zone of the parabola - ! with zone averages a_(j-1) and a_j, a_(j+1) - ! --> a_j^n - deltaj = dzj/(dzjm1+dzj+dzjp1)* & - ( & - (2._WP*dzjm1+dzj )/(dzjp1+dzj)*(ttf(nz+1,n)-ttf(nz ,n)) + & - (dzj +2._WP*dzjp1)/(dzjm1+dzj)*(ttf(nz ,n)-ttf(nz-1,n)) & - ) - ! --> a_(j+1)^n - deltajp1 = dzjp1/(dzj+dzjp1+dzjp2)* & - ( & - (2._WP*dzj+dzjp1 )/(dzjp2+dzjp1)*(ttf(nz+2,n)-ttf(nz+1,n)) + & - (dzjp1+2._WP*dzjp2)/(dzj +dzjp1)*(ttf(nz+1,n)-ttf(nz ,n)) & - ) - !___________________________________________________________________ - ! condition (1.8) - ! --> This modification leads to a somewhat steeper representation of - ! discontinuities in the solution. It also guarantees that a_(j+0.5) - ! lies in the range of values defined by a_j; and a_(j+1); - if ( (ttf(nz+1,n)-ttf(nz ,n))*(ttf(nz ,n)-ttf(nz-1,n)) > 0._WP ) then - deltaj = min( abs(deltaj), & - 2._WP*abs(ttf(nz+1,n)-ttf(nz ,n)),& - 2._WP*abs(ttf(nz ,n)-ttf(nz-1,n)) & - )*sign(1.0_WP,deltaj) - else - deltaj = 0.0_WP - endif - if ( (ttf(nz+2,n)-ttf(nz+1,n))*(ttf(nz+1,n)-ttf(nz ,n)) > 0._WP ) then - deltajp1 = min( abs(deltajp1),& - 2._WP*abs(ttf(nz+2,n)-ttf(nz+1,n)),& - 2._WP*abs(ttf(nz+1,n)-ttf(nz,n)) & - )*sign(1.0_WP,deltajp1) - else - deltajp1 = 0.0_WP - endif - !___________________________________________________________________ - ! equation (1.6) - ! --> calcualte a_(j+0.5) - ! nz+1 is the interface betweel layers (segments) nz and nz+1 - tv(nz+1)= ttf(nz,n) & - + dzj/(dzj+dzjp1)*(ttf(nz+1,n)-ttf(nz,n)) & - + 1._WP/(dzjm1+dzj+dzjp1+dzjp2) * & - ( & - (2._WP*dzjp1*dzj)/(dzj+dzjp1)* & - ((dzjm1+dzj)/(2._WP*dzj+dzjp1) - (dzjp2+dzjp1)/(2._WP*dzjp1+dzj))*(ttf(nz+1,n)-ttf(nz,n)) & - - dzj*(dzjm1+dzj)/(2._WP*dzj+dzjp1)*deltajp1 & - + dzjp1*(dzjp1+dzjp2)/(dzj+2._WP*dzjp1)*deltaj & - ) - !tv(nz+1)=max(min(ttf(nz, n), ttf(nz+1, n)), min(max(ttf(nz, n), ttf(nz+1, n)), tv(nz+1))) - end do ! --> do nz=2,nzmax-3 - - tvert(1:nzmax)=0._WP - ! loop over layers (segments) - do nz=nzmin, nzmax-1 - if ((W(nz,n)<=0._WP) .AND. (W(nz+1,n)>=0._WP)) CYCLE - counter=counter+1 - aL=tv(nz) - aR=tv(nz+1) - if ((aR-ttf(nz, n))*(ttf(nz, n)-aL)<=0._WP) then - ! write(*,*) aL, ttf(nz, n), aR - overshoot_counter=overshoot_counter+1 - aL =ttf(nz, n) - aR =ttf(nz, n) - end if - if ((aR-aL)*(ttf(nz, n)-0.5_WP*(aL+aR))> (aR-aL)**2/6._WP) then - aL =3._WP*ttf(nz, n)-2._WP*aR - end if - if ((aR-aL)*(ttf(nz, n)-0.5_WP*(aR+aL))<-(aR-aL)**2/6._WP) then - aR =3._WP*ttf(nz, n)-2._WP*aL - end if - - dzj = hnode(nz,n) - aj=6.0_WP*(ttf(nz, n)-0.5_WP*(aL+aR)) - - if (W(nz,n)>0._WP) then - x=min(W(nz,n)*dt/dzj, 1._WP) - tvert(nz )=(-aL-0.5_WP*x*(aR-aL+(1._WP-2._WP/3._WP*x)*aj)) - tvert(nz )=tvert(nz) ! compute 2nd moment for DVD - tvert(nz )=tvert(nz)*area(nz,n)*W(nz,n) - end if - - if (W(nz+1,n)<0._WP) then - x=min(-W(nz+1,n)*dt/dzj, 1._WP) - tvert(nz+1)=(-aR+0.5_WP*x*(aR-aL-(1._WP-2._WP/3._WP*x)*aj)) - tvert(nz+1)=tvert(nz+1) ! compute 2nd moment for DVD - tvert(nz+1)=tvert(nz+1)*area(nz+1,n)*W(nz+1,n) - end if - end do - - !_______________________________________________________________________ - ! Surface flux - tvert(nzmin)= -tv(nzmin)*W(nzmin,n)*area(nzmin,n) - ! Zero bottom flux - tvert(nzmax)=0.0_WP - flux(nzmin:nzmax, n)=tvert(nzmin:nzmax)-flux(nzmin:nzmax, n) - end do ! --> do n=1, myDim_nod2D -! if (mype==0) write(*,*) 'PPM overshoot statistics:', real(overshoot_counter)/real(counter) -end subroutine adv_tra_vert_ppm -! -! -!=============================================================================== -subroutine adv_tra_ver_cdiff(w, ttf, partit, mesh, flux, init_zero) - use MOD_MESH - use MOD_TRACER - use MOD_PARTIT - use g_comm_auto - implicit none - type(t_partit),intent(in), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP), intent(in) :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: W (mesh%nl, partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(inout) :: flux(mesh%nl, partit%myDim_nod2D) - logical, optional :: init_zero - integer :: n, nz, nzmax, nzmin - real(kind=WP) :: tvert(mesh%nl), tv -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - - if (present(init_zero))then - if (init_zero) flux=0.0_WP - else - flux=0.0_WP - end if - - do n=1, myDim_nod2D - !_______________________________________________________________________ - nzmax=nlevels_nod2D(n)-1 - nzmin=ulevels_nod2D(n) - - !_______________________________________________________________________ - ! Surface flux - tvert(nzmin)= -W(nzmin,n)*ttf(nzmin,n)*area(nzmin,n) - - !_______________________________________________________________________ - ! Zero bottom flux - tvert(nzmax+1)=0.0_WP - - !_______________________________________________________________________ - ! Other levels - do nz=nzmin+1, nzmax - tv=0.5_WP*(ttf(nz-1,n)+ttf(nz,n)) - tvert(nz)= -tv*W(nz,n)*area(nz,n) - end do - - !_______________________________________________________________________ - flux(nzmin:nzmax, n)=tvert(nzmin:nzmax)-flux(nzmin:nzmax, n) - end do ! --> do n=1, myDim_nod2D -end subroutine adv_tra_ver_cdiff diff --git a/src/temp/oce_modules.F90 b/src/temp/oce_modules.F90 deleted file mode 100755 index 2da93f461..000000000 --- a/src/temp/oce_modules.F90 +++ /dev/null @@ -1,267 +0,0 @@ - -! Modules of cell-vertex ocean model -! S. Danilov, 2012 (sergey.danilov@awi.de) -! SI units are used - -!========================================================== -MODULE o_PARAM -integer, parameter :: WP=8 ! Working precision -integer, parameter :: MAX_PATH=4096 ! Maximum file path length -integer :: mstep -real(kind=WP), parameter :: pi=3.14159265358979 -real(kind=WP), parameter :: rad=pi/180.0_WP -real(kind=WP), parameter :: density_0=1030.0_WP -real(kind=WP), parameter :: density_0_r=1.0_WP/density_0 ! [m^3/kg] -real(kind=WP), parameter :: g=9.81_WP -real(kind=WP), parameter :: r_earth=6367500.0_WP -real(kind=WP), parameter :: omega=2*pi/(3600.0_WP*24.0_WP) -real(kind=WP), parameter :: vcpw=4.2e6 ![J/m^3/K] water heat cap -real(kind=WP), parameter :: inv_vcpw = 1._WP / vcpw ! inverse, to replace divide by multiply -real(kind=WP), parameter :: small=1.0e-8 !small value -integer :: state_equation = 1 !1 - full equation of state, 0 - linear equation of state - -real(kind=WP) :: C_d= 0.0025_WP ! Bottom drag coefficient -real(kind=WP) :: kappa=0.4 !von Karman's constant -real(kind=WP) :: mix_coeff_PP=0.01_WP ! mixing coef for PP scheme -real(kind=WP) :: gamma0=0.01! [m/s], gamma0*len*dt is the background viscosity -real(kind=WP) :: gamma1=0.1! [non dim.], or computation of the flow aware viscosity -real(kind=WP) :: gamma2=10.! [s/m], is only used in easy backscatter option -real(kind=WP) :: Div_c =1.0_WP !modified Leith viscosity weight -real(kind=WP) :: Leith_c=1.0_WP !Leith viscosity weight. It needs vorticity! -real(kind=WP) :: easy_bs_return=1.0 !backscatter option only (how much to return) -real(kind=WP) :: A_ver=0.001_WP ! Vertical harm. visc. -integer :: visc_option=5 -logical :: uke_scaling=.true. -real(kind=WP) :: uke_scaling_factor=1._WP -real(kind=WP) :: rosb_dis=1._WP -integer :: smooth_back=2 -integer :: smooth_dis=2 -integer :: smooth_back_tend=4 -real(kind=WP) :: K_back=600._WP -real(kind=WP) :: c_back=0.1_8 -real(kind=WP) :: K_hor=10._WP -real(kind=WP) :: K_ver=0.00001_WP -real(kind=WP) :: scale_area=2.0e8 -real(kind=WP) :: surf_relax_T= 0.0_WP -real(kind=WP) :: surf_relax_S= 10.0_WP/(60*3600.0_WP*24) -logical :: balance_salt_water =.true. -real(kind=WP) :: clim_relax= 1.0_WP/(10*3600.0_WP*24) -real(kind=WP) :: clim_decay, clim_growth - ! set to 0.0 if no relaxation -logical :: ref_sss_local=.false. -real(kind=WP) :: ref_sss=34.7 -logical :: Fer_GM =.false. !flag for Ferrari et al. (2010) GM scheme -real(kind=WP) :: K_GM_max = 3000. -real(kind=WP) :: K_GM_min = 2.0 -integer :: K_GM_bvref = 2 ! 0...surface, 1...bottom mixlay, 2...mean over mixlay -real(kind=WP) :: K_GM_resscalorder = 2.0 -real(kind=WP) :: K_GM_rampmax = 40.0 ! Resol >K_GM_rampmax[km] GM full -real(kind=WP) :: K_GM_rampmin = 30.0 ! Resol replace string by int comparison -real(KIND=WP) :: Ricr = 0.3_WP ! critical bulk Richardson Number -real(KIND=WP) :: concv = 1.6_WP ! constant for pure convection (eqn. 23) (Large 1.5-1.6; MOM default 1.8) - -logical :: hbl_diag =.false. ! writen boundary layer depth -logical :: use_global_tides=.false. ! tidal potential will be computed and used in the SSH gradient computation -! Time stepping -! real(kind=WP) :: alpha=1.0_WP, theta=1.0_WP ! implicitness for -real(kind=WP) :: alpha=1.0_WP, theta=1.0_WP ! implicitness for - ! elevation and divergence -real(kind=WP) :: epsilon=0.1_WP ! AB2 offset -! Tracers -logical :: i_vert_visc= .true. -logical :: w_split =.false. -real(kind=WP) :: w_max_cfl=1.e-5_WP - -logical :: SPP=.false. - -TYPE tracer_source3d_type - integer :: locID - integer :: ID - integer, allocatable, dimension(:) :: ind2 -END TYPE tracer_source3d_type - -type(tracer_source3d_type), & - allocatable, dimension(:) :: ptracers_restore -integer :: ptracers_restore_total=0 - - -! Momentum -logical :: free_slip=.false. - ! false=no slip -integer :: mom_adv=2 - ! 1 vector control volumes, p1 velocities - ! 2 scalar control volumes - ! 3 vector invariant - -logical :: open_b=.false. ! Reserved - -!_______________________________________________________________________________ -!--> mixing enhancement than can be applied via subroutine mo_convect(mesh) -! additionally to every mixing scheme i.e. KPP, PP, cvmix_KPP, cvmix_PP, cvmix_TKE - -! Switch for Monin-Obukov TB04 mixing --> can be additionally applied for all mixing schemes -! --> definetely recommented for KPP -logical :: use_momix = .true. !.false. !Monin-Obukhov -> TB04 mixing on/off -real(kind=WP) :: momix_lat = -50.0_WP ! latitudinal treshhold to apply mo_on enhanced -! convection -logical :: use_instabmix = .true. -real(kind=WP) :: instabmix_kv = 0.1 - -! Switch for enhanced wind mixing --> nasty trick from pp mixing in FESOM1.4 -logical :: use_windmix = .false. -real(kind=WP) :: windmix_kv = 1.e-3 -integer :: windmix_nl = 2 - -!_______________________________________________________________________________ -! use non-constant reference density if .false. density_ref=density_0 -logical :: use_density_ref = .false. -real(kind=WP) :: density_ref_T = 2.0_WP -real(kind=WP) :: density_ref_S = 34.0_WP - -!_______________________________________________________________________________ -! use k-profile nonlocal fluxes -logical :: use_kpp_nonlclflx = .false. - -!_______________________________________________________________________________ -! *** active tracer cutoff -logical :: limit_salinity=.true. !set an allowed range for salinity -real(kind=WP) :: salinity_min=5.0 !minimal salinity -real(kind=WP) :: coeff_limit_salinity=0.0023 !m/s, coefficient to restore s to s_min - - namelist /tracer_cutoff/ limit_salinity, salinity_min, coeff_limit_salinity - -! *** others *** - real(kind=WP) :: time_sum=0.0 ! for runtime estimate - -!___________________________________________ -! Pressure Gradient Force calculation (pgf) -! calculation of pgf either: -! only linfs: -! > 'nemo' ... like NEMO (interpolate to elemental depth, inter-/extrapolation) -! linfs, zlevel, zstar: -! > 'shchepetkin' ... based on density jacobian -! > 'cubicspline' ... like in FESOM1.4 -! > 'easypgf' ... interpolate pressure on elemental depth -character(20) :: which_pgf='shchepetkin' - - - NAMELIST /oce_dyn/ state_equation, C_d, A_ver, gamma0, gamma1, gamma2, Leith_c, Div_c, easy_bs_return, & - scale_area, mom_adv, free_slip, i_vert_visc, w_split, w_max_cfl, SPP,& - Fer_GM, K_GM_max, K_GM_min, K_GM_bvref, K_GM_resscalorder, K_GM_rampmax, K_GM_rampmin, & - scaling_Ferreira, scaling_Rossby, scaling_resolution, scaling_FESOM14, & - Redi, visc_sh_limit, mix_scheme, Ricr, concv, which_pgf, visc_option, alpha, theta, use_density_ref, & - K_back, c_back, uke_scaling, uke_scaling_factor, smooth_back, smooth_dis, & - smooth_back_tend, rosb_dis - - NAMELIST /tracer_phys/ diff_sh_limit, Kv0_const, double_diffusion, K_ver, K_hor, surf_relax_T, surf_relax_S, & - balance_salt_water, clim_relax, ref_sss_local, ref_sss, & - use_momix, momix_lat, momix_kv, & - use_instabmix, instabmix_kv, & - use_windmix, windmix_kv, windmix_nl, & - use_kpp_nonlclflx - -END MODULE o_PARAM -!========================================================== -MODULE o_ARRAYS -USE o_PARAM -IMPLICIT NONE -! Arrays are described in subroutine array_setup -real(kind=WP), allocatable, target :: Wvel(:,:), Wvel_e(:,:), Wvel_i(:,:) -real(kind=WP), allocatable :: UV(:,:,:) -real(kind=WP), allocatable :: UV_rhs(:,:,:), UV_rhsAB(:,:,:) -real(kind=WP), allocatable :: uke(:,:), v_back(:,:), uke_back(:,:), uke_dis(:,:), uke_dif(:,:) -real(kind=WP), allocatable :: uke_rhs(:,:), uke_rhs_old(:,:) -real(kind=WP), allocatable :: UV_dis_tend(:,:,:), UV_back_tend(:,:,:), UV_total_tend(:,:,:), UV_dis_tend_node(:,:,:) -real(kind=WP), allocatable :: UV_dis_posdef_b2(:,:), UV_dis_posdef(:,:), UV_back_posdef(:,:) -real(kind=WP), allocatable :: eta_n(:), d_eta(:) -real(kind=WP), allocatable :: ssh_rhs(:), hpressure(:,:) -real(kind=WP), allocatable :: CFL_z(:,:) -real(kind=WP), allocatable :: stress_surf(:,:) -real(kind=WP), allocatable :: stress_node_surf(:,:) -REAL(kind=WP), ALLOCATABLE :: stress_atmoce_x(:) -REAL(kind=WP), ALLOCATABLE :: stress_atmoce_y(:) -real(kind=WP), allocatable :: heat_flux(:), Tsurf(:) -real(kind=WP), allocatable :: heat_flux_in(:) !to keep the unmodified (by SW penetration etc.) heat flux -real(kind=WP), allocatable :: water_flux(:), Ssurf(:) -real(kind=WP), allocatable :: virtual_salt(:), relax_salt(:) -real(kind=WP), allocatable :: Tclim(:,:), Sclim(:,:) -real(kind=WP), allocatable :: Visc(:,:) -real(kind=WP), allocatable :: Tsurf_t(:,:), Ssurf_t(:,:) -real(kind=WP), allocatable :: tau_x_t(:,:), tau_y_t(:,:) -real(kind=WP), allocatable :: heat_flux_t(:,:), heat_rel_t(:,:), heat_rel(:) -real(kind=WP), allocatable :: coriolis(:), coriolis_node(:) -real(kind=WP), allocatable :: relax2clim(:) -real(kind=WP), allocatable :: MLD1(:), MLD2(:), MLD3(:) -integer, allocatable :: MLD1_ind(:), MLD2_ind(:), MLD3_ind(:) -real(kind=WP), allocatable :: ssh_gp(:) -!Tracer gradients&RHS -real(kind=WP), allocatable :: ttrhs(:,:) -real(kind=WP), allocatable :: tr_xy(:,:,:) -real(kind=WP), allocatable :: tr_z(:,:) - -! Auxiliary arrays for vector-invariant form of momentum advection -real(kind=WP), allocatable,dimension(:,:) :: vorticity - -!Viscosity and diff coefs -real(kind=WP), allocatable,dimension(:,:) :: Av,Kv -real(kind=WP), allocatable,dimension(:,:,:) :: Kv_double -real(kind=WP), allocatable,dimension(:) :: Kv0 -!Velocities interpolated to nodes -real(kind=WP), allocatable,dimension(:,:,:) :: Unode - -! Auxiliary arrays to store Redi-GM fields -real(kind=WP), allocatable,dimension(:,:,:) :: neutral_slope -real(kind=WP), allocatable,dimension(:,:,:) :: slope_tapered -real(kind=WP), allocatable,dimension(:,:,:) :: sigma_xy -real(kind=WP), allocatable,dimension(:,:) :: sw_beta, sw_alpha -real(kind=WP), allocatable,dimension(:) :: dens_flux -!real(kind=WP), allocatable,dimension(:,:,:) :: tsh, tsv, tsh_nodes -!real(kind=WP), allocatable,dimension(:,:) :: hd_flux,vd_flux -!Isoneutral diffusivities (or xy diffusivities if Redi=.false) -real(kind=WP), allocatable :: Ki(:,:) - -! --> auxiliary array to store an intermediate part of the rhs computations. -real(kind=WP), allocatable,dimension(:) :: ssh_rhs_old !, ssh_rhs_old2 !PS -real(kind=WP) :: is_nonlinfs - -!_______________________________________________________________________________ -! Arrays added for pressure gradient force calculation -real(kind=WP), allocatable,dimension(:,:) :: density_m_rho0 -real(kind=WP), allocatable,dimension(:,:) :: density_m_rho0_slev -real(kind=WP), allocatable,dimension(:,:) :: density_ref -real(kind=WP), allocatable,dimension(:,:) :: density_dmoc -real(kind=WP), allocatable,dimension(:,:) :: pgf_x, pgf_y - -!_______________________________________________________________________________ -!!PS ! dummy arrays -real(kind=WP), allocatable,dimension(:,:) :: dum_3d_n !, dum_3d_e -!!PS real(kind=WP), allocatable,dimension(:) :: dum_2d_n, dum_2d_e - -!_______________________________________________________________________________ -!Monin-Obukhov correction -real(kind=WP),allocatable :: mo(:,:),mixlength(:) -!GM_stuff -real(kind=WP),allocatable :: bvfreq(:,:),mixlay_dep(:),bv_ref(:) - -real(kind=WP), allocatable :: fer_UV(:,:,:), fer_wvel(:,:) -real(kind=WP), target, allocatable :: fer_c(:), fer_scal(:), fer_K(:,:), fer_gamma(:,:,:) - -real(kind=WP), allocatable :: ice_rejected_salt(:) -END MODULE o_ARRAYS -!========================================================== diff --git a/test_pgierz/actual_ugrid_ncdump.txt b/test/ugrid_complicance/actual_ugrid_ncdump.txt similarity index 100% rename from test_pgierz/actual_ugrid_ncdump.txt rename to test/ugrid_complicance/actual_ugrid_ncdump.txt diff --git a/test_pgierz/example_ugrid_ncdump.txt b/test/ugrid_complicance/example_ugrid_ncdump.txt similarity index 100% rename from test_pgierz/example_ugrid_ncdump.txt rename to test/ugrid_complicance/example_ugrid_ncdump.txt diff --git a/test_pgierz/test_one_timestep.sh b/test_pgierz/test_one_timestep.sh deleted file mode 100755 index e28dc5724..000000000 --- a/test_pgierz/test_one_timestep.sh +++ /dev/null @@ -1,32 +0,0 @@ -#!/bin/bash -# Run simples FESOM2 test in a container. -# -# With singularity on ollie -# -# module load singularity/3.9.1 -# cd fesom2 -# singularity exec /home/ollie/nkolduno/SINGULARITY/fesom_refactoring.sif ./test.sh -# -# With docker on Linux/Mac -# docker run -it -v "$(pwd)"/fesom2:/fesom/fesom2 koldunovn/fesom2_test:refactoring /bin/bash -# cd fesom2 -# ./test.sh -# - -set -e - -machine="docker" -tests="test_pi" - -for test in $tests; do - - ./configure.sh ubuntu - echo $test - mkrun pi $test -m $machine - cd work_pi_one_timestep - chmod +x job_docker_new - ./job_docker_new_one_timestep - fcheck . - cd ../ - -done