diff --git a/mesh_part/CMakeLists.txt b/mesh_part/CMakeLists.txt index 05bd5dfde..7938ae9d1 100644 --- a/mesh_part/CMakeLists.txt +++ b/mesh_part/CMakeLists.txt @@ -46,6 +46,7 @@ set(sources_Fortran ${src_home}/gen_modules_config.F90 ${src_home}/gen_modules_partitioning.F90 ${src_home}/gen_modules_rotate_grid.F90 + ${src_home}/mod_mesh_utils.F90 ${src_home}/fvom_init.F90 ${src_home}/oce_local.F90 ${src_home}/gen_comm.F90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 17ec8405c..2d2c0ea46 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -177,7 +177,7 @@ add_custom_target(fesom_version_info-generated.F90 ALL #if(FESOM_STANDALONE) # list(REMOVE_ITEM sources_Fortran ${src_home}/cpl_driver.F90) #endif() -list(REMOVE_ITEM sources_Fortran ${src_home}/fvom_init.F90 ${src_home}/oce_local.F90 ${src_home}/gen_comm.F90) +list(REMOVE_ITEM sources_Fortran ${src_home}/fvom_init.F90 ${src_home}/mod_mesh_utils.F90 ${src_home}/oce_local.F90 ${src_home}/gen_comm.F90) list(REMOVE_ITEM sources_Fortran ${src_home}/fesom_main.F90) find_package( NETCDF REQUIRED COMPONENTS C Fortran ) diff --git a/src/MOD_MESH.F90 b/src/MOD_MESH.F90 index 04cff77f9..7d9b080ce 100644 --- a/src/MOD_MESH.F90 +++ b/src/MOD_MESH.F90 @@ -1,6 +1,6 @@ !========================================================== MODULE MOD_MESH -USE O_PARAM +USE O_PARAM, only : WP USE MOD_WRITE_BINARY_ARRAYS USE MOD_READ_BINARY_ARRAYS USE, intrinsic :: ISO_FORTRAN_ENV, only : int32 diff --git a/src/MOD_PARTIT.F90 b/src/MOD_PARTIT.F90 index 2e8330a4b..1459cd863 100644 --- a/src/MOD_PARTIT.F90 +++ b/src/MOD_PARTIT.F90 @@ -1,7 +1,6 @@ !========================================================== ! Variables to organize parallel work module MOD_PARTIT -USE O_PARAM USE, intrinsic :: ISO_FORTRAN_ENV, only : int32 USE MOD_WRITE_BINARY_ARRAYS USE MOD_READ_BINARY_ARRAYS diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index d5ab69c69..06e17fdf8 100644 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -14,8 +14,9 @@ module cpl_driver ! use mod_oasis ! oasis module use g_config, only : dt, use_icebergs, lwiso, compute_oasis_corners - use o_param, only : rad + use o_param, only : rad, WP USE MOD_PARTIT + use oce_mesh_module, only : elem_center, edge_center use mpi implicit none save @@ -393,7 +394,7 @@ subroutine cpl_oasis3mct_define_unstr(partit, mesh) use mod_oasis, only: oasis_write_area, oasis_write_mask implicit none save - type(t_mesh), intent(in), target :: mesh + type(t_mesh), intent(in), target :: mesh type(t_partit), intent(inout), target :: partit !------------------------------------------------------------------- ! Definition of grid and field information for ocean diff --git a/src/fesom_meshdiag.F90 b/src/fesom_meshdiag.F90 index fc11a6fb3..7a0d0e7be 100644 --- a/src/fesom_meshdiag.F90 +++ b/src/fesom_meshdiag.F90 @@ -19,6 +19,8 @@ program fesom_meshdiag use g_config use g_comm_auto use io_mesh_info + use oce_setup_step_module, only: tracer_init, ocean_setup, dynamics_init, arrays_init + use oce_mesh_module, only: mesh_setup, check_mesh_consistency use, intrinsic :: iso_fortran_env, only : real32 implicit none diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 2ed97d36e..2c577d95c 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -19,15 +19,15 @@ module fesom_main_storage_module use io_mesh_info use diagnostics use mo_tidal - use tracer_init_interface - use ocean_setup_interface + use oce_setup_step_module, only: tracer_init, ocean_setup, dynamics_init, arrays_init use ice_setup_interface use ocean2ice_interface use oce_fluxes_interface - use update_atm_forcing_interface - use before_oce_step_interface - use oce_timestep_ale_interface - use read_mesh_interface + use forcing_coupling_interfaces, only: update_atm_forcing + use oce_setup_step_module, only: before_oce_step + use oce_ale_module, only: oce_timestep_ale, restart_thickness_ale + use oce_mesh_module, only: mesh_setup, check_mesh_consistency + use oce_dyn_module, only: compute_vel_nodes use fesom_version_info_module use command_line_options_module use, intrinsic :: iso_fortran_env, only : real32 diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 76a8891c0..ba5044cf7 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -12,82 +12,20 @@ !> Main driver routine for initialization program MAIN - use o_PARAM - use MOD_MESH - use MOD_PARTIT - use g_CONFIG - use g_rotate_grid + use o_PARAM, only: MAX_PATH + use MOD_MESH, only: t_mesh + use MOD_PARTIT, only: t_partit + use g_CONFIG, only: paths, geometry, run_config, machine, cyclic_length, rad, & + alphaEuler, betaEuler, gammaEuler, use_cavity + use g_rotate_grid, only: set_mesh_transform_matrix use iso_fortran_env, only: error_unit + use mod_mesh_utils, only: read_mesh_ini, read_mesh_cavity, test_tri_ini, & + find_edges_ini, find_elem_neighbors_ini, find_levels, & + find_levels_cavity, stiff_mat_ini, set_par_support_ini, & + communication_ini implicit none -interface - subroutine read_mesh_ini(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine read_mesh_ini -end interface - -interface - subroutine test_tri_ini(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine test_tri_ini -end interface -interface - subroutine find_edges_ini(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine find_edges_ini -end interface -interface - subroutine find_elem_neighbors_ini(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine find_elem_neighbors_ini -end interface -interface - subroutine find_levels(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine find_levels -end interface -interface - subroutine stiff_mat_ini(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine stiff_mat_ini -end interface -interface - subroutine set_par_support_ini(partit, mesh) - use mod_mesh - use mod_partit - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine set_par_support_ini -end interface -interface - subroutine communication_ini(partit, mesh) - use mod_mesh - use mod_partit - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine communication_ini -end interface - -interface - subroutine read_mesh_cavity(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine read_mesh_cavity -end interface - -interface - subroutine find_levels_cavity(mesh) - use mod_mesh - type(t_mesh), intent(inout) , target :: mesh - end subroutine find_levels_cavity -end interface character(len=MAX_PATH) :: nmlfile !> name of configuration namelist file integer :: start_t, interm_t, finish_t, rate_t @@ -157,1760 +95,3 @@ end subroutine find_levels_cavity print '("**** Total time = ",f12.3," seconds. ****")', & real(finish_t-start_t)/real(rate_t) end program MAIN -!============================================================================= -!> @brief -!> Reads mesh files -subroutine read_mesh_ini(mesh) -USE MOD_MESH -USE o_PARAM -use g_CONFIG -use g_rotate_grid -! -IMPLICIT NONE -! -type(t_mesh), intent(inout), target :: mesh -INTEGER :: nq -INTEGER :: n1,n2,n3 -INTEGER :: n, nz, exit_flag -REAL(kind=WP) :: x1, x2, gx1, gx2 -INTEGER :: tag -INTEGER, allocatable :: elem_data(:) -INTEGER :: i_error -#include "associate_mesh_ini.h" -! =================== -! Surface mesh -! =================== - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: read elem2d.out & nod2d.out '//achar(27)//'[0m' - - open (20,file=trim(meshpath)//'nod2d.out', status='old') - open (21,file=trim(meshpath)//'elem2d.out', status='old') - READ(20,*) mesh%nod2D - ALLOCATE(mesh%coord_nod2D(2,mesh%nod2D)) - coord_nod2D => mesh%coord_nod2D !required after the allocation, otherwise the pointer remains undefined - - do n=1, mesh%nod2D - read(20,*) nq, x1, x2, tag - x1=x1*rad - x2=x2*rad - if (force_rotation) then - gx1=x1 - gx2=x2 - call g2r(gx1, gx2, x1, x2) - end if - mesh%coord_nod2D(1,n)=x1 - mesh%coord_nod2D(2,n)=x2 - end do - CLOSE(20) - READ(21,*) mesh%elem2D - ALLOCATE(mesh%elem2D_nodes(4,mesh%elem2D)) - elem2D_nodes => mesh%elem2D_nodes !required after the allocation, otherwise the pointer remains undefined - ALLOCATE(elem_data(4*mesh%elem2D)) - elem_data(:)=-1 - - ! meshes with quads have 4 columns, but TsunAWI grids may be - ! purely triangular, with 3 columns each. Test, how many - ! columns there are! - read(21,*,iostat=i_error) elem_data(1:4*mesh%elem2D) - if (i_error == 0) then ! There is a fourth column => quad or mixed mesh (not working yet!) - mesh%elem2D_nodes = reshape(elem_data, shape(mesh%elem2D_nodes)) - else ! No fourth column => triangles only - mesh%elem2D_nodes(1:3,:) = reshape(elem_data, shape(mesh%elem2D_nodes(1:3,:))) - mesh%elem2D_nodes(4,:) = mesh%elem2D_nodes(1,:) - end if - - deallocate(elem_data) - CLOSE(21) - -write(*,*) '=========================' -write(*,*) 'Mesh is read' -write(*,*) '=========================' -END SUBROUTINE read_mesh_ini -!============================================================================= -!> @brief -!> Reads mesh files -subroutine read_mesh_cavity(mesh) - use mod_mesh - use o_PARAM - use g_CONFIG - implicit none - - type(t_mesh), intent(inout), target :: mesh - integer :: node, auxi - character(len=MAX_PATH) :: fname - logical :: file_exist=.False. -#include "associate_mesh_ini.h" - - !___________________________________________________________________________ - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: read cavity depth '//achar(27)//'[0m' - - !___________________________________________________________________________ - ! read depth of cavity-ocean boundary - if (use_cavityonelem) then - fname = trim(meshpath)//'cavity_depth@elem.out' - else - fname = trim(meshpath)//'cavity_depth@node.out' - end if - file_exist=.False. - inquire(file=trim(fname),exist=file_exist) - if (file_exist) then - open (21,file=fname, status='old') - if (use_cavityonelem) then - allocate(mesh%cavity_depth(mesh%elem2d)) - else - allocate(mesh%cavity_depth(mesh%nod2D)) - end if - cavity_depth => mesh%cavity_depth - else - write(*,*) '____________________________________________________________________' - write(*,*) ' ERROR: could not find cavity file:', fname - write(*,*) ' --> stop partitioning here !' - write(*,*) '____________________________________________________________________' - stop - end if - - !___________________________________________________________________________ - auxi=mesh%nod2D - if (use_cavityonelem) auxi=mesh%elem2d - do node=1, auxi - read(21,*) mesh%cavity_depth(node) - end do - - !___________________________________________________________________________ - close(21) - -end subroutine read_mesh_cavity - -!======================================================================= -!> @brief -!> Check the order of nodes in triangles; correct it if necessary to make -!! it same sense (clockwise) -SUBROUTINE test_tri_ini(mesh) -USE MOD_MESH -USE o_PARAM -USE g_CONFIG -IMPLICIT NONE - -real(kind=WP) :: a(2), b(2), c(2), r -integer :: n, nx, elnodes(3) -type(t_mesh), intent(inout), target :: mesh -#include "associate_mesh_ini.h" - - DO n=1, elem2D - elnodes=elem2D_nodes(1:3,n) - - a=coord_nod2D(:,elnodes(1)) - b=coord_nod2D(:,elnodes(2))-a - c=coord_nod2D(:,elnodes(3))-a - - if(b(1)>cyclic_length/2.) b(1)=b(1)-cyclic_length - if(b(1)<-cyclic_length/2.) b(1)=b(1)+cyclic_length - if(c(1)>cyclic_length/2.) c(1)=c(1)-cyclic_length - if(c(1)<-cyclic_length/2.) c(1)=c(1)+cyclic_length - - - r=b(1)*c(2)-b(2)*c(1) - if (r>0) then - ! Vector b is to right of c - ! Exchange second and third nodes: - - nx=elnodes(2) - elnodes(2)=elnodes(3) - elnodes(3)=nx - elem2D_nodes(1:3,n)=elnodes - end if - END DO -END SUBROUTINE test_tri_ini -!========================================================================= -!> @brief -!> Finds edges. Creates 3 files: edgenum.out, edges.out, edge_tri.out -SUBROUTINE find_edges_ini(mesh) -USE MOD_MESH -USE o_PARAM -USE g_CONFIG -use g_rotate_grid -IMPLICIT NONE - -interface - subroutine elem_center(elem, x, y, mesh) - USE MOD_MESH - USE g_CONFIG - integer, intent(in) :: elem - real(kind=WP), intent(out) :: x, y - type(t_mesh), intent(in), target :: mesh - end subroutine elem_center -end interface - -integer, allocatable :: aux1(:), ne_num(:), ne_pos(:,:), nn_num(:), nn_pos(:,:) -integer :: counter, counter_in, n, k, q -integer :: elem, elem1, elems(2), q1, q2 -integer :: elnodes(4), ed(2), flag, eledges(4) -integer :: temp(100), node -real(kind=WP) :: xc(2), xe(2), ax(3), amin -type(t_mesh), intent(inout), target :: mesh -#include "associate_mesh_ini.h" -! ==================== -! (a) find edges. To make the procedure fast -! one needs neighbourhood arrays -! ==================== -print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' -print *, achar(27)//'[7;1m' //' -->: compute edge connectivity '//achar(27)//'[0m' - -allocate(ne_num(nod2d), ne_pos(MAX_ADJACENT, nod2D), nn_num(nod2D)) -ne_num=0 -DO n=1,elem2D - elnodes=elem2D_nodes(:,n) - q1=4 - if(elnodes(1)==elnodes(4)) q1=3 - DO q=1,q1 - ne_num(elnodes(q))=ne_num(elnodes(q))+1 - if (ne_num(elnodes(q)) > MAX_ADJACENT ) then - print *,'Recompile with larger value for MAX_ADJACENT.' - stop - else - ne_pos(ne_num(elnodes(q)),elnodes(q))=n - endif - END Do -END DO ! neighbor elements are found - -! count neighbour nodes -! In quads we should count the nodes that are -! connected by edges! -allocate(aux1(nod2D)) -aux1=0 - -DO n=1, nod2D - counter=0 - DO k=1, ne_num(n) - elem=ne_pos(k,n) - elnodes=elem2D_nodes(:,elem) - if(elnodes(1)==elnodes(4)) then - DO q=1,3 - if(elnodes(q)==n) CYCLE - if(aux1(elnodes(q)).ne.1) then - counter=counter+1 - aux1(elnodes(q))=1 - temp(counter)=elnodes(q) - end if - END DO - else - ! Find the position of n in elnodes: - if(elnodes(1)==n .or. elnodes(3)==n) then - ed(1)=elnodes(2) - ed(2)=elnodes(4) - else - ed(1)=elnodes(1) - ed(2)=elnodes(3) - end if - DO q=1,2 - if(aux1(ed(q)).ne.1) then - counter=counter+1 - aux1(ed(q))=1 - temp(counter)=ed(q) - end if - END DO - end if - END DO - nn_num(n)=counter - aux1(temp(1:counter))=0 -END DO - -allocate(nn_pos(maxval(nn_num)+1,nod2D)) -nn_pos = -1 -aux1=0 - -DO n=1, nod2D - counter=0 - DO k=1, ne_num(n) - elem=ne_pos(k,n) - elnodes=elem2D_nodes(:,elem) - if(elnodes(1)==elnodes(4)) then - DO q=1,3 - if(elnodes(q)==n) CYCLE - if(aux1(elnodes(q)).ne.1) then - counter=counter+1 - aux1(elnodes(q))=1 - temp(counter)=elnodes(q) - end if - END DO - else - ! Find the position of n in elnodes: - if(elnodes(1)==n .or. elnodes(3)==n) then - ed(1)=elnodes(2) - ed(2)=elnodes(4) - else - ed(1)=elnodes(1) - ed(2)=elnodes(3) - end if - DO q=1,2 - if(aux1(ed(q)).ne.1) then - counter=counter+1 - aux1(ed(q))=1 - temp(counter)=ed(q) - end if - END DO - end if - END DO - nn_num(n)=counter+1 - aux1(temp(1:counter))=0 - nn_pos(1,n)=n - nn_pos(2:counter+1,n)=temp(1:counter) -END DO -deallocate(aux1) -! neighboring nodes are found. First in the list is the node itself - -! ==================== -! (b) Find edges and elements containing them. -! Write information to auxiliary file -! ==================== -! open(10, file='edges.out') - - ! Count edges: - ! ==================== - ! form edges with n by cycling over neighbor - ! nodes (if edges are not accounted yet). - ! New edges are added only if neighbor>n - ! ==================== - counter = 0 - DO n=1,nod2D - counter = counter + count(nn_pos(2:nn_num(n),n)>nn_pos(1,n)) - end do - edge2D=counter - - allocate(mesh%edges (2, edge2D)) - allocate(mesh%edge_tri(2, edge2D)) - edges => mesh%edges !required after the allocation, otherwise the pointer remains undefined - edge_tri => mesh%edge_tri !required after the allocation, otherwise the pointer remains undefined - counter_in=0 - - DO n=1,nod2D - DO q=2,nn_num(n) - node=nn_pos(q,n) - if(node0.0_WP) then - ! Vector drawn to the center of the first triangle is to the right - ! of the edge vector. Triangles have to be exchanged: - elem=edge_tri(1,n) - elem1=edge_tri(2,n) - if(elem1>0) then !TODO - edge_tri(1,n)=elem1 - edge_tri(2,n)=elem - else - if (elem<=0) write(*,*) '2 neighbouring elems are on the ground.' - elem=edges(2,n) ! change the order of nodes - edges(2,n)=edges(1,n) - edges(1,n)=elem - end if - end if - END DO - - ! ==================== - ! (e) We need an array inverse to edge_tri listing edges - ! of a given triangle - ! ==================== - allocate(mesh%elem_edges(4,elem2D)) - elem_edges => mesh%elem_edges !required after the allocation, otherwise the pointer remains undefined - allocate(aux1(elem2D)) - aux1=0 - DO n=1, edge2D - DO k=1,2 - q=edge_tri(k,n) ! triangle number - if (q>0) then - aux1(q)=aux1(q)+1 - elem_edges(aux1(q),q)=n - end if - END DO - END DO - deallocate(aux1) - ! We order the edges in this list so that they - ! are listed in the same rotation sense as nodes. - ! First is the edge formed by elnodes(1:2), and so on - DO elem=1,elem2D - elnodes=elem2D_nodes(:,elem) - q1=4 - if(elnodes(1)==elnodes(4)) q1=3 - eledges=elem_edges(:,elem) - DO q=1,q1-1 - DO k=1,q1 - if(((edges(1,eledges(k))==elnodes(q)).and. & - (edges(2,eledges(k))==elnodes(q+1))).or. & - ((edges(1,eledges(k))==elnodes(q+1)).and. & - (edges(2,eledges(k))==elnodes(q)))) then - elem_edges(q,elem)=eledges(k) - exit - end if - END DO - END DO - DO k=1,q1 - if(((edges(1,eledges(k))==elnodes(q1)).and. & - (edges(2,eledges(k))==elnodes(1))).or. & - ((edges(1,eledges(k))==elnodes(1)).and. & - (edges(2,eledges(k))==elnodes(q1)))) then - elem_edges(q1,elem)=eledges(k) - exit - end if - END DO - if(q1==3) elem_edges(4,elem)=elem_edges(1,elem) - END DO - - !> The edge and elem lists agree in the sense that edge 1 does not - !! contain node 1 and so on - open(11, file=trim(meshpath)//'edgenum.out') - write(11,*) edge2D - write(11,*) edge2D_in - close(11) - open(10, file=trim(meshpath)//'edges.out') - open(12, file=trim(meshpath)//'edge_tri.out') - do n=1,edge2D - write(10,*) edges(:,n) - write(12,*) edge_tri(:,n) - end do - close(10) - close(12) - deallocate(ne_num, ne_pos) -END SUBROUTINE find_edges_ini -!========================================================================= -!> @brief -!> Finds elemental and nodal levels. -!> Does some thresholding: if (depth>zbar(4)) x=zbar(4) -!> Fixes rough topography, by converting some oceans cells to ground cell(reflected by changing levels arrays) -!> Creates 2 files: elvls.out, nlvls.out -subroutine find_levels(mesh) - use g_config - use mod_mesh - implicit none - INTEGER :: nodes(3), elems(3), eledges(3) - integer :: elem1, j, n, nneighb, q, node, i, nz, auxi - integer :: count_iter, count_neighb_open, exit_flag, fileID=111 - real(kind=WP) :: x, dmean - logical :: file_exist - integer :: max_iter=1000 - character(MAX_PATH) :: file_name - type(t_mesh), intent(inout), target :: mesh - -#include "associate_mesh_ini.h" - - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: read bottom depth '//achar(27)//'[0m' - - !___________________________________________________________________________ - ! allocate depth - if (use_depthonelem) then - allocate(mesh%depth(elem2D)) - else - allocate(mesh%depth(nod2D)) - end if - depth => mesh%depth !required after the allocation, otherwise the pointer remains undefined - - !______________________________________________________________________________ - ! read depth from aux3d.out - if (trim(use_depthfile)=='aux3d') then - ! check if aux3d.out file does exist - file_exist=.False. - file_name=trim(meshpath)//'aux3d.out' - inquire(file=trim(file_name),exist=file_exist) - !_______________________________________________________________________ - if (file_exist) then - write(*," (A, A)") ' read file:',trim(file_name) - !___________________________________________________________________ - ! load fesom2.0 aux3d.out file - open(fileID, file=file_name) - - ! read the number of levels - read(fileID,*) nl - allocate(mesh%zbar(nl)) ! their standard depths - - ! read full depth levels - zbar => mesh%zbar !required after the allocation, otherwise the pointer remains undefined - read(fileID,*) zbar - if(zbar(2)>0) zbar=-zbar ! zbar is negative - - ! compute mid depth levels - allocate(mesh%Z(nl-1)) - Z => mesh%Z !required after the allocation, otherwise the pointer remains undefined - Z=zbar(1:nl-1)+zbar(2:nl) ! mid-depths of cells - Z=0.5_WP*Z - else - write(*,*) '____________________________________________________________________' - write(*,*) ' ERROR: You want to use aux3d.out file to define your depth, but ' - write(*,*) ' the file seems not to exist' - write(*,*) ' --> check in namelist.config, the flag use_depthfile must be' - write(*,*) ' use_depthfile= "aux3d" or "depth@" ' - write(*,*) ' --> model stops here' - write(*,*) '____________________________________________________________________' - stop - end if - !___________________________________________________________________________ - ! read depth from depth@node.out or depth@elem.out - elseif (trim(use_depthfile)=='depth@') then - !_______________________________________________________________________ - ! load file depth_zlev.out --> contains number of model levels and full depth - ! levels - file_exist=.False. - file_name=trim(meshpath)//'depth_zlev.out' - inquire(file=trim(file_name),exist=file_exist) - if (file_exist) then - write(*," (A, A)") ' read file:',trim(file_name) - !___________________________________________________________________ - ! load fesom2.0 aux3d.out file - open(fileID, file=file_name) - - ! read the number of levels - read(fileID,*) nl - allocate(mesh%zbar(nl)) ! their standard depths - - ! read full depth levels - zbar => mesh%zbar !required after the allocation, otherwise the pointer remains undefined - read(fileID,*) zbar - if(zbar(2)>0) zbar=-zbar ! zbar is negative - - ! compute mid depth levels - allocate(mesh%Z(nl-1)) - Z => mesh%Z !required after the allocation, otherwise the pointer remains undefined - Z=zbar(1:nl-1)+zbar(2:nl) ! mid-depths of cells - Z=0.5_WP*Z - - close(fileID) - else - write(*,*) '____________________________________________________________________' - write(*,*) ' ERROR: You want to use depth@elem.out or depth@node.out file, therefore' - write(*,*) ' you also need the file depth_zlev.out which contains the model ' - write(*,*) ' number of layers and the depth of model levels. This file seems ' - write(*,*) ' not to exist' - write(*,*) ' --> check in namelist.config, the flag use_depthfile must be' - write(*,*) ' use_depthfile= "aux3d" or "depth@" and your meshfolder' - write(*,*) ' --> model stops here' - write(*,*) '____________________________________________________________________' - stop - endif - - !_______________________________________________________________________ - ! load file depth@elem.out or depth@node.out contains topography either at - ! nodes or elements - if (use_depthonelem) then - file_name=trim(meshpath)//'depth@elem.out' - else - file_name=trim(meshpath)//'depth@node.out' - end if - inquire(file=trim(file_name),exist=file_exist) - if (file_exist) then - write(*," (A, A)") ' read file:',trim(file_name) - open(fileID, file=file_name) - else - write(*,*) '____________________________________________________________________' - write(*,*) ' ERROR: You want to use depth@elem.out or depth@node.out file to ' - write(*,*) ' define your depth, but the file seems not to exist' - write(*,*) ' --> check in namelist.config, the flag use_depthfile must be' - write(*,*) ' use_depthfile= "aux3d" or "depth@" and your meshfolder ' - write(*,*) ' --> model stops here' - write(*,*) '____________________________________________________________________' - stop - end if - end if - - !___________________________________________________________________________ - ! read topography from file - auxi = nod2d - if (use_depthonelem) auxi = elem2d -! write(*,*) ' use_depthonelem = ',use_depthonelem -! write(*,*) ' auxi =',auxi - DO n = 1, auxi - read(fileID,*) x - if (x>0) x=-x - if (x>zbar(thers_zbar_lev)) x=zbar(thers_zbar_lev) !TODO KK thresholding for depth - depth(n)=x - END DO - close(fileID) - if(depth(2)>0) depth=-depth ! depth is negative - - !___________________________________________________________________________ - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: compute elem, vertice bottom depth index '//achar(27)//'[0m' - - allocate(mesh%nlevels(elem2D)) - nlevels => mesh%nlevels !required after the allocation, otherwise the pointer remains undefined - allocate(mesh%nlevels_nod2D(nod2D)) - nlevels_nod2D => mesh%nlevels_nod2D !required after the allocation, otherwise the pointer remains undefined - - !___________________________________________________________________________ - ! Compute the initial number number of elementa levels, based on the vertice - ! depth information - do n=1, elem2D - - !_______________________________________________________________________ - if (use_depthonelem) then - dmean = depth(n) ! depth is already defined on elements - else - nodes=elem2D_nodes(1:3,n) - !___________________________________________________________________ - ! depth of element is shallowest depth of sorounding vertices - if (trim(which_depth_n2e) .eq. 'min') then ; dmean=maxval(depth(nodes)) - ! depth of element is deepest depth of sorounding vertices - elseif (trim(which_depth_n2e) .eq. 'max') then ; dmean=minval(depth(nodes)) - ! DEFAULT: depth of element is mean depth of sorounding vertices - elseif (trim(which_depth_n2e) .eq. 'mean') then; dmean=sum(depth(nodes))/3.0 - end if - end if - - !_______________________________________________________________________ - exit_flag=0 - do nz=1,nl-1 - if(Z(nz)=0) nlevels(n)=thers_zbar_lev - - ! set minimum number of levels to --> thers_lev=5 - if(nlevels(n) do n=1, elem2D - - !___________________________________________________________________________ - ! write out vertical level indices before iterative geometric adaption to - ! exclude isolated cells - !_______________________________________________________________________ - file_name=trim(meshpath)//'elvls_raw.out' - open(fileID, file=file_name) - do n=1,elem2D - write(fileID,*) nlevels(n) - end do - close(fileID) - - !___________________________________________________________________________ - ! check for isolated cells (cells with at least two boundary faces or three - ! boundary vertices) and eliminate them --> FESOM2.0 doesn't like these kind - ! of cells - do nz=thers_zbar_lev+1,nl - exit_flag=0 - count_iter=0 - - !_______________________________________________________________________ - ! iteration loop within each layer - do while((exit_flag==0).and.(count_iter if elem2D_nodes(1,n) == elem2D_nodes(4,n): True --> q=3 --> triangular mesh - ! --> if elem2D_nodes(1,n) == elem2D_nodes(4,n): False --> q=4 --> quad mesh - nneighb = merge(3,4,elem2D_nodes(1,n) == elem2D_nodes(4,n)) - ! - ! +---isolated bottom cell - ! ._______________ | _______________________. - ! |###|###|###|###|___ | ___|###|###|###|###|###|###| - ! |###|###|###|###|###| | ___|###|###|###|###|###|###|###| - ! |###|###|###|###|###| | |###|###|###|###| BOTTOM |###| - ! |###|###|###|###|###|_v_|###|###|###|###|###|###|###|###| - ! |###|###|###|###|###|###|###|###|###|###|###|###|###|###| - ! - if (nlevels(n)>=nz) then - count_neighb_open=0 - elems=elem_neighbors(1:3,n) - !___________________________________________________________ - ! loop over neighbouring triangles - do i=1,nneighb - if (elems(i)>0) then - if (nlevels(elems(i))>=nz) then - !count neighbours - count_neighb_open=count_neighb_open+1 - endif - endif - enddo - - !___________________________________________________________ - ! check how many open faces to neighboring triangles the cell - ! has, if there are less than 2 its isolated (a cell should - ! have at least 2 valid neighbours) - if (count_neighb_open<2) then - ! if cell is "isolated", and the one levels shallower bottom - ! cell would be shallower than the minimum vertical level - ! treshhold (thers_lev). --> in this make sorrounding elements - ! one level deeper to reconnect the isolated cell - if (nz-10) then - nlevels(elems(i)) = max(nlevels(elems(i)),nz) - end if - end do - - !if cell is "isolated" convert to one level shallower bottom cell - else - nlevels(n)=nz-1 - end if - !force recheck for all current ocean cells - exit_flag=0 - - end if - end if ! --> if (nlevels(n)>=nz) then - end do ! --> do n=1,elem2D - end do ! --> do while((exit_flag==0).and.(count1<1000)) - write(*,"(A, I5, A, i5, A, I3)") ' -[iter ]->: nlevel, iter/maxiter=',count_iter,'/',max_iter,', nz=',nz - end do ! --> do nz=4,nl - - !___________________________________________________________________________ - ! vertical vertice level index of ocean bottom boundary - write(*,"(A)" ) ' -[compu]->: nlevels_nod2D ' - nlevels_nod2D=0 - do n=1,elem2D - q = merge(3,4,elem2D_nodes(1,n) == elem2D_nodes(4,n)) - do j=1,q - node=elem2D_nodes(j,n) - if(nlevels_nod2D(node): compute elem, vertice cavity depth index '//achar(27)//'[0m' - - !___________________________________________________________________________ - allocate(mesh%ulevels(elem2D)) - ulevels => mesh%ulevels - allocate(mesh%ulevels_nod2D(nod2D)) - ulevels_nod2D => mesh%ulevels_nod2D - - !___________________________________________________________________________ - ! Compute level position of ocean-cavity boundary - cavity_maxlev=0 - do elem=1, elem2D - - !_______________________________________________________________________ - if (use_cavityonelem) then - dmean = cavity_depth(elem) - else - nodes=elem2D_nodes(1:3,elem) - !_______________________________________________________________________ - ! depth of element is shallowest depth of sorounding vertices - if (trim(which_depth_n2e) .eq. 'min') then ; dmean=maxval(cavity_depth(nodes)) - ! depth of element is deepest depth of sorounding vertices - elseif (trim(which_depth_n2e) .eq. 'max') then ; dmean=minval(cavity_depth(nodes)) - ! DEFAULT: depth of element is mean depth of sorounding vertices - elseif (trim(which_depth_n2e) .eq. 'mean') then ; dmean=sum(cavity_depth(nodes))/3.0 - end if - end if - - !_______________________________________________________________________ - ! vertical elem level index of cavity-ocean boundary - ulevels(elem) = 1 - if (dmean<0.0_WP) ulevels(elem) = 2 - - do nz=1,nlevels(elem)-1 - if (Z(nz) should not be - ! possible in FESOM2.0 - ! loop over all cavity levels - allocate(elemreducelvl(elem2d),elemfixlvl(elem2d)) - allocate(numelemtonode(nl)) - - !___________________________________________________________________________ - ! outer iteration loop - count_iter2 = 0 - exit_flag2 = 0 - elemfixlvl = .False. - do while((exit_flag2==0) .and. (count_iter2 tri mesh, nneighb = 4 --> quad mesh - nneighb = merge(3,4,elem2D_nodes(1,elem) == elem2D_nodes(4,elem)) - elems = elem_neighbors(1:3,elem) - ! - ! .___________________________.~~~~~~~~~~~~~~~~~~~~~~~~~~ - ! |###|###|###|###|###|###|###| - ! |# CAVITY |###| . |###|###| OCEAN - ! |###|###|###|###|/|\|###| - ! |###|###| | - ! |###| +-- Not good can lead to isolated cells - ! - - !___________________________________________________________ - ! (1st) Ask the Question: is nz at element elem an here an - ! valid layer in the ocean - if ( nz >= ulevels(elem) .and. nz0) then ! if its a valid boundary triangle, 0=missing value - ! check for isolated cell - if ( ulevels(elems(j))<= nz .and. & - nlevels(elems(j))> nz ) then - !count the open faces to neighboring cells - count_neighb_open=count_neighb_open+1 - endif - end if - end do ! --> do i = 1, nneighb - - !_______________________________________________________ - ! (3rd) check how many open faces to neighboring triangles the cell - ! has, if there are less than 2 its isolated (a cell should - ! have at least 2 valid neighbours) - ! --> in this case shift cavity-ocean interface one level down - if (count_neighb_open<=1) then - count_isoelem = count_isoelem+1 - - ! (4th.1 ) if cell elem is isolated convert it to a deeper ocean level - ! except when this levels would remain less than 3 valid - ! bottom levels --> in case make the levels of all sorounding - ! triangles shallower - if ( (nlevels(elem)-(nz+1))>=3 .and. & - (elemreducelvl(elem) .eqv. .False.) .and. & - (elemfixlvl( elem) .eqv. .False.) & - ) then - ulevels(elem)=nz+1 - - ! (4th.2) can not increase depth anymore to eleminate - ! isolated cell, otherwise lessthan 3 valid layers - ! therefor reduce depth of ONE!!! of the neighbouring - ! triangles. Choose triangle whos depth is already closest - ! to nz - else - !PS replace this with do j=1,3... because of - !PS indice -999 conflict in elems, ulevels(-999) - !PS not possible - !PS idx = minloc(ulevels(elems)-nz, 1, MASK=( (elems>0) .and. ((ulevels(elems)-nz)>0) ) ) - val=nl - do j = 1, 3 - if (elems(j)>0) then ! if its a valid boundary triangle, 0=missing value - if (ulevels(elems(j))-nz>0 .and. ulevels(elems(j))-nz do i = 1, nneighb - - ulevels( elems(idx)) = nz - elemreducelvl(elems(idx)) = .True. - end if - - ! force recheck for all current ocean cells - exit_flag1=0 - end if ! --> if (count_neighb_open<2) then - - end if ! --> if ( nz >= ulevels(elem) .and. nz do elem=1,elem2D - - end do ! --> do while((exit_flag==0).and.(count_iter<1000)) - write(*,"(A, I5, A, i5, A, I3)") ' -[iter ]->: ulevel, iter/maxiter=',count_iter,'/',max_iter,', nz=',nz - end do ! --> do nz=1,cavity_maxlev - - !_______________________________________________________________________ - ! compute vertical vertice level index of cavity_ocean boundary - write(*,"(A)" ) ' -[compu]->: ulevels_nod2D ' - ulevels_nod2D = nl - do elem=1,elem2D - nneighb = merge(3,4,elem2D_nodes(1,elem) == elem2D_nodes(4,elem)) - !___________________________________________________________________ - ! loop over neighbouring triangles - do j=1,nneighb - node=elem2D_nodes(j,elem) - ulevels_nod2D(node)=min(ulevels_nod2D(node),ulevels(elem)) - end do - end do ! --> do elem=1,elem2D - - !_______________________________________________________________________ - ! check if all constrains for ulevel and ulevels_nod2D is fullfilled - exit_flag2 = 1 - do elem=1,elem2D - if (ulevels(elem)>=nlevels(elem)) then - write(*,*) ' -[check]->: elem cavity depth deeper or equal bottom depth, elem=',elem - exit_flag2 = 0 - end if - - if (nlevels(elem)-ulevels(elem)<3) then - write(*,*) ' -[check]->: less than three valid elem ocean layers, elem=',elem - exit_flag2 = 0 - end if - end do ! --> do elem=1,elem2D - - do node=1,nod2D - if (ulevels_nod2D(node)>=nlevels_nod2D(node)) then - write(*,*) ' -[check]->: vertice cavity depth deeper or equal bottom depth, node=', node - exit_flag2 = 0 - end if - - if (nlevels_nod2D(node)-ulevels_nod2D(node)<3) then - write(*,*) ' -[check]->: less than three valid vertice ocean layers, node=', node - exit_flag2 = 0 - end if - end do ! --> do node=1,nod2D - - do elem=1,elem2D - if (ulevels(elem)< maxval(ulevels_nod2D(elem2D_nodes(:,elem))) ) then - write(*,*) ' -[check]->: found elem cavity shallower than its valid maximum cavity vertice depths, elem=', elem2d - exit_flag2 = 0 - end if - end do ! --> do elem=1,elem2D - - !_______________________________________________________________________ - ! compute how many triangle elements contribute to every vertice in every - ! layer - ! - ! --> What can happen is that a node point in the middle of the vertical - ! domain can become isolated due to the cavity constrains. The model - ! would not be able to deal with this kind of situation. So we must - ! prevent it by adapting ulevels! - ! O node - ! _._ - ! _/ | \_ - ! _/ | \_ - ! _/ | \_ - ! elem(1) elem(2) elem(3)... <--elem=nod_in_elem2D(j,node) - ! ._______. ulevel(elem2)=30 - ! |_______| - ! |_______| - ! |_______| - ! |_______| - ! |_______| nlevel(elem2)=38 - ! - ! In this possible gap node points - ! would have no neighboring elements - ! - ! ulevel(elem1)=42 ._______. ._______. ulevel(elem3)=42 - ! |_______| |_______| - ! |_______| |_______| - ! |_______| |_______| - ! |_______| |_______| - ! nlevel(elem1)=46 |_______| |_______| - ! |_______| nlevel(elem3)=48 - ! - ! --> Problem here is we want to keep nlevels fixed so what we can do is - ! to set ulevels(elem1) and ulevels(elem3) towards nlevel(elem2) - count_iter=0 - do node=1, nod2D - !___________________________________________________________________ - ! check if there is a possible gap as described above that would - ! allow for node points without neighbors - min_nlvl = nl - max_ulvl = 1 - do j=1, nod_in_elem2D_num(node) - elem=nod_in_elem2D(j, node) - min_nlvl = min(min_nlvl, nlevels(elem)) - max_ulvl = max(max_ulvl, ulevels(elem)) - end do - - ! found a potential gap - if (min_nlvl < max_ulvl) then - - !_______________________________________________________________ - ! compute how many triangle elements contribute to vertice in - ! every layer check if there are layers where a node point has - ! only one or even zero neighboring triangles. - numelemtonode=0 - do j=1, nod_in_elem2D_num(node) - elem=nod_in_elem2D(j, node) - do nz=ulevels(elem), nlevels(elem)-1 - numelemtonode(nz) = numelemtonode(nz) + 1 - end do - end do - - !_______________________________________________________________ - ! check in which depth level is an isolated node - nzloop: do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 - !___________________________________________________________ - ! nodes has zero neighbouring triangles and is completely - ! isolated need to adapt ulevels --> inflicts another - ! outher iteration loop (exit_flag2=0). It needs at least - ! two neighboring elements so the node is considered as - ! connected. Search the index of the two elements where ulevels>nz - ! but that are closest to nz - if (numelemtonode(nz)==0) then - count_iter = count_iter+1 - write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz) ,' neighb. triangle: n=', node, ', nz=',nz - !_______________________________________________________ - allocate(aux_arr(nod_in_elem2D_num(node)), aux_idx(nod_in_elem2D_num(node))) - aux_arr(:) = ulevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) - aux_arr(:) = aux_arr(:) - nz - ! fill array with index of element - do j=1, nod_in_elem2D_num(node) - aux_idx(j) = j - end do - ! index of closest elem to nz where ulevel>nz - idx = minloc(aux_arr, 1, MASK=((aux_arr>0)) ) - ! index of second closest elem to nz where ulevel>nz - idx2 = minloc(aux_arr, 1, MASK=((aux_arr>0) .and. (aux_idx/=idx)) ) - deallocate(aux_arr, aux_idx) - ulevels( nod_in_elem2D(idx ,node)) = nz - ulevels( nod_in_elem2D(idx2,node)) = nz - elemfixlvl(nod_in_elem2D(idx ,node)) = .True. - elemfixlvl(nod_in_elem2D(idx2,node)) = .True. - - !_______________________________________________________ - ! inflict another outer iteration loop - exit_flag2 = 0 - - !_______________________________________________________ - ! if the upper most isolated layer is fixed all layers below should be fixed as well - ! --> exit loop do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 - exit nzloop - - !___________________________________________________________ - ! nodes has one neighbouring triangles it needs at least - ! another neighboring elements so the node is considered as - ! connected - elseif (numelemtonode(nz)==1) then - count_iter = count_iter+1 - write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz) ,'neighb. triangle: n=', node, ', nz=',nz - !_______________________________________________________ - allocate(aux_arr(nod_in_elem2D_num(node)), aux_idx(nod_in_elem2D_num(node))) - aux_arr(:) = ulevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) - aux_arr(:) = aux_arr(:) - nz - ! fill array with index of element - do j=1, nod_in_elem2D_num(node) - aux_idx(j) = j - end do - ! index of closest elem to nz where ulevel>nz - idx = minloc(aux_arr, 1, MASK=((aux_arr>0)) ) - deallocate(aux_arr, aux_idx) - ulevels( nod_in_elem2D(idx,node)) = nz - elemfixlvl(nod_in_elem2D(idx,node)) = .True. - - !_______________________________________________________ - ! inflict another outer iteration loop - exit_flag2 = 0 - - !_______________________________________________________ - ! if the upper most isolated layer is fixed all layers below should be fixed as well - ! --> exit loop do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 - exit nzloop - - end if - end do nzloop ! --> do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 - end if ! --> if (min_nlvl < max_ulvl) then - end do ! --> do node=1, nod2D - - !_______________________________________________________________________ - ! check if cavity geometry did converge - if (exit_flag2 == 0) then - print *, achar(27)//'[33m' //'____________________________________________________________'//achar(27)//'[0m' - print *, ' -['//achar(27)//'[33m'//'WARN'//achar(27)//'[0m'//']->: Cavity geom. not converged yet, do further outer iteration' - write(*,"(A, I3, A, I3)") ' iter step ', count_iter2,'/', max_iter2 - write(*,*) - end if - - !_______________________________________________________________________ - end do ! --> do while((exit_flag2==0) .and. (count_iter2: Cavity geometry constrains did not converge !!! *\(>︿<)/*'//achar(27)//'[0m' - write(*,*) - stop - else - write(*,*) - print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' - print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge !!! *\(^o^)/*' - end if - - !___________________________________________________________________________ - ! write out cavity mesh files for vertice and elemental position of - ! vertical cavity-ocean boundary - ! write out elemental cavity-ocean boundary level - file_name=trim(meshpath)//'cavity_elvls.out' - open(20, file=file_name) - do elem=1,elem2D - write(20,*) ulevels(elem) - enddo - close(20) - - ! write out vertice cavity-ocean boundary level + yes/no cavity flag - file_name=trim(meshpath)//'cavity_nlvls.out' - open(20, file=file_name) - do node=1,nod2D - write(20,*) ulevels_nod2D(node) - enddo - close(20) - -end subroutine find_levels_cavity - - -!=================================================================== - -subroutine edge_center(n1, n2, x, y, mesh) -USE MOD_MESH -USE g_CONFIG -! -! Returns coordinates of edge center in x and y -! -implicit none -integer :: n1, n2 ! nodes of the edge -real(kind=WP) :: x, y, a(2), b(2) -type(t_mesh), intent(inout), target :: mesh -#include "associate_mesh_ini.h" - -a=coord_nod2D(:,n1) -b=coord_nod2D(:,n2) -if(a(1)-b(1)>cyclic_length/2.0) a(1)=a(1)-cyclic_length -if(a(1)-b(1)<-cyclic_length/2.0) b(1)=b(1)-cyclic_length -x=0.5_WP*(a(1)+b(1)) -y=0.5_WP*(a(2)+b(2)) -end subroutine edge_center -!==================================================================== -subroutine elem_center(elem, x, y, mesh) -! -! Returns coordinates of elem center in x and y -! -USE MOD_MESH -USE g_CONFIG -implicit none -integer, intent(in) :: elem -integer :: elnodes(3), k -real(kind=WP), intent(out) :: x, y -real(kind=WP) :: ax(3), amin -type(t_mesh), intent(in), target :: mesh -#include "associate_mesh_ini.h" - - elnodes=elem2D_nodes(1:3,elem) - ax=coord_nod2D(1, elnodes) - amin=minval(ax) - DO k=1,3 - if(ax(k)-amin>cyclic_length/2.0) ax(k)=ax(k)-cyclic_length - END DO - x=sum(ax)/3.0_WP - y=sum(coord_nod2D(2,elnodes))/3.0_WP - -end subroutine elem_center -!======================================================================= -SUBROUTINE find_elem_neighbors_ini(mesh) -! For each element three its element neighbors are found -USE MOD_MESH -implicit none -integer :: elem, eledges(3), elem1, j, n, elnodes(3) -type(t_mesh), intent(inout), target :: mesh -#include "associate_mesh_ini.h" - -allocate(mesh%elem_neighbors(4,elem2D)) -elem_neighbors => mesh%elem_neighbors !required after the allocation, otherwise the pointer remains undefined -elem_neighbors=0 -DO elem=1,elem2D - - eledges=elem_edges(1:3,elem) - DO j=1,3 - elem1=edge_tri(1,eledges(j)) - if(elem1==elem) elem1=edge_tri(2,eledges(j)) - elem_neighbors(j,elem)=elem1 - END DO - -END DO - ! Among elem_neighbors there can be negative numbers. These correspond to - ! boundary elements for which neighbours are absent. However, an element - ! should have at least two valid neighbors - - ! Test that there are at least two neighbors at the surface: - -DO elem=1,elem2D - elem1=0 - DO j=1,3 - if(elem_neighbors(j,elem)>0) elem1=elem1+1 - END DO - if (elem1<2) then - write(*,*) 'find_elem_neighbors_ini:Insufficient number of neighbors ',elem - write(*,*) 'find_elem_neighbors_ini:Elem neighbors ',elem_neighbors(:,elem) - write(*,*) '____________________________________________________________________' - write(*,*) ' ERROR: The mesh you want to partitioning contains triangles that' - write(*,*) ' have just one neighbor, this was OK for FESOM1.4 but not' - write(*,*) ' for FESOM2.0. ' - write(*,*) ' ' - write(*,*) ' ######################################### ' - write(*,*) ' ################### o ################### ' - write(*,*) ' ################# ./|\. ################# ' - write(*,*) ' Land ### ./|||||\. ### Land ' - write(*,*) ' ############## /|||||||||\ ############## ' - write(*,*) ' --o-----------o-----------o-----------o-- ' - write(*,*) ' ./ \. ./ \. ./ \. ./ \. ' - write(*,*) ' \. ./ \. ./ \. ./ ' - write(*,*) ' \ / \ / \ / ' - write(*,*) ' ------o-----------o-----------o------- ' - write(*,*) ' ./ \. ./ \. ./ \. ' - write(*,*) ' ' - write(*,*) ' Take a programm of your choice (Python, Matlab ...) and ' - write(*,*) ' eliminate these triangles and the corresponding ' - write(*,*) ' unconnected vertice and try to re-partitioning again ' - write(*,*) '____________________________________________________________________' - STOP - end if -END DO - - ! The rotation sense: corresponds to edges, and edges correspond - ! to nodes - - ! ============= - ! To facilitate computations the neibourhood - ! information is assembled - ! ============= - allocate(mesh%nod_in_elem2D_num(nod2D)) - nod_in_elem2D_num => mesh%nod_in_elem2D_num !required after the allocation, otherwise the pointer remains undefined - nod_in_elem2D_num=0 - do n=1,elem2D - elnodes=elem2D_nodes(1:3,n) - nod_in_elem2D_num(elnodes)=nod_in_elem2D_num(elnodes)+1 - end do - allocate(mesh%nod_in_elem2D(maxval(nod_in_elem2D_num),nod2D)) - nod_in_elem2D => mesh%nod_in_elem2D - nod_in_elem2D=0 - - nod_in_elem2D_num=0 - do n=1,elem2D - elnodes=elem2D_nodes(1:3,n) - nod_in_elem2D_num(elnodes)=nod_in_elem2D_num(elnodes)+1 - do j=1, 3 - nod_in_elem2D(nod_in_elem2D_num(elnodes(j)),elnodes(j))=n - end do - end do -END SUBROUTINE find_elem_neighbors_ini -!=================================================================== -! Stiffness matrix for the elevation -subroutine stiff_mat_ini(mesh) - use MOD_MESH - - ! - implicit none - integer :: i, j, n, q, el, elem_nodes_max, nod(4) - integer, allocatable :: num_ne(:), ne(:,:) - ! - type(t_mesh), intent(inout), target :: mesh -#include "associate_mesh_ini.h" - - ssh_stiff%dim = nod2D - allocate(mesh%ssh_stiff%rowptr(nod2D+1)) - ssh_stiff => mesh%ssh_stiff !required after the allocation, otherwise the pointer remains undefined - - allocate(num_ne(nod2D), ne(MAX_ADJACENT,nod2D)) - num_ne(:) = 0 - ne(:,:) = -1 - - ! Check the maximum number of nodes in an element (FESOM triangular meshes = 3, Hybrid meshes = 4) - elem_nodes_max = size(elem2D_nodes, 1) - - ! Fill node adjacency info - ! all nodes in an element are adjacent in the sense of being halo nodes - ! (also the opposite nodes of a quad: no edge, but the indirect connection - ! should be taken into account by metis domain decomposition) - do el=1,elem2D - nod(1:elem_nodes_max) = elem2D_nodes(1:elem_nodes_max,el) ! Fortran-numbering - q = elem_nodes_max - if (nod(1) == nod(elem_nodes_max)) q = q-1 ! triangle - - do i=2,q - do j=1,i-1 - if (all(ne(:,nod(i)) /= nod(j))) then - num_ne(nod(i)) = num_ne(nod(i)) + 1 - num_ne(nod(j)) = num_ne(nod(j)) + 1 - - if (max(num_ne(nod(i)), num_ne(nod(j))) > MAX_ADJACENT ) then - print *,'Recompile with larger value for MAX_ADJACENT.' - stop - else - ne(num_ne(nod(i)), nod(i)) = nod(j) - ne(num_ne(nod(j)), nod(j)) = nod(i) - endif - endif - end do - end do - end do - -! copy adjacency matrix to CSR-format - ssh_stiff%rowptr(1) = 1 - do n=1,nod2D - ssh_stiff%rowptr(n+1) = ssh_stiff%rowptr(n)+num_ne(n) - end do - - allocate(ssh_stiff%colind(ssh_stiff%rowptr(nod2D+1)-1)) - ssh_stiff => mesh%ssh_stiff !required after the allocation, otherwise the pointer remains undefined - - !required after the allocation, otherwise the pointer remains undefined - do n=1,nod2D - ssh_stiff%colind(ssh_stiff%rowptr(n):ssh_stiff%rowptr(n+1)-1) = ne(1:num_ne(n),n) - end do - - deallocate(num_ne, ne) - -end subroutine stiff_mat_ini - -!=================================================================== -! Setup of communication arrays -subroutine communication_ini(partit, mesh) - use MOD_MESH - USE g_CONFIG - USE MOD_PARTIT - use omp_lib - implicit none - - integer :: n - character*10 :: npes_string - character(MAX_PATH) :: dist_mesh_dir - LOGICAL :: L_EXISTS - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit -#include "associate_part_def.h" -#include "associate_mesh_ini.h" -#include "associate_part_ass.h" !only my - - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: compute communication arrays '//achar(27)//'[0m' - - ! Create the distributed mesh subdirectory - write(npes_string,"(I10)") npes - dist_mesh_dir=trim(meshpath)//'dist_'//trim(ADJUSTL(npes_string))//'/' - INQUIRE(file=trim(dist_mesh_dir), EXIST=L_EXISTS) - if (.not. L_EXISTS) call system('mkdir '//trim(dist_mesh_dir)) - -#ifdef OMP_MAX_THREADS -!$OMP PARALLEL NUM_THREADS(OMP_MAX_THREADS) - if (OMP_GET_THREAD_NUM() == 0) then - write(*,*) 'Setting up communication arrays using ', OMP_GET_NUM_THREADS(), ' threads' - endif -#else -!$OMP PARALLEL NUM_THREADS(1) - write(*,*) 'Setting up communication arrays using 1 thread (serially)' -#endif - -!$OMP DO - do n = 0, npes-1 - mype = n ! mype is threadprivate and must not be iterator - call communication_nodn(partit, mesh) - call communication_elemn(partit, mesh) - call save_dist_mesh(partit, mesh) ! Write out communication file com_infoxxxxx.out - end do -!$OMP END DO -!$OMP END PARALLEL - - deallocate(mesh%elem_neighbors) - deallocate(mesh%elem_edges) - deallocate(partit%part) - write(*,*) 'Communication arrays have been set up' -end subroutine communication_ini -!================================================================= -subroutine set_par_support_ini(partit, mesh) - use iso_c_binding, only: idx_t=>C_INT32_T - use MOD_MESH - use MOD_PARTIT - use g_config - implicit none -interface - subroutine check_partitioning(partit, mesh) - use MOD_MESH - use MOD_PARTIT - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine check_partitioning -end interface - - integer :: n, j, k, nini, nend, ierr - integer(idx_t) :: np(10) - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - interface - subroutine do_partit(n,ptr,adj,wgt,np,part) bind(C) - use iso_c_binding, only: idx_t=>C_INT32_T - integer(idx_t), intent(in) :: n, ptr(*), adj(*), wgt(*), np(*) - integer(idx_t), intent(out) :: part(*) - end subroutine do_partit - end interface - -#include "associate_part_def.h" -#include "associate_mesh_ini.h" -#include "associate_part_ass.h" - - if (mype==0) then - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: compute partitioning '//achar(27)//'[0m' - end if - - ! Construct partitioning vector - if (n_levels<1 .OR. n_levels>10) then - print *,'Number of hierarchic partition levels is out of range [1-10]! Aborting...' - stop - end if - - np(:) = n_part(:) ! Number of partitions on each hierarchy level - if (n_part(1) == 0) then ! Backward compatibility case: Take the number of - np(1) = npes ! partitions from the number of MPI processes - n_levels = 1 - end if - if (n_levels < 10) then ! 0 is an indicator of the last hierarchy level - np(n_levels+1) = 0 - end if - - allocate(partit%part(nod2D)) - part=>partit%part - part=0 - - npes = PRODUCT(np(1:n_levels)) - if(npes<2) then - print *,'Total number of parallel partitions is less than one! Aborting...' - stop - end if - - write(*,*) 'Calling partit for npes=', np - call do_partit(ssh_stiff%dim, ssh_stiff%rowptr, ssh_stiff%colind, & - nlevels_nod2D, np, part) - -write(*,*) 'DONE' -write(*,*) size(partit%part) - call check_partitioning(partit, mesh) - - write(*,*) 'Partitioning is done.' - -! The stiffness matrix is no longer needed. - deallocate(mesh%ssh_stiff%rowptr) - deallocate(mesh%ssh_stiff%colind) - - !NR No longer needed - last use was as weight for partitioning - deallocate(mesh%nlevels_nod2D) -end subroutine set_par_support_ini -!======================================================================= -subroutine check_partitioning(partit, mesh) - - ! In general, METIS 5 has several advantages compared to METIS 4, e.g., - ! * neighbouring tasks get neighbouring partitions (important for multicore computers!) - ! * lower maximum of weights per partition (better load balancing) - ! * lower memory demand - ! - ! BUT: there might be outliers, single nodes connected to their partition by - ! only one edge or even completely isolated. This spoils everything :-( - ! - ! This routine checks for isolated nodes and moves them to an adjacent partition, - ! trying not to spoil the load balance. - - use MOD_MESH - use MOD_PARTIT - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(inout), target :: mesh - integer :: i, j, k, n, n_iso, n_iter, is, ie, kmax, np, cnt - integer :: nod_per_partition(2,0:partit%npes-1) - integer :: max_nod_per_part(2), min_nod_per_part(2) - integer :: average_nod_per_part(2), node_neighb_part(100) - logical :: already_counted, found_part - - integer :: max_adjacent_nodes - integer, allocatable :: ne_part(:), ne_part_num(:), ne_part_load(:,:) -#include "associate_part_def.h" -#include "associate_mesh_ini.h" -#include "associate_part_ass.h" !just for partit%part - - if (mype==0) then - print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;1m' //' -->: check partitioning '//achar(27)//'[0m' - end if - - ! Check load balancing - do i=0, npes-1 - nod_per_partition(1,i) = count(part(:) == i) - nod_per_partition(2,i) = sum(nlevels_nod2D, part(:) == i) - enddo - - min_nod_per_part(1) = minval( nod_per_partition(1,:)) - min_nod_per_part(2) = minval( nod_per_partition(2,:)) - - max_nod_per_part(1) = maxval( nod_per_partition(1,:)) - max_nod_per_part(2) = maxval( nod_per_partition(2,:)) - - average_nod_per_part(1) = nod2D / npes - average_nod_per_part(2) = sum(nlevels_nod2D(:)) / npes - - ! Now check for isolated nodes (connect by one or even no edge to other - ! nodes of its partition) and repair, if possible - - max_adjacent_nodes = maxval(ssh_stiff%rowptr(2:nod2D+1) - ssh_stiff%rowptr(1:nod2D)) - allocate(ne_part(max_adjacent_nodes), ne_part_num(max_adjacent_nodes), & - ne_part_load(2,max_adjacent_nodes)) - - print *,' ' - print *,'Check for isolated nodes ========' - n_iso = 0 - checkloop: do n=1,nod2D - is = ssh_stiff%rowptr(n) - ie = ssh_stiff%rowptr(n+1) -1 - cnt = ie-is+1 - - node_neighb_part(1:cnt) = part(ssh_stiff%colind(is:ie)) - if (count(node_neighb_part(1:cnt) == part(n)) <= 1) then - - n_iso = n_iso+1 - print *,'Isolated node',n, 'in partition', part(n) - print *,'Neighbouring nodes are in partitions', node_neighb_part(1:cnt) - - ! count the adjacent nodes of the other PEs - - np=1 - ne_part(1) = node_neighb_part(1) - ne_part_num(1) = 1 - ne_part_load(1,1) = nod_per_partition(1,ne_part(1)) + 1 - ne_part_load(2,1) = nod_per_partition(2,ne_part(1)) + nlevels_nod2D(n) - - do i=1,cnt - if (node_neighb_part(i)==part(n)) cycle - already_counted = .false. - do k=1,np - if (node_neighb_part(i) == ne_part(k)) then - ne_part_num(k) = ne_part_num(k) + 1 - already_counted = .true. - exit - endif - enddo - if (.not. already_counted) then - np = np+1 - ne_part(np) = node_neighb_part(i) - ne_part_num(np) = 1 - ne_part_load(1,np) = nod_per_partition(1,ne_part(np)) + 1 - ne_part_load(2,np) = nod_per_partition(2,ne_part(np)) + nlevels_nod2D(n) - endif - enddo - - ! Now, check for two things: The load balance, and if - ! there is more than one node of that partition. - ! Otherwise, it would become isolated again. - - ! Best choice would be the partition with most adjacent nodes (edgecut!) - ! Choose, if it does not decrease the load balance. - ! (There might be two partitions with the same number of adjacent - ! nodes. Don't care about this here) - - kmax = maxloc(ne_part_num(1:np),1) - - if (ne_part_num(kmax) <= 1) then - ! No chance - this is probably a boundary - ! node that has only two neighbors. - cycle checkloop - endif - - if (ne_part_load(1,kmax) <= max_nod_per_part(1) .and. & - ne_part_load(2,kmax) <= max_nod_per_part(2) ) then - k = kmax - else - ! Don't make it too compicated. Reject partitions that have only one - ! adjacent node. Take the next not violating the load balance. - found_part = .false. - do k=1,np - if (ne_part_num(k)==1 .or. k==kmax) cycle - - if (ne_part_load(1,k) <= max_nod_per_part(1) .and. & - ne_part_load(2,k) <= max_nod_per_part(2) ) then - - found_part = .true. - exit - endif - enddo - - if (.not. found_part) then - ! Ok, don't think to much. Simply go for minimized edge cut. - k = kmax - endif - endif - - ! Adjust the load balancing - - nod_per_partition(1,ne_part(k)) = nod_per_partition(1,ne_part(k)) + 1 - nod_per_partition(2,ne_part(k)) = nod_per_partition(2,ne_part(k)) + nlevels_nod2D(n) - nod_per_partition(1,part(n)) = nod_per_partition(1,part(n)) - 1 - nod_per_partition(2,part(n)) = nod_per_partition(2,part(n)) - nlevels_nod2D(n) - - ! And, finally, move nod n to other partition - part(n) = ne_part(k) - print *,'Node',n,'is moved to part',part(n) - endif - enddo checkloop - - deallocate(ne_part, ne_part_num, ne_part_load) - - print *,'=== LOAD BALANCING ===' - print *,'2D nodes: min, aver, max per part',min_nod_per_part(1), & - average_nod_per_part(1),max_nod_per_part(1) - - write(*,"('2D nodes: percent min, aver, max ',f8.3,'%, 100%, ',f8.3,'%')") & - 100.*real(min_nod_per_part(1)) / real(average_nod_per_part(1)), & - 100.*real(max_nod_per_part(1)) / real(average_nod_per_part(1)) - - print *,'3D nodes: Min, aver, max per part',min_nod_per_part(2), & - average_nod_per_part(2),max_nod_per_part(2) - write(*,"('3D nodes: percent min, aver, max ',f8.3,'%, 100%, ',f8.3,'%')") & - 100.*real(min_nod_per_part(2)) / real(average_nod_per_part(2)), & - 100.*real(max_nod_per_part(2)) / real(average_nod_per_part(2)) - -end subroutine check_partitioning - diff --git a/src/gen_comm.F90 b/src/gen_comm.F90 index 26d318008..c47c0a3a3 100755 --- a/src/gen_comm.F90 +++ b/src/gen_comm.F90 @@ -1,14 +1,26 @@ -! Cell-vertex finite-volume version -! Contains: Routines that support parallelization -! set_par_support_ini run in the initialization phase. -! The communication rules are saved. -! set_par_support in the main phase just allocates memory for buffer -! arrays, the rest is read together with mesh from saved files. +!> @file mod_gen_comm.F90 +!! @brief Communication routines for FESOM mesh partitioner +!! @details Cell-vertex finite-volume version +!! Contains: Routines that support parallelization +!! set_par_support_ini run in the initialization phase. +!! The communication rules are saved. +!! set_par_support in the main phase just allocates memory for buffer +!! arrays, the rest is read together with mesh from saved files. + +module mod_gen_comm + + use MOD_MESH + use MOD_PARTIT + use MOD_PARSUP + + implicit none + + public :: communication_nodn, communication_elemn, mymesh + +contains + !======================================================================= subroutine communication_nodn(partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP implicit none type(t_mesh), intent(in), target :: mesh type(t_partit), intent(inout), target :: partit @@ -19,7 +31,8 @@ subroutine communication_nodn(partit, mesh) integer :: IERR #include "associate_part_def.h" #include "associate_mesh_ini.h" -#include "associate_part_ass.h" !part only +#include "associate_part_ass.h" + !part only ! Assume we have 2D partitioning vector in part. Find communication rules ! Reduce allocation: find all neighboring PE nd_count = count(part(1:nod2d) == mype) @@ -221,9 +234,6 @@ end subroutine communication_nodn !========================================================================== subroutine communication_elemn(partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP implicit none type(t_mesh), intent(in), target :: mesh @@ -236,7 +246,8 @@ subroutine communication_elemn(partit, mesh) integer :: IERR #include "associate_part_def.h" #include "associate_mesh_ini.h" -#include "associate_part_ass.h" !part only +#include "associate_part_ass.h" + !part only ! Assume we have 2D partitioning vector in part. Find communication ! rules. An elem is external to element n if neither of its nodes ! belongs to PE, but it is among the neighbors. Element n belongs to PE if @@ -527,9 +538,6 @@ subroutine communication_elemn(partit, mesh) end subroutine communication_elemn !========================================================================== subroutine mymesh(partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP implicit none type(t_mesh), intent(in), target :: mesh @@ -656,4 +664,5 @@ subroutine mymesh(partit, mesh) ! myList_edge2D(1:myDim_edge2D) contains owned edges + ! shared edges which mype updates end subroutine mymesh -!================================================================= + +end module mod_gen_comm diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 3e7e85d4a..86f912e10 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -1,77 +1,14 @@ -module force_flux_consv_interface - interface - subroutine force_flux_consv(field2d, mask, n, h, do_stats, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in), target :: mesh - type(t_partit), intent(inout), target :: partit - real(kind=WP), intent (inout) :: field2d(partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent (in) :: mask(partit%myDim_nod2D+partit%eDim_nod2D) - integer, intent (in) :: n, h - logical, intent (in) :: do_stats - end subroutine force_flux_consv - end interface -end module force_flux_consv_interface -module compute_residual_interface - interface - subroutine compute_residual(field2d, mask, n, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in), target :: mesh - type(t_partit), intent(inout), target :: partit - real(kind=WP), intent (in) :: field2d(partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent (in) :: mask(partit%myDim_nod2D+partit%eDim_nod2D) - integer, intent (in) :: n - end subroutine compute_residual - end interface -end module compute_residual_interface -module integrate_2D_interface - interface - subroutine integrate_2D(flux_global, flux_local, eff_vol, field2d, mask, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in), target :: mesh - type(t_partit), intent(in), target :: partit - real(kind=WP), intent (out) :: flux_global(2), flux_local(2) - real(kind=WP), intent (out) :: eff_vol(2) - real(kind=WP), intent (in) :: field2d(partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent (in) :: mask(partit%myDim_nod2D +partit%eDim_nod2D) - end subroutine integrate_2D - end interface -end module integrate_2D_interface - -module update_atm_forcing_interface - interface - subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) - USE MOD_TRACER - USE MOD_ICE - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_MESH - USE MOD_DYN - integer, intent(in) :: istep - type(t_ice), intent(inout), target :: ice - type(t_tracer), intent(in), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - type(t_dyn) , intent(in), target :: dynamics - end subroutine update_atm_forcing - end interface -end module update_atm_forcing_interface - -module net_rec_from_atm_interface - interface - subroutine net_rec_from_atm(action, partit) - USE MOD_PARTIT - USE MOD_PARSUP - logical, intent(in) :: action - type(t_partit), intent(inout), target :: partit - end subroutine net_rec_from_atm - end interface -end module net_rec_from_atm_interface +module forcing_coupling_interfaces + use mod_mesh + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_TRACER + USE MOD_ICE + USE MOD_DYN + + implicit none + +contains ! Routines for updating ocean surface forcing fields !------------------------------------------------------------------------- subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) @@ -89,7 +26,6 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) use g_config use g_comm_auto use g_rotate_grid - use net_rec_from_atm_interface use g_sbf, only: sbc_do use g_sbf, only: atmdata, i_totfl, i_xwind, i_ywind, i_xstre, i_ystre, i_humi, i_qsr, i_qlw, i_tair, i_prec, i_mslp, i_cloud, i_snow, & l_xwind, l_ywind, l_xstre, l_ystre, l_humi, l_qsr, l_qlw, l_tair, l_prec, l_mslp, l_cloud, l_snow @@ -97,7 +33,6 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) use cpl_driver #endif use gen_bulk - use force_flux_consv_interface implicit none integer, intent(in) :: istep @@ -623,8 +558,6 @@ SUBROUTINE force_flux_consv(field2d, mask, n, h, do_stats, partit, mesh) USE MOD_PARSUP use cpl_driver, only : nrecv, cpl_recv, a2o_fcorr_stat use o_PARAM, only : mstep, WP - use compute_residual_interface - use integrate_2D_interface IMPLICIT NONE type(t_mesh), intent(in), target :: mesh type(t_partit), intent(inout), target :: partit @@ -743,7 +676,6 @@ SUBROUTINE compute_residual(field2d, mask, n, partit, mesh) use MOD_MESH USE MOD_PARTIT USE MOD_PARSUP - use integrate_2D_interface IMPLICIT NONE type(t_mesh), intent(in), target :: mesh @@ -875,3 +807,5 @@ SUBROUTINE net_rec_from_atm(action, partit) end if END SUBROUTINE net_rec_from_atm #endif + +end module forcing_coupling_interfaces diff --git a/src/gen_ic3d.F90 b/src/gen_ic3d.F90 index 4af214cd6..3d3a5d243 100644 --- a/src/gen_ic3d.F90 +++ b/src/gen_ic3d.F90 @@ -20,6 +20,7 @@ MODULE g_ic3d USE g_comm_auto USE g_support USE g_config, only: dummy, ClimateDataPath, use_cavity + use oce_ale_pressure_bv_module, only: insitu2pot IMPLICIT NONE @@ -496,7 +497,6 @@ SUBROUTINE do_ic3d(tracers, partit, mesh) !! !! ** Purpose : read 3D initial conditions for tracers from netcdf and interpolate on model grid !!---------------------------------------------------------------------- - USE insitu2pot_interface IMPLICIT NONE type(t_mesh), intent(in), target :: mesh type(t_partit), intent(inout), target :: partit diff --git a/src/gen_modules_cvmix_kpp.F90 b/src/gen_modules_cvmix_kpp.F90 index b72706d6d..a70a59e1f 100644 --- a/src/gen_modules_cvmix_kpp.F90 +++ b/src/gen_modules_cvmix_kpp.F90 @@ -33,6 +33,7 @@ module g_cvmix_kpp use g_forcing_arrays use g_support use o_mixing_KPP_mod + use oce_ale_pressure_bv_module, only: densityJM_components implicit none !___Parameter for the init of KPP___________________________________________ @@ -431,8 +432,8 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! oce_ale_pressure_bf.F90 --> subroutine pressure_bv ! --> dbsfc(nz,node) !!PS call densityJM_components(temp(1,node), salt(1,node), sfc_bulk_0, sfc_bulk_pz, sfc_bulk_pz2, sfc_rhopot, mesh) - call densityJM_components(temp(nun,node), salt(nun,node), sfc_bulk_0, sfc_bulk_pz, sfc_bulk_pz2, sfc_rhopot, mesh) - call densityJM_components(temp(nz, node), salt(nz, node), bulk_0, bulk_pz, bulk_pz2, rhopot, mesh) + call densityJM_components(temp(nun,node), salt(nun,node), sfc_bulk_0, sfc_bulk_pz, sfc_bulk_pz2, sfc_rhopot) + call densityJM_components(temp(nz, node), salt(nz, node), bulk_0, bulk_pz, bulk_pz2, rhopot) rho_nz = bulk_0 + Z_3d_n(nz,node)*(bulk_pz + Z_3d_n(nz,node)*bulk_pz2) rho_nz = rho_nz*rhopot/(rho_nz+0.1_WP*Z_3d_n(nz,node))-density_0 rho_sfc = sfc_bulk_0 + Z_3d_n(nz,node)*(sfc_bulk_pz + Z_3d_n(nz,node)*sfc_bulk_pz2) @@ -489,8 +490,8 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! --> bring density of surface point adiabatically to the same ! depth level as the deep point --> than calculate bouyancy ! difference - call densityJM_components(sfc_temp, sfc_salt, sfc_bulk_0, sfc_bulk_pz, sfc_bulk_pz2, sfc_rhopot, mesh) - call densityJM_components(temp(nz,node), salt(nz,node), bulk_0, bulk_pz, bulk_pz2, rhopot, mesh) + call densityJM_components(sfc_temp, sfc_salt, sfc_bulk_0, sfc_bulk_pz, sfc_bulk_pz2, sfc_rhopot) + call densityJM_components(temp(nz,node), salt(nz,node), bulk_0, bulk_pz, bulk_pz2, rhopot) rho_nz = bulk_0 + Z_3d_n(nz,node)*(bulk_pz + Z_3d_n(nz,node)*bulk_pz2) rho_nz = rho_nz*rhopot/(rho_nz+0.1_WP*Z_3d_n(nz,node))-density_0 rho_sfc = sfc_bulk_0 + Z_3d_n(nz,node)*(sfc_bulk_pz + Z_3d_n(nz,node)*sfc_bulk_pz2) diff --git a/src/mod_mesh_utils.F90 b/src/mod_mesh_utils.F90 new file mode 100644 index 000000000..0577c8dd5 --- /dev/null +++ b/src/mod_mesh_utils.F90 @@ -0,0 +1,1750 @@ +!> @file mod_mesh_utils.F90 +!! @brief Mesh utility subroutines for FESOM mesh partitioner +!! @author Extracted from fvom_init.F90 +!! @details Contains utility subroutines used by the FESOM mesh partitioner tool. +!! These subroutines were originally part of fvom_init.F90 but have been separated +!! into this module for better organization and to avoid compilation issues. + +module mod_mesh_utils + + use MOD_MESH + use o_PARAM + use g_CONFIG + + implicit none + + public :: read_mesh_ini, read_mesh_cavity, test_tri_ini, find_edges_ini + public :: find_elem_neighbors_ini, find_levels, find_levels_cavity + public :: edge_center, elem_center, stiff_mat_ini, communication_ini + public :: set_par_support_ini, check_partitioning + +contains + +!============================================================================= +!> @brief +!> Reads mesh files +subroutine read_mesh_ini(mesh) +use g_rotate_grid +! +IMPLICIT NONE +! +type(t_mesh), intent(inout), target :: mesh +INTEGER :: nq +INTEGER :: n1,n2,n3 +INTEGER :: n, nz, exit_flag +REAL(kind=WP) :: x1, x2, gx1, gx2 +INTEGER :: tag +INTEGER, allocatable :: elem_data(:) +INTEGER :: i_error +#include "associate_mesh_ini.h" +! =================== +! Surface mesh +! =================== + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: read elem2d.out & nod2d.out '//achar(27)//'[0m' + + open (20,file=trim(meshpath)//'nod2d.out', status='old') + open (21,file=trim(meshpath)//'elem2d.out', status='old') + READ(20,*) mesh%nod2D + ALLOCATE(mesh%coord_nod2D(2,mesh%nod2D)) + coord_nod2D => mesh%coord_nod2D !required after the allocation, otherwise the pointer remains undefined + + do n=1, mesh%nod2D + read(20,*) nq, x1, x2, tag + x1=x1*rad + x2=x2*rad + if (force_rotation) then + gx1=x1 + gx2=x2 + call g2r(gx1, gx2, x1, x2) + end if + mesh%coord_nod2D(1,n)=x1 + mesh%coord_nod2D(2,n)=x2 + end do + CLOSE(20) + READ(21,*) mesh%elem2D + ALLOCATE(mesh%elem2D_nodes(4,mesh%elem2D)) + elem2D_nodes => mesh%elem2D_nodes !required after the allocation, otherwise the pointer remains undefined + ALLOCATE(elem_data(4*mesh%elem2D)) + elem_data(:)=-1 + + ! meshes with quads have 4 columns, but TsunAWI grids may be + ! purely triangular, with 3 columns each. Test, how many + ! columns there are! + read(21,*,iostat=i_error) elem_data(1:4*mesh%elem2D) + if (i_error == 0) then ! There is a fourth column => quad or mixed mesh (not working yet!) + mesh%elem2D_nodes = reshape(elem_data, shape(mesh%elem2D_nodes)) + else ! No fourth column => triangles only + mesh%elem2D_nodes(1:3,:) = reshape(elem_data, shape(mesh%elem2D_nodes(1:3,:))) + mesh%elem2D_nodes(4,:) = mesh%elem2D_nodes(1,:) + end if + + deallocate(elem_data) + CLOSE(21) + +write(*,*) '=========================' +write(*,*) 'Mesh is read' +write(*,*) '=========================' +END SUBROUTINE read_mesh_ini + +!============================================================================= +!> @brief +!> Reads mesh cavity files +subroutine read_mesh_cavity(mesh) + implicit none + + type(t_mesh), intent(inout), target :: mesh + integer :: node, auxi + character(len=MAX_PATH) :: fname + logical :: file_exist=.False. +#include "associate_mesh_ini.h" + + !___________________________________________________________________________ + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: read cavity depth '//achar(27)//'[0m' + + !___________________________________________________________________________ + ! read depth of cavity-ocean boundary + if (use_cavityonelem) then + fname = trim(meshpath)//'cavity_depth@elem.out' + else + fname = trim(meshpath)//'cavity_depth@node.out' + end if + file_exist=.False. + inquire(file=trim(fname),exist=file_exist) + if (file_exist) then + open (21,file=fname, status='old') + if (use_cavityonelem) then + allocate(mesh%cavity_depth(mesh%elem2d)) + else + allocate(mesh%cavity_depth(mesh%nod2D)) + end if + cavity_depth => mesh%cavity_depth + else + write(*,*) '____________________________________________________________________' + write(*,*) ' ERROR: could not find cavity file:', fname + write(*,*) ' --> stop partitioning here !' + write(*,*) '____________________________________________________________________' + stop + end if + + !___________________________________________________________________________ + auxi=mesh%nod2D + if (use_cavityonelem) auxi=mesh%elem2d + do node=1, auxi + read(21,*) mesh%cavity_depth(node) + end do + + !___________________________________________________________________________ + close(21) + +end subroutine read_mesh_cavity + +!======================================================================= +!> @brief +!> Check the order of nodes in triangles; correct it if necessary to make +!! it same sense (clockwise) +SUBROUTINE test_tri_ini(mesh) +IMPLICIT NONE + +real(kind=WP) :: a(2), b(2), c(2), r +integer :: n, nx, elnodes(3) +type(t_mesh), intent(inout), target :: mesh +#include "associate_mesh_ini.h" + + DO n=1, elem2D + elnodes=elem2D_nodes(1:3,n) + + a=coord_nod2D(:,elnodes(1)) + b=coord_nod2D(:,elnodes(2))-a + c=coord_nod2D(:,elnodes(3))-a + + if(b(1)>cyclic_length/2.) b(1)=b(1)-cyclic_length + if(b(1)<-cyclic_length/2.) b(1)=b(1)+cyclic_length + if(c(1)>cyclic_length/2.) c(1)=c(1)-cyclic_length + if(c(1)<-cyclic_length/2.) c(1)=c(1)+cyclic_length + + + r=b(1)*c(2)-b(2)*c(1) + if (r>0) then + ! Vector b is to right of c + ! Exchange second and third nodes: + + nx=elnodes(2) + elnodes(2)=elnodes(3) + elnodes(3)=nx + elem2D_nodes(1:3,n)=elnodes + end if + END DO +END SUBROUTINE test_tri_ini + +!========================================================================= +!> @brief +!> Finds edges. Creates 3 files: edgenum.out, edges.out, edge_tri.out +SUBROUTINE find_edges_ini(mesh) +use g_rotate_grid +IMPLICIT NONE + +integer, allocatable :: aux1(:), ne_num(:), ne_pos(:,:), nn_num(:), nn_pos(:,:) +integer :: counter, counter_in, n, k, q +integer :: elem, elem1, elems(2), q1, q2 +integer :: elnodes(4), ed(2), flag, eledges(4) +integer :: temp(100), node +real(kind=WP) :: xc(2), xe(2), ax(3), amin +type(t_mesh), intent(inout), target :: mesh +#include "associate_mesh_ini.h" +! ==================== +! (a) find edges. To make the procedure fast +! one needs neighbourhood arrays +! ==================== +print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' +print *, achar(27)//'[7;1m' //' -->: compute edge connectivity '//achar(27)//'[0m' + +allocate(ne_num(nod2d), ne_pos(MAX_ADJACENT, nod2D), nn_num(nod2D)) +ne_num=0 +DO n=1,elem2D + elnodes=elem2D_nodes(:,n) + q1=4 + if(elnodes(1)==elnodes(4)) q1=3 + DO q=1,q1 + ne_num(elnodes(q))=ne_num(elnodes(q))+1 + if (ne_num(elnodes(q)) > MAX_ADJACENT ) then + print *,'Recompile with larger value for MAX_ADJACENT.' + stop + else + ne_pos(ne_num(elnodes(q)),elnodes(q))=n + endif + END Do +END DO ! neighbor elements are found + +! count neighbour nodes +! In quads we should count the nodes that are +! connected by edges! +allocate(aux1(nod2D)) +aux1=0 + +DO n=1, nod2D + counter=0 + DO k=1, ne_num(n) + elem=ne_pos(k,n) + elnodes=elem2D_nodes(:,elem) + if(elnodes(1)==elnodes(4)) then + DO q=1,3 + if(elnodes(q)==n) CYCLE + if(aux1(elnodes(q)).ne.1) then + counter=counter+1 + aux1(elnodes(q))=1 + temp(counter)=elnodes(q) + end if + END DO + else + ! Find the position of n in elnodes: + if(elnodes(1)==n .or. elnodes(3)==n) then + ed(1)=elnodes(2) + ed(2)=elnodes(4) + else + ed(1)=elnodes(1) + ed(2)=elnodes(3) + end if + DO q=1,2 + if(aux1(ed(q)).ne.1) then + counter=counter+1 + aux1(ed(q))=1 + temp(counter)=ed(q) + end if + END DO + end if + END DO + nn_num(n)=counter + aux1(temp(1:counter))=0 +END DO + +allocate(nn_pos(maxval(nn_num)+1,nod2D)) +nn_pos = -1 +aux1=0 + +DO n=1, nod2D + counter=0 + DO k=1, ne_num(n) + elem=ne_pos(k,n) + elnodes=elem2D_nodes(:,elem) + if(elnodes(1)==elnodes(4)) then + DO q=1,3 + if(elnodes(q)==n) CYCLE + if(aux1(elnodes(q)).ne.1) then + counter=counter+1 + aux1(elnodes(q))=1 + temp(counter)=elnodes(q) + end if + END DO + else + ! Find the position of n in elnodes: + if(elnodes(1)==n .or. elnodes(3)==n) then + ed(1)=elnodes(2) + ed(2)=elnodes(4) + else + ed(1)=elnodes(1) + ed(2)=elnodes(3) + end if + DO q=1,2 + if(aux1(ed(q)).ne.1) then + counter=counter+1 + aux1(ed(q))=1 + temp(counter)=ed(q) + end if + END DO + end if + END DO + nn_num(n)=counter+1 + aux1(temp(1:counter))=0 + nn_pos(1,n)=n + nn_pos(2:counter+1,n)=temp(1:counter) +END DO +deallocate(aux1) +! neighboring nodes are found. First in the list is the node itself + +! ==================== +! (b) Find edges and elements containing them. +! Write information to auxiliary file +! ==================== +! open(10, file='edges.out') + + ! Count edges: + ! ==================== + ! form edges with n by cycling over neighbor + ! nodes (if edges are not accounted yet). + ! New edges are added only if neighbor>n + ! ==================== + counter = 0 + DO n=1,nod2D + counter = counter + count(nn_pos(2:nn_num(n),n)>nn_pos(1,n)) + end do + edge2D=counter + + allocate(mesh%edges (2, edge2D)) + allocate(mesh%edge_tri(2, edge2D)) + edges => mesh%edges !required after the allocation, otherwise the pointer remains undefined + edge_tri => mesh%edge_tri !required after the allocation, otherwise the pointer remains undefined + counter_in=0 + + DO n=1,nod2D + DO q=2,nn_num(n) + node=nn_pos(q,n) + if(node0.0_WP) then + ! Vector drawn to the center of the first triangle is to the right + ! of the edge vector. Triangles have to be exchanged: + elem=edge_tri(1,n) + elem1=edge_tri(2,n) + if(elem1>0) then !TODO + edge_tri(1,n)=elem1 + edge_tri(2,n)=elem + else + if (elem<=0) write(*,*) '2 neighbouring elems are on the ground.' + elem=edges(2,n) ! change the order of nodes + edges(2,n)=edges(1,n) + edges(1,n)=elem + end if + end if + END DO + + ! ==================== + ! (e) We need an array inverse to edge_tri listing edges + ! of a given triangle + ! ==================== + allocate(mesh%elem_edges(4,elem2D)) + elem_edges => mesh%elem_edges !required after the allocation, otherwise the pointer remains undefined + allocate(aux1(elem2D)) + aux1=0 + DO n=1, edge2D + DO k=1,2 + q=edge_tri(k,n) ! triangle number + if (q>0) then + aux1(q)=aux1(q)+1 + elem_edges(aux1(q),q)=n + end if + END DO + END DO + deallocate(aux1) + ! We order the edges in this list so that they + ! are listed in the same rotation sense as nodes. + ! First is the edge formed by elnodes(1:2), and so on + DO elem=1,elem2D + elnodes=elem2D_nodes(:,elem) + q1=4 + if(elnodes(1)==elnodes(4)) q1=3 + eledges=elem_edges(:,elem) + DO q=1,q1-1 + DO k=1,q1 + if(((edges(1,eledges(k))==elnodes(q)).and. & + (edges(2,eledges(k))==elnodes(q+1))).or. & + ((edges(1,eledges(k))==elnodes(q+1)).and. & + (edges(2,eledges(k))==elnodes(q)))) then + elem_edges(q,elem)=eledges(k) + exit + end if + END DO + END DO + DO k=1,q1 + if(((edges(1,eledges(k))==elnodes(q1)).and. & + (edges(2,eledges(k))==elnodes(1))).or. & + ((edges(1,eledges(k))==elnodes(1)).and. & + (edges(2,eledges(k))==elnodes(q1)))) then + elem_edges(q1,elem)=eledges(k) + exit + end if + END DO + if(q1==3) elem_edges(4,elem)=elem_edges(1,elem) + END DO + + !> The edge and elem lists agree in the sense that edge 1 does not + !! contain node 1 and so on + open(11, file=trim(meshpath)//'edgenum.out') + write(11,*) edge2D + write(11,*) edge2D_in + close(11) + open(10, file=trim(meshpath)//'edges.out') + open(12, file=trim(meshpath)//'edge_tri.out') + do n=1,edge2D + write(10,*) edges(:,n) + write(12,*) edge_tri(:,n) + end do + close(10) + close(12) + deallocate(ne_num, ne_pos) +END SUBROUTINE find_edges_ini + +!======================================================================= +SUBROUTINE find_elem_neighbors_ini(mesh) +! For each element three its element neighbors are found +implicit none +integer :: elem, eledges(3), elem1, j, n, elnodes(3) +type(t_mesh), intent(inout), target :: mesh +#include "associate_mesh_ini.h" + +allocate(mesh%elem_neighbors(4,elem2D)) +elem_neighbors => mesh%elem_neighbors !required after the allocation, otherwise the pointer remains undefined +elem_neighbors=0 +DO elem=1,elem2D + + eledges=elem_edges(1:3,elem) + DO j=1,3 + elem1=edge_tri(1,eledges(j)) + if(elem1==elem) elem1=edge_tri(2,eledges(j)) + elem_neighbors(j,elem)=elem1 + END DO + +END DO + ! Among elem_neighbors there can be negative numbers. These correspond to + ! boundary elements for which neighbours are absent. However, an element + ! should have at least two valid neighbors + + ! Test that there are at least two neighbors at the surface: + +DO elem=1,elem2D + elem1=0 + DO j=1,3 + if(elem_neighbors(j,elem)>0) elem1=elem1+1 + END DO + if (elem1<2) then + write(*,*) 'find_elem_neighbors_ini:Insufficient number of neighbors ',elem + write(*,*) 'find_elem_neighbors_ini:Elem neighbors ',elem_neighbors(:,elem) + write(*,*) '____________________________________________________________________' + write(*,*) ' ERROR: The mesh you want to partitioning contains triangles that' + write(*,*) ' have just one neighbor, this was OK for FESOM1.4 but not' + write(*,*) ' for FESOM2.0. ' + write(*,*) ' ' + write(*,*) ' ######################################### ' + write(*,*) ' ################### o ################### ' + write(*,*) ' ################# ./|\. ################# ' + write(*,*) ' Land ### ./|||||\. ### Land ' + write(*,*) ' ############## /|||||||||\ ############## ' + write(*,*) ' --o-----------o-----------o-----------o-- ' + write(*,*) ' ./ \. ./ \. ./ \. ./ \. ' + write(*,*) ' \. ./ \. ./ \. ./ ' + write(*,*) ' \ / \ / \ / ' + write(*,*) ' ------o-----------o-----------o------- ' + write(*,*) ' ./ \. ./ \. ./ \. ' + write(*,*) ' ' + write(*,*) ' Take a programm of your choice (Python, Matlab ...) and ' + write(*,*) ' eliminate these triangles and the corresponding ' + write(*,*) ' unconnected vertice and try to re-partitioning again ' + write(*,*) '____________________________________________________________________' + STOP + end if +END DO + + ! The rotation sense: corresponds to edges, and edges correspond + ! to nodes + + ! ============= + ! To facilitate computations the neibourhood + ! information is assembled + ! ============= + allocate(mesh%nod_in_elem2D_num(nod2D)) + nod_in_elem2D_num => mesh%nod_in_elem2D_num !required after the allocation, otherwise the pointer remains undefined + nod_in_elem2D_num=0 + do n=1,elem2D + elnodes=elem2D_nodes(1:3,n) + nod_in_elem2D_num(elnodes)=nod_in_elem2D_num(elnodes)+1 + end do + allocate(mesh%nod_in_elem2D(maxval(nod_in_elem2D_num),nod2D)) + nod_in_elem2D => mesh%nod_in_elem2D + nod_in_elem2D=0 + + nod_in_elem2D_num=0 + do n=1,elem2D + elnodes=elem2D_nodes(1:3,n) + nod_in_elem2D_num(elnodes)=nod_in_elem2D_num(elnodes)+1 + do j=1, 3 + nod_in_elem2D(nod_in_elem2D_num(elnodes(j)),elnodes(j))=n + end do + end do +END SUBROUTINE find_elem_neighbors_ini + +!========================================================================= +!> @brief +!> Finds elemental and nodal levels. +!> Does some thresholding: if (depth>zbar(4)) x=zbar(4) +!> Fixes rough topography, by converting some oceans cells to ground cell(reflected by changing levels arrays) +!> Creates 2 files: elvls.out, nlvls.out +subroutine find_levels(mesh) + implicit none + INTEGER :: nodes(3), elems(3), eledges(3) + integer :: elem1, j, n, nneighb, q, node, i, nz, auxi + integer :: count_iter, count_neighb_open, exit_flag, fileID=111 + real(kind=WP) :: x, dmean + logical :: file_exist + integer :: max_iter=1000 + character(MAX_PATH) :: file_name + type(t_mesh), intent(inout), target :: mesh + +#include "associate_mesh_ini.h" + + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: read bottom depth '//achar(27)//'[0m' + + !___________________________________________________________________________ + ! allocate depth + if (use_depthonelem) then + allocate(mesh%depth(elem2D)) + else + allocate(mesh%depth(nod2D)) + end if + depth => mesh%depth !required after the allocation, otherwise the pointer remains undefined + + !______________________________________________________________________________ + ! read depth from aux3d.out + if (trim(use_depthfile)=='aux3d') then + ! check if aux3d.out file does exist + file_exist=.False. + file_name=trim(meshpath)//'aux3d.out' + inquire(file=trim(file_name),exist=file_exist) + !_______________________________________________________________________ + if (file_exist) then + write(*," (A, A)") ' read file:',trim(file_name) + !___________________________________________________________________ + ! load fesom2.0 aux3d.out file + open(fileID, file=file_name) + + ! read the number of levels + read(fileID,*) nl + allocate(mesh%zbar(nl)) ! their standard depths + + ! read full depth levels + zbar => mesh%zbar !required after the allocation, otherwise the pointer remains undefined + read(fileID,*) zbar + if(zbar(2)>0) zbar=-zbar ! zbar is negative + + ! compute mid depth levels + allocate(mesh%Z(nl-1)) + Z => mesh%Z !required after the allocation, otherwise the pointer remains undefined + Z=zbar(1:nl-1)+zbar(2:nl) ! mid-depths of cells + Z=0.5_WP*Z + else + write(*,*) '____________________________________________________________________' + write(*,*) ' ERROR: You want to use aux3d.out file to define your depth, but ' + write(*,*) ' the file seems not to exist' + write(*,*) ' --> check in namelist.config, the flag use_depthfile must be' + write(*,*) ' use_depthfile= "aux3d" or "depth@" ' + write(*,*) ' --> model stops here' + write(*,*) '____________________________________________________________________' + stop + end if + !___________________________________________________________________________ + ! read depth from depth@node.out or depth@elem.out + elseif (trim(use_depthfile)=='depth@') then + !_______________________________________________________________________ + ! load file depth_zlev.out --> contains number of model levels and full depth + ! levels + file_exist=.False. + file_name=trim(meshpath)//'depth_zlev.out' + inquire(file=trim(file_name),exist=file_exist) + if (file_exist) then + write(*," (A, A)") ' read file:',trim(file_name) + !___________________________________________________________________ + ! load fesom2.0 aux3d.out file + open(fileID, file=file_name) + + ! read the number of levels + read(fileID,*) nl + allocate(mesh%zbar(nl)) ! their standard depths + + ! read full depth levels + zbar => mesh%zbar !required after the allocation, otherwise the pointer remains undefined + read(fileID,*) zbar + if(zbar(2)>0) zbar=-zbar ! zbar is negative + + ! compute mid depth levels + allocate(mesh%Z(nl-1)) + Z => mesh%Z !required after the allocation, otherwise the pointer remains undefined + Z=zbar(1:nl-1)+zbar(2:nl) ! mid-depths of cells + Z=0.5_WP*Z + + close(fileID) + else + write(*,*) '____________________________________________________________________' + write(*,*) ' ERROR: You want to use depth@elem.out or depth@node.out file, therefore' + write(*,*) ' you also need the file depth_zlev.out which contains the model ' + write(*,*) ' number of layers and the depth of model levels. This file seems ' + write(*,*) ' not to exist' + write(*,*) ' --> check in namelist.config, the flag use_depthfile must be' + write(*,*) ' use_depthfile= "aux3d" or "depth@" and your meshfolder' + write(*,*) ' --> model stops here' + write(*,*) '____________________________________________________________________' + stop + endif + + !_______________________________________________________________________ + ! load file depth@elem.out or depth@node.out contains topography either at + ! nodes or elements + if (use_depthonelem) then + file_name=trim(meshpath)//'depth@elem.out' + else + file_name=trim(meshpath)//'depth@node.out' + end if + inquire(file=trim(file_name),exist=file_exist) + if (file_exist) then + write(*," (A, A)") ' read file:',trim(file_name) + open(fileID, file=file_name) + else + write(*,*) '____________________________________________________________________' + write(*,*) ' ERROR: You want to use depth@elem.out or depth@node.out file to ' + write(*,*) ' define your depth, but the file seems not to exist' + write(*,*) ' --> check in namelist.config, the flag use_depthfile must be' + write(*,*) ' use_depthfile= "aux3d" or "depth@" and your meshfolder ' + write(*,*) ' --> model stops here' + write(*,*) '____________________________________________________________________' + stop + end if + end if + + !___________________________________________________________________________ + ! read topography from file + auxi = nod2d + if (use_depthonelem) auxi = elem2d +! write(*,*) ' use_depthonelem = ',use_depthonelem +! write(*,*) ' auxi =',auxi + DO n = 1, auxi + read(fileID,*) x + if (x>0) x=-x + if (x>zbar(thers_zbar_lev)) x=zbar(thers_zbar_lev) !TODO KK thresholding for depth + depth(n)=x + END DO + close(fileID) + if(depth(2)>0) depth=-depth ! depth is negative + + !___________________________________________________________________________ + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute elem, vertice bottom depth index '//achar(27)//'[0m' + + allocate(mesh%nlevels(elem2D)) + nlevels => mesh%nlevels !required after the allocation, otherwise the pointer remains undefined + allocate(mesh%nlevels_nod2D(nod2D)) + nlevels_nod2D => mesh%nlevels_nod2D !required after the allocation, otherwise the pointer remains undefined + + !___________________________________________________________________________ + ! Compute the initial number number of elementa levels, based on the vertice + ! depth information + do n=1, elem2D + + !_______________________________________________________________________ + if (use_depthonelem) then + dmean = depth(n) ! depth is already defined on elements + else + nodes=elem2D_nodes(1:3,n) + !___________________________________________________________________ + ! depth of element is shallowest depth of sorounding vertices + if (trim(which_depth_n2e) .eq. 'min') then ; dmean=maxval(depth(nodes)) + ! depth of element is deepest depth of sorounding vertices + elseif (trim(which_depth_n2e) .eq. 'max') then ; dmean=minval(depth(nodes)) + ! DEFAULT: depth of element is mean depth of sorounding vertices + elseif (trim(which_depth_n2e) .eq. 'mean') then; dmean=sum(depth(nodes))/3.0 + end if + end if + + !_______________________________________________________________________ + exit_flag=0 + do nz=1,nl-1 + if(Z(nz)=0) nlevels(n)=thers_zbar_lev + + ! set minimum number of levels to --> thers_lev=5 + if(nlevels(n) do n=1, elem2D + + !___________________________________________________________________________ + ! write out vertical level indices before iterative geometric adaption to + ! exclude isolated cells + !_______________________________________________________________________ + file_name=trim(meshpath)//'elvls_raw.out' + open(fileID, file=file_name) + do n=1,elem2D + write(fileID,*) nlevels(n) + end do + close(fileID) + + !___________________________________________________________________________ + ! check for isolated cells (cells with at least two boundary faces or three + ! boundary vertices) and eliminate them --> FESOM2.0 doesn't like these kind + ! of cells + do nz=thers_zbar_lev+1,nl + exit_flag=0 + count_iter=0 + + !_______________________________________________________________________ + ! iteration loop within each layer + do while((exit_flag==0).and.(count_iter if elem2D_nodes(1,n) == elem2D_nodes(4,n): True --> q=3 --> triangular mesh + ! --> if elem2D_nodes(1,n) == elem2D_nodes(4,n): False --> q=4 --> quad mesh + nneighb = merge(3,4,elem2D_nodes(1,n) == elem2D_nodes(4,n)) + ! + ! +---isolated bottom cell + ! ._______________ | _______________________. + ! |###|###|###|###|___ | ___|###|###|###|###|###|###| + ! |###|###|###|###|###| | ___|###|###|###|###|###|###|###| + ! |###|###|###|###|###| | |###|###|###|###| BOTTOM |###| + ! |###|###|###|###|###|_v_|###|###|###|###|###|###|###|###| + ! |###|###|###|###|###|###|###|###|###|###|###|###|###|###| + ! + if (nlevels(n)>=nz) then + count_neighb_open=0 + elems=elem_neighbors(1:3,n) + !___________________________________________________________ + ! loop over neighbouring triangles + do i=1,nneighb + if (elems(i)>0) then + if (nlevels(elems(i))>=nz) then + !count neighbours + count_neighb_open=count_neighb_open+1 + endif + endif + enddo + + !___________________________________________________________ + ! check how many open faces to neighboring triangles the cell + ! has, if there are less than 2 its isolated (a cell should + ! have at least 2 valid neighbours) + if (count_neighb_open<2) then + ! if cell is "isolated", and the one levels shallower bottom + ! cell would be shallower than the minimum vertical level + ! treshhold (thers_lev). --> in this make sorrounding elements + ! one level deeper to reconnect the isolated cell + if (nz-10) then + nlevels(elems(i)) = max(nlevels(elems(i)),nz) + end if + end do + + !if cell is "isolated" convert to one level shallower bottom cell + else + nlevels(n)=nz-1 + end if + !force recheck for all current ocean cells + exit_flag=0 + + end if + end if ! --> if (nlevels(n)>=nz) then + end do ! --> do n=1,elem2D + end do ! --> do while((exit_flag==0).and.(count1<1000)) + write(*,"(A, I5, A, i5, A, I3)") ' -[iter ]->: nlevel, iter/maxiter=',count_iter,'/',max_iter,', nz=',nz + end do ! --> do nz=4,nl + + !___________________________________________________________________________ + ! vertical vertice level index of ocean bottom boundary + write(*,"(A)" ) ' -[compu]->: nlevels_nod2D ' + nlevels_nod2D=0 + do n=1,elem2D + q = merge(3,4,elem2D_nodes(1,n) == elem2D_nodes(4,n)) + do j=1,q + node=elem2D_nodes(j,n) + if(nlevels_nod2D(node): compute elem, vertice cavity depth index '//achar(27)//'[0m' + + !___________________________________________________________________________ + allocate(mesh%ulevels(elem2D)) + ulevels => mesh%ulevels + allocate(mesh%ulevels_nod2D(nod2D)) + ulevels_nod2D => mesh%ulevels_nod2D + + !___________________________________________________________________________ + ! Compute level position of ocean-cavity boundary + cavity_maxlev=0 + do elem=1, elem2D + + !_______________________________________________________________________ + if (use_cavityonelem) then + dmean = cavity_depth(elem) + else + nodes=elem2D_nodes(1:3,elem) + !_______________________________________________________________________ + ! depth of element is shallowest depth of sorounding vertices + if (trim(which_depth_n2e) .eq. 'min') then ; dmean=maxval(cavity_depth(nodes)) + ! depth of element is deepest depth of sorounding vertices + elseif (trim(which_depth_n2e) .eq. 'max') then ; dmean=minval(cavity_depth(nodes)) + ! DEFAULT: depth of element is mean depth of sorounding vertices + elseif (trim(which_depth_n2e) .eq. 'mean') then ; dmean=sum(cavity_depth(nodes))/3.0 + end if + end if + + !_______________________________________________________________________ + ! vertical elem level index of cavity-ocean boundary + ulevels(elem) = 1 + if (dmean<0.0_WP) ulevels(elem) = 2 + + do nz=1,nlevels(elem)-1 + if (Z(nz) should not be + ! possible in FESOM2.0 + ! loop over all cavity levels + allocate(elemreducelvl(elem2d),elemfixlvl(elem2d)) + allocate(numelemtonode(nl)) + + !___________________________________________________________________________ + ! outer iteration loop + count_iter2 = 0 + exit_flag2 = 0 + elemfixlvl = .False. + do while((exit_flag2==0) .and. (count_iter2 tri mesh, nneighb = 4 --> quad mesh + nneighb = merge(3,4,elem2D_nodes(1,elem) == elem2D_nodes(4,elem)) + elems = elem_neighbors(1:3,elem) + ! + ! .___________________________.~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! |###|###|###|###|###|###|###| + ! |# CAVITY |###| . |###|###| OCEAN + ! |###|###|###|###|/|\|###| + ! |###|###| | + ! |###| +-- Not good can lead to isolated cells + ! + + !___________________________________________________________ + ! (1st) Ask the Question: is nz at element elem an here an + ! valid layer in the ocean + if ( nz >= ulevels(elem) .and. nz0) then ! if its a valid boundary triangle, 0=missing value + ! check for isolated cell + if ( ulevels(elems(j))<= nz .and. & + nlevels(elems(j))> nz ) then + !count the open faces to neighboring cells + count_neighb_open=count_neighb_open+1 + endif + end if + end do ! --> do i = 1, nneighb + + !_______________________________________________________ + ! (3rd) check how many open faces to neighboring triangles the cell + ! has, if there are less than 2 its isolated (a cell should + ! have at least 2 valid neighbours) + ! --> in this case shift cavity-ocean interface one level down + if (count_neighb_open<=1) then + count_isoelem = count_isoelem+1 + + ! (4th.1 ) if cell elem is isolated convert it to a deeper ocean level + ! except when this levels would remain less than 3 valid + ! bottom levels --> in case make the levels of all sorounding + ! triangles shallower + if ( (nlevels(elem)-(nz+1))>=3 .and. & + (elemreducelvl(elem) .eqv. .False.) .and. & + (elemfixlvl( elem) .eqv. .False.) & + ) then + ulevels(elem)=nz+1 + + ! (4th.2) can not increase depth anymore to eleminate + ! isolated cell, otherwise lessthan 3 valid layers + ! therefor reduce depth of ONE!!! of the neighbouring + ! triangles. Choose triangle whos depth is already closest + ! to nz + else + !PS replace this with do j=1,3... because of + !PS indice -999 conflict in elems, ulevels(-999) + !PS not possible + !PS idx = minloc(ulevels(elems)-nz, 1, MASK=( (elems>0) .and. ((ulevels(elems)-nz)>0) ) ) + val=nl + do j = 1, 3 + if (elems(j)>0) then ! if its a valid boundary triangle, 0=missing value + if (ulevels(elems(j))-nz>0 .and. ulevels(elems(j))-nz do i = 1, nneighb + + ulevels( elems(idx)) = nz + elemreducelvl(elems(idx)) = .True. + end if + + ! force recheck for all current ocean cells + exit_flag1=0 + end if ! --> if (count_neighb_open<2) then + + end if ! --> if ( nz >= ulevels(elem) .and. nz do elem=1,elem2D + + end do ! --> do while((exit_flag==0).and.(count_iter<1000)) + write(*,"(A, I5, A, i5, A, I3)") ' -[iter ]->: ulevel, iter/maxiter=',count_iter,'/',max_iter,', nz=',nz + end do ! --> do nz=1,cavity_maxlev + + !_______________________________________________________________________ + ! compute vertical vertice level index of cavity_ocean boundary + write(*,"(A)" ) ' -[compu]->: ulevels_nod2D ' + ulevels_nod2D = nl + do elem=1,elem2D + nneighb = merge(3,4,elem2D_nodes(1,elem) == elem2D_nodes(4,elem)) + !___________________________________________________________________ + ! loop over neighbouring triangles + do j=1,nneighb + node=elem2D_nodes(j,elem) + ulevels_nod2D(node)=min(ulevels_nod2D(node),ulevels(elem)) + end do + end do ! --> do elem=1,elem2D + + !_______________________________________________________________________ + ! check if all constrains for ulevel and ulevels_nod2D is fullfilled + exit_flag2 = 1 + do elem=1,elem2D + if (ulevels(elem)>=nlevels(elem)) then + write(*,*) ' -[check]->: elem cavity depth deeper or equal bottom depth, elem=',elem + exit_flag2 = 0 + end if + + if (nlevels(elem)-ulevels(elem)<3) then + write(*,*) ' -[check]->: less than three valid elem ocean layers, elem=',elem + exit_flag2 = 0 + end if + end do ! --> do elem=1,elem2D + + do node=1,nod2D + if (ulevels_nod2D(node)>=nlevels_nod2D(node)) then + write(*,*) ' -[check]->: vertice cavity depth deeper or equal bottom depth, node=', node + exit_flag2 = 0 + end if + + if (nlevels_nod2D(node)-ulevels_nod2D(node)<3) then + write(*,*) ' -[check]->: less than three valid vertice ocean layers, node=', node + exit_flag2 = 0 + end if + end do ! --> do node=1,nod2D + + do elem=1,elem2D + if (ulevels(elem)< maxval(ulevels_nod2D(elem2D_nodes(:,elem))) ) then + write(*,*) ' -[check]->: found elem cavity shallower than its valid maximum cavity vertice depths, elem=', elem2d + exit_flag2 = 0 + end if + end do ! --> do elem=1,elem2D + + !_______________________________________________________________________ + ! compute how many triangle elements contribute to every vertice in every + ! layer + ! + ! --> What can happen is that a node point in the middle of the vertical + ! domain can become isolated due to the cavity constrains. The model + ! would not be able to deal with this kind of situation. So we must + ! prevent it by adapting ulevels! + ! O node + ! _._ + ! _/ | \_ + ! _/ | \_ + ! _/ | \_ + ! elem(1) elem(2) elem(3)... <--elem=nod_in_elem2D(j,node) + ! ._______. ulevel(elem2)=30 + ! |_______| + ! |_______| + ! |_______| + ! |_______| + ! |_______| nlevel(elem2)=38 + ! + ! In this possible gap node points + ! would have no neighboring elements + ! + ! ulevel(elem1)=42 ._______. ._______. ulevel(elem3)=42 + ! |_______| |_______| + ! |_______| |_______| + ! |_______| |_______| + ! |_______| |_______| + ! nlevel(elem1)=46 |_______| |_______| + ! |_______| nlevel(elem3)=48 + ! + ! --> Problem here is we want to keep nlevels fixed so what we can do is + ! to set ulevels(elem1) and ulevels(elem3) towards nlevel(elem2) + count_iter=0 + do node=1, nod2D + !___________________________________________________________________ + ! check if there is a possible gap as described above that would + ! allow for node points without neighbors + min_nlvl = nl + max_ulvl = 1 + do j=1, nod_in_elem2D_num(node) + elem=nod_in_elem2D(j, node) + min_nlvl = min(min_nlvl, nlevels(elem)) + max_ulvl = max(max_ulvl, ulevels(elem)) + end do + + ! found a potential gap + if (min_nlvl < max_ulvl) then + + !_______________________________________________________________ + ! compute how many triangle elements contribute to vertice in + ! every layer check if there are layers where a node point has + ! only one or even zero neighboring triangles. + numelemtonode=0 + do j=1, nod_in_elem2D_num(node) + elem=nod_in_elem2D(j, node) + do nz=ulevels(elem), nlevels(elem)-1 + numelemtonode(nz) = numelemtonode(nz) + 1 + end do + end do + + !_______________________________________________________________ + ! check in which depth level is an isolated node + nzloop: do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 + !___________________________________________________________ + ! nodes has zero neighbouring triangles and is completely + ! isolated need to adapt ulevels --> inflicts another + ! outher iteration loop (exit_flag2=0). It needs at least + ! two neighboring elements so the node is considered as + ! connected. Search the index of the two elements where ulevels>nz + ! but that are closest to nz + if (numelemtonode(nz)==0) then + count_iter = count_iter+1 + write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz) ,' neighb. triangle: n=', node, ', nz=',nz + !_______________________________________________________ + allocate(aux_arr(nod_in_elem2D_num(node)), aux_idx(nod_in_elem2D_num(node))) + aux_arr(:) = ulevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) + aux_arr(:) = aux_arr(:) - nz + ! fill array with index of element + do j=1, nod_in_elem2D_num(node) + aux_idx(j) = j + end do + ! index of closest elem to nz where ulevel>nz + idx = minloc(aux_arr, 1, MASK=((aux_arr>0)) ) + ! index of second closest elem to nz where ulevel>nz + idx2 = minloc(aux_arr, 1, MASK=((aux_arr>0) .and. (aux_idx/=idx)) ) + deallocate(aux_arr, aux_idx) + ulevels( nod_in_elem2D(idx ,node)) = nz + ulevels( nod_in_elem2D(idx2,node)) = nz + elemfixlvl(nod_in_elem2D(idx ,node)) = .True. + elemfixlvl(nod_in_elem2D(idx2,node)) = .True. + + !_______________________________________________________ + ! inflict another outer iteration loop + exit_flag2 = 0 + + !_______________________________________________________ + ! if the upper most isolated layer is fixed all layers below should be fixed as well + ! --> exit loop do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 + exit nzloop + + !___________________________________________________________ + ! nodes has one neighbouring triangles it needs at least + ! another neighboring elements so the node is considered as + ! connected + elseif (numelemtonode(nz)==1) then + count_iter = count_iter+1 + write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz) ,'neighb. triangle: n=', node, ', nz=',nz + !_______________________________________________________ + allocate(aux_arr(nod_in_elem2D_num(node)), aux_idx(nod_in_elem2D_num(node))) + aux_arr(:) = ulevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) + aux_arr(:) = aux_arr(:) - nz + ! fill array with index of element + do j=1, nod_in_elem2D_num(node) + aux_idx(j) = j + end do + ! index of closest elem to nz where ulevel>nz + idx = minloc(aux_arr, 1, MASK=((aux_arr>0)) ) + deallocate(aux_arr, aux_idx) + ulevels( nod_in_elem2D(idx,node)) = nz + elemfixlvl(nod_in_elem2D(idx,node)) = .True. + + !_______________________________________________________ + ! inflict another outer iteration loop + exit_flag2 = 0 + + !_______________________________________________________ + ! if the upper most isolated layer is fixed all layers below should be fixed as well + ! --> exit loop do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 + exit nzloop + + end if + end do nzloop ! --> do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 + end if ! --> if (min_nlvl < max_ulvl) then + end do ! --> do node=1, nod2D + + !_______________________________________________________________________ + ! check if cavity geometry did converge + if (exit_flag2 == 0) then + print *, achar(27)//'[33m' //'____________________________________________________________'//achar(27)//'[0m' + print *, ' -['//achar(27)//'[33m'//'WARN'//achar(27)//'[0m'//']->: Cavity geom. not converged yet, do further outer iteration' + write(*,"(A, I3, A, I3)") ' iter step ', count_iter2,'/', max_iter2 + write(*,*) + end if + + !_______________________________________________________________________ + end do ! --> do while((exit_flag2==0) .and. (count_iter2: Cavity geometry constrains did not converge !!! *\(>︿<)/*'//achar(27)//'[0m' + write(*,*) + stop + else + write(*,*) + print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' + print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge !!! *\(^o^)/*' + end if + + !___________________________________________________________________________ + ! write out cavity mesh files for vertice and elemental position of + ! vertical cavity-ocean boundary + ! write out elemental cavity-ocean boundary level + file_name=trim(meshpath)//'cavity_elvls.out' + open(20, file=file_name) + do elem=1,elem2D + write(20,*) ulevels(elem) + enddo + close(20) + + ! write out vertice cavity-ocean boundary level + yes/no cavity flag + file_name=trim(meshpath)//'cavity_nlvls.out' + open(20, file=file_name) + do node=1,nod2D + write(20,*) ulevels_nod2D(node) + enddo + close(20) + +end subroutine find_levels_cavity + +!=================================================================== +subroutine edge_center(n1, n2, x, y, mesh) +! +! Returns coordinates of edge center in x and y +! +implicit none +integer :: n1, n2 ! nodes of the edge +real(kind=WP) :: x, y, a(2), b(2) +type(t_mesh), intent(in), target :: mesh +#include "associate_mesh_ini.h" + +a=coord_nod2D(:,n1) +b=coord_nod2D(:,n2) +if(a(1)-b(1)>cyclic_length/2.0) a(1)=a(1)-cyclic_length +if(a(1)-b(1)<-cyclic_length/2.0) b(1)=b(1)-cyclic_length +x=0.5_WP*(a(1)+b(1)) +y=0.5_WP*(a(2)+b(2)) +end subroutine edge_center + +!==================================================================== +subroutine elem_center(elem, x, y, mesh) +! +! Returns coordinates of elem center in x and y +! +implicit none +integer, intent(in) :: elem +integer :: elnodes(3), k +real(kind=WP), intent(out) :: x, y +real(kind=WP) :: ax(3), amin +type(t_mesh), intent(in), target :: mesh +#include "associate_mesh_ini.h" + + elnodes=elem2D_nodes(1:3,elem) + ax=coord_nod2D(1, elnodes) + amin=minval(ax) + DO k=1,3 + if(ax(k)-amin>cyclic_length/2.0) ax(k)=ax(k)-cyclic_length + END DO + x=sum(ax)/3.0_WP + y=sum(coord_nod2D(2,elnodes))/3.0_WP + +end subroutine elem_center + +!=================================================================== +! Stiffness matrix for the elevation +subroutine stiff_mat_ini(mesh) + + ! + implicit none + integer :: i, j, n, q, el, elem_nodes_max, nod(4) + integer, allocatable :: num_ne(:), ne(:,:) + ! + type(t_mesh), intent(inout), target :: mesh +#include "associate_mesh_ini.h" + + ssh_stiff%dim = nod2D + allocate(mesh%ssh_stiff%rowptr(nod2D+1)) + ssh_stiff => mesh%ssh_stiff !required after the allocation, otherwise the pointer remains undefined + + allocate(num_ne(nod2D), ne(MAX_ADJACENT,nod2D)) + num_ne(:) = 0 + ne(:,:) = -1 + + ! Check the maximum number of nodes in an element (FESOM triangular meshes = 3, Hybrid meshes = 4) + elem_nodes_max = size(elem2D_nodes, 1) + + ! Fill node adjacency info + ! all nodes in an element are adjacent in the sense of being halo nodes + ! (also the opposite nodes of a quad: no edge, but the indirect connection + ! should be taken into account by metis domain decomposition) + do el=1,elem2D + nod(1:elem_nodes_max) = elem2D_nodes(1:elem_nodes_max,el) ! Fortran-numbering + q = elem_nodes_max + if (nod(1) == nod(elem_nodes_max)) q = q-1 ! triangle + + do i=2,q + do j=1,i-1 + if (all(ne(:,nod(i)) /= nod(j))) then + num_ne(nod(i)) = num_ne(nod(i)) + 1 + num_ne(nod(j)) = num_ne(nod(j)) + 1 + + if (max(num_ne(nod(i)), num_ne(nod(j))) > MAX_ADJACENT ) then + print *,'Recompile with larger value for MAX_ADJACENT.' + stop + else + ne(num_ne(nod(i)), nod(i)) = nod(j) + ne(num_ne(nod(j)), nod(j)) = nod(i) + endif + endif + end do + end do + end do + +! copy adjacency matrix to CSR-format + ssh_stiff%rowptr(1) = 1 + do n=1,nod2D + ssh_stiff%rowptr(n+1) = ssh_stiff%rowptr(n)+num_ne(n) + end do + + allocate(ssh_stiff%colind(ssh_stiff%rowptr(nod2D+1)-1)) + ssh_stiff => mesh%ssh_stiff !required after the allocation, otherwise the pointer remains undefined + + !required after the allocation, otherwise the pointer remains undefined + do n=1,nod2D + ssh_stiff%colind(ssh_stiff%rowptr(n):ssh_stiff%rowptr(n+1)-1) = ne(1:num_ne(n),n) + end do + + deallocate(num_ne, ne) + +end subroutine stiff_mat_ini + +!=================================================================== +! Setup of communication arrays +subroutine communication_ini(partit, mesh) + use MOD_PARTIT + use mod_gen_comm, only: communication_nodn, communication_elemn + use omp_lib + implicit none + + interface + subroutine save_dist_mesh(partit, mesh) + use MOD_MESH + use MOD_PARTIT + type(t_mesh), intent(in), target :: mesh + type(t_partit), intent(in), target :: partit + end subroutine save_dist_mesh + end interface + + integer :: n + character*10 :: npes_string + character(MAX_PATH) :: dist_mesh_dir + LOGICAL :: L_EXISTS + type(t_mesh), intent(inout), target :: mesh + type(t_partit), intent(inout), target :: partit +#include "associate_part_def.h" +#include "associate_mesh_ini.h" +#include "associate_part_ass.h" + + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute communication arrays '//achar(27)//'[0m' + + ! Create the distributed mesh subdirectory + write(npes_string,"(I10)") npes + dist_mesh_dir=trim(meshpath)//'dist_'//trim(ADJUSTL(npes_string))//'/' + INQUIRE(file=trim(dist_mesh_dir), EXIST=L_EXISTS) + if (.not. L_EXISTS) call system('mkdir '//trim(dist_mesh_dir)) + +#ifdef OMP_MAX_THREADS +!$OMP PARALLEL NUM_THREADS(OMP_MAX_THREADS) + if (OMP_GET_THREAD_NUM() == 0) then + write(*,*) 'Setting up communication arrays using ', OMP_GET_NUM_THREADS(), ' threads' + endif +#else +!$OMP PARALLEL NUM_THREADS(1) + write(*,*) 'Setting up communication arrays using 1 thread (serially)' +#endif + +!$OMP DO + do n = 0, npes-1 + mype = n ! mype is threadprivate and must not be iterator + call communication_nodn(partit, mesh) + call communication_elemn(partit, mesh) + call save_dist_mesh(partit, mesh) ! Write out communication file com_infoxxxxx.out + end do +!$OMP END DO +!$OMP END PARALLEL + + deallocate(mesh%elem_neighbors) + deallocate(mesh%elem_edges) + deallocate(partit%part) + write(*,*) 'Communication arrays have been set up' +end subroutine communication_ini + +!================================================================= +subroutine set_par_support_ini(partit, mesh) + use iso_c_binding, only: idx_t=>C_INT32_T + use MOD_PARTIT + implicit none + + integer :: n, j, k, nini, nend, ierr + integer(idx_t) :: np(10) + type(t_mesh), intent(inout), target :: mesh + type(t_partit), intent(inout), target :: partit + interface + subroutine do_partit(n,ptr,adj,wgt,np,part) bind(C) + use iso_c_binding, only: idx_t=>C_INT32_T + integer(idx_t), intent(in) :: n, ptr(*), adj(*), wgt(*), np(*) + integer(idx_t), intent(out) :: part(*) + end subroutine do_partit + end interface + +#include "associate_part_def.h" +#include "associate_mesh_ini.h" +#include "associate_part_ass.h" + + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute partitioning '//achar(27)//'[0m' + end if + + ! Construct partitioning vector + if (n_levels<1 .OR. n_levels>10) then + print *,'Number of hierarchic partition levels is out of range [1-10]! Aborting...' + stop + end if + + np(:) = n_part(:) ! Number of partitions on each hierarchy level + if (n_part(1) == 0) then ! Backward compatibility case: Take the number of + np(1) = npes ! partitions from the number of MPI processes + n_levels = 1 + end if + if (n_levels < 10) then ! 0 is an indicator of the last hierarchy level + np(n_levels+1) = 0 + end if + + allocate(partit%part(nod2D)) + part=>partit%part + part=0 + + npes = PRODUCT(np(1:n_levels)) + if(npes<2) then + print *,'Total number of parallel partitions is less than one! Aborting...' + stop + end if + + write(*,*) 'Calling partit for npes=', np + call do_partit(ssh_stiff%dim, ssh_stiff%rowptr, ssh_stiff%colind, & + nlevels_nod2D, np, part) + +write(*,*) 'DONE' +write(*,*) size(partit%part) + call check_partitioning(partit, mesh) + + write(*,*) 'Partitioning is done.' + +! The stiffness matrix is no longer needed. + deallocate(mesh%ssh_stiff%rowptr) + deallocate(mesh%ssh_stiff%colind) + + !NR No longer needed - last use was as weight for partitioning + deallocate(mesh%nlevels_nod2D) +end subroutine set_par_support_ini + +!======================================================================= +subroutine check_partitioning(partit, mesh) + + ! In general, METIS 5 has several advantages compared to METIS 4, e.g., + ! * neighbouring tasks get neighbouring partitions (important for multicore computers!) + ! * lower maximum of weights per partition (better load balancing) + ! * lower memory demand + ! + ! BUT: there might be outliers, single nodes connected to their partition by + ! only one edge or even completely isolated. This spoils everything :-( + ! + ! This routine checks for isolated nodes and moves them to an adjacent partition, + ! trying not to spoil the load balance. + + use MOD_PARTIT + type(t_partit), intent(inout), target :: partit + type(t_mesh), intent(inout), target :: mesh + integer :: i, j, k, n, n_iso, n_iter, is, ie, kmax, np, cnt + integer :: nod_per_partition(2,0:partit%npes-1) + integer :: max_nod_per_part(2), min_nod_per_part(2) + integer :: average_nod_per_part(2), node_neighb_part(100) + logical :: already_counted, found_part + + integer :: max_adjacent_nodes + integer, allocatable :: ne_part(:), ne_part_num(:), ne_part_load(:,:) +#include "associate_part_def.h" +#include "associate_mesh_ini.h" +#include "associate_part_ass.h" +!just for partit%part + + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: check partitioning '//achar(27)//'[0m' + end if + + ! Check load balancing + do i=0, npes-1 + nod_per_partition(1,i) = count(part(:) == i) + nod_per_partition(2,i) = sum(nlevels_nod2D, part(:) == i) + enddo + + min_nod_per_part(1) = minval( nod_per_partition(1,:)) + min_nod_per_part(2) = minval( nod_per_partition(2,:)) + + max_nod_per_part(1) = maxval( nod_per_partition(1,:)) + max_nod_per_part(2) = maxval( nod_per_partition(2,:)) + + average_nod_per_part(1) = nod2D / npes + average_nod_per_part(2) = sum(nlevels_nod2D(:)) / npes + + ! Now check for isolated nodes (connect by one or even no edge to other + ! nodes of its partition) and repair, if possible + + max_adjacent_nodes = maxval(ssh_stiff%rowptr(2:nod2D+1) - ssh_stiff%rowptr(1:nod2D)) + allocate(ne_part(max_adjacent_nodes), ne_part_num(max_adjacent_nodes), & + ne_part_load(2,max_adjacent_nodes)) + + print *,' ' + print *,'Check for isolated nodes ========' + n_iso = 0 + checkloop: do n=1,nod2D + is = ssh_stiff%rowptr(n) + ie = ssh_stiff%rowptr(n+1) -1 + cnt = ie-is+1 + + node_neighb_part(1:cnt) = part(ssh_stiff%colind(is:ie)) + if (count(node_neighb_part(1:cnt) == part(n)) <= 1) then + + n_iso = n_iso+1 + print *,'Isolated node',n, 'in partition', part(n) + print *,'Neighbouring nodes are in partitions', node_neighb_part(1:cnt) + + ! count the adjacent nodes of the other PEs + + np=1 + ne_part(1) = node_neighb_part(1) + ne_part_num(1) = 1 + ne_part_load(1,1) = nod_per_partition(1,ne_part(1)) + 1 + ne_part_load(2,1) = nod_per_partition(2,ne_part(1)) + nlevels_nod2D(n) + + do i=1,cnt + if (node_neighb_part(i)==part(n)) cycle + already_counted = .false. + do k=1,np + if (node_neighb_part(i) == ne_part(k)) then + ne_part_num(k) = ne_part_num(k) + 1 + already_counted = .true. + exit + endif + enddo + if (.not. already_counted) then + np = np+1 + ne_part(np) = node_neighb_part(i) + ne_part_num(np) = 1 + ne_part_load(1,np) = nod_per_partition(1,ne_part(np)) + 1 + ne_part_load(2,np) = nod_per_partition(2,ne_part(np)) + nlevels_nod2D(n) + endif + enddo + + ! Now, check for two things: The load balance, and if + ! there is more than one node of that partition. + ! Otherwise, it would become isolated again. + + ! Best choice would be the partition with most adjacent nodes (edgecut!) + ! Choose, if it does not decrease the load balance. + ! (There might be two partitions with the same number of adjacent + ! nodes. Don't care about this here) + + kmax = maxloc(ne_part_num(1:np),1) + + if (ne_part_num(kmax) <= 1) then + ! No chance - this is probably a boundary + ! node that has only two neighbors. + cycle checkloop + endif + + if (ne_part_load(1,kmax) <= max_nod_per_part(1) .and. & + ne_part_load(2,kmax) <= max_nod_per_part(2) ) then + k = kmax + else + ! Don't make it too compicated. Reject partitions that have only one + ! adjacent node. Take the next not violating the load balance. + found_part = .false. + do k=1,np + if (ne_part_num(k)==1 .or. k==kmax) cycle + + if (ne_part_load(1,k) <= max_nod_per_part(1) .and. & + ne_part_load(2,k) <= max_nod_per_part(2) ) then + + found_part = .true. + exit + endif + enddo + + if (.not. found_part) then + ! Ok, don't think to much. Simply go for minimized edge cut. + k = kmax + endif + endif + + ! Adjust the load balancing + + nod_per_partition(1,ne_part(k)) = nod_per_partition(1,ne_part(k)) + 1 + nod_per_partition(2,ne_part(k)) = nod_per_partition(2,ne_part(k)) + nlevels_nod2D(n) + nod_per_partition(1,part(n)) = nod_per_partition(1,part(n)) - 1 + nod_per_partition(2,part(n)) = nod_per_partition(2,part(n)) - nlevels_nod2D(n) + + ! And, finally, move nod n to other partition + part(n) = ne_part(k) + print *,'Node',n,'is moved to part',part(n) + endif + enddo checkloop + + deallocate(ne_part, ne_part_num, ne_part_load) + + print *,'=== LOAD BALANCING ===' + print *,'2D nodes: min, aver, max per part',min_nod_per_part(1), & + average_nod_per_part(1),max_nod_per_part(1) + + write(*,"('2D nodes: percent min, aver, max ',f8.3,'%, 100%, ',f8.3,'%')") & + 100.*real(min_nod_per_part(1)) / real(average_nod_per_part(1)), & + 100.*real(max_nod_per_part(1)) / real(average_nod_per_part(1)) + + print *,'3D nodes: Min, aver, max per part',min_nod_per_part(2), & + average_nod_per_part(2),max_nod_per_part(2) + write(*,"('3D nodes: percent min, aver, max ',f8.3,'%, 100%, ',f8.3,'%')") & + 100.*real(min_nod_per_part(2)) / real(average_nod_per_part(2)), & + 100.*real(max_nod_per_part(2)) / real(average_nod_per_part(2)) + +end subroutine check_partitioning + +end module mod_mesh_utils diff --git a/src/oce_adv_tra_driver.F90 b/src/oce_adv_tra_driver.F90 index 9b9606ea8..5def28377 100644 --- a/src/oce_adv_tra_driver.F90 +++ b/src/oce_adv_tra_driver.F90 @@ -1,45 +1,22 @@ -module oce_adv_tra_driver_interfaces - interface - subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, mesh) - use MOD_MESH - use MOD_TRACER - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - 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 - type(t_dyn) , intent(inout), target :: dynamics - 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 do_oce_adv_tra - end interface -end module oce_adv_tra_driver_interfaces +module oce_adv_tra_driver_module + use MOD_MESH + use MOD_TRACER + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + use g_comm_auto + use diagnostics + use oce_adv_tra_hor_interfaces + use oce_adv_tra_ver_interfaces + use oce_adv_tra_fct_module, only: oce_tra_adv_fct + + implicit none + + private + public :: do_oce_adv_tra, oce_tra_adv_flux2dtracer + +contains -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 - USE MOD_PARSUP - 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 oce_tra_adv_flux2dtracer - end interface -end module oce_tra_adv_flux2dtracer_interface ! ! !=============================================================================== @@ -53,8 +30,8 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, use diagnostics, only: ldiag_DVD use oce_adv_tra_hor_interfaces use oce_adv_tra_ver_interfaces - use oce_adv_tra_fct_interfaces - use oce_tra_adv_flux2dtracer_interface + use oce_adv_tra_fct_module, only: oce_tra_adv_fct + ! oce_tra_adv_flux2dtracer is now in the same module implicit none real(kind=WP), intent(in), target :: dt integer, intent(in) :: tr_num @@ -573,3 +550,5 @@ subroutine oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, flux_h, flux_v, partit, #endif end subroutine oce_tra_adv_flux2dtracer + +end module oce_adv_tra_driver_module diff --git a/src/oce_adv_tra_fct.F90 b/src/oce_adv_tra_fct.F90 index 1083843aa..b239cf52b 100644 --- a/src/oce_adv_tra_fct.F90 +++ b/src/oce_adv_tra_fct.F90 @@ -1,36 +1,17 @@ -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 - USE MOD_PARSUP - type(t_mesh), intent(in), target :: mesh - type(t_partit),intent(inout), target :: partit - type(t_tracer_work), intent(inout), target :: twork - end subroutine oce_adv_tra_fct_init +module oce_adv_tra_fct_module + use MOD_MESH + use MOD_TRACER + USE MOD_PARTIT + USE MOD_PARSUP + USE g_comm_auto + + implicit none + + private + public :: oce_adv_tra_fct_init, oce_tra_adv_fct + +contains - 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 - USE MOD_PARSUP - 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_nod2D) - real(kind=WP), intent(inout) :: fct_minus(mesh%nl, partit%myDim_nod2D) - real(kind=WP), intent(inout) :: AUX(:,:,:) !a large auxuary array - end subroutine oce_tra_adv_fct - end interface -end module oce_adv_tra_fct_interfaces -! -! !=============================================================================== subroutine oce_adv_tra_fct_init(twork, partit, mesh) use MOD_MESH @@ -498,3 +479,5 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, !$ACC END DATA #endif end subroutine oce_tra_adv_fct + +end module oce_adv_tra_fct_module diff --git a/src/oce_adv_tra_hor.F90 b/src/oce_adv_tra_hor.F90 index 3586f683c..18086a8c8 100644 --- a/src/oce_adv_tra_hor.F90 +++ b/src/oce_adv_tra_hor.F90 @@ -63,6 +63,7 @@ end module oce_adv_tra_hor_interfaces !=============================================================================== subroutine adv_tra_hor_upw1(vel, ttf, partit, mesh, flux, o_init_zero) use MOD_MESH + use O_PARAM, only: r_earth USE MOD_PARTIT USE MOD_PARSUP use g_comm_auto diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 11f7601e5..319b7d52c 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1,195 +1,28 @@ - -module compute_CFLz_interface - interface - subroutine compute_CFLz(dynamics, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine compute_CFLz - end interface -end module compute_CFLz_interface - -module compute_Wvel_split_interface - interface - subroutine compute_Wvel_split(dynamics, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine compute_Wvel_split - end interface -end module compute_Wvel_split_interface - -module compute_vert_vel_transpv_interface - interface - subroutine compute_vert_vel_transpv(dynamics, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine compute_vert_vel_transpv - end interface -end module compute_vert_vel_transpv_interface - -module oce_ale_interfaces - interface - subroutine init_bottom_elem_thickness(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine init_bottom_elem_thickness - - subroutine init_bottom_node_thickness(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine init_bottom_node_thickness - - subroutine init_surface_elem_depth(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine init_surface_elem_depth - - subroutine init_surface_node_depth(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine init_surface_node_depth - - subroutine impl_vert_visc_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine impl_vert_visc_ale - - subroutine update_stiff_mat_ale(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine update_stiff_mat_ale - - subroutine compute_ssh_rhs_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine compute_ssh_rhs_ale - - subroutine solve_ssh_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine solve_ssh_ale - - subroutine compute_hbar_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine compute_hbar_ale - - subroutine vert_vel_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine vert_vel_ale - - subroutine update_thickness_ale(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine update_thickness_ale - end interface -end module oce_ale_interfaces - -module init_ale_interface - interface - subroutine init_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine init_ale - end interface -end module init_ale_interface - -module init_thickness_ale_interface - interface - subroutine init_thickness_ale(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine init_thickness_ale - end interface -end module init_thickness_ale_interface - -module oce_timestep_ale_interface - interface - subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use MOD_DYN - use MOD_ICE - integer , intent(in) :: n - type(t_dyn) , intent(inout), target :: dynamics - type(t_ice), intent(inout), target :: ice - type(t_tracer), intent(inout), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout), target :: mesh - end subroutine oce_timestep_ale - end interface -end module oce_timestep_ale_interface - +module oce_ale_module + USE o_PARAM + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + USE o_ARRAYS + USE g_config + USE g_forcing_param, only: use_virt_salt + use solver_module, only: ssh_solve_preconditioner, ssh_solve_cg + use oce_dyn_module, only: viscosity_filter, update_vel, compute_apegen, check_viscopt, compute_ke_wrho + use oce_ale_pressure_bv_module, only: pressure_bv, pressure_force_4_linfs, pressure_force_4_zxxxx, & + sw_alpha_beta, compute_neutral_slope, compute_sigma_xy + + implicit none + + private + public :: init_ale, init_bottom_elem_thickness, init_bottom_node_thickness, & + init_surface_elem_depth, init_surface_node_depth, init_thickness_ale, & + update_thickness_ale, impl_vert_visc_ale, update_stiff_mat_ale, & + compute_ssh_rhs_ale, solve_ssh_ale, compute_hbar_ale, vert_vel_ale, & + oce_timestep_ale, compute_CFLz, compute_Wvel_split, compute_vert_vel_transpv, & + init_stiff_mat_ale, restart_thickness_ale + +contains ! CONTENT: ! ------------ @@ -226,7 +59,6 @@ subroutine init_ale(dynamics, partit, mesh) USE g_config, only: which_ale, use_cavity, use_partial_cell, ib_async_mode USE g_forcing_param, only: use_virt_salt - use oce_ale_interfaces Implicit NONE ! kh 18.03.21 @@ -2115,8 +1947,6 @@ subroutine vert_vel_ale(dynamics, partit, mesh) use g_comm_auto use io_RESTART !!PS use g_forcing_arrays !!PS - use compute_Wvel_split_interface - use compute_CFLz_interface implicit none type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit @@ -2685,8 +2515,6 @@ subroutine compute_vert_vel_transpv(dynamics, partit, mesh) use o_ARRAYS, only: water_flux use g_config, only: dt, which_ale use g_comm_auto - use compute_Wvel_split_interface - use compute_CFLz_interface implicit none !___________________________________________________________________________ type(t_dyn) , intent(inout), target :: dynamics @@ -3062,8 +2890,6 @@ subroutine solve_ssh_ale(dynamics, partit, mesh) USE MOD_DYN use g_comm_auto use g_config, only: which_ale - use ssh_solve_preconditioner_interface - use ssh_solve_cg_interface implicit none type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit @@ -3327,14 +3153,9 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) use g_cvmix_kpp use g_cvmix_tidal use Toy_Channel_Soufflet - use oce_ale_interfaces - use compute_vert_vel_transpv_interface use compute_ssh_split_explicit_interface - use pressure_bv_interface - use pressure_force_4_linfs_interface - use pressure_force_4_zxxxx_interface - use compute_vel_rhs_interface - use solve_tracers_ale_interface + use oce_ale_vel_rhs_module, only: compute_vel_rhs + use oce_ale_tracer_module, only: solve_tracers_ale use write_step_info_interface use check_blowup_interface use fer_solve_interface @@ -3853,4 +3674,4 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) write(*,*) end if end subroutine oce_timestep_ale - +end module diff --git a/src/oce_ale_pressure_bv.F90 b/src/oce_ale_pressure_bv.F90 index 5be3b0c46..d8c730b2d 100644 --- a/src/oce_ale_pressure_bv.F90 +++ b/src/oce_ale_pressure_bv.F90 @@ -1,193 +1,107 @@ -module densityJM_components_interface - interface - subroutine densityJM_components(t, s, bulk_0, bulk_pz, bulk_pz2, rhopot) - USE MOD_PARSUP - USE o_param - real(kind=WP), intent(IN) :: t,s - real(kind=WP), intent(OUT) :: bulk_0, bulk_pz, bulk_pz2, rhopot - end subroutine densityJM_components - end interface -end module densityJM_components_interface - -module density_linear_interface - interface - subroutine density_linear(t, s, bulk_0, bulk_pz, bulk_pz2, rho_out) - USE MOD_PARSUP - USE o_param - real(kind=WP), intent(IN) :: t,s - real(kind=WP), intent(OUT) :: bulk_0, bulk_pz, bulk_pz2, rho_out - end subroutine density_linear - end interface -end module density_linear_interface - -module pressure_force_4_linfs_fullcell_interface - interface - subroutine pressure_force_4_linfs_fullcell(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine pressure_force_4_linfs_fullcell - end interface -end module pressure_force_4_linfs_fullcell_interface -module pressure_force_4_linfs_nemo_interface - interface - subroutine pressure_force_4_linfs_nemo(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer), intent(in), target :: tracers - end subroutine pressure_force_4_linfs_nemo - end interface -end module pressure_force_4_linfs_nemo_interface -module pressure_force_4_linfs_shchepetkin_interface - interface - subroutine pressure_force_4_linfs_shchepetkin(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine pressure_force_4_linfs_shchepetkin - end interface -end module pressure_force_4_linfs_shchepetkin_interface -module pressure_force_4_linfs_easypgf_interface - interface - subroutine pressure_force_4_linfs_easypgf(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_tracer), intent(in), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - end subroutine pressure_force_4_linfs_easypgf - end interface -end module pressure_force_4_linfs_easypgf_interface -module pressure_force_4_linfs_cubicspline_interface - interface - subroutine pressure_force_4_linfs_cubicspline(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine pressure_force_4_linfs_cubicspline - end interface -end module pressure_force_4_linfs_cubicspline_interface -module pressure_force_4_linfs_cavity_interface - interface - subroutine pressure_force_4_linfs_cavity(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine pressure_force_4_linfs_cavity - end interface -end module pressure_force_4_linfs_cavity_interface -module pressure_force_4_zxxxx_shchepetkin_interface - interface - subroutine pressure_force_4_zxxxx_shchepetkin(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine pressure_force_4_zxxxx_shchepetkin - end interface -end module pressure_force_4_zxxxx_shchepetkin_interface -module pressure_force_4_zxxxx_easypgf_interface - interface - subroutine pressure_force_4_zxxxx_easypgf(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer), intent(in), target :: tracers - end subroutine pressure_force_4_zxxxx_easypgf - end interface -end module pressure_force_4_zxxxx_easypgf_interface -module pressure_force_4_zxxxx_cubicspline_interface - interface - subroutine pressure_force_4_zxxxx_cubicspline(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine pressure_force_4_zxxxx_cubicspline - end interface -end module pressure_force_4_zxxxx_cubicspline_interface -module init_ref_density_interface - interface - subroutine init_ref_density(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine init_ref_density - end interface -end module init_ref_density_interface -module insitu2pot_interface - interface - subroutine insitu2pot(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer), intent(in), target :: tracers - end subroutine insitu2pot - end interface -end module insitu2pot_interface -module pressure_bv_interface - interface - subroutine pressure_bv(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer), intent(in), target :: tracers - end subroutine pressure_bv - end interface -end module pressure_bv_interface -module pressure_force_4_linfs_interface - interface - subroutine pressure_force_4_linfs(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer), intent(in), target :: tracers - end subroutine pressure_force_4_linfs - end interface -end module pressure_force_4_linfs_interface -module pressure_force_4_zxxxx_interface - interface - subroutine pressure_force_4_zxxxx(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer), intent(in), target :: tracers - end subroutine pressure_force_4_zxxxx - end interface -end module pressure_force_4_zxxxx_interface +module oce_ale_pressure_bv_module + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_TRACER + USE o_PARAM + USE o_ARRAYS + use g_config + use g_comm_auto + ! USE diagnostics, only: ldiag_dMOC ! Temporarily commented out due to compilation issues + + implicit none + + private + public :: pressure_bv, pressure_force_4_linfs, pressure_force_4_zxxxx + public :: pressure_force_4_linfs_fullcell, pressure_force_4_linfs_nemo + public :: pressure_force_4_linfs_shchepetkin, pressure_force_4_linfs_easypgf + public :: pressure_force_4_linfs_cubicspline, pressure_force_4_linfs_cavity + public :: pressure_force_4_zxxxx_shchepetkin, pressure_force_4_zxxxx_easypgf + public :: pressure_force_4_zxxxx_cubicspline, init_ref_density, insitu2pot + public :: densityJM_components, density_linear, densityJM_local + public :: sw_alpha_beta, compute_sigma_xy, compute_neutral_slope + public :: ptheta + public :: atg + +contains + +! +! +!=============================================================================== +function ptheta(s,t,p,pr) result(ptheta_result) + ! Compute local potential temperature at pr + ! using bryden 1973 polynomial for adiabatic lapse rate + ! and runge-kutta 4-th order integration algorithm. + ! ref: bryden,h.,1973,deep-sea res.,20,401-408 + ! fofonoff,n.,1977,deep-sea res.,24,489-491 + ! units: + ! pressure p decibars + ! temperature t deg celsius (ipts-68) + ! salinity s (ipss-78) + ! reference prs pr decibars + ! potential tmp. theta deg celsius + ! checkvalue: theta= 36.89073 c,s=40 (ipss-78),t=40 deg c, + ! p=10000 decibars,pr=0 decibars + ! + ! Coded by ?? + ! Reviewed by ?? + !-------------------------------------------------------- + + implicit none + real(kind=WP), intent(in) :: s, t, p, pr + real(kind=WP) :: ptheta_result + real(kind=WP) :: h, xk, q + real(kind=WP) :: t_local, p_local + + t_local = t + p_local = p + h = pr - p_local + xk = h*atg(s,t_local,p_local) + t_local = t_local + 0.5_WP*xk + q = xk + p_local = p_local + 0.5_WP*h + xk = h*atg(s,t_local,p_local) + t_local = t_local + 0.29289322_WP*(xk-q) + q = 0.58578644_WP*xk + 0.121320344_WP*q + xk = h*atg(s,t_local,p_local) + t_local = t_local + 1.707106781_WP*(xk-q) + q = 3.414213562_WP*xk - 4.121320344_WP*q + p_local = p_local + 0.5_WP*h + xk = h*atg(s,t_local,p_local) + ptheta_result = t_local + (xk-2.0_WP*q)/6.0_WP +end function ptheta +! +! +! +!=============================================================================== +function atg(s,t,p) result(atg_result) + ! adiabatic temperature gradient deg c per decibar + ! ref: bryden,h.,1973,deep-sea res.,20,401-408 + ! units: + ! pressure p decibars + ! temperature t deg celsius (ipts-68) + ! salinity s (ipss-78) + ! adiabatic atg deg. c/decibar + ! checkvalue: atg=3.255976e-4 c/dbar for s=40 (ipss-78), + ! t=40 deg c,p0=10000 decibars + ! + ! Coded by ?? + ! Reviewed by ?? + !-------------------------------------------------------- + + implicit none + real(kind=WP), intent(in) :: s, t, p + real(kind=WP) :: atg_result + real(kind=WP) :: ds + + ds = s - 35.0_WP + atg_result = (((-2.1687e-16_WP*t+1.8676e-14_WP)*t-4.6206e-13_WP)*p & + +((2.7759e-12_WP*t-1.1351e-10_WP)*ds+((-5.4481e-14_WP*t & + +8.733e-12_WP)*t-6.7795e-10_WP)*t+1.8741e-8_WP))*p & + +(-4.2393e-8_WP*t+1.8932e-6_WP)*ds & + +((6.6228e-10_WP*t-6.836e-8_WP)*t+8.5258e-6_WP)*t+3.5803e-5_WP +end function atg +! ! ! !=============================================================================== @@ -204,9 +118,7 @@ subroutine pressure_bv(tracers, partit, mesh) USE o_ARRAYS USE g_support USE o_mixing_KPP_mod, only: dbsfc - USE diagnostics, only: ldiag_dMOC - use densityJM_components_interface - use density_linear_interface + ! densityJM_components and density_linear are now in the same module IMPLICIT NONE type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -299,7 +211,7 @@ subroutine pressure_bv(tracers, partit, mesh) !NR and did not vectorize the full loop. !_______________________________________________________________________ ! calculate density for MOC - if (ldiag_dMOC) then + if (.false.) then ! Temporarily disabled: ldiag_dMOC !!PS do nz=1, nl1 do nz=nzmin, nzmax-1 rho(nz) = bulk_0(nz) - 2000._WP*(bulk_pz(nz) -2000._WP*bulk_pz2(nz)) @@ -506,12 +418,7 @@ subroutine pressure_force_4_linfs(tracers, partit, mesh) USE MOD_PARTIT USE MOD_PARSUP use mod_tracer - use pressure_force_4_linfs_fullcell_interface - use pressure_force_4_linfs_nemo_interface - use pressure_force_4_linfs_shchepetkin_interface - use pressure_force_4_linfs_cubicspline_interface - use pressure_force_4_linfs_cavity_interface - use pressure_force_4_linfs_easypgf_interface + ! All pressure_force subroutines are now in the same module implicit none type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -628,8 +535,6 @@ subroutine pressure_force_4_linfs_nemo(tracers, partit, mesh) use MOD_TRACER use o_ARRAYS use g_config - use densityJM_components_interface - use density_linear_interface implicit none type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -1074,8 +979,6 @@ subroutine pressure_force_4_linfs_easypgf(tracers, partit, mesh) use MOD_TRACER use o_ARRAYS use g_config - use densityJM_components_interface - use density_linear_interface implicit none type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -1872,9 +1775,7 @@ subroutine pressure_force_4_zxxxx(tracers, partit, mesh) USE MOD_PARSUP use mod_tracer use g_config - use pressure_force_4_zxxxx_shchepetkin_interface - use pressure_force_4_zxxxx_cubicspline_interface - use pressure_force_4_zxxxx_easypgf_interface + ! All pressure_force subroutines are now in the same module implicit none type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -2104,8 +2005,6 @@ subroutine pressure_force_4_zxxxx_shchepetkin(partit, mesh) USE MOD_PARSUP use o_ARRAYS use g_config - use densityJM_components_interface - use density_linear_interface implicit none type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -2352,8 +2251,6 @@ subroutine pressure_force_4_zxxxx_easypgf(tracers, partit, mesh) use MOD_TRACER use o_ARRAYS use g_config - use densityJM_components_interface - use density_linear_interface implicit none type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -2565,7 +2462,6 @@ SUBROUTINE densityJM_local(t, s, pz, rho_out, partit, mesh) USE MOD_PARSUP !, only: par_ex,pe_status USE o_ARRAYS USE o_PARAM -use densityJM_components_interface IMPLICIT NONE ! @@ -2667,83 +2563,6 @@ end subroutine densityJM_components ! ! !=============================================================================== -function ptheta(s,t,p,pr) - ! Compute local potential temperature at pr - ! using bryden 1973 polynomial for adiabatic lapse rate - ! and runge-kutta 4-th order integration algorithm. - ! ref: bryden,h.,1973,deep-sea res.,20,401-408 - ! fofonoff,n.,1977,deep-sea res.,24,489-491 - ! units: - ! pressure p decibars - ! temperature t deg celsius (ipts-68) - ! salinity s (ipss-78) - ! reference prs pr decibars - ! potential tmp. theta deg celsius - ! checkvalue: theta= 36.89073 c,s=40 (ipss-78),t=40 deg c, - ! p=10000 decibars,pr=0 decibars - ! - ! Coded by ?? - ! Reviewed by ?? - !-------------------------------------------------------- - - use o_param, only: WP - implicit none - real(kind=WP) :: ptheta, s, t, p, pr - real(kind=WP) :: h, xk, q - real(kind=WP), external :: atg - - h = pr - p - xk = h*atg(s,t,p) - t = t + 0.5_WP*xk - q = xk - p = p + 0.5_WP*h - xk = h*atg(s,t,p) - t = t + 0.29289322_WP*(xk-q) - q = 0.58578644_WP*xk + 0.121320344_WP*q - xk = h*atg(s,t,p) - t = t + 1.707106781_WP*(xk-q) - q = 3.414213562_WP*xk - 4.121320344_WP*q - p = p + 0.5_WP*h - xk = h*atg(s,t,p) - ptheta = t + (xk-2.0_WP*q)/6.0_WP - return -end function ptheta -! -! -! -!=============================================================================== -function atg(s,t,p) - ! adiabatic temperature gradient deg c per decibar - ! ref: bryden,h.,1973,deep-sea res.,20,401-408 - ! units: - ! pressure p decibars - ! temperature t deg celsius (ipts-68) - ! salinity s (ipss-78) - ! adiabatic atg deg. c/decibar - ! checkvalue: atg=3.255976e-4 c/dbar for s=40 (ipss-78), - ! t=40 deg c,p0=10000 decibars - ! - ! Coded by ?? - ! Reviewed by ?? - !-------------------------------------------------------- - - use o_param, only: WP - implicit none - real(kind=WP) atg, s, t, p, ds - - ds = s - 35.0_WP - atg = (((-2.1687e-16_WP*t+1.8676e-14_WP)*t-4.6206e-13_WP)*p & - +((2.7759e-12_WP*t-1.1351e-10_WP)*ds+((-5.4481e-14_WP*t & - +8.733e-12_WP)*t-6.7795e-10_WP)*t+1.8741e-8_WP))*p & - +(-4.2393e-8_WP*t+1.8932e-6_WP)*ds & - +((6.6228e-10_WP*t-6.836e-8_WP)*t+8.5258e-6_WP)*t+3.5803e-5_WP - - return -end function atg -! -! -! -!=============================================================================== subroutine sw_alpha_beta(TF1,SF1, partit, mesh) ! DESCRIPTION: ! A function to calculate the thermal expansion coefficient @@ -3011,7 +2830,7 @@ subroutine insitu2pot(tracers, partit, mesh) type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit type(t_tracer), intent(in), target :: tracers - real(kind=WP), external :: ptheta + ! ptheta is now a module function, no need for external declaration real(kind=WP) :: pp, pr, tt, ss integer :: n, nz, nzmin, nzmax real(kind=WP), dimension(:,:), pointer :: temp, salt @@ -3089,7 +2908,6 @@ subroutine init_ref_density(partit, mesh) USE MOD_PARSUP use o_PARAM use o_ARRAYS - use densityJM_components_interface implicit none !___________________________________________________________________________ @@ -3126,3 +2944,5 @@ subroutine init_ref_density(partit, mesh) !$OMP END PARALLEL DO if(mype==0) write(*,*) ' --> compute reference density' end subroutine init_ref_density + +end module oce_ale_pressure_bv_module diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 568f4fe05..a4e1f2840 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -1,133 +1,28 @@ -module diff_part_hor_redi_interface - interface - subroutine diff_part_hor_redi(tracer, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - type(t_tracer), intent(inout), target :: tracer - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine diff_part_hor_redi - end interface -end module diff_part_hor_redi_interface - -module diff_ver_part_expl_ale_interface - interface - subroutine diff_ver_part_expl_ale(tr_num, tracer, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - integer , intent(in) , target :: tr_num - type(t_tracer), intent(inout), target :: tracer - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine diff_ver_part_expl_ale - end interface -end module diff_ver_part_expl_ale_interface - -module diff_ver_part_redi_expl_interface - interface - subroutine diff_ver_part_redi_expl(tracer, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - type(t_tracer), intent(inout), target :: tracer - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine diff_ver_part_redi_expl - end interface -end module diff_ver_part_redi_expl_interface - -module diff_ver_part_impl_ale_interface - interface - subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracer, ice, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use MOD_DYN - use mod_ice - integer , intent(in) , target :: tr_num - type(t_dyn) , intent(inout), target :: dynamics - type(t_tracer), intent(inout), target :: tracer - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - type(t_ice) , intent(in) , target :: ice - end subroutine diff_ver_part_impl_ale - end interface -end module diff_ver_part_impl_ale_interface - -module diff_tracers_ale_interface - interface - subroutine diff_tracers_ale(tr_num, dynamics, tracer, ice, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use mod_ice - use MOD_DYN - integer , intent(in), target :: tr_num - type(t_dyn) , intent(inout), target :: dynamics - type(t_tracer), intent(inout), target :: tracer - type(t_ice), intent(in), target :: ice - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine diff_tracers_ale - end interface -end module diff_tracers_ale_interface - -module bc_surface_interface - interface - function bc_surface(n, id, sval, nzmin, partit, mesh, sst, sss, a_ice) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - integer , intent(in) :: n, id, nzmin - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - real(kind=WP) :: bc_surface - real(kind=WP), intent(in) :: sval, sst, sss, a_ice - end function bc_surface - end interface -end module bc_surface_interface - -module diff_part_bh_interface - interface - subroutine diff_part_bh(tr_num, dynamics, tracer, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use MOD_DYN - integer , intent(in) , target :: tr_num - type(t_dyn) , intent(inout), target :: dynamics - type(t_tracer), intent(inout), target :: tracer - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine diff_part_bh - end interface -end module diff_part_bh_interface - -module solve_tracers_ale_interface - interface - subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use MOD_DYN - USE MOD_ICE - type(t_ice) , intent(in), target :: ice - type(t_dyn) , intent(inout), target :: dynamics - type(t_tracer), intent(inout), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine solve_tracers_ale - end interface -end module solve_tracers_ale_interface +module oce_ale_tracer_module + use mod_mesh + USE MOD_PARTIT + USE MOD_PARSUP + use mod_tracer + USE MOD_DYN + use mod_ice + USE o_PARAM + USE o_ARRAYS + USE g_CONFIG + use g_comm_auto + use g_forcing_arrays + use g_forcing_param + use diagnostics + use mod_transit + + implicit none + + private + public :: solve_tracers_ale, diff_tracers_ale, diff_ver_part_expl_ale, & + diff_ver_part_impl_ale, diff_ver_part_redi_expl, diff_part_hor_redi, & + diff_part_bh, bc_surface + +contains + ! ! !=============================================================================== @@ -147,8 +42,8 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) use Toy_Channel_Dbgyre use o_ARRAYS, only: heat_flux use g_forcing_arrays, only: sw_3d - use diff_tracers_ale_interface - use oce_adv_tra_driver_interfaces + ! diff_tracers_ale is now in the same module + use oce_adv_tra_driver_module, only: do_oce_adv_tra #if defined(__recom) use recom_glovar use recom_config @@ -334,11 +229,7 @@ subroutine diff_tracers_ale(tr_num, dynamics, tracers, ice, partit, mesh) use MOD_DYN use o_arrays use o_tracers - use diff_part_hor_redi_interface - use diff_ver_part_expl_ale_interface - use diff_ver_part_redi_expl_interface - use diff_ver_part_impl_ale_interface - use diff_part_bh_interface + ! All these subroutines are now in the same module #if defined(__recom) use ver_sinking_recom_interface use diff_ver_recom_expl_interface @@ -567,7 +458,7 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) use g_forcing_arrays use o_mixing_KPP_mod !for ghats _GO_ use g_cvmix_kpp, only: kpp_nonlcltranspT, kpp_nonlcltranspS, kpp_oblmixc - use bc_surface_interface + ! bc_surface is now in the same module use mod_ice use iceberg_params implicit none @@ -1668,3 +1559,4 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit, mesh, sst, sss, aice) RETURN end function bc_surface +end module oce_ale_tracer_module diff --git a/src/oce_ale_vel_rhs.F90 b/src/oce_ale_vel_rhs.F90 index 2cf0ab2b0..ece5c2506 100644 --- a/src/oce_ale_vel_rhs.F90 +++ b/src/oce_ale_vel_rhs.F90 @@ -1,33 +1,23 @@ +module oce_ale_vel_rhs_module + USE MOD_ICE + USE MOD_DYN + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_MESH + use o_ARRAYS + use o_PARAM + use g_CONFIG + use g_forcing_param + use g_forcing_arrays + use g_comm_auto + use g_sbf + + implicit none + + private + public :: compute_vel_rhs, momentum_adv_scalar -module compute_vel_rhs_interface - interface - subroutine compute_vel_rhs(ice, dynamics, partit, mesh) - USE MOD_ICE - USE MOD_DYN - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_MESH - type(t_ice) , intent(inout), target :: ice - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine compute_vel_rhs - end interface -end module compute_vel_rhs_interface - -module momentum_adv_scalar_interface - interface - subroutine momentum_adv_scalar(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine momentum_adv_scalar - end interface -end module momentum_adv_scalar_interface +contains ! ! @@ -45,7 +35,6 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) use g_forcing_arrays, only: press_air use g_comm_auto use g_sbf, only: l_mslp - use momentum_adv_scalar_interface use momentum_adv_scalar_transpv_interface implicit none type(t_ice) , intent(inout), target :: ice @@ -588,6 +577,5 @@ subroutine momentum_adv_scalar(dynamics, partit, mesh) !$OMP END PARALLEL end subroutine momentum_adv_scalar - -! =================================================================== +end module oce_ale_vel_rhs_module diff --git a/src/oce_dyn.F90 b/src/oce_dyn.F90 index 0d20fbd80..40ac6d362 100755 --- a/src/oce_dyn.F90 +++ b/src/oce_dyn.F90 @@ -1,3 +1,20 @@ +module oce_dyn_module + use mod_mesh + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + USE o_PARAM + USE g_CONFIG + use g_comm_auto + + implicit none + + private + public :: update_vel, compute_vel_nodes, viscosity_filter, visc_filt_bcksct, & + visc_filt_bilapl, visc_filt_bidiff, compute_ke_wrho, compute_apegen, & + check_viscopt, check_validviscopt_5 + +contains ! A set of routines for computing the horizonlal viscosity ! the control parameters (their default values) are: @@ -9,76 +26,12 @@ ! We however, try to keep dynamics%visc_gamma1<0.1 ! 3. dynamics%visc_gamma2 is dimensional (1/velocity). If it is 10, then the respective term dominates starting from |u|=0.1 m/s an so on. It is only used in: ! (5) visc_filt_bcksct, (6) visc_filt_bilapl, (7) visc_filt_bidiff -module visc_filt_bcksct_interface - interface - subroutine visc_filt_bcksct(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - - end subroutine visc_filt_bcksct - end interface -end module visc_filt_bcksct_interface - -module visc_filt_bilapl_interface - interface - subroutine visc_filt_bilapl(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - - end subroutine visc_filt_bilapl - end interface -end module visc_filt_bilapl_interface - -module visc_filt_bidiff_interface - interface - subroutine visc_filt_bidiff(dynamics, partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - - end subroutine visc_filt_bidiff - end interface -end module visc_filt_bidiff_interface - -module check_validviscopt_interface - interface - subroutine check_validviscopt_5(partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine check_validviscopt_5 - end interface -end module check_validviscopt_interface ! ! Contains routines needed for computations of dynamics. ! includes: update_vel, compute_vel_nodes !_______________________________________________________________________________ SUBROUTINE update_vel(dynamics, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - USE o_PARAM - USE g_CONFIG - use g_comm_auto - IMPLICIT NONE type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit type(t_mesh) , intent(in) , target :: mesh @@ -175,13 +128,6 @@ end subroutine update_vel ! !_______________________________________________________________________________ subroutine compute_vel_nodes(dynamics, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - USE o_PARAM - use g_comm_auto - IMPLICIT NONE type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit type(t_mesh) , intent(in) , target :: mesh @@ -233,9 +179,7 @@ subroutine viscosity_filter(option, dynamics, partit, mesh) USE MOD_PARTIT USE MOD_PARSUP use MOD_DYN - use visc_filt_bcksct_interface - use visc_filt_bilapl_interface - use visc_filt_bidiff_interface + ! visc_filt subroutines are now in the same module use g_backscatter IMPLICIT NONE integer :: option @@ -859,7 +803,7 @@ subroutine check_viscopt(dynamics, partit, mesh) USE MOD_PARTIT USE MOD_PARSUP USE MOD_MESH - USE check_validviscopt_interface + ! check_validviscopt_5 is now in the same module IMPLICIT NONE type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit @@ -991,3 +935,5 @@ subroutine check_validviscopt_5(partit, mesh) end if end subroutine check_validviscopt_5 + +end module oce_dyn_module diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 7c55b3dbc..a1e0b316e 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1,137 +1,19 @@ -module read_mesh_interface - interface - subroutine read_mesh(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine read_mesh - end interface -end module read_mesh_interface -module find_levels_interface - interface - subroutine find_levels(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine find_levels - end interface -end module find_levels_interface -module find_levels_cavity_interface - interface - subroutine find_levels_cavity(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine find_levels_cavity - end interface -end module find_levels_cavity_interface -module test_tri_interface - interface - subroutine test_tri(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine test_tri - end interface -end module test_tri_interface -module load_edges_interface - interface - subroutine load_edges(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine load_edges - end interface -end module load_edges_interface -module find_neighbors_interface - interface - subroutine find_neighbors(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine find_neighbors - end interface -end module find_neighbors_interface -module mesh_areas_interface - interface - subroutine mesh_areas(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine mesh_areas - end interface -end module mesh_areas_interface -module elem_center_interface - interface - subroutine elem_center(elem, x, y, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - integer :: elem - real(kind=WP), intent(inout) :: x, y - type(t_mesh), intent(inout), target :: mesh - end subroutine elem_center - end interface -end module elem_center_interface -module edge_center_interface - interface - subroutine edge_center(n1, n2, x, y, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - integer :: n1, n2 - real(kind=WP), intent(inout):: x, y - type(t_mesh), intent(inout), target :: mesh - end subroutine edge_center - end interface -end module edge_center_interface -module mesh_auxiliary_arrays_interface - interface - subroutine mesh_auxiliary_arrays(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine mesh_auxiliary_arrays - end interface -end module mesh_auxiliary_arrays_interface -module find_levels_min_e2n_interface - interface - subroutine find_levels_min_e2n(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine find_levels_min_e2n - end interface -end module find_levels_min_e2n_interface -module check_total_volume_interface - interface - subroutine check_total_volume(partit, mesh) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - type(t_mesh), intent(inout), target :: mesh - type(t_partit), intent(inout), target :: partit - end subroutine check_total_volume - end interface -end module check_total_volume_interface +module oce_mesh_module + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE g_config, only: flag_debug + USE g_ROTATE_grid + use par_support_interfaces + + implicit none + + private + public :: mesh_setup, read_mesh, find_levels, find_levels_cavity, test_tri, load_edges, & + find_neighbors, mesh_areas, elem_center, edge_center, mesh_auxiliary_arrays, & + find_levels_min_e2n, check_total_volume, check_mesh_consistency + +contains ! Driving routine. The distributed mesh information and mesh proper ! are read from files. @@ -140,22 +22,6 @@ end module check_total_volume_interface ! Array sizes vary (sometimes we need only myDim, yet sometimes more)! ! S. Danilov, 2012 SUBROUTINE mesh_setup(partit, mesh) -USE MOD_MESH -USE MOD_PARTIT -USE MOD_PARSUP -USE g_config, only: flag_debug -USE g_ROTATE_grid -use read_mesh_interface -use find_levels_interface -use find_levels_cavity_interface -use mesh_auxiliary_arrays_interface -use test_tri_interface -use load_edges_interface -use find_levels_min_e2n_interface -use find_neighbors_interface -use mesh_areas_interface -use par_support_interfaces -IMPLICIT NONE type(t_mesh), intent(inout) :: mesh type(t_partit), intent(inout), target :: partit @@ -1976,7 +1842,7 @@ SUBROUTINE find_neighbors(partit, mesh) USE MOD_PARSUP USE g_ROTATE_grid use g_comm_auto -use elem_center_interface +! elem_center is now in the same module implicit none type(t_mesh), intent(inout), target :: mesh type(t_partit), intent(inout), target :: partit @@ -2120,7 +1986,7 @@ subroutine edge_center(n1, n2, x, y, mesh) implicit none integer :: n1, n2 ! nodes of the edge real(kind=WP), intent(inout) :: x, y -type(t_mesh), intent(inout), target :: mesh +type(t_mesh), intent(in), target :: mesh real(kind=WP) :: a(2), b(2) a=mesh%coord_nod2D(:,n1) @@ -2138,7 +2004,7 @@ subroutine elem_center(elem, x, y, mesh) USE g_CONFIG implicit none real(kind=WP), intent(inout) :: x, y -type(t_mesh), intent(inout), target :: mesh +type(t_mesh), intent(in), target :: mesh integer :: elem, elnodes(3), k real(kind=WP) :: ax(3), amin @@ -2438,8 +2304,7 @@ SUBROUTINE mesh_auxiliary_arrays(partit, mesh) USE g_CONFIG, only: rotated_grid, force_rotation USE g_ROTATE_grid use g_comm_auto -use elem_center_interface -use edge_center_interface +! elem_center and edge_center are now in the same module IMPLICIT NONE integer :: n,j,q, elnodes(3), ed(2), elem, el(2), elnodes_(3),node @@ -2892,6 +2757,5 @@ subroutine check_total_volume(partit, mesh) end if end subroutine check_total_volume -! -! -!_______________________________________________________________________________ + +end module oce_mesh_module diff --git a/src/oce_muscl_adv.F90 b/src/oce_muscl_adv.F90 index cb3c0638a..5da469eeb 100755 --- a/src/oce_muscl_adv.F90 +++ b/src/oce_muscl_adv.F90 @@ -1,16 +1,19 @@ -module find_up_downwind_triangles_interface - interface - subroutine find_up_downwind_triangles(twork, partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - use MOD_TRACER - type(t_mesh), intent(in) , target :: mesh - type(t_partit), intent(inout), target :: partit - type(t_tracer_work), intent(inout), target :: twork - end subroutine find_up_downwind_triangles - end interface -end module find_up_downwind_triangles_interface +module oce_muscl_adv_module + use MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + use MOD_TRACER + use o_ARRAYS + use o_PARAM + use g_comm_auto + use g_config + + implicit none + + private + public :: muscl_adv_init, find_up_downwind_triangles, fill_up_dn_grad + +contains ! A set of routines to implement MUSCL-type of advection ! For description, see Abalakin, I., Dervieux, A., Kozubskaya, T., 2002. A @@ -38,7 +41,7 @@ subroutine muscl_adv_init(twork, partit, mesh) use o_PARAM use g_comm_auto use g_config - use find_up_downwind_triangles_interface + ! find_up_downwind_triangles is now in the same module IMPLICIT NONE integer :: n, k, n1, n2 @@ -523,3 +526,5 @@ SUBROUTINE fill_up_dn_grad(twork, partit, mesh) !$OMP END DO !$OMP END PARALLEL END SUBROUTINE fill_up_dn_grad + +end module oce_muscl_adv_module diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index ec4f93d72..9e7cb7c15 100755 --- a/src/oce_setup_step.F90 +++ b/src/oce_setup_step.F90 @@ -1,80 +1,4 @@ -module oce_initial_state_interface - interface - subroutine oce_initial_state(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - type(t_tracer), intent(inout), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in) , target :: mesh - end subroutine oce_initial_state - end interface -end module oce_initial_state_interface - -module tracer_init_interface - interface - subroutine tracer_init(tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - type(t_tracer), intent(inout), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(in) , target :: mesh - end subroutine tracer_init - end interface -end module tracer_init_interface - -module dynamics_init_interface - interface - subroutine dynamics_init(dynamics, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - use MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine dynamics_init - end interface -end module dynamics_init_interface - -module ocean_setup_interface - interface - subroutine ocean_setup(dynamics, tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_tracer), intent(inout), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(inout) , target :: mesh - end subroutine ocean_setup - end interface -end module ocean_setup_interface - -module before_oce_step_interface - interface - subroutine before_oce_step(dynamics, tracers, partit, mesh) - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - use mod_tracer - use MOD_DYN - type(t_dyn) , intent(inout), target :: dynamics - type(t_tracer), intent(inout), target :: tracers - type(t_partit), intent(inout), target :: partit - type(t_mesh) , intent(in) , target :: mesh - end subroutine before_oce_step - end interface -end module before_oce_step_interface -! -! -!_______________________________________________________________________________ -subroutine ocean_setup(dynamics, tracers, partit, mesh) +module oce_setup_step_module USE MOD_MESH USE MOD_PARTIT USE MOD_PARSUP @@ -92,11 +16,22 @@ subroutine ocean_setup(dynamics, tracers, partit, mesh) use g_backscatter use Toy_Channel_Soufflet use Toy_Channel_Dbgyre - use oce_initial_state_interface - use oce_adv_tra_fct_interfaces - use init_ale_interface - use init_thickness_ale_interface - IMPLICIT NONE + use oce_ale_module, only: init_ale, init_thickness_ale, init_stiff_mat_ale + use oce_adv_tra_fct_module, only: oce_adv_tra_fct_init + use oce_muscl_adv_module, only: muscl_adv_init + use oce_ale_pressure_bv_module, only: init_ref_density + + implicit none + + private + public :: oce_initial_state, tracer_init, dynamics_init, ocean_setup, before_oce_step, arrays_init + +contains + +!=============================================================================== +! Ocean setup subroutine +!=============================================================================== +subroutine ocean_setup(dynamics, tracers, partit, mesh) type(t_dyn) , intent(inout), target :: dynamics type(t_tracer), intent(inout), target :: tracers type(t_partit), intent(inout), target :: partit @@ -1399,3 +1334,5 @@ SUBROUTINE before_oce_step(dynamics, tracers, partit, mesh) END SELECT end if END SUBROUTINE before_oce_step + +end module oce_setup_step_module diff --git a/src/oce_tracer_mod.F90 b/src/oce_tracer_mod.F90 index 2fb963fe0..2b78dbc26 100755 --- a/src/oce_tracer_mod.F90 +++ b/src/oce_tracer_mod.F90 @@ -4,6 +4,7 @@ MODULE o_tracers USE MOD_TRACER USE MOD_PARTIT USE MOD_PARSUP +use oce_muscl_adv_module, only: fill_up_dn_grad IMPLICIT NONE CONTAINS diff --git a/src/solver.F90 b/src/solver.F90 index 2a1a95af1..eb9e6ee4c 100644 --- a/src/solver.F90 +++ b/src/solver.F90 @@ -1,32 +1,17 @@ -module ssh_solve_preconditioner_interface - interface - subroutine ssh_solve_preconditioner(solverinfo, partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_solverinfo), intent(inout), target :: solverinfo - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(inout), target :: mesh - end subroutine ssh_solve_preconditioner - end interface -end module ssh_solve_preconditioner_interface +module solver_module + use MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + USE g_comm_auto + + implicit none + + private + public :: ssh_solve_preconditioner, ssh_solve_cg + +contains -module ssh_solve_cg_interface - interface - subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - type(t_solverinfo), intent(inout), target :: solverinfo - type(t_partit), intent(inout), target :: partit - type(t_mesh), intent(inout), target :: mesh - real(kind=WP), intent(inout) :: x(partit%myDim_nod2D+partit%eDim_nod2D) - real(kind=WP), intent(in) :: rhs(partit%myDim_nod2D+partit%eDim_nod2D) - end subroutine ssh_solve_cg - end interface -end module ssh_solve_cg_interface !========================================================================= subroutine ssh_solve_preconditioner(solverinfo, partit, mesh) ! Preconditioner follows MITgcm (JGR, 102,5753-5766, 1997) @@ -40,11 +25,6 @@ subroutine ssh_solve_preconditioner(solverinfo, partit, mesh) ! paper cited) is, in reality, one iteration of the ! Jacobi method, with symmetrization. We need symmetrization to be able to use ! the conjugate gradient method. - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - USE g_comm_auto IMPLICIT NONE type(t_solverinfo), intent(inout), target :: solverinfo type(t_partit), intent(inout), target :: partit @@ -102,11 +82,6 @@ subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) ! ! I tried first to follow the MITgcm paper, but I have doubts about ! their computations of beta. The variant below -- see Wikipedia. - USE MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE MOD_DYN - USE g_comm_auto IMPLICIT NONE type(t_solverinfo), intent(inout), target :: solverinfo type(t_partit), intent(inout), target :: partit @@ -278,7 +253,7 @@ subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) ! At the end: The result is in X, but it needs a halo exchange. call exchange_nod(x, partit) !$OMP BARRIER -end subroutine ssh_solve_cg +end subroutine ssh_solve_cg -! =================================================================== +end module solver_module