From e741b38c961767ca80067fe5785b929ce5ef9ee1 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 29 Jan 2021 14:20:54 +0100 Subject: [PATCH 01/40] try to fix cavity related bug when computing the scalar area when cavity is used --- src/MOD_MESH.F90 | 2 +- src/associate_mesh.h | 4 +++- src/cavity_param.F90 | 17 ++++++++++++++-- src/gen_modules_diag.F90 | 3 ++- src/gen_support.F90 | 3 ++- src/ice_fct.F90 | 13 ++++++++++-- src/ice_oce_coupling.F90 | 21 +++++++++++++++----- src/oce_ale.F90 | 16 ++++++++++----- src/oce_ale_tracer.F90 | 34 +++++++++++++++---------------- src/oce_mesh.F90 | 43 ++++++++++++++++++++++++++++++++-------- src/oce_modules.F90 | 1 + src/oce_setup_step.F90 | 8 ++++++-- 12 files changed, 120 insertions(+), 45 deletions(-) diff --git a/src/MOD_MESH.F90 b/src/MOD_MESH.F90 index 53f45e1bd..795242a8d 100644 --- a/src/MOD_MESH.F90 +++ b/src/MOD_MESH.F90 @@ -18,7 +18,7 @@ MODULE MOD_MESH TYPE T_MESH integer :: nod2D ! the number of 2D nodes -real(kind=WP) :: ocean_area +real(kind=WP) :: ocean_area, ocean_areawithcav real(kind=WP), allocatable, dimension(:,:) :: coord_nod2D, geo_coord_nod2D integer :: edge2D ! the number of 2D edges integer :: edge2D_in ! the number of internal 2D edges diff --git a/src/associate_mesh.h b/src/associate_mesh.h index 554819e88..30611d44b 100644 --- a/src/associate_mesh.h +++ b/src/associate_mesh.h @@ -3,6 +3,7 @@ integer , pointer :: elem2D integer , pointer :: edge2D integer , pointer :: edge2D_in real(kind=WP) , pointer :: ocean_area +real(kind=WP) , pointer :: ocean_areawithcav integer , pointer :: nl real(kind=WP), dimension(:,:), pointer :: coord_nod2D, geo_coord_nod2D integer, dimension(:,:) , pointer :: elem2D_nodes @@ -35,7 +36,8 @@ nod2D => mesh%nod2D elem2D => mesh%elem2D edge2D => mesh%edge2D edge2D_in => mesh%edge2D_in -ocean_area => mesh%ocean_area +ocean_area => mesh%ocean_area +ocean_areawithcav => mesh%ocean_areawithcav nl => mesh%nl !!$coord_nod2D => mesh%coord_nod2D diff --git a/src/cavity_param.F90 b/src/cavity_param.F90 index ba6e2270c..68f1735e9 100644 --- a/src/cavity_param.F90 +++ b/src/cavity_param.F90 @@ -356,14 +356,14 @@ end subroutine cavity_heat_water_fluxes_2eq subroutine cavity_momentum_fluxes(mesh) use MOD_MESH use o_PARAM , only: density_0, C_d, WP - use o_ARRAYS, only: UV, stress_surf + use o_ARRAYS, only: UV, Unode, stress_surf, stress_node_surf use i_ARRAYS, only: u_w, v_w use g_PARSUP implicit none !___________________________________________________________________________ type(t_mesh), intent(inout) , target :: mesh - integer :: elem, elnodes(3), nzmin + integer :: elem, elnodes(3), nzmin, node real(kind=WP) :: aux #include "associate_mesh.h" @@ -381,6 +381,19 @@ subroutine cavity_momentum_fluxes(mesh) stress_surf(1,elem)=-aux*UV(1,nzmin,elem) stress_surf(2,elem)=-aux*UV(2,nzmin,elem) end do + + !___________________________________________________________________________ + do node=1,myDim_nod2D+eDim_nod2D + nzmin = ulevels_nod2d(node) + if(nzmin==1) cycle + + ! momentum stress: + ! need to check the sensitivity to the drag coefficient + ! here I use the bottom stress coefficient, which is 3e-3, for this FO2 work. + aux=sqrt(Unode(1,nzmin,node)**2+Unode(2,nzmin,node)**2)*density_0*C_d + stress_node_surf(1,node)=-aux*Unode(1,nzmin,node) + stress_node_surf(2,node)=-aux*Unode(2,nzmin,node) + end do end subroutine cavity_momentum_fluxes ! ! diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index 63fd0c98e..c1e1ec3d8 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -134,7 +134,8 @@ subroutine diag_curl_stress_surf(mode, mesh) end if END DO DO n=1, myDim_nod2D+eDim_nod2D - curl_stress_surf(n)=curl_stress_surf(n)/area(1,n) + !!PS curl_stress_surf(n)=curl_stress_surf(n)/area(1,n) + curl_stress_surf(n)=curl_stress_surf(n)/area(ulevels_nod2D(n),n) END DO end subroutine diag_curl_stress_surf ! ============================================================== diff --git a/src/gen_support.F90 b/src/gen_support.F90 index df1af60bc..eee45c814 100644 --- a/src/gen_support.F90 +++ b/src/gen_support.F90 @@ -272,7 +272,8 @@ subroutine integrate_nod_2D(data, int2D, mesh) #include "associate_mesh.h" lval=0.0_WP do row=1, myDim_nod2D - lval=lval+data(row)*area(1,row) + !!PS lval=lval+data(row)*area(1,row) + lval=lval+data(row)*area(ulevels_nod2D(row),row) end do int2D=0.0_WP diff --git a/src/ice_fct.F90 b/src/ice_fct.F90 index dcee2c15c..e82f2e11e 100755 --- a/src/ice_fct.F90 +++ b/src/ice_fct.F90 @@ -663,6 +663,10 @@ SUBROUTINE ice_mass_matrix_fill(mesh) do n=1,3 row=elnodes(n) if(row>myDim_nod2D) cycle + !___________________________________________________________________ + ! if node is cavity cycle over + if(ulevels_nod2d(row)>1) cycle + ! Global-to-local neighbourhood correspondence DO q=1,nn_num(row) col_pos(nn_pos(q,row))=q @@ -670,6 +674,10 @@ SUBROUTINE ice_mass_matrix_fill(mesh) offset=ssh_stiff%rowptr(row)-ssh_stiff%rowptr(1) DO q=1,3 col=elnodes(q) + !___________________________________________________________________ + ! if node is cavity cycle over + if(ulevels_nod2d(col)>1) cycle + ipos=offset+col_pos(col) mass_matrix(ipos)=mass_matrix(ipos)+elem_area(elem)/12.0_WP if(q==n) then @@ -689,7 +697,8 @@ SUBROUTINE ice_mass_matrix_fill(mesh) offset=ssh_stiff%rowptr(q)-ssh_stiff%rowptr(1)+1 n=ssh_stiff%rowptr(q+1)-ssh_stiff%rowptr(1) aa=sum(mass_matrix(offset:n)) - if(abs(area(1,q)-aa)>.1_WP) then + !!PS if(abs(area(1,q)-aa)>.1_WP) then + if(abs(area(ulevels_nod2d(q),q)-aa)>.1_WP) then iflag=q flag=1 endif @@ -698,7 +707,7 @@ SUBROUTINE ice_mass_matrix_fill(mesh) offset=ssh_stiff%rowptr(iflag)-ssh_stiff%rowptr(1)+1 n=ssh_stiff%rowptr(iflag+1)-ssh_stiff%rowptr(1) aa=sum(mass_matrix(offset:n)) - write(*,*) '#### MASS MATRIX PROBLEM', mype, iflag, aa, area(1,iflag) + write(*,*) '#### MASS MATRIX PROBLEM', mype, iflag, aa, area(1,iflag), ulevels_nod2D(iflag) endif deallocate(col_pos) END SUBROUTINE ice_mass_matrix_fill diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 14c1fe29d..15db83c06 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -38,6 +38,10 @@ subroutine oce_fluxes_mom(mesh) stress_iceoce_x(n)=0.0_WP stress_iceoce_y(n)=0.0_WP end if + + ! total surface stress (iceoce+atmoce) on nodes + stress_node_surf(1,n) = stress_iceoce_x(n)*a_ice(n) + stress_atmoce_x(n)*(1.0_WP-a_ice(n)) + stress_node_surf(2,n) = stress_iceoce_y(n)*a_ice(n) + stress_atmoce_y(n)*(1.0_WP-a_ice(n)) end do !___________________________________________________________________________ @@ -48,10 +52,12 @@ subroutine oce_fluxes_mom(mesh) !_______________________________________________________________________ elnodes=elem2D_nodes(:,elem) - stress_surf(1,elem)=sum(stress_iceoce_x(elnodes)*a_ice(elnodes) + & - stress_atmoce_x(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP - stress_surf(2,elem)=sum(stress_iceoce_y(elnodes)*a_ice(elnodes) + & - stress_atmoce_y(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP + !!PS stress_surf(1,elem)=sum(stress_iceoce_x(elnodes)*a_ice(elnodes) + & + !!PS stress_atmoce_x(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP + !!PS stress_surf(2,elem)=sum(stress_iceoce_y(elnodes)*a_ice(elnodes) + & + !!PS stress_atmoce_y(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP + stress_surf(1,elem)=sum(stress_node_surf(1,elnodes))/3.0_WP + stress_surf(2,elem)=sum(stress_node_surf(2,elnodes))/3.0_WP END DO !___________________________________________________________________________ @@ -177,6 +183,9 @@ subroutine oce_fluxes(mesh) if (use_cavity) call cavity_heat_water_fluxes_3eq(mesh) !!PS if (use_cavity) call cavity_heat_water_fluxes_2eq(mesh) +!!PS where(ulevels_nod2D>1) heat_flux=0.0_WP +!!PS where(ulevels_nod2D>1) water_flux=0.0_WP + !___________________________________________________________________________ call exchange_nod(heat_flux, water_flux) @@ -268,7 +277,9 @@ subroutine oce_fluxes(mesh) ! here the + sign must be used because we switched up the sign of the ! water_flux with water_flux = -fresh_wa_flux, but evap, prec_... and runoff still ! have there original sign - water_flux=water_flux+net/ocean_area + ! if use_cavity=.false. --> ocean_area == ocean_areawithcav + !! water_flux=water_flux+net/ocean_area + water_flux=water_flux+net/ocean_areawithcav !___________________________________________________________________________ if (use_sw_pene) call cal_shortwave_rad(mesh) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 44cd3ea86..174697209 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -572,7 +572,10 @@ subroutine init_thickness_ale(mesh) end if !___________________________________________________________________________ ! Fill in ssh_rhs_old - ssh_rhs_old=(hbar-hbar_old)*area(1,:)/dt + !!PS ssh_rhs_old=(hbar-hbar_old)*area(1,:)/dt + do n=1,myDim_nod2D+eDim_nod2D + ssh_rhs_old(n)=(hbar(n)-hbar_old(n))*area(ulevels_nod2D(n),n)/dt + end do ! -->see equation (14) FESOM2:from finite elements to finie volume eta_n=alpha*hbar_old+(1.0_WP-alpha)*hbar @@ -1258,7 +1261,8 @@ subroutine init_stiff_mat_ale(mesh) ! Mass matrix part do row=1, myDim_nod2D offset = ssh_stiff%rowptr(row) - SSH_stiff%values(offset) = SSH_stiff%values(offset)+ area(1,row)/dt + !!PS SSH_stiff%values(offset) = SSH_stiff%values(offset)+ area(1,row)/dt + SSH_stiff%values(offset) = SSH_stiff%values(offset)+ area(ulevels_nod2D(row),row)/dt end do deallocate(n_pos,n_num) @@ -2550,9 +2554,11 @@ subroutine oce_timestep_ale(n, mesh) t0=MPI_Wtime() -!!PS water_flux = 0.0_WP -!!PS heat_flux = 0.0_WP -!!PS stress_surf= 0.0_WP + water_flux = 0.0_WP + heat_flux = 0.0_WP + stress_surf= 0.0_WP + stress_node_surf= 0.0_WP + !___________________________________________________________________________ ! calculate equation of state, density, pressure and mixed layer depths if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call pressure_bv'//achar(27)//'[0m' diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 34e2ff27d..24f77d628 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -613,13 +613,13 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) !_______________________________________________________________________ ! case of activated shortwave penetration into the ocean, ad 3d contribution - if (use_sw_pene .and. tr_num==1) then - !!PS do nz=1, nzmax-1 - do nz=nzmin, nzmax-1 - zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! - tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n)*area(nz+1,n)/area(nz,n))*zinv - end do - end if +!!PS if (use_sw_pene .and. tr_num==1) then +!!PS !!PS do nz=1, nzmax-1 +!!PS do nz=nzmin, nzmax-1 +!!PS zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! +!!PS tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n)*area(nz+1,n)/area(nz,n))*zinv +!!PS end do +!!PS end if !_______________________________________________________________________ ! The first row contains also the boundary condition from heatflux, @@ -634,7 +634,7 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! v (+) v (+) ! !!PS tr(1)= tr(1)+bc_surface(n, tracer_id(tr_num)) - tr(nzmin)= tr(nzmin)+bc_surface(n, tracer_id(tr_num),mesh) +!!PS tr(nzmin)= tr(nzmin)+bc_surface(n, tracer_id(tr_num),mesh) !_______________________________________________________________________ ! The forward sweep algorithm to solve the three-diagonal matrix @@ -719,15 +719,15 @@ subroutine diff_ver_part_redi_expl(mesh) Tx=0.0_WP Ty=0.0_WP do k=1, nod_in_elem2D_num(n) - elem=nod_in_elem2D(k,n) - !!PS if(nz.LE.(nlevels(elem)-1)) then - if( nz.LE.(nlevels(elem)-1) .and. nz.GE.(ulevels(elem))) then - Tx=Tx+tr_xy(1,nz,elem)*elem_area(elem) - Ty=Ty+tr_xy(2,nz,elem)*elem_area(elem) - endif - end do - tr_xynodes(1,nz,n)=tx/3.0_WP/area(nz,n) - tr_xynodes(2,nz,n)=ty/3.0_WP/area(nz,n) + elem=nod_in_elem2D(k,n) + !!PS if(nz.LE.(nlevels(elem)-1)) then + if( nz.LE.(nlevels(elem)-1) .and. nz.GE.(ulevels(elem))) then + Tx=Tx+tr_xy(1,nz,elem)*elem_area(elem) + Ty=Ty+tr_xy(2,nz,elem)*elem_area(elem) + endif + end do + tr_xynodes(1,nz,n)=tx/3.0_WP/area(nz,n) + tr_xynodes(2,nz,n)=ty/3.0_WP/area(nz,n) end do end do diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index e641b4de3..763c31020 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1811,8 +1811,8 @@ SUBROUTINE mesh_areas(mesh) ! area(nl, myDim_nod2D) -integer :: n,j,q, elnodes(3), ed(2), elem, nz -real(kind=WP) :: a(2), b(2), ax, ay, lon, lat, vol +integer :: n,j,q, elnodes(3), ed(2), elem, nz,nzmin +real(kind=WP) :: a(2), b(2), ax, ay, lon, lat, vol, vol2 real(kind=WP), allocatable,dimension(:) :: work_array real(kind=WP) :: t0, t1 type(t_mesh), intent(inout), target :: mesh @@ -1855,13 +1855,12 @@ SUBROUTINE mesh_areas(mesh) DO n=1, myDim_nod2D DO j=1,mesh%nod_in_elem2D_num(n) elem=mesh%nod_in_elem2D(j,n) - !!PS DO nz=mesh%ulevels(elem),mesh%nlevels(elem)-1 - DO nz=1,mesh%nlevels(elem)-1 + DO nz=mesh%ulevels(elem),mesh%nlevels(elem)-1 + !!PS DO nz=1,mesh%nlevels(elem)-1 mesh%area(nz,n)=mesh%area(nz,n)+mesh%elem_area(elem)/3.0_WP END DO END DO END DO - ! Only areas through which there is exchange are counted ! =========== @@ -1870,7 +1869,17 @@ SUBROUTINE mesh_areas(mesh) mesh%elem_area=mesh%elem_area*r_earth*r_earth mesh%area=mesh%area*r_earth*r_earth - call exchange_nod(mesh%area) +call exchange_nod(mesh%area) + + +!!PS do n=1, myDim_nod2D +!!PS nzmin = mesh%ulevels_nod2d(n) +!!PS if (nzmin>1) then +!!PS write(*,*) ' --> mesh area:', mype, n, nzmin, mesh%area(nzmin,n),mesh%area(nzmin+1,n),mesh%area(nzmin+2,n) +!!PS end if +!!PS end do + + do n=1,myDim_nod2d+eDim_nod2D do nz=1,mesh%nl @@ -1884,9 +1893,19 @@ SUBROUTINE mesh_areas(mesh) ! coordinates are in radians, edge_dxdy are in meters, ! and areas are in m^2 +!!PS do n=1,myDim_nod2d+eDim_nod2D +!!PS mesh%area_inv(1:mesh%ulevels_nod2D(n)-1,n) = 0.0_WP +!!PS mesh%area(1:mesh%ulevels_nod2D(n)-1,n) = 0.0_WP +!!PS end do + + allocate(work_array(myDim_nod2D)) - mesh%mesh_resolution=sqrt(mesh%area(1, :)/pi)*2._WP + !!PS mesh%mesh_resolution=sqrt(mesh%area(1, :)/pi)*2._WP + do n=1,myDim_nod2d+eDim_nod2D + mesh%mesh_resolution(n)=sqrt(mesh%area(mesh%ulevels_nod2D(n), n)/pi)*2._WP + end do + DO q=1, 3 !apply mass matrix N times to smooth the field DO n=1, myDim_nod2D vol=0._WP @@ -1907,12 +1926,18 @@ SUBROUTINE mesh_areas(mesh) deallocate(work_array) vol=0.0_WP + vol2=0.0_WP do n=1, myDim_nod2D + vol2=vol2+mesh%area(mesh%ulevels_nod2D(n), n) + if (mesh%ulevels_nod2D(n)>1) cycle vol=vol+mesh%area(1, n) end do mesh%ocean_area=0.0 + mesh%ocean_areawithcav=0.0 call MPI_AllREDUCE(vol, mesh%ocean_area, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(vol2, mesh%ocean_areawithcav, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + MPI_COMM_FESOM, MPIerr) if (mype==0) then write(*,*) mype, 'Mesh statistics:' @@ -1921,10 +1946,12 @@ SUBROUTINE mesh_areas(mesh) ' MinScArea ', minval(mesh%area(1,:)) write(*,*) mype, 'Edges: ', mesh%edge2D, ' internal ', mesh%edge2D_in if (mype==0) then - write(*,*) 'Total ocean area is: ', mesh%ocean_area, ' m^2' + write(*,*) 'Total ocean surface area is : ', mesh%ocean_area, ' m^2' + write(*,*) 'Total ocean surface area wth cavity is: ', mesh%ocean_areawithcav, ' m^2' end if endif + t1=MPI_Wtime() if (mype==0) then write(*,*) 'mesh_areas finished in ', t1-t0, ' seconds' diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index 5ab678381..b1b646aba 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -205,6 +205,7 @@ MODULE o_ARRAYS real(kind=WP), allocatable :: ssh_rhs(:), hpressure(:,:) real(kind=WP), allocatable :: CFL_z(:,:) real(kind=WP), allocatable :: stress_surf(:,:) +real(kind=WP), allocatable :: stress_node_surf(:,:) REAL(kind=WP), ALLOCATABLE :: stress_atmoce_x(:) REAL(kind=WP), ALLOCATABLE :: stress_atmoce_y(:) real(kind=WP), allocatable :: T_rhs(:,:) diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index 236878576..fe4273765 100755 --- a/src/oce_setup_step.F90 +++ b/src/oce_setup_step.F90 @@ -247,6 +247,7 @@ SUBROUTINE array_setup(mesh) ! ================ allocate(Tclim(nl-1,node_size), Sclim(nl-1, node_size)) allocate(stress_surf(2,myDim_elem2D)) !!! Attention, it is shorter !!! +allocate(stress_node_surf(2,node_size)) allocate(stress_atmoce_x(node_size), stress_atmoce_y(node_size)) allocate(relax2clim(node_size)) allocate(heat_flux(node_size), Tsurf(node_size)) @@ -371,8 +372,11 @@ SUBROUTINE array_setup(mesh) Ssurf_old=0.0_WP !PS real_salt_flux=0.0_WP - stress_atmoce_x=0. - stress_atmoce_y=0. + + stress_surf =0.0_WP + stress_node_surf =0.0_WP + stress_atmoce_x =0.0_WP + stress_atmoce_y =0.0_WP tr_arr=0.0_WP tr_arr_old=0.0_WP From 215bec169c5530b87b61b9f2f0801a6914142d2d Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 3 Feb 2021 14:00:28 +0100 Subject: [PATCH 02/40] fix bug mesh partitioning, that caused in cavity case to have isolated elements --- src/fvom_init.F90 | 98 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 83 insertions(+), 15 deletions(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 810899d93..2a2c73972 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -831,7 +831,9 @@ subroutine find_levels_cavity(mesh) integer :: exit_flag, count_iter, count_neighb_open, nneighb, cavity_maxlev real(kind=WP) :: dmean character*200 :: file_name - type(t_mesh), intent(inout), target :: mesh + integer, allocatable, dimension(:,:) :: numelemtonode, idxelemtonode + integer, allocatable, dimension(:) :: elemreducelvl + type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" !___________________________________________________________________________ @@ -839,8 +841,7 @@ subroutine find_levels_cavity(mesh) ulevels => mesh%ulevels allocate(mesh%ulevels_nod2D(nod2D)) ulevels_nod2D => mesh%ulevels_nod2D -!!PS allocate(mesh%cavity_flag_n(nod2D)) -!!PS cavity_flag_n => mesh%cavity_flag_n + !___________________________________________________________________________ ! Compute level position of ocean-cavity boundary @@ -862,9 +863,9 @@ subroutine find_levels_cavity(mesh) ulevels(elem) = 1 do nz=1,nlevels(elem)-1 !!PS if(Z(nz) I need 4 valid full depth, 3 valid mid-depth levels + ulevels(elem)=nz ! to compute shechpetkin PGF exit end if end do @@ -872,14 +873,16 @@ subroutine find_levels_cavity(mesh) cavity_maxlev = max(cavity_maxlev,ulevels(elem)) end do + !___________________________________________________________________________ ! Eliminate cells that have two cavity boundary faces --> should not be ! possible in FESOM2.0 ! loop over all cavity levels + allocate(elemreducelvl(elem2d)) do nz=1,cavity_maxlev exit_flag=0 count_iter=0 - + elemreducelvl=0 !_______________________________________________________________________ ! iteration loop within each layer do while((exit_flag==0).and.(count_iter<1000)) @@ -915,6 +918,22 @@ subroutine find_levels_cavity(mesh) end if end do ! --> do i = 1, nneighb +!!PS if (elem==133438) then +!!PS write(*,*) +!!PS write(*,*) 'nz =', nz +!!PS write(*,*) 'elem =', elem +!!PS write(*,*) 'ulevels(elem) =', ulevels(elem) +!!PS write(*,*) 'nlevels(elem) =', nlevels(elem) +!!PS write(*,*) 'elemreducelvl(elem)=',elemreducelvl(elem) +!!PS write(*,*) 'elems =', elems +!!PS write(*,*) 'ulevels(elems) =', ulevels(elems) +!!PS write(*,*) 'nlevels(elems) =', nlevels(elems) +!!PS write(*,*) 'nlvl-ulvl =', nlevels(elems)-ulevels(elems) +!!PS write(*,*) 'elemreducelvl(elems)=',elemreducelvl(elems) +!!PS write(*,*) 'count_neighb_open =',count_neighb_open +!!PS write(*,*) +!!PS end if + !___________________________________________________________ ! check how many open faces to neighboring triangles the cell ! has, if there are less than 2 its isolated (a cell should @@ -925,13 +944,33 @@ subroutine find_levels_cavity(mesh) ! except when this levels would remain less than 3 valid ! bottom levels --> in case make the levels of all sorounding ! one level shallower - if (nlevels(elem)-(nz+1)<=3) then + if (nlevels(elem)-(nz+1)<3) then do j = 1, nneighb - if (elems(j)>0 .and. ulevels(elems(j))>1 ) ulevels(elems(j)) = min(ulevels(elems(j)),nz) + if (elems(j)>0) then + if (ulevels(elems(j))>1 .and. ulevels(elems(j))>ulevels(elem) ) then + ulevels(elems(j)) = min(ulevels(elems(j)),nz) + elemreducelvl(elems(j))=1 + end if + end if end do - else + + ! if ulevel of element was already made once shalower, do + ! not allow ti be deepen again + elseif (elemreducelvl(elem)==0) then ulevels(elem)=nz+1 + ! if ulevel of element still hasent enough neighbors, can not + ! be deepend any more (elemreducelvl(elem)==1), so allow + ! neighbouring elemnts to become shallower + else + do j = 1, nneighb + if (elems(j)>0) then + if (ulevels(elems(j))>1 .and. ulevels(elems(j))>ulevels(elem) ) then + ulevels(elems(j)) = min(ulevels(elems(j)),nz) + elemreducelvl(elems(j))=1 + end if + end if + end do end if !force recheck for all current ocean cells @@ -942,6 +981,7 @@ subroutine find_levels_cavity(mesh) end do ! --> do elem=1,elem2D end do ! --> do while((exit_flag==0).and.(count_iter<1000)) end do ! --> do nz=1,cavity_maxlev + deallocate(elemreducelvl) !___________________________________________________________________________ ! vertical vertice level index of cavity_ocean boundary @@ -952,11 +992,7 @@ subroutine find_levels_cavity(mesh) ! loop over neighbouring triangles do j=1,nneighb node=elem2D_nodes(j,elem) -!!PS if(ulevels_nod2D(node)<=ulevels(elem)) then -!!PS ulevels_nod2D(node)=ulevels(elem) -!!PS end if ulevels_nod2D(node)=min(ulevels_nod2D(node),ulevels(elem)) - end do end do @@ -972,6 +1008,7 @@ subroutine find_levels_cavity(mesh) write(*,*) ' ulevels,nlevels = ',ulevels(elem), nlevels(elem) write(*,*) ' ulevels(neighb) = ',ulevels(elem_neighbors(1:3,elem)) write(*,*) ' nlevels(neighb) = ',nlevels(elem_neighbors(1:3,elem)) + call par_ex(0) end if end do @@ -996,7 +1033,38 @@ subroutine find_levels_cavity(mesh) call par_ex(0) end if end do -!!PS !___________________________________________________________________________ + + !___________________________________________________________________________ + ! check how many triangle elements contribute to every vertice in every layer + ! every vertice in every layer should be connected to at least two triangle + ! elements ! + allocate(numelemtonode(nl,nod2D),idxelemtonode(nl,nod2D)) + numelemtonode=0 + idxelemtonode=0 + do node=1, nod2D + do j=1,nod_in_elem2D_num(node) + elem=nod_in_elem2D(j,node) + do nz=ulevels(elem),nlevels(elem)-1 + numelemtonode(nz,node) = numelemtonode(nz,node) + 1 + idxelemtonode(nz,node) = elem + end do + end do + end do + + exit_flag = 0 + do node=1, nod2D + do nz=1,nl + if (numelemtonode(nz,node)== 1) then + write(*,*) 'ERROR: found vertice with just one triangle:', mype, nz, 'node=',node, ulevels_nod2d(node), nlevels_nod2D(node), & + 'elem=', idxelemtonode(nz,node), ulevels(idxelemtonode(nz,node)), nlevels(idxelemtonode(nz,node)) + exit_flag = 1 + end if + end do + end do + deallocate(numelemtonode,idxelemtonode) + if (exit_flag == 1) call par_ex(0) + + !___________________________________________________________________________ !!PS ! compute nodal cavity flag: 1 yes cavity/ 0 no cavity !!PS cavity_flag = 0 !!PS do node=1,nod2D From 4f6b8a357a296737f861d57ebcec8dc9027d8ae9 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 17 Feb 2021 17:22:29 +0100 Subject: [PATCH 03/40] introduce new variable for the scalar cell area, make here the difference between the mid-scalar cell area and area of the upper scalar cell edge. In case on no cavity both are identical but with cavity they are different from each other. Since below the cavity, the scalar cell area is defined by the area of the lower edge of the scalar control volume --- src/MOD_MESH.F90 | 2 +- src/associate_mesh.h | 14 +- src/oce_mesh.F90 | 442 +++++++++++++++++++++++++++---------------- src/oce_modules.F90 | 2 +- 4 files changed, 294 insertions(+), 166 deletions(-) diff --git a/src/MOD_MESH.F90 b/src/MOD_MESH.F90 index 795242a8d..0e275067a 100644 --- a/src/MOD_MESH.F90 +++ b/src/MOD_MESH.F90 @@ -63,7 +63,7 @@ MODULE MOD_MESH ! ! !___horizontal mesh info________________________________________________________ -real(kind=WP), allocatable, dimension(:,:) :: area, area_inv +real(kind=WP), allocatable, dimension(:,:) :: area, area_inv, areasvol, areasvol_inv real(kind=WP), allocatable, dimension(:) :: mesh_resolution ! diff --git a/src/associate_mesh.h b/src/associate_mesh.h index 30611d44b..26c197ec6 100644 --- a/src/associate_mesh.h +++ b/src/associate_mesh.h @@ -23,7 +23,7 @@ real(kind=WP), dimension(:,:), pointer :: gradient_sca integer, dimension(:) , pointer :: bc_index_nod2D real(kind=WP), dimension(:) , pointer :: zbar, Z, elem_depth integer, dimension(:) , pointer :: nlevels, nlevels_nod2D, nlevels_nod2D_min -real(kind=WP), dimension(:,:), pointer :: area, area_inv +real(kind=WP), dimension(:,:), pointer :: area, area_inv, areasvol, areasvol_inv real(kind=WP), dimension(:) , pointer :: mesh_resolution real(kind=WP), dimension(:) , pointer :: lump2d_north, lump2d_south type(sparse_matrix) , pointer :: ssh_stiff @@ -47,6 +47,7 @@ nl => mesh%nl !!$edge_tri => mesh%edge_tri !!$elem_edges => mesh%elem_edges !!$elem_area => mesh%elem_area +!!$node_area => mesh%node_area !!$edge_dxdy => mesh%edge_dxdy !!$edge_cross_dxdy => mesh%edge_cross_dxdy !!$elem_cos => mesh%elem_cos @@ -66,7 +67,8 @@ nl => mesh%nl !!$nlevels => mesh%nlevels !!$nlevels_nod2D => mesh%nlevels_nod2D !!$nlevels_nod2D_min => mesh%nlevels_nod2D_min -!!$area => mesh%area +!!$area => mesh%area +!!$area2 => mesh%area2 !!$area_inv => mesh%area_inv !!$mesh_resolution => mesh%mesh_resolution !!$ssh_stiff => mesh%ssh_stiff @@ -84,7 +86,7 @@ elem2D_nodes(1:3, 1:myDim_elem2D+eDim_elem2D+eXDim_elem2D) => mesh%elem2D_nodes edges(1:2,1:myDim_edge2D+eDim_edge2D) => mesh%edges edge_tri(1:2,1:myDim_edge2D+eDim_edge2D) => mesh%edge_tri elem_edges(1:3,1:myDim_elem2D) => mesh%elem_edges -elem_area(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D) => mesh%elem_area +elem_area(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D) => mesh%elem_area edge_dxdy(1:2,1:myDim_edge2D+eDim_edge2D) => mesh%edge_dxdy edge_cross_dxdy(1:4,1:myDim_edge2D+eDim_edge2D) => mesh%edge_cross_dxdy elem_cos(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D) => mesh%elem_cos @@ -104,8 +106,10 @@ elem_depth => mesh%elem_depth ! never used, not even allocated nlevels(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D) => mesh%nlevels nlevels_nod2D(1:myDim_nod2D+eDim_nod2D) => mesh%nlevels_nod2D nlevels_nod2D_min(1:myDim_nod2D+eDim_nod2D) => mesh%nlevels_nod2D_min -area(1:mesh%nl,1:myDim_nod2d+eDim_nod2D) => mesh%area -area_inv(1:mesh%nl,1:myDim_nod2d+eDim_nod2D) => mesh%area_inv +area(1:mesh%nl,1:myDim_nod2d+eDim_nod2D) => mesh%area +areasvol(1:mesh%nl,1:myDim_nod2d+eDim_nod2D) => mesh%areasvol +area_inv(1:mesh%nl,1:myDim_nod2d+eDim_nod2D) => mesh%area_inv +areasvol_inv(1:mesh%nl,1:myDim_nod2d+eDim_nod2D) => mesh%areasvol_inv mesh_resolution(1:myDim_nod2d+eDim_nod2D) => mesh%mesh_resolution ssh_stiff => mesh%ssh_stiff lump2d_north(1:myDim_nod2d) => mesh%lump2d_north diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 763c31020..2a6d5980a 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -118,13 +118,17 @@ SUBROUTINE mesh_setup(mesh) call set_mesh_transform_matrix !(rotated grid) call read_mesh(mesh) call set_par_support(mesh) - call find_levels(mesh) - - if (use_cavity) call find_levels_cavity(mesh) - +!!PS call find_levels(mesh) +!!PS +!!PS if (use_cavity) call find_levels_cavity(mesh) +!!PS call test_tri(mesh) call load_edges(mesh) call find_neighbors(mesh) + + call find_levels(mesh) + if (use_cavity) call find_levels_cavity(mesh) + call find_levels_min_e2n(mesh) call mesh_areas(mesh) call mesh_auxiliary_arrays(mesh) @@ -904,7 +908,8 @@ subroutine find_levels_cavity(mesh) integer, allocatable, dimension(:) :: ibuff real(kind=WP) :: t0, t1 logical :: file_exist=.False. - integer :: elem, elnodes(3), ule, uln(3) + integer :: elem, elnodes(3), ule, uln(3), node, j, nz + integer, allocatable, dimension(:,:) :: numelemtonode !NR Cannot include the pointers before the targets are allocated... !NR #include "associate_mesh.h" @@ -1283,6 +1288,30 @@ subroutine find_levels_cavity(mesh) end if end do + + !___________________________________________________________________________ + allocate(numelemtonode(mesh%nl,myDim_nod2d+eDim_nod2D)) + numelemtonode=0 + do node=1, myDim_nod2D+eDim_nod2D + do j=1,mesh%nod_in_elem2D_num(node) + elem=mesh%nod_in_elem2D(j,node) + do nz=mesh%ulevels(elem),mesh%nlevels(elem)-1 + numelemtonode(nz,node) = numelemtonode(nz,node) + 1 + end do + end do + end do + + ! check how many triangle elements contribute to every vertice in every layer + ! every vertice in every layer should be connected to at least two triangle + ! elements ! + do node=1, myDim_nod2D+eDim_nod2D + do nz=1,mesh%nl + if (numelemtonode(nz,node)== 1) then + write(*,*) 'ERROR A: found vertice with just one triangle:', mype, node, nz + end if + end do + end do + end subroutine find_levels_cavity ! ! @@ -1799,164 +1828,257 @@ subroutine elem_center(elem, x, y, mesh) end subroutine elem_center !========================================================================== SUBROUTINE mesh_areas(mesh) -USE MOD_MESH -USE o_PARAM -USE g_PARSUP -USE g_ROTATE_GRID -use g_comm_auto -IMPLICIT NONE -! Collects auxilliary information on the mesh -! Allocated and filled in are: -! elem_area(myDim_elem2D) -! area(nl, myDim_nod2D) - - -integer :: n,j,q, elnodes(3), ed(2), elem, nz,nzmin -real(kind=WP) :: a(2), b(2), ax, ay, lon, lat, vol, vol2 -real(kind=WP), allocatable,dimension(:) :: work_array -real(kind=WP) :: t0, t1 -type(t_mesh), intent(inout), target :: mesh - -!NR Cannot include the pointers before the targets are allocated... -!NR #include "associate_mesh.h" - -t0=MPI_Wtime() - - allocate(mesh%elem_area(myDim_elem2D+eDim_elem2D+eXDim_elem2D)) - !allocate(elem_area(myDim_elem2D)) - allocate(mesh%area(mesh%nl,myDim_nod2d+eDim_nod2D)) !! Extra size just for simplicity - !! in some further routines - allocate(mesh%area_inv(mesh%nl,myDim_nod2d+eDim_nod2D)) - allocate(mesh%mesh_resolution(myDim_nod2d+eDim_nod2D)) - ! ============ - ! The areas of triangles: - ! ============ - DO n=1, myDim_elem2D - !DO n=1, myDim_elem2D+eDim_elem2D+eXDim_elem2D - elnodes=mesh%elem2D_nodes(:,n) - ay=sum(mesh%coord_nod2D(2,elnodes))/3.0_WP - ay=cos(ay) - if (cartesian) ay=1.0_WP - a = mesh%coord_nod2D(:,elnodes(2))-mesh%coord_nod2D(:,elnodes(1)) - b = mesh%coord_nod2D(:,elnodes(3))-mesh%coord_nod2D(:,elnodes(1)) - call trim_cyclic(a(1)) - call trim_cyclic(b(1)) - a(1)=a(1)*ay - b(1)=b(1)*ay - mesh%elem_area(n)=0.5_WP*abs(a(1)*b(2)-b(1)*a(2)) - END DO - call exchange_elem(mesh%elem_area) - ! ============= - ! Scalar element - ! areas at different levels (there can be partly land) - ! ============= - - mesh%area=0.0_WP - DO n=1, myDim_nod2D - DO j=1,mesh%nod_in_elem2D_num(n) - elem=mesh%nod_in_elem2D(j,n) - DO nz=mesh%ulevels(elem),mesh%nlevels(elem)-1 - !!PS DO nz=1,mesh%nlevels(elem)-1 - mesh%area(nz,n)=mesh%area(nz,n)+mesh%elem_area(elem)/3.0_WP - END DO - END DO - END DO - ! Only areas through which there is exchange are counted - - ! =========== - ! Update to proper dimension - ! =========== - mesh%elem_area=mesh%elem_area*r_earth*r_earth - mesh%area=mesh%area*r_earth*r_earth - -call exchange_nod(mesh%area) - - -!!PS do n=1, myDim_nod2D -!!PS nzmin = mesh%ulevels_nod2d(n) -!!PS if (nzmin>1) then -!!PS write(*,*) ' --> mesh area:', mype, n, nzmin, mesh%area(nzmin,n),mesh%area(nzmin+1,n),mesh%area(nzmin+2,n) -!!PS end if -!!PS end do - - + USE MOD_MESH + USE o_PARAM + USE o_arrays, only: dum_3d_n + USE g_PARSUP + USE g_ROTATE_GRID + use g_comm_auto + IMPLICIT NONE + ! Collects auxilliary information on the mesh + ! Allocated and filled in are: + ! elem_area(myDim_elem2D) + ! area(nl, myDim_nod2D) + + integer :: n,j,q, elnodes(3), ed(2), elem, nz,nzmin, nzmax + real(kind=WP) :: a(2), b(2), ax, ay, lon, lat, vol, vol2 + real(kind=WP), allocatable,dimension(:) :: work_array + integer, allocatable,dimension(:,:) :: cavity_contribut + real(kind=WP) :: t0, t1 + type(t_mesh), intent(inout), target :: mesh + + !NR Cannot include the pointers before the targets are allocated... + !NR #include "associate_mesh.h" -do n=1,myDim_nod2d+eDim_nod2D - do nz=1,mesh%nl - if (mesh%area(nz,n) > 0._WP) then - mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) - else - mesh%area_inv(nz,n) = 0._WP - end if - end do -end do - ! coordinates are in radians, edge_dxdy are in meters, - ! and areas are in m^2 - -!!PS do n=1,myDim_nod2d+eDim_nod2D -!!PS mesh%area_inv(1:mesh%ulevels_nod2D(n)-1,n) = 0.0_WP -!!PS mesh%area(1:mesh%ulevels_nod2D(n)-1,n) = 0.0_WP -!!PS end do + t0=MPI_Wtime() + + ! area of triangles + allocate(mesh%elem_area(myDim_elem2D+eDim_elem2D+eXDim_elem2D)) + + ! area of upper edge and lower edge of scalar cell: size nl x node + allocate(mesh%area(mesh%nl,myDim_nod2d+eDim_nod2D)) + + ! "mid" area of scalar cell in case of cavity area \= areasvol, size: nl-1 x node + allocate(mesh%areasvol(mesh%nl-1,myDim_nod2d+eDim_nod2D)) + + ! area inverse + allocate(mesh%area_inv(mesh%nl,myDim_nod2d+eDim_nod2D)) + allocate(mesh%areasvol_inv(mesh%nl-1,myDim_nod2d+eDim_nod2D)) + + ! resolution at nodes + allocate(mesh%mesh_resolution(myDim_nod2d+eDim_nod2D)) + + !___compute triangle areas__________________________________________________ + do n=1, myDim_elem2D + elnodes=mesh%elem2D_nodes(:,n) + ay=sum(mesh%coord_nod2D(2,elnodes))/3.0_WP + ay=cos(ay) + if (cartesian) ay=1.0_WP + a = mesh%coord_nod2D(:,elnodes(2))-mesh%coord_nod2D(:,elnodes(1)) + b = mesh%coord_nod2D(:,elnodes(3))-mesh%coord_nod2D(:,elnodes(1)) + call trim_cyclic(a(1)) + call trim_cyclic(b(1)) + a(1)=a(1)*ay + b(1)=b(1)*ay + mesh%elem_area(n)=0.5_WP*abs(a(1)*b(2)-b(1)*a(2)) + end do + call exchange_elem(mesh%elem_area) + + !___compute areas of upper/lower scalar cell edge___________________________ + ! areas at different levels (there can be partly land) + ! --> only areas through which there is exchange are counted + ! + !-----------------------------~+~~~~~~~+~~~ + ! ############################ | | + ! ############################ | | layer k-3 + ! #################### ._______|_______|___area_k-2 + ! ## CAVITY ######## | / / / | | + ! #################### |/ /°/ /| | layer k-2 --> Transport: w_k-2*A_k-1 + ! ############ ._______|_/_/_/_|_______|___area_k-1 -> A_k-1 lower prisma area defines + ! ############ | | | | scalar area under the cavity + ! ############ | ° | | | layer k-1 + !______________|_______|_______|_______|___area_k + ! | | / / / | | | + ! | |/ /°/ /| | | layer k --> Transport: w_k*A_k + !______|_______|_/_/_/_|_______|_______|___area_k+1 -> A_k upper prisma face area defines + ! | | | | | scalar area of cell + ! | | ° | | | layer k+1 + !______|_______|_______|_______|_______|___area_k+2 + ! #############| | | | + ! #############| ° | | | layer k+2 + ! #############|_______|_______|_______|___area_k+3 + ! #####################| | | + ! #####################| | | layer k+3 + ! ## BOTTOM #########|_______|_______|___area_k+4 + ! #############################| | + ! #############################| | : + ! #############################|_______|___area_k+5 + ! ######################################### + mesh%area = 0.0_WP + mesh%areasvol = 0.0_WP + + if (use_cavity) then + allocate(cavity_contribut(mesh%nl,myDim_nod2d+eDim_nod2D)) + cavity_contribut = 0 + end if + + do n=1, myDim_nod2D+eDim_nod2D + do j=1,mesh%nod_in_elem2D_num(n) + elem=mesh%nod_in_elem2D(j,n) + + !___________________________________________________________________ + ! compute scalar area of prisms at different depth layers. In normal + ! case without cavity the area of the scalar cell corresponds to the + ! area of the upper edge of the prism --> if there is cavity its + ! different. Directly under the cavity the area of scalar cell + ! corresponds to the area of the lower edge + nzmin = mesh%ulevels(elem) + nzmax = mesh%nlevels(elem) + do nz=nzmin,nzmax + mesh%area(nz,n)=mesh%area(nz,n)+mesh%elem_area(elem)/3.0_WP + end do + + !___________________________________________________________________ + ! how many ocean-cavity triangles contribute to an upper edge of a + ! scalar area + if (use_cavity) then + do nz=1,nzmin-1 + cavity_contribut(nz,n)=cavity_contribut(nz,n)+1 + end do + end if + end do + end do + + !___compute "mid" scalar cell area__________________________________________ + ! for cavity case: redefine scalar cell area from upper edge of prism to + ! lower edge of prism if a cavity triangle is present at the upper scalar + ! cell edge + if (use_cavity) then + do n = 1, myDim_nod2D+eDim_nod2D + nzmin = mesh%ulevels_nod2d(n) + nzmax = mesh%nlevels_nod2d(n)-1 + do nz=nzmin,nzmax + if (cavity_contribut(nz,n)>0) then + mesh%areasvol(nz,n) = mesh%area(min(nz+1,nzmax),n) + end if + end do + end do + deallocate(cavity_contribut) + ! for non cavity case: the "mid" area of the scalar cell always corresponds to + ! the area of the upper scalar cell edge + else + do n = 1, myDim_nod2D+eDim_nod2D + nzmin = mesh%ulevels_nod2d(n) + nzmax = mesh%nlevels_nod2d(n)-1 + do nz=nzmin,nzmax + mesh%areasvol(nz,n) = mesh%area(nz,n) + end do + end do + end if + + ! update to proper dimension + ! coordinates are in radians, edge_dxdy are in meters, + ! and areas are in m^2 + mesh%elem_area = mesh%elem_area*r_earth*r_earth + mesh%area = mesh%area *r_earth*r_earth + mesh%areasvol = mesh%areasvol *r_earth*r_earth + call exchange_nod(mesh%area) + call exchange_nod(mesh%areasvol) + + !___compute inverse area____________________________________________________ + mesh%area_inv = 0.0_WP + do n=1,myDim_nod2d+eDim_nod2D + nzmin = mesh%ulevels_nod2d(n) + nzmax = mesh%nlevels_nod2d(n) + do nz=nzmin,nzmax + mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) +!!PS if (mesh%area(nz,n) > 0._WP) then +!!PS mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) +!!PS else +!!PS mesh%area_inv(nz,n) = 0._WP +!!PS end if + end do + end do + + if (use_cavity) then + mesh%areasvol_inv = 0.0_WP + do n=1,myDim_nod2d+eDim_nod2D + nzmin = mesh%ulevels_nod2d(n) + nzmax = mesh%nlevels_nod2d(n)-1 + do nz=nzmin,nzmax + mesh%areasvol_inv(nz,n) = 1._WP/mesh%areasvol(nz,n) +!!PS if (mesh%areasvol(nz,n) > 0._WP) then +!!PS mesh%areasvol_inv(nz,n) = 1._WP/mesh%areasvol(nz,n) +!!PS else +!!PS mesh%areasvol_inv(nz,n) = 0._WP +!!PS end if + end do + end do + else + mesh%areasvol_inv = mesh%area_inv + endif - - allocate(work_array(myDim_nod2D)) - !!PS mesh%mesh_resolution=sqrt(mesh%area(1, :)/pi)*2._WP - do n=1,myDim_nod2d+eDim_nod2D - mesh%mesh_resolution(n)=sqrt(mesh%area(mesh%ulevels_nod2D(n), n)/pi)*2._WP - end do + !___compute scalar cell resolution__________________________________________ + allocate(work_array(myDim_nod2D)) + !!PS mesh%mesh_resolution=sqrt(mesh%area(1, :)/pi)*2._WP + do n=1,myDim_nod2d+eDim_nod2D + mesh%mesh_resolution(n)=sqrt(mesh%areasvol(mesh%ulevels_nod2d(n),n)/pi)*2._WP + end do - DO q=1, 3 !apply mass matrix N times to smooth the field - DO n=1, myDim_nod2D - vol=0._WP - work_array(n)=0._WP - DO j=1, mesh%nod_in_elem2D_num(n) - elem=mesh%nod_in_elem2D(j, n) - elnodes=mesh%elem2D_nodes(:,elem) - work_array(n)=work_array(n)+sum(mesh%mesh_resolution(elnodes))/3._WP*mesh%elem_area(elem) - vol=vol+mesh%elem_area(elem) - END DO - work_array(n)=work_array(n)/vol + ! smooth resolution + DO q=1, 3 !apply mass matrix N times to smooth the field + DO n=1, myDim_nod2D + vol=0._WP + work_array(n)=0._WP + DO j=1, mesh%nod_in_elem2D_num(n) + elem=mesh%nod_in_elem2D(j, n) + elnodes=mesh%elem2D_nodes(:,elem) + work_array(n)=work_array(n)+sum(mesh%mesh_resolution(elnodes))/3._WP*mesh%elem_area(elem) + vol=vol+mesh%elem_area(elem) + END DO + work_array(n)=work_array(n)/vol + END DO + DO n=1,myDim_nod2D + mesh%mesh_resolution(n)=work_array(n) + ENDDO + call exchange_nod(mesh%mesh_resolution) END DO - DO n=1,myDim_nod2D - mesh%mesh_resolution(n)=work_array(n) - ENDDO - call exchange_nod(mesh%mesh_resolution) - END DO - deallocate(work_array) - - vol=0.0_WP - vol2=0.0_WP - do n=1, myDim_nod2D - vol2=vol2+mesh%area(mesh%ulevels_nod2D(n), n) - if (mesh%ulevels_nod2D(n)>1) cycle - vol=vol+mesh%area(1, n) - end do - mesh%ocean_area=0.0 - mesh%ocean_areawithcav=0.0 - call MPI_AllREDUCE(vol, mesh%ocean_area, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & - MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(vol2, mesh%ocean_areawithcav, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & - MPI_COMM_FESOM, MPIerr) - -if (mype==0) then - write(*,*) mype, 'Mesh statistics:' - write(*,*) mype, 'maxArea ',maxval(mesh%elem_area), ' MinArea ', minval(mesh%elem_area) - write(*,*) mype, 'maxScArea ',maxval(mesh%area(1,:)), & - ' MinScArea ', minval(mesh%area(1,:)) - write(*,*) mype, 'Edges: ', mesh%edge2D, ' internal ', mesh%edge2D_in - if (mype==0) then - write(*,*) 'Total ocean surface area is : ', mesh%ocean_area, ' m^2' - write(*,*) 'Total ocean surface area wth cavity is: ', mesh%ocean_areawithcav, ' m^2' - end if -endif - + deallocate(work_array) + + !___compute total ocean areas with/without cavity___________________________ + vol = 0.0_WP + vol2= 0.0_WP + do n=1, myDim_nod2D + vol2=vol2+mesh%area(mesh%ulevels_nod2D(n), n) ! area also under cavity + if (mesh%ulevels_nod2D(n)>1) cycle + vol=vol+mesh%area(1, n) ! area only surface + end do + mesh%ocean_area=0.0 + mesh%ocean_areawithcav=0.0 + call MPI_AllREDUCE(vol, mesh%ocean_area, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(vol2, mesh%ocean_areawithcav, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + MPI_COMM_FESOM, MPIerr) + + !___write mesh statistics___________________________________________________ + if (mype==0) then + write(*,*) '____________________________________________________________________' + write(*,*) ' --> mesh statistics:', mype + write(*,*) mype, 'maxArea ',maxval(mesh%elem_area), ' MinArea ', minval(mesh%elem_area) + write(*,*) mype, 'maxScArea ',maxval(mesh%area(1,:)), & + ' MinScArea ', minval(mesh%area(1,:)) + write(*,*) mype, 'Edges: ', mesh%edge2D, ' internal ', mesh%edge2D_in + if (mype==0) then + write(*,*) ' > Total ocean surface area is : ', mesh%ocean_area, ' m^2' + write(*,*) ' > Total ocean surface area wth cavity is: ', mesh%ocean_areawithcav, ' m^2' + end if + endif -t1=MPI_Wtime() -if (mype==0) then - write(*,*) 'mesh_areas finished in ', t1-t0, ' seconds' - write(*,*) '=========================' -endif + t1=MPI_Wtime() + if (mype==0) then + write(*,*) ' > mesh_areas finished in ', t1-t0, ' seconds' + endif END SUBROUTINE mesh_areas !=================================================================== @@ -2288,10 +2410,12 @@ SUBROUTINE mesh_auxiliary_arrays(mesh) do i=1, myDim_nod2D if (mesh%geo_coord_nod2D(2, i) > 0) then nn=nn+1 - mesh%lump2d_north(i)=mesh%area(1, i) +!!PS mesh%lump2d_north(i)=mesh%area(1, i)! --> TEST_cavity + mesh%lump2d_north(i)=mesh%node_area(i) else ns=ns+1 - mesh%lump2d_south(i)=mesh%area(1, i) +!!PS mesh%lump2d_south(i)=mesh%area(1, i)! --> TEST_cavity + mesh%lump2d_south(i)=mesh%node_area(i) end if end do diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index b1b646aba..460828137 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -304,7 +304,7 @@ MODULE o_ARRAYS !_______________________________________________________________________________ !!PS ! dummy arrays -!!PS real(kind=WP), allocatable,dimension(:,:) :: dum_3d_n, dum_3d_e +real(kind=WP), allocatable,dimension(:,:) :: dum_3d_n !, dum_3d_e !!PS real(kind=WP), allocatable,dimension(:) :: dum_2d_n, dum_2d_e !_______________________________________________________________________________ From 08d0d24d486f0c5b1006d1a07b5b0f7c63315056 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 17 Feb 2021 17:40:29 +0100 Subject: [PATCH 04/40] improve computation of mid scalar cell area in oce_mesh.F90 --- src/oce_mesh.F90 | 36 +++++++++++++++++++----------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 2a6d5980a..cc5a9e599 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1950,7 +1950,7 @@ SUBROUTINE mesh_areas(mesh) end do !___compute "mid" scalar cell area__________________________________________ - ! for cavity case: redefine scalar cell area from upper edge of prism to + ! for cavity case: redefine "mid" scalar cell area from upper edge of prism to ! lower edge of prism if a cavity triangle is present at the upper scalar ! cell edge if (use_cavity) then @@ -1960,6 +1960,8 @@ SUBROUTINE mesh_areas(mesh) do nz=nzmin,nzmax if (cavity_contribut(nz,n)>0) then mesh%areasvol(nz,n) = mesh%area(min(nz+1,nzmax),n) + else + mesh%areasvol(nz,n) = mesh%area(nz,n) end if end do end do @@ -2027,23 +2029,23 @@ SUBROUTINE mesh_areas(mesh) end do ! smooth resolution - DO q=1, 3 !apply mass matrix N times to smooth the field - DO n=1, myDim_nod2D - vol=0._WP - work_array(n)=0._WP - DO j=1, mesh%nod_in_elem2D_num(n) - elem=mesh%nod_in_elem2D(j, n) - elnodes=mesh%elem2D_nodes(:,elem) - work_array(n)=work_array(n)+sum(mesh%mesh_resolution(elnodes))/3._WP*mesh%elem_area(elem) - vol=vol+mesh%elem_area(elem) - END DO - work_array(n)=work_array(n)/vol - END DO - DO n=1,myDim_nod2D - mesh%mesh_resolution(n)=work_array(n) - ENDDO + do q=1, 3 !apply mass matrix N times to smooth the field + do n=1, myDim_nod2D + vol=0._WP + work_array(n)=0._WP + do j=1, mesh%nod_in_elem2D_num(n) + elem=mesh%nod_in_elem2D(j, n) + elnodes=mesh%elem2D_nodes(:,elem) + work_array(n)=work_array(n)+sum(mesh%mesh_resolution(elnodes))/3._WP*mesh%elem_area(elem) + vol=vol+mesh%elem_area(elem) + end do + work_array(n)=work_array(n)/vol + end do + do n=1,myDim_nod2D + mesh%mesh_resolution(n)=work_array(n) + end do call exchange_nod(mesh%mesh_resolution) - END DO + end do deallocate(work_array) !___compute total ocean areas with/without cavity___________________________ From d746b16c94d55b24d23e1d4f8dc50b855628d482 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 17 Feb 2021 17:44:02 +0100 Subject: [PATCH 05/40] exchange area_inv with areasvol_inv in subroutine momentum_adv_scalar --- src/oce_ale_vel_rhs.F90 | 312 ++++++++++++++++++++++------------------ 1 file changed, 174 insertions(+), 138 deletions(-) diff --git a/src/oce_ale_vel_rhs.F90 b/src/oce_ale_vel_rhs.F90 index a362bd1fe..210c7f28d 100644 --- a/src/oce_ale_vel_rhs.F90 +++ b/src/oce_ale_vel_rhs.F90 @@ -168,150 +168,186 @@ subroutine momentum_adv_scalar(mesh) #include "associate_mesh.h" -!_______________________________________________________________________________ -do n=1,myDim_nod2d - nl1 = nlevels_nod2D(n)-1 - ul1 = ulevels_nod2D(n) - wu(1:nl1+1) = 0._WP - wv(1:nl1+1) = 0._WP - - do k=1,nod_in_elem2D_num(n) - el = nod_in_elem2D(k,n) - nle = nlevels(el)-1 - ule = ulevels(el) - !___________________________________________________________________________ - ! The vertical part for each element is collected - !!PS wu(1) = wu(1) + UV(1,1,el)*elem_area(el) - !!PS wv(1) = wv(1) + UV(2,1,el)*elem_area(el) - wu(ule) = wu(ule) + UV(1,ule,el)*elem_area(el) - wv(ule) = wv(ule) + UV(2,ule,el)*elem_area(el) - - !!PS wu(2:nle) = wu(2:nle) + 0.5_WP*(UV(1,2:nle,el)+UV(1,1:nle-1,el))*elem_area(el) - !!PS wv(2:nle) = wv(2:nle) + 0.5_WP*(UV(2,2:nle,el)+UV(2,1:nle-1,el))*elem_area(el) - wu(ule+1:nle) = wu(ule+1:nle) + 0.5_WP*(UV(1,ule+1:nle,el)+UV(1,ule:nle-1,el))*elem_area(el) - wv(ule+1:nle) = wv(ule+1:nle) + 0.5_WP*(UV(2,ule+1:nle,el)+UV(2,ule:nle-1,el))*elem_area(el) - enddo - - !!PS wu(1:nl1) = wu(1:nl1)*Wvel_e(1:nl1,n) - !!PS wv(1:nl1) = wv(1:nl1)*Wvel_e(1:nl1,n) - wu(ul1:nl1) = wu(ul1:nl1)*Wvel_e(ul1:nl1,n) - wv(ul1:nl1) = wv(ul1:nl1)*Wvel_e(ul1:nl1,n) - - !!PS do nz=1,nl1 - do nz=ul1,nl1 - ! Here 1/3 because 1/3 of the area is related to the node - Unode_rhs(1,nz,n) = - (wu(nz) - wu(nz+1) ) / (3._WP*hnode(nz,n)) - Unode_rhs(2,nz,n) = - (wv(nz) - wv(nz+1) ) / (3._WP*hnode(nz,n)) - - enddo - - ! To get a clean checksum, set the remaining values to zero - Unode_rhs(1:2,nl1+1:nl-1,n) = 0._WP - Unode_rhs(1:2,1:ul1-1 ,n) = 0._WP -end do - - -!_______________________________________________________________________________ -DO ed=1, myDim_edge2D - nod = edges(:,ed) - el1 = edge_tri(1,ed) - el2 = edge_tri(2,ed) - nl1 = nlevels(el1)-1 - ul1 = ulevels(el1) - - !___________________________________________________________________________ - ! The horizontal part - !!PS un1(1:nl1) = UV(2,1:nl1,el1)*edge_cross_dxdy(1,ed) & - !!PS - UV(1,1:nl1,el1)*edge_cross_dxdy(2,ed) - un1(ul1:nl1) = UV(2,ul1:nl1,el1)*edge_cross_dxdy(1,ed) & - - UV(1,ul1:nl1,el1)*edge_cross_dxdy(2,ed) - !___________________________________________________________________________ - if (el2>0) then - nl2 = nlevels(el2)-1 - ul2 = ulevels(el2) - - !!PS un2(1:nl2) = - UV(2,1:nl2,el2)*edge_cross_dxdy(3,ed) & - !!PS + UV(1,1:nl2,el2)*edge_cross_dxdy(4,ed) - un2(ul2:nl2) = - UV(2,ul2:nl2,el2)*edge_cross_dxdy(3,ed) & - + UV(1,ul2:nl2,el2)*edge_cross_dxdy(4,ed) - - ! fill with zeros to combine the loops - ! Usually, no or only a very few levels have to be filled. In this case, - ! computing "zeros" is cheaper than the loop overhead. - un1(nl1+1:max(nl1,nl2)) = 0._WP - un2(nl2+1:max(nl1,nl2)) = 0._WP - un1(1:ul1-1) = 0._WP - un2(1:ul2-1) = 0._WP - - ! first edge node - ! Do not calculate on Halo nodes, as the result will not be used. - ! The "if" is cheaper than the avoided computiations. - if (nod(1) <= myDim_nod2d) then - !!PS do nz=1, max(nl1,nl2) - do nz=min(ul1,ul2), max(nl1,nl2) - Unode_rhs(1,nz,nod(1)) = Unode_rhs(1,nz,nod(1)) + un1(nz)*UV(1,nz,el1) + un2(nz)*UV(1,nz,el2) - Unode_rhs(2,nz,nod(1)) = Unode_rhs(2,nz,nod(1)) + un1(nz)*UV(2,nz,el1) + un2(nz)*UV(2,nz,el2) - end do - endif - - if (nod(2) <= myDim_nod2d) then - !!PS do nz=1, max(nl1,nl2) - do nz=min(ul1,ul2), max(nl1,nl2) - Unode_rhs(1,nz,nod(2)) = Unode_rhs(1,nz,nod(2)) - un1(nz)*UV(1,nz,el1) - un2(nz)*UV(1,nz,el2) - Unode_rhs(2,nz,nod(2)) = Unode_rhs(2,nz,nod(2)) - un1(nz)*UV(2,nz,el1) - un2(nz)*UV(2,nz,el2) - end do - endif + !___________________________________________________________________________ + ! 1st. compute vertical momentum advection component: w * du/dz, w*dv/dz + do n=1,myDim_nod2d + nl1 = nlevels_nod2D(n)-1 + ul1 = ulevels_nod2D(n) + wu(1:nl1+1) = 0._WP + wv(1:nl1+1) = 0._WP + + !_______________________________________________________________________ + ! loop over adjacent elements of vertice n + do k=1,nod_in_elem2D_num(n) + el = nod_in_elem2D(k,n) + !___________________________________________________________________ + nle = nlevels(el)-1 + ule = ulevels(el) + + !___________________________________________________________________ + ! accumulate horizontal velocities at full depth levels (top and + ! bottom faces of prism) + ! account here also for boundary condition below cavity --> + ! horizontal velocity at cavity-ocean interce ule (if ule>1) must be + ! zero ??? + if (ule==1) then + wu(ule) = wu(ule) + UV(1,ule,el)*elem_area(el) + wv(ule) = wv(ule) + UV(2,ule,el)*elem_area(el) + end if + + ! interpolate horizontal velocity from mid-depth levels to full + ! depth levels of upper and lower prism faces and average over adjacent + ! elements of vertice n + wu(ule+1:nle) = wu(ule+1:nle) + 0.5_WP*(UV(1,ule+1:nle,el)+UV(1,ule:nle-1,el))*elem_area(el) + wv(ule+1:nle) = wv(ule+1:nle) + 0.5_WP*(UV(2,ule+1:nle,el)+UV(2,ule:nle-1,el))*elem_area(el) + enddo + + !_______________________________________________________________________ + ! multiply w*du and w*dv + wu(ul1:nl1) = wu(ul1:nl1)*Wvel_e(ul1:nl1,n) + wv(ul1:nl1) = wv(ul1:nl1)*Wvel_e(ul1:nl1,n) + + !_______________________________________________________________________ + ! compute w*du/dz, w*dv/dz + do nz=ul1,nl1 +!!PS if (ul1>1) write(*,*) mype, wu(ul1:nl1) + ! Here 1/3 because 1/3 of the area is related to the node --> comes from + ! averaging the elemental velocities + Unode_rhs(1,nz,n) = - (wu(nz) - wu(nz+1) ) / (3._WP*hnode(nz,n)) + Unode_rhs(2,nz,n) = - (wv(nz) - wv(nz+1) ) / (3._WP*hnode(nz,n)) + + enddo + + !_______________________________________________________________________ + ! To get a clean checksum, set the remaining values to zero + Unode_rhs(1:2,nl1+1:nl-1,n) = 0._WP + Unode_rhs(1:2,1:ul1-1 ,n) = 0._WP + end do - else ! ed is a boundary edge, there is only the contribution from el1 - if (nod(1) <= myDim_nod2d) then - !!PS do nz=1, nl1 - do nz=ul1, nl1 + !___________________________________________________________________________ + ! 2nd. compute horizontal advection component: u*du/dx, u*dv/dx & v*du/dy, v*dv/dy + ! loop over triangle edges + do ed=1, myDim_edge2D + nod = edges(:,ed) + el1 = edge_tri(1,ed) + el2 = edge_tri(2,ed) + nl1 = nlevels(el1)-1 + ul1 = ulevels(el1) + + !_______________________________________________________________________ + ! compute horizontal normal velocity with respect to the edge from triangle + ! centroid towards triangel edge mid-pointe for element el1 + ! .o. + ! ./ \. + ! ./ el1 \. + ! ./ x \. + ! ./ |-------\.-----------------edge_cross_dxdy(1:2,ed) --> (dx,dy) + ! / |->n_vec \ + ! nod(1) o----------O----------o nod(2) + ! \. |->n_vec ./ + ! \. |------./------------------edge_cross_dxdy(3:4,ed) --> (dx,dy) + ! \. x ./ + ! \. el2 ./ + ! \. ./ + ! ° + un1(ul1:nl1) = UV(2,ul1:nl1,el1)*edge_cross_dxdy(1,ed) & + - UV(1,ul1:nl1,el1)*edge_cross_dxdy(2,ed) + + !_______________________________________________________________________ + ! compute horizontal normal velocity with respect to the edge from triangle + ! centroid towards triangel edge mid-pointe for element el2 when it is valid + ! --> if its a boundary triangle el2 will be not valid + if (el2>0) then ! --> el2 is valid element + nl2 = nlevels(el2)-1 + ul2 = ulevels(el2) - Unode_rhs(1,nz,nod(1)) = Unode_rhs(1,nz,nod(1)) + un1(nz)*UV(1,nz,el1) - Unode_rhs(2,nz,nod(1)) = Unode_rhs(2,nz,nod(1)) + un1(nz)*UV(2,nz,el1) - end do - endif - ! second edge node - if (nod(2) <= myDim_nod2d) then - !!PS do nz=1, nl1 - do nz=ul1, nl1 - Unode_rhs(1,nz,nod(2)) = Unode_rhs(1,nz,nod(2)) - un1(nz)*UV(1,nz,el1) - Unode_rhs(2,nz,nod(2)) = Unode_rhs(2,nz,nod(2)) - un1(nz)*UV(2,nz,el1) - end do - endif - endif - -end do + un2(ul2:nl2) = - UV(2,ul2:nl2,el2)*edge_cross_dxdy(3,ed) & + + UV(1,ul2:nl2,el2)*edge_cross_dxdy(4,ed) + + ! fill with zeros to combine the loops + ! Usually, no or only a very few levels have to be filled. In this case, + ! computing "zeros" is cheaper than the loop overhead. + un1(nl1+1:max(nl1,nl2)) = 0._WP + un2(nl2+1:max(nl1,nl2)) = 0._WP + un1(1:ul1-1) = 0._WP + un2(1:ul2-1) = 0._WP + + ! first edge node + ! Do not calculate on Halo nodes, as the result will not be used. + ! The "if" is cheaper than the avoided computiations. + if (nod(1) <= myDim_nod2d) then + do nz=min(ul1,ul2), max(nl1,nl2) + ! add w*du/dz+(u*du/dx+v*du/dy) & w*dv/dz+(u*dv/dx+v*dv/dy) + Unode_rhs(1,nz,nod(1)) = Unode_rhs(1,nz,nod(1)) + un1(nz)*UV(1,nz,el1) + un2(nz)*UV(1,nz,el2) + Unode_rhs(2,nz,nod(1)) = Unode_rhs(2,nz,nod(1)) + un1(nz)*UV(2,nz,el1) + un2(nz)*UV(2,nz,el2) + end do + endif + + ! second edge node + if (nod(2) <= myDim_nod2d) then + do nz=min(ul1,ul2), max(nl1,nl2) + ! add w*du/dz+(u*du/dx+v*du/dy) & w*dv/dz+(u*dv/dx+v*dv/dy) + Unode_rhs(1,nz,nod(2)) = Unode_rhs(1,nz,nod(2)) - un1(nz)*UV(1,nz,el1) - un2(nz)*UV(1,nz,el2) + Unode_rhs(2,nz,nod(2)) = Unode_rhs(2,nz,nod(2)) - un1(nz)*UV(2,nz,el1) - un2(nz)*UV(2,nz,el2) + end do + endif + + else ! el2 is not a valid element --> ed is a boundary edge, there is only the contribution from el1 + ! first edge node + if (nod(1) <= myDim_nod2d) then + do nz=ul1, nl1 + ! add w*du/dz+(u*du/dx+v*du/dy) & w*dv/dz+(u*dv/dx+v*dv/dy) + Unode_rhs(1,nz,nod(1)) = Unode_rhs(1,nz,nod(1)) + un1(nz)*UV(1,nz,el1) + Unode_rhs(2,nz,nod(1)) = Unode_rhs(2,nz,nod(1)) + un1(nz)*UV(2,nz,el1) + end do ! --> do nz=ul1, nl1 + endif + + ! second edge node + if (nod(2) <= myDim_nod2d) then + !!PS do nz=1, nl1 + do nz=ul1, nl1 + ! add w*du/dz+(u*du/dx+v*du/dy) & w*dv/dz+(u*dv/dx+v*dv/dy) + Unode_rhs(1,nz,nod(2)) = Unode_rhs(1,nz,nod(2)) - un1(nz)*UV(1,nz,el1) + Unode_rhs(2,nz,nod(2)) = Unode_rhs(2,nz,nod(2)) - un1(nz)*UV(2,nz,el1) + end do ! --> do nz=ul1, nl1 + endif + endif ! --> if (el2>0) then + end do ! --> do ed=1, myDim_edge2D -!_______________________________________________________________________________ -do n=1,myDim_nod2d - nl1 = nlevels_nod2D(n)-1 - ul1 = ulevels_nod2D(n) - - !!PS Unode_rhs(1,1:nl1,n) = Unode_rhs(1,1:nl1,n) *area_inv(1:nl1,n) - !!PS Unode_rhs(2,1:nl1,n) = Unode_rhs(2,1:nl1,n) *area_inv(1:nl1,n) - Unode_rhs(1,ul1:nl1,n) = Unode_rhs(1,ul1:nl1,n) *area_inv(ul1:nl1,n) - Unode_rhs(2,ul1:nl1,n) = Unode_rhs(2,ul1:nl1,n) *area_inv(ul1:nl1,n) -end do + !___________________________________________________________________________ + ! divide total nodal advection by scalar area + do n=1,myDim_nod2d + nl1 = nlevels_nod2D(n)-1 + ul1 = ulevels_nod2D(n) +!!PS Unode_rhs(1,ul1:nl1,n) = Unode_rhs(1,ul1:nl1,n) *area_inv(ul1:nl1,n) ! --> TEST_cavity +!!PS Unode_rhs(2,ul1:nl1,n) = Unode_rhs(2,ul1:nl1,n) *area_inv(ul1:nl1,n) ! --> TEST_cavity + Unode_rhs(1,ul1:nl1,n) = Unode_rhs(1,ul1:nl1,n) *areasvol_inv(ul1:nl1,n) + Unode_rhs(2,ul1:nl1,n) = Unode_rhs(2,ul1:nl1,n) *areasvol_inv(ul1:nl1,n) + + IF (ANY(areasvol_inv(ul1:nl1,n) < 1.e-15)) THEN + WRITE(*,*) "BLA, BLA" + write(*,*) areasvol_inv(ul1:nl1,n) + write(*,*) area(ul1:nl1,n) + CALL PAR_EX + STOP + END IF + end do !-->do n=1,myDim_nod2d -!_______________________________________________________________________________ -call exchange_nod(Unode_rhs) + !___________________________________________________________________________ + call exchange_nod(Unode_rhs) -!_______________________________________________________________________________ -do el=1, myDim_elem2D - nl1 = nlevels(el)-1 - ul1 = ulevels(el) - !!PS UV_rhsAB(1:2,1:nl1,el) = UV_rhsAB(1:2,1:nl1,el) & - !!PS + elem_area(el)*(Unode_rhs(1:2,1:nl1,elem2D_nodes(1,el)) & - !!PS + Unode_rhs(1:2,1:nl1,elem2D_nodes(2,el)) & - !!PS + Unode_rhs(1:2,1:nl1,elem2D_nodes(3,el))) / 3.0_WP - UV_rhsAB(1:2,ul1:nl1,el) = UV_rhsAB(1:2,ul1:nl1,el) & - + elem_area(el)*(Unode_rhs(1:2,ul1:nl1,elem2D_nodes(1,el)) & - + Unode_rhs(1:2,ul1:nl1,elem2D_nodes(2,el)) & - + Unode_rhs(1:2,ul1:nl1,elem2D_nodes(3,el))) / 3.0_WP - -end do + !___________________________________________________________________________ + ! convert total nodal advection from vertice --> elements + do el=1, myDim_elem2D + nl1 = nlevels(el)-1 + ul1 = ulevels(el) + UV_rhsAB(1:2,ul1:nl1,el) = UV_rhsAB(1:2,ul1:nl1,el) & + + elem_area(el)*(Unode_rhs(1:2,ul1:nl1,elem2D_nodes(1,el)) & + + Unode_rhs(1:2,ul1:nl1,elem2D_nodes(2,el)) & + + Unode_rhs(1:2,ul1:nl1,elem2D_nodes(3,el))) / 3.0_WP + + end do ! --> do el=1, myDim_elem2D end subroutine momentum_adv_scalar From 9991f35b5242ef1353eeeda481e19c86f94cc5d6 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 17 Feb 2021 17:53:21 +0100 Subject: [PATCH 06/40] improve computation of mid scalar cell area in oce_mesh.F90 --- src/oce_mesh.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index cc5a9e599..8e0fea2c1 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -2412,12 +2412,10 @@ SUBROUTINE mesh_auxiliary_arrays(mesh) do i=1, myDim_nod2D if (mesh%geo_coord_nod2D(2, i) > 0) then nn=nn+1 -!!PS mesh%lump2d_north(i)=mesh%area(1, i)! --> TEST_cavity - mesh%lump2d_north(i)=mesh%node_area(i) + mesh%lump2d_north(i)=mesh%areasvol(mesh%ulevels_nod2d(i), i) else ns=ns+1 -!!PS mesh%lump2d_south(i)=mesh%area(1, i)! --> TEST_cavity - mesh%lump2d_south(i)=mesh%node_area(i) + mesh%lump2d_south(i)=mesh%area(mesh%ulevels_nod2d(i), i) end if end do @@ -2468,7 +2466,7 @@ SUBROUTINE check_mesh_consistency(mesh) aux=0._WP do n=1, myDim_nod2D do nz=mesh%ulevels_nod2D(n), mesh%nlevels_nod2D(n)-1 - aux(nz)=aux(nz)+mesh%area(nz, n) + aux(nz)=aux(nz)+mesh%areasvol(nz, n) end do end do call MPI_AllREDUCE(aux, vol_n, mesh%nl, MPI_DOUBLE_PRECISION, MPI_SUM, & From 8eca26748a514c4819b5398081a4ad46ba6a6690 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 17 Feb 2021 17:54:26 +0100 Subject: [PATCH 07/40] exchange area_inv with areasvol_inv in subroutine momentum_adv_scalar --- src/oce_ale_vel_rhs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oce_ale_vel_rhs.F90 b/src/oce_ale_vel_rhs.F90 index 210c7f28d..1b54398ac 100644 --- a/src/oce_ale_vel_rhs.F90 +++ b/src/oce_ale_vel_rhs.F90 @@ -328,7 +328,7 @@ subroutine momentum_adv_scalar(mesh) IF (ANY(areasvol_inv(ul1:nl1,n) < 1.e-15)) THEN WRITE(*,*) "BLA, BLA" write(*,*) areasvol_inv(ul1:nl1,n) - write(*,*) area(ul1:nl1,n) + write(*,*) areasvol(ul1:nl1,n) CALL PAR_EX STOP END IF From 4ff09f40272a8dad3fe1c7ffdb0400c70ea02aa4 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 12 Mar 2021 15:52:00 +0100 Subject: [PATCH 08/40] area definition for the scalar area and for the area of the upper edge of the scalar prism. If there is no cavity these two areas are identical but in presence of cavity they can be different --- src/oce_mesh.F90 | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 8e0fea2c1..2d8324ac3 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1859,11 +1859,11 @@ SUBROUTINE mesh_areas(mesh) allocate(mesh%area(mesh%nl,myDim_nod2d+eDim_nod2D)) ! "mid" area of scalar cell in case of cavity area \= areasvol, size: nl-1 x node - allocate(mesh%areasvol(mesh%nl-1,myDim_nod2d+eDim_nod2D)) + allocate(mesh%areasvol(mesh%nl,myDim_nod2d+eDim_nod2D)) ! area inverse allocate(mesh%area_inv(mesh%nl,myDim_nod2d+eDim_nod2D)) - allocate(mesh%areasvol_inv(mesh%nl-1,myDim_nod2d+eDim_nod2D)) + allocate(mesh%areasvol_inv(mesh%nl,myDim_nod2d+eDim_nod2D)) ! resolution at nodes allocate(mesh%mesh_resolution(myDim_nod2d+eDim_nod2D)) @@ -1914,14 +1914,12 @@ SUBROUTINE mesh_areas(mesh) ! #############################| | : ! #############################|_______|___area_k+5 ! ######################################### - mesh%area = 0.0_WP - mesh%areasvol = 0.0_WP - if (use_cavity) then allocate(cavity_contribut(mesh%nl,myDim_nod2d+eDim_nod2D)) cavity_contribut = 0 end if + mesh%area = 0.0_WP do n=1, myDim_nod2D+eDim_nod2D do j=1,mesh%nod_in_elem2D_num(n) elem=mesh%nod_in_elem2D(j,n) @@ -1933,7 +1931,7 @@ SUBROUTINE mesh_areas(mesh) ! different. Directly under the cavity the area of scalar cell ! corresponds to the area of the lower edge nzmin = mesh%ulevels(elem) - nzmax = mesh%nlevels(elem) + nzmax = mesh%nlevels(elem)-1 do nz=nzmin,nzmax mesh%area(nz,n)=mesh%area(nz,n)+mesh%elem_area(elem)/3.0_WP end do @@ -1953,6 +1951,7 @@ SUBROUTINE mesh_areas(mesh) ! for cavity case: redefine "mid" scalar cell area from upper edge of prism to ! lower edge of prism if a cavity triangle is present at the upper scalar ! cell edge + mesh%areasvol = 0.0_WP if (use_cavity) then do n = 1, myDim_nod2D+eDim_nod2D nzmin = mesh%ulevels_nod2d(n) @@ -1994,12 +1993,12 @@ SUBROUTINE mesh_areas(mesh) nzmin = mesh%ulevels_nod2d(n) nzmax = mesh%nlevels_nod2d(n) do nz=nzmin,nzmax - mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) -!!PS if (mesh%area(nz,n) > 0._WP) then -!!PS mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) -!!PS else -!!PS mesh%area_inv(nz,n) = 0._WP -!!PS end if +!!PS mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) + if (mesh%area(nz,n) > 0._WP) then + mesh%area_inv(nz,n) = 1._WP/mesh%area(nz,n) + else + mesh%area_inv(nz,n) = 0._WP + end if end do end do @@ -2009,12 +2008,12 @@ SUBROUTINE mesh_areas(mesh) nzmin = mesh%ulevels_nod2d(n) nzmax = mesh%nlevels_nod2d(n)-1 do nz=nzmin,nzmax - mesh%areasvol_inv(nz,n) = 1._WP/mesh%areasvol(nz,n) -!!PS if (mesh%areasvol(nz,n) > 0._WP) then -!!PS mesh%areasvol_inv(nz,n) = 1._WP/mesh%areasvol(nz,n) -!!PS else -!!PS mesh%areasvol_inv(nz,n) = 0._WP -!!PS end if +!!PS mesh%areasvol_inv(nz,n) = 1._WP/mesh%areasvol(nz,n) + if (mesh%areasvol(nz,n) > 0._WP) then + mesh%areasvol_inv(nz,n) = 1._WP/mesh%areasvol(nz,n) + else + mesh%areasvol_inv(nz,n) = 0._WP + end if end do end do else From 83887097c8e6153b09400cb90af73ef785560ca5 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 12 Mar 2021 15:57:56 +0100 Subject: [PATCH 09/40] exchange area with areasvol where it is neccessary --- src/oce_ale.F90 | 33 ++++++++--------- src/oce_ale_tracer.F90 | 77 ++++++++++++++++++++++------------------ src/oce_ale_vel_rhs.F90 | 8 ----- src/oce_vel_rhs_vinv.F90 | 4 +-- 4 files changed, 62 insertions(+), 60 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 174697209..967c74dc2 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -574,7 +574,7 @@ subroutine init_thickness_ale(mesh) ! Fill in ssh_rhs_old !!PS ssh_rhs_old=(hbar-hbar_old)*area(1,:)/dt do n=1,myDim_nod2D+eDim_nod2D - ssh_rhs_old(n)=(hbar(n)-hbar_old(n))*area(ulevels_nod2D(n),n)/dt + ssh_rhs_old(n)=(hbar(n)-hbar_old(n))*areasvol(ulevels_nod2D(n),n)/dt ! --> TEST_cavity end do ! -->see equation (14) FESOM2:from finite elements to finie volume @@ -1260,9 +1260,12 @@ subroutine init_stiff_mat_ale(mesh) ! 2nd do first term of lhs od equation (18) of "FESOM2 from finite element to finite volumes" ! Mass matrix part do row=1, myDim_nod2D + ! if cavity no time derivative for eta in case of rigid lid approximation at + ! thee cavity-ocean interface, which means cavity-ocean interface is not allowed + ! to move vertically. + if (ulevels_nod2D(row)>1) cycle offset = ssh_stiff%rowptr(row) - !!PS SSH_stiff%values(offset) = SSH_stiff%values(offset)+ area(1,row)/dt - SSH_stiff%values(offset) = SSH_stiff%values(offset)+ area(ulevels_nod2D(row),row)/dt + SSH_stiff%values(offset) = SSH_stiff%values(offset)+ areasvol(ulevels_nod2D(row),row)/dt end do deallocate(n_pos,n_num) @@ -1551,8 +1554,7 @@ subroutine compute_ssh_rhs_ale(mesh) if ( .not. trim(which_ALE)=='linfs') then do n=1,myDim_nod2D nzmin = ulevels_nod2D(n) - ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*area(nzmin,n)+(1.0_WP-alpha)*ssh_rhs_old(n) - !!PS ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*area(1,n)+(1.0_WP-alpha)*ssh_rhs_old(n) + ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n)+(1.0_WP-alpha)*ssh_rhs_old(n) end do else do n=1,myDim_nod2D @@ -1646,7 +1648,7 @@ subroutine compute_hbar_ale(mesh) !!PS end if if (.not. trim(which_ALE)=='linfs') then do n=1,myDim_nod2D - ssh_rhs_old(n)=ssh_rhs_old(n)-water_flux(n)*area(ulevels_nod2D(n),n) + ssh_rhs_old(n)=ssh_rhs_old(n)-water_flux(n)*areasvol(ulevels_nod2D(n),n) end do call exchange_nod(ssh_rhs_old) end if @@ -1658,7 +1660,7 @@ subroutine compute_hbar_ale(mesh) !!PS call exchange_nod(hbar) hbar_old=hbar do n=1,myDim_nod2D - hbar(n)=hbar_old(n)+ssh_rhs_old(n)*dt/area(ulevels_nod2D(n),n) + hbar(n)=hbar_old(n)+ssh_rhs_old(n)*dt/areasvol(ulevels_nod2D(n),n) end do call exchange_nod(hbar) @@ -1702,7 +1704,7 @@ subroutine vert_vel_ale(mesh) use g_forcing_arrays !!PS implicit none - integer :: el(2), enodes(2), n, nz, ed, nzmin, nzmax + integer :: el(2), enodes(2), n, nz, ed, nzmin, nzmax, uln1, uln2, nln1, nln2 real(kind=WP) :: c1, c2, deltaX1, deltaY1, deltaX2, deltaY2, dd, dd1, dddt, cflmax !_______________________________ @@ -1738,7 +1740,7 @@ subroutine vert_vel_ale(mesh) ! do it with gauss-law: int( div(u_vec)*dV) = int( u_vec * n_vec * dS ) nzmin = ulevels(el(1)) nzmax = nlevels(el(1))-1 - !!PS do nz=nlevels(el(1))-1,1,-1 + do nz = nzmax, nzmin, -1 ! --> h * u_vec * n_vec ! --> e_vec = (dx,dy), n_vec = (-dy,dx); @@ -1754,7 +1756,6 @@ subroutine vert_vel_ale(mesh) fer_Wvel(nz,enodes(1))=fer_Wvel(nz,enodes(1))+c1 fer_Wvel(nz,enodes(2))=fer_Wvel(nz,enodes(2))-c1 end if - end do !_______________________________________________________________________ @@ -1764,8 +1765,8 @@ subroutine vert_vel_ale(mesh) deltaX2=edge_cross_dxdy(3,ed) deltaY2=edge_cross_dxdy(4,ed) nzmin = ulevels(el(2)) - nzmax = nlevels(el(2))-1 - !!PS do nz=nlevels(el(2))-1,1,-1 + nzmax = nlevels(el(2))-1 + do nz = nzmax, nzmin, -1 c2=-(UV(2,nz,el(2))*deltaX2 - UV(1,nz,el(2))*deltaY2)*helem(nz,el(2)) Wvel(nz,enodes(1))=Wvel(nz,enodes(1))+c2 @@ -1790,7 +1791,7 @@ subroutine vert_vel_ale(mesh) do n=1, myDim_nod2D nzmin = ulevels_nod2D(n) nzmax = nlevels_nod2d(n)-1 - !!PS do nz=nl-1,1,-1 + do nz=nzmax,nzmin,-1 Wvel(nz,n)=Wvel(nz,n)+Wvel(nz+1,n) if (Fer_GM) then @@ -1805,11 +1806,11 @@ subroutine vert_vel_ale(mesh) do n=1, myDim_nod2D nzmin = ulevels_nod2D(n) nzmax = nlevels_nod2d(n)-1 - !!PS do nz=1,nlevels_nod2D(n)-1 + do nz=nzmin,nzmax Wvel(nz,n)=Wvel(nz,n)/area(nz,n) if (Fer_GM) then - fer_Wvel(nz,n)=fer_Wvel(nz,n)/area(nz,n) + fer_Wvel(nz,n)=fer_Wvel(nz,n)/area(nz,n) end if end do end do @@ -2744,7 +2745,7 @@ subroutine oce_timestep_ale(n, mesh) !___________________________________________________________________________ ! solve tracer equation if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call solve_tracers_ale'//achar(27)//'[0m' - call solve_tracers_ale(mesh) +!!PS call solve_tracers_ale(mesh) t8=MPI_Wtime() !___________________________________________________________________________ diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 24f77d628..7e57cfbc8 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -364,9 +364,9 @@ subroutine diff_ver_part_expl_ale(tr_num, mesh) !_______________________________________________________________________ !!PS do nz=1,nl1-1 do nz=ul1,nl1-1 - del_ttf(nz,n) = del_ttf(nz,n) + (vd_flux(nz) - vd_flux(nz+1))/(zbar_3d_n(nz,n)-zbar_3d_n(nz+1,n))*dt/area(nz,n) + del_ttf(nz,n) = del_ttf(nz,n) + (vd_flux(nz) - vd_flux(nz+1))/(zbar_3d_n(nz,n)-zbar_3d_n(nz+1,n))*dt/areasvol(nz,n) end do - del_ttf(nl1,n) = del_ttf(nl1,n) + (vd_flux(nl1)/(zbar_3d_n(nl1,n)-zbar_3d_n(nl1+1,n)))*dt/area(nl1,n) + del_ttf(nl1,n) = del_ttf(nl1,n) + (vd_flux(nl1)/(zbar_3d_n(nl1,n)-zbar_3d_n(nl1+1,n)))*dt/areasvol(nl1,n) end do ! --> do n=1, myDim_nod2D end subroutine diff_ver_part_expl_ale @@ -413,7 +413,7 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! --> h^(n+0.5)* (dTnew) = dt*(K_33*d/dz*dTnew) + RHS ! --> solve for dT_new ! - ! ----------- zbar_1, V_1 (Volume eq. to Area) + ! ----------- zbar_1, V_1 (Skalar Volume), A_1 (Area of edge), no Cavity A1==V1, yes Cavity A1 !=V1 ! Z_1 o T_1 ! ----------- zbar_2, V_2 ! Z_2 o T_2 @@ -422,8 +422,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! ----------- zbar_4 ! : ! --> Difference Quotient at Volume _2: ddTnew_2/dt + d/dz*K_33 d/dz*dTnew_2 = 0 --> homogene solution - ! V2*dTnew_2 *h^(n+0.5) = -dt * [ (dTnew_1-dTnew_2)/(Z_1-Z_2)*V_2 + (dTnew_2-dTnew_3)/(Z_2-Z_3)*V_3 ] + RHS - ! dTnew_2 *h^(n+0.5) = -dt * [ (dTnew_1-dTnew_2)/(Z_1-Z_2)*V_2 + (dTnew_2-dTnew_3)/(Z_2-Z_3)*V_3/V_2 ] + RHS + ! V2*dTnew_2 *h^(n+0.5) = -dt * [ (dTnew_1-dTnew_2)/(Z_1-Z_2)*A_2 + (dTnew_2-dTnew_3)/(Z_2-Z_3)*A_3 ] + RHS + ! dTnew_2 *h^(n+0.5) = -dt * [ (dTnew_1-dTnew_2)/(Z_1-Z_2)*A_2/V_2 + (dTnew_2-dTnew_3)/(Z_2-Z_3)*A_3/V_2 ] + RHS ! | | ! v v ! diffusive flux towards diffusive flux towards @@ -432,11 +432,11 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! --> solve coefficents for homogene part ! dTnew_2 *h^(n+0.5) = -dt * [ a*dTnew_1 + b*dTnew_2 + c*dTnew_3 ] ! - ! --> a = -dt*K_33/(Z_1-Z_2) + ! --> a = -dt*K_33/(Z_1-Z_2)*A_2/V_2 ! - ! --> c = -dt*K_33/(Z_2-Z_3)*V_3/V_2 + ! --> c = -dt*K_33/(Z_2-Z_3)*A_3/V_2 ! - ! --> b = h^(n+0.5) -[ dt*K_33/(Z_1-Z_2) + dt*K_33/(Z_2-Z_3)*V_3/V_2 ] = -(a+c) + h^(n+0.5) + ! --> b = h^(n+0.5) -[ dt*K_33/(Z_1-Z_2)*A_2/V_2 + dt*K_33/(Z_2-Z_3)*A_3/V_2 ] = -(a+c) + h^(n+0.5) !___________________________________________________________________________ ! loop over local nodes @@ -490,15 +490,21 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! layer dependent coefficients for for solving dT(1)/dt+d/dz*K_33*d/dz*T(1) = ... a(nz)=0.0_WP !!PS c(nz)=-(Kv(2,n)+Ty1)*zinv2*zinv*area(nz+1,n)/area(nz,n) - c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/area(nz,n) + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=-c(nz)+hnode_new(nz,n) ! ale ! update from the vertical advection --> comes from splitting of vert ! velocity into explicite and implicite contribution if (do_wimpl) then - v_adv=zinv*area(nz+1,n)/area(nz,n) - b(nz)=b(nz)+Wvel_i(nz, n)*zinv-min(0._WP, Wvel_i(nz+1, n))*v_adv - c(nz)=c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv + !!PS v_adv =zinv*area(nz+1,n)/areasvol(nz,n) + !!PS b(nz) =b(nz)+Wvel_i(nz, n)*zinv-min(0._WP, Wvel_i(nz+1, n))*v_adv + !!PS c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv + v_adv =zinv*area(nz ,n)/areasvol(nz,n) + b(nz) =b(nz)+Wvel_i(nz, n)*v_adv + + v_adv =zinv*area(nz+1,n)/areasvol(nz,n) + b(nz) =b(nz)-min(0._WP, Wvel_i(nz+1, n))*v_adv + c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv end if ! backup zinv2 for next depth level zinv1=zinv2 @@ -518,8 +524,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) Ty =Ty *isredi Ty1=Ty1*isredi ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... - a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv - c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/area(nz,n) + a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv*area(nz ,n)/areasvol(nz,n) + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=-a(nz)-c(nz)+hnode_new(nz,n) ! backup zinv2 for next depth level @@ -527,10 +533,12 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! update from the vertical advection if (do_wimpl) then - v_adv=zinv + !!PS v_adv=zinv + v_adv=zinv*area(nz ,n)/areasvol(nz,n) a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv - v_adv=v_adv*area(nz+1,n)/area(nz,n) + !!PS v_adv=v_adv*areasvol(nz+1,n)/areasvol(nz,n) + v_adv=zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=b(nz)-min(0._WP, Wvel_i(nz+1, n))*v_adv c(nz)=c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv end if @@ -547,13 +555,14 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) (zbar_n(nz)-Z_n(nz)) *zinv1 *slope_tapered(3,nz,n)**2 *Ki(nz,n) Ty =Ty *isredi ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... - a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv + a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv*area(nz ,n)/areasvol(nz,n) c(nz)=0.0_WP b(nz)=-a(nz)+hnode_new(nz,n) ! update from the vertical advection if (do_wimpl) then - v_adv=zinv + !!PS v_adv=zinv + v_adv=zinv*area(nz ,n)/areasvol(nz,n) a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv end if @@ -613,13 +622,13 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) !_______________________________________________________________________ ! case of activated shortwave penetration into the ocean, ad 3d contribution -!!PS if (use_sw_pene .and. tr_num==1) then -!!PS !!PS do nz=1, nzmax-1 -!!PS do nz=nzmin, nzmax-1 -!!PS zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! -!!PS tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n)*area(nz+1,n)/area(nz,n))*zinv -!!PS end do -!!PS end if + if (use_sw_pene .and. tr_num==1) then + !!PS do nz=1, nzmax-1 + do nz=nzmin, nzmax-1 + zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! + tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n)*area(nz+1,n)/areasvol(nz,n))*zinv + end do + end if !_______________________________________________________________________ ! The first row contains also the boundary condition from heatflux, @@ -634,7 +643,7 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! v (+) v (+) ! !!PS tr(1)= tr(1)+bc_surface(n, tracer_id(tr_num)) -!!PS tr(nzmin)= tr(nzmin)+bc_surface(n, tracer_id(tr_num),mesh) + tr(nzmin)= tr(nzmin)+bc_surface(n, tracer_id(tr_num),mesh) !_______________________________________________________________________ ! The forward sweep algorithm to solve the three-diagonal matrix @@ -726,8 +735,8 @@ subroutine diff_ver_part_redi_expl(mesh) Ty=Ty+tr_xy(2,nz,elem)*elem_area(elem) endif end do - tr_xynodes(1,nz,n)=tx/3.0_WP/area(nz,n) - tr_xynodes(2,nz,n)=ty/3.0_WP/area(nz,n) + tr_xynodes(1,nz,n)=tx/3.0_WP/areasvol(nz,n) + tr_xynodes(2,nz,n)=ty/3.0_WP/areasvol(nz,n) end do end do @@ -756,13 +765,13 @@ subroutine diff_ver_part_redi_expl(mesh) !!PS do nz=2,nl1 do nz=ul1+1,nl1 vd_flux(nz)=(Z_n(nz-1)-zbar_n(nz))*(slope_tapered(1,nz-1,n)*tr_xynodes(1,nz-1,n)+slope_tapered(2,nz-1,n)*tr_xynodes(2,nz-1,n))*Ki(nz-1,n) - vd_flux(nz)=vd_flux(nz)+& - (zbar_n(nz)-Z_n(nz)) *(slope_tapered(1,nz,n) *tr_xynodes(1,nz,n) +slope_tapered(2,nz,n) *tr_xynodes(2,nz,n)) *Ki(nz,n) - vd_flux(nz)=vd_flux(nz)/(Z_n(nz-1)-Z_n(nz))*area(nz,n) + vd_flux(nz)=vd_flux(nz)+& + (zbar_n(nz)-Z_n(nz)) *(slope_tapered(1,nz,n) *tr_xynodes(1,nz,n) +slope_tapered(2,nz,n) *tr_xynodes(2,nz,n)) *Ki(nz,n) + vd_flux(nz)=vd_flux(nz)/(Z_n(nz-1)-Z_n(nz))*area(nz,n) enddo !!PS do nz=1,nl1 do nz=ul1,nl1 - del_ttf(nz,n) = del_ttf(nz,n)+(vd_flux(nz) - vd_flux(nz+1))*dt/area(nz,n) + del_ttf(nz,n) = del_ttf(nz,n)+(vd_flux(nz) - vd_flux(nz+1))*dt/areasvol(nz,n) enddo end do end subroutine diff_ver_part_redi_expl! @@ -908,8 +917,8 @@ subroutine diff_part_hor_redi(mesh) if (ul2>0) ul12=min(ul1,ul2) !!PS del_ttf(1:nl12,enodes(1))=del_ttf(1:nl12,enodes(1))+rhs1(1:nl12)*dt/area(1:nl12,enodes(1)) !!PS del_ttf(1:nl12,enodes(2))=del_ttf(1:nl12,enodes(2))+rhs2(1:nl12)*dt/area(1:nl12,enodes(2)) - del_ttf(ul12:nl12,enodes(1))=del_ttf(ul12:nl12,enodes(1))+rhs1(ul12:nl12)*dt/area(ul12:nl12,enodes(1)) - del_ttf(ul12:nl12,enodes(2))=del_ttf(ul12:nl12,enodes(2))+rhs2(ul12:nl12)*dt/area(ul12:nl12,enodes(2)) + del_ttf(ul12:nl12,enodes(1))=del_ttf(ul12:nl12,enodes(1))+rhs1(ul12:nl12)*dt/areasvol(ul12:nl12,enodes(1)) + del_ttf(ul12:nl12,enodes(2))=del_ttf(ul12:nl12,enodes(2))+rhs2(ul12:nl12)*dt/areasvol(ul12:nl12,enodes(2)) end do end subroutine diff_part_hor_redi diff --git a/src/oce_ale_vel_rhs.F90 b/src/oce_ale_vel_rhs.F90 index 1b54398ac..0f5b0ac6e 100644 --- a/src/oce_ale_vel_rhs.F90 +++ b/src/oce_ale_vel_rhs.F90 @@ -324,14 +324,6 @@ subroutine momentum_adv_scalar(mesh) !!PS Unode_rhs(2,ul1:nl1,n) = Unode_rhs(2,ul1:nl1,n) *area_inv(ul1:nl1,n) ! --> TEST_cavity Unode_rhs(1,ul1:nl1,n) = Unode_rhs(1,ul1:nl1,n) *areasvol_inv(ul1:nl1,n) Unode_rhs(2,ul1:nl1,n) = Unode_rhs(2,ul1:nl1,n) *areasvol_inv(ul1:nl1,n) - - IF (ANY(areasvol_inv(ul1:nl1,n) < 1.e-15)) THEN - WRITE(*,*) "BLA, BLA" - write(*,*) areasvol_inv(ul1:nl1,n) - write(*,*) areasvol(ul1:nl1,n) - CALL PAR_EX - STOP - END IF end do !-->do n=1,myDim_nod2d !___________________________________________________________________________ diff --git a/src/oce_vel_rhs_vinv.F90 b/src/oce_vel_rhs_vinv.F90 index a09e6658e..849b5aea9 100755 --- a/src/oce_vel_rhs_vinv.F90 +++ b/src/oce_vel_rhs_vinv.F90 @@ -92,7 +92,7 @@ subroutine relative_vorticity(mesh) nl1 = nlevels_nod2D(n) !!PS DO nz=1,nlevels_nod2D(n)-1 DO nz=ul1,nl1-1 - vorticity(nz,n)=vorticity(nz,n)/area(nz,n) + vorticity(nz,n)=vorticity(nz,n)/areasvol(nz,n) END DO END DO @@ -151,7 +151,7 @@ subroutine compute_vel_rhs_vinv(mesh) !vector invariant !!PS DO nz=1, nlevels_nod2D(n)-1 DO nz=nzmin, nzmax-1 !DO nz=1, nl-1 - KE_node(nz,n)=KE_node(nz,n)/(6._WP*area(nz,n)) !NR divide by 6 here + KE_node(nz,n)=KE_node(nz,n)/(6._WP*areasvol(nz,n)) !NR divide by 6 here END DO END DO From 4152b2d1b6801eb42c2035c261dd57f4f93db0c3 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 12 Mar 2021 15:58:59 +0100 Subject: [PATCH 10/40] exchange area with areasvol where it is neccessary --- src/oce_adv_tra_driver.F90 | 8 ++++---- src/oce_adv_tra_fct.F90 | 4 ++-- src/oce_adv_tra_ver.F90 | 27 ++++++++++++++++++++------- 3 files changed, 26 insertions(+), 13 deletions(-) diff --git a/src/oce_adv_tra_driver.F90 b/src/oce_adv_tra_driver.F90 index 9454a007b..988b6a090 100644 --- a/src/oce_adv_tra_driver.F90 +++ b/src/oce_adv_tra_driver.F90 @@ -117,7 +117,7 @@ subroutine do_oce_adv_tra(ttf, ttfAB, vel, w, wi, we, do_Xmoment, dttf_h, dttf_v nl1 = nlevels_nod2D(n) !!PS do nz=1, nlevels_nod2D(n)-1 do nz= nu1, nl1-1 - fct_LO(nz,n)=(ttf(nz,n)*hnode(nz,n)+(fct_LO(nz,n)+(adv_flux_ver(nz, n)-adv_flux_ver(nz+1, n)))*dt/area(nz,n))/hnode_new(nz,n) + fct_LO(nz,n)=(ttf(nz,n)*hnode(nz,n)+(fct_LO(nz,n)+(adv_flux_ver(nz, n)-adv_flux_ver(nz+1, n)))*dt/areasvol(nz,n))/hnode_new(nz,n) end do end do @@ -237,7 +237,7 @@ subroutine oce_tra_adv_flux2dtracer(dttf_h, dttf_v, flux_h, flux_v, mesh, use_lo nu1 = ulevels_nod2D(n) nl1 = nlevels_nod2D(n) do nz=nu1,nl1-1 - dttf_v(nz,n)=dttf_v(nz,n) + (flux_v(nz,n)-flux_v(nz+1,n))*dt/area(nz,n) + dttf_v(nz,n)=dttf_v(nz,n) + (flux_v(nz,n)-flux_v(nz+1,n))*dt/areasvol(nz,n) end do end do @@ -262,8 +262,8 @@ subroutine oce_tra_adv_flux2dtracer(dttf_h, dttf_v, flux_h, flux_v, mesh, use_lo !!PS do nz=1, max(nl1, nl2) do nz=nu12, nl12 - dttf_h(nz,enodes(1))=dttf_h(nz,enodes(1))+flux_h(nz,edge)*dt/area(nz,enodes(1)) - dttf_h(nz,enodes(2))=dttf_h(nz,enodes(2))-flux_h(nz,edge)*dt/area(nz,enodes(2)) + dttf_h(nz,enodes(1))=dttf_h(nz,enodes(1))+flux_h(nz,edge)*dt/areasvol(nz,enodes(1)) + dttf_h(nz,enodes(2))=dttf_h(nz,enodes(2))-flux_h(nz,edge)*dt/areasvol(nz,enodes(2)) end do end do end subroutine oce_tra_adv_flux2dtracer diff --git a/src/oce_adv_tra_fct.F90 b/src/oce_adv_tra_fct.F90 index 4af76fdf7..62db9a7d3 100644 --- a/src/oce_adv_tra_fct.F90 +++ b/src/oce_adv_tra_fct.F90 @@ -268,9 +268,9 @@ subroutine oce_tra_adv_fct(dttf_h, dttf_v, ttf, lo, adf_h, adf_v, mesh) nu1=ulevels_nod2D(n) nl1=nlevels_nod2D(n) do nz=nu1,nl1-1 - flux=fct_plus(nz,n)*dt/area(nz,n)+flux_eps + flux=fct_plus(nz,n)*dt/areasvol(nz,n)+flux_eps fct_plus(nz,n)=min(1.0_WP,fct_ttf_max(nz,n)/flux) - flux=fct_minus(nz,n)*dt/area(nz,n)-flux_eps + flux=fct_minus(nz,n)*dt/areasvol(nz,n)-flux_eps fct_minus(nz,n)=min(1.0_WP,fct_ttf_min(nz,n)/flux) end do end do diff --git a/src/oce_adv_tra_ver.F90 b/src/oce_adv_tra_ver.F90 index 8bc872575..7cdf4acbc 100644 --- a/src/oce_adv_tra_ver.F90 +++ b/src/oce_adv_tra_ver.F90 @@ -146,19 +146,28 @@ subroutine adv_tra_vert_impl(ttf, w, mesh) ! 1/dz(nz) zinv=1.0_WP*dt ! no .../(zbar(1)-zbar(2)) because of ALE + !!PS a(nz)=0.0_WP + !!PS v_adv=zinv*areasvol(nz+1,n)/areasvol(nz,n) + !!PS b(nz)= hnode_new(nz,n)+W(nz, n)*zinv-min(0._WP, W(nz+1, n))*v_adv + !!PS c(nz)=-max(0._WP, W(nz+1, n))*v_adv + a(nz)=0.0_WP - v_adv=zinv*area(nz+1,n)/area(nz,n) - b(nz)= hnode_new(nz,n)+W(nz, n)*zinv-min(0._WP, W(nz+1, n))*v_adv + v_adv=zinv*area(nz ,n)/areasvol(nz,n) + b(nz)= hnode_new(nz,n)+W(nz, n)*v_adv + + v_adv=zinv*area(nz+1,n)/areasvol(nz,n) + b(nz)= b(nz)-min(0._WP, W(nz+1, n))*v_adv c(nz)=-max(0._WP, W(nz+1, n))*v_adv !_______________________________________________________________________ ! Regular part of coefficients: --> 2nd...nl-2 layer do nz=nzmin+1, nzmax-2 ! update from the vertical advection - a(nz)=min(0._WP, W(nz, n))*zinv - b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*zinv + v_adv=zinv*area(nz ,n)/areasvol(nz,n) + a(nz)=min(0._WP, W(nz, n))*v_adv + b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*v_adv - v_adv=zinv*area(nz+1,n)/area(nz,n) + v_adv=zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=b(nz)-min(0._WP, W(nz+1, n))*v_adv c(nz)= -max(0._WP, W(nz+1, n))*v_adv end do ! --> do nz=2, nzmax-2 @@ -167,8 +176,12 @@ subroutine adv_tra_vert_impl(ttf, w, mesh) ! Regular part of coefficients: --> nl-1 layer nz=nzmax-1 ! update from the vertical advection - a(nz)= min(0._WP, W(nz, n))*zinv - b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*zinv + !!PS a(nz)= min(0._WP, W(nz, n))*zinv + !!PS b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*zinv + !!PS c(nz)=0.0_WP + v_adv=zinv*area(nz ,n)/areasvol(nz,n) + a(nz)= min(0._WP, W(nz, n))*v_adv + b(nz)=hnode_new(nz,n)+max(0._WP, W(nz, n))*v_adv c(nz)=0.0_WP !_______________________________________________________________________ From 067c582f832669a85b00862359d1bac4834dbbb8 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 12 Mar 2021 15:59:38 +0100 Subject: [PATCH 11/40] exchange area with areasvol where it is neccessary --- src/gen_support.F90 | 4 +++- src/oce_spp.F90 | 7 ++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/gen_support.F90 b/src/gen_support.F90 index eee45c814..7a5a2752e 100644 --- a/src/gen_support.F90 +++ b/src/gen_support.F90 @@ -274,6 +274,7 @@ subroutine integrate_nod_2D(data, int2D, mesh) do row=1, myDim_nod2D !!PS lval=lval+data(row)*area(1,row) lval=lval+data(row)*area(ulevels_nod2D(row),row) +!!PS lval=lval+data(row)*node_area(row) ! --> TEST_cavity end do int2D=0.0_WP @@ -300,7 +301,8 @@ subroutine integrate_nod_3D(data, int3D, mesh) do row=1, myDim_nod2D !!PS do k=1, nlevels_nod2D(row)-1 do k=ulevels_nod2D(row), nlevels_nod2D(row)-1 - lval=lval+data(k, row)*area(k,row)*hnode_new(k,row) + lval=lval+data(k, row)*area(k,row)*hnode_new(k,row) ! --> TEST_cavity +!!PS lval=lval+data(k, row)*area2(k,row)*hnode_new(k,row) end do end do int3D=0.0_WP diff --git a/src/oce_spp.F90 b/src/oce_spp.F90 index 83cc17e36..26de0ef8e 100644 --- a/src/oce_spp.F90 +++ b/src/oce_spp.F90 @@ -28,7 +28,7 @@ subroutine cal_rejected_salt(mesh) aux=rhoice/rhowat*dt do row=1, myDim_nod2d +eDim_nod2D! myDim is sufficient - if (thdgr(row)>0.0_WP) then + if (thdgr(row)>0.0_WP .and. ulevels_nod2D(row)==1) then ice_rejected_salt(row)= & (S_oc_array(row)-Sice)*thdgr(row)*aux*area(1, row) !unit: psu m3 @@ -63,6 +63,7 @@ subroutine app_rejected_salt(mesh) #include "associate_mesh.h" do row=1,myDim_nod2d+eDim_nod2D ! myDim is sufficient + if (ulevels_nod2D(row)>1) cycle if (ice_rejected_salt(row)<=0.0_WP) cycle ! do not parameterize brine rejection in regions with low salinity ! 1. it leads to further decrease of SSS @@ -92,10 +93,10 @@ subroutine app_rejected_salt(mesh) !!PS end do !!PS endif if (kml>nzmin) then - tr_arr(nzmin,row,2)=tr_arr(nzmin,row,2)-ice_rejected_salt(row)/area(1,row)/hnode(1,row) + tr_arr(nzmin,row,2)=tr_arr(nzmin,row,2)-ice_rejected_salt(row)/areasvol(1,row)/hnode(1,row) spar(nzmin+1:kml)=spar(nzmin+1:kml)/sum(spar(nzmin+1:kml)) do k=nzmin+1,kml - tr_arr(k,row,2)=tr_arr(k,row,2)+ice_rejected_salt(row)*spar(k)/area(k,row)/hnode(k,row) + tr_arr(k,row,2)=tr_arr(k,row,2)+ice_rejected_salt(row)*spar(k)/areasvol(k,row)/hnode(k,row) end do endif endif From 69f71684ad643d335e745689b798c41c864204bf Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 12 Mar 2021 18:14:37 +0100 Subject: [PATCH 12/40] exchange area with areasvol where it is neccessary --- src/gen_modules_cvmix_idemix.F90 | 26 ++++++++++++++------------ src/gen_modules_diag.F90 | 4 ++-- src/gen_support.F90 | 6 ++---- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/gen_modules_cvmix_idemix.F90 b/src/gen_modules_cvmix_idemix.F90 index 98929506d..e8203868b 100644 --- a/src/gen_modules_cvmix_idemix.F90 +++ b/src/gen_modules_cvmix_idemix.F90 @@ -272,7 +272,7 @@ subroutine calc_cvmix_idemix(mesh) implicit none type(t_mesh), intent(in), target :: mesh integer :: node, elem, edge, node_size - integer :: nz, nln, nl1, nl2, nl12, nu1, nu2, nu12 + integer :: nz, nln, nl1, nl2, nl12, nu1, nu2, nu12, uln integer :: elnodes1(3), elnodes2(3), el(2), ednodes(2) real(kind=WP) :: dz_trr(mesh%nl), dz_trr2(mesh%nl), bvfreq2(mesh%nl), vflux, dz_el, aux, cflfac real(kind=WP) :: grad_v0Eiw(2), deltaX1, deltaY1, deltaX2, deltaY2 @@ -286,19 +286,20 @@ subroutine calc_cvmix_idemix(mesh) node_size = myDim_nod2D do node = 1,node_size nln = nlevels_nod2D(node)-1 + uln = ulevels_nod2D(node) !___________________________________________________________________ ! calculate for TKE square of Brünt-Väisälä frequency, be aware that ! bvfreq contains already the squared brünt Väisälä frequency ... bvfreq2 = 0.0_WP - bvfreq2(2:nln) = bvfreq(2:nln,node) + bvfreq2(uln:nln) = bvfreq(uln:nln,node) !___________________________________________________________________ ! dz_trr distance between tracer points, surface and bottom dz_trr is half ! the layerthickness ... dz_trr = 0.0_WP - dz_trr(2:nln) = abs(Z_3d_n(1:nln-1,node)-Z_3d_n(2:nln,node)) - dz_trr(1) = hnode(1,node)/2.0_WP + dz_trr(uln+1:nln) = abs(Z_3d_n(uln:nln-1,node)-Z_3d_n(uln+1:nln,node)) + dz_trr(uln) = hnode(uln,node)/2.0_WP dz_trr(nln+1) = hnode(nln,node)/2.0_WP !___________________________________________________________________ @@ -389,23 +390,24 @@ subroutine calc_cvmix_idemix(mesh) ! number of above bottom levels at node nln = nlevels_nod2D(node)-1 + uln = ulevels_nod2D(node) ! thickness of mid-level to mid-level interface at node dz_trr = 0.0_WP - dz_trr(2:nln) = Z_3d_n(1:nln-1,node)-Z_3d_n(2:nln,node) - dz_trr(1) = hnode(1,node)/2.0_WP + dz_trr(uln+1:nln) = Z_3d_n(uln:nln-1,node)-Z_3d_n(uln+1:nln,node) + dz_trr(uln) = hnode(uln,node)/2.0_WP dz_trr(nln+1) = hnode(nln,node)/2.0_WP ! surface cell - vol_wcelli(1,node) = 1/(area(1,node)*dz_trr(1)) - aux = sqrt(cflfac*(area(1,node)/pi*4.0_WP)/(idemix_tau_h*dt/idemix_n_hor_iwe_prop_iter)) - iwe_v0(1,node) = min(iwe_v0(1,node),aux) + vol_wcelli(uln,node) = 1/(areasvol(uln,node)*dz_trr(uln)) + aux = sqrt(cflfac*(area(uln,node)/pi*4.0_WP)/(idemix_tau_h*dt/idemix_n_hor_iwe_prop_iter)) + iwe_v0(uln,node) = min(iwe_v0(uln,node),aux) ! bulk cells !!PS do nz=2,nln - do nz=ulevels_nod2D(node)+1,nln + do nz=uln+1,nln ! inverse volumne - vol_wcelli(nz,node) = 1/(area(nz-1,node)*dz_trr(nz)) + vol_wcelli(nz,node) = 1/(areasvol(nz-1,node)*dz_trr(nz)) ! restrict iwe_v0 aux = sqrt(cflfac*(area(nz-1,node)/pi*4.0_WP)/(idemix_tau_h*dt/idemix_n_hor_iwe_prop_iter)) @@ -415,7 +417,7 @@ subroutine calc_cvmix_idemix(mesh) end do ! bottom cell - vol_wcelli(nln+1,node) = 1/(area(nln,node)*dz_trr(nln+1)) + vol_wcelli(nln+1,node) = 1/(areasvol(nln,node)*dz_trr(nln+1)) aux = sqrt(cflfac*(area(nln,node)/pi*4.0_WP)/(idemix_tau_h*dt/idemix_n_hor_iwe_prop_iter)) iwe_v0(nln+1,node) = min(iwe_v0(nln+1,node),aux) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index c1e1ec3d8..a8477094e 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -135,7 +135,7 @@ subroutine diag_curl_stress_surf(mode, mesh) END DO DO n=1, myDim_nod2D+eDim_nod2D !!PS curl_stress_surf(n)=curl_stress_surf(n)/area(1,n) - curl_stress_surf(n)=curl_stress_surf(n)/area(ulevels_nod2D(n),n) + curl_stress_surf(n)=curl_stress_surf(n)/areasvol(ulevels_nod2D(n),n) END DO end subroutine diag_curl_stress_surf ! ============================================================== @@ -210,7 +210,7 @@ subroutine diag_curl_vel3(mode, mesh) DO n=1, myDim_nod2D !!PS DO nz=1, nlevels_nod2D(n)-1 DO nz=ulevels_nod2D(n), nlevels_nod2D(n)-1 - curl_vel3(nz,n)=curl_vel3(nz,n)/area(nz,n) + curl_vel3(nz,n)=curl_vel3(nz,n)/areasvol(nz,n) END DO END DO end subroutine diag_curl_vel3 diff --git a/src/gen_support.F90 b/src/gen_support.F90 index 7a5a2752e..bb9e24ad6 100644 --- a/src/gen_support.F90 +++ b/src/gen_support.F90 @@ -273,8 +273,7 @@ subroutine integrate_nod_2D(data, int2D, mesh) lval=0.0_WP do row=1, myDim_nod2D !!PS lval=lval+data(row)*area(1,row) - lval=lval+data(row)*area(ulevels_nod2D(row),row) -!!PS lval=lval+data(row)*node_area(row) ! --> TEST_cavity + lval=lval+data(row)*areasvol(ulevels_nod2D(row),row) end do int2D=0.0_WP @@ -301,8 +300,7 @@ subroutine integrate_nod_3D(data, int3D, mesh) do row=1, myDim_nod2D !!PS do k=1, nlevels_nod2D(row)-1 do k=ulevels_nod2D(row), nlevels_nod2D(row)-1 - lval=lval+data(k, row)*area(k,row)*hnode_new(k,row) ! --> TEST_cavity -!!PS lval=lval+data(k, row)*area2(k,row)*hnode_new(k,row) + lval=lval+data(k, row)*areasvol(k,row)*hnode_new(k,row) ! --> TEST_cavity end do end do int3D=0.0_WP From 452f597b7313649151671a8e46cb94366a5674ca Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 12 Mar 2021 18:16:26 +0100 Subject: [PATCH 13/40] only print info where is cavity in write_step_info.F90 --- src/write_step_info.F90 | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/write_step_info.F90 b/src/write_step_info.F90 index d3d7d47fd..3be6a08cd 100644 --- a/src/write_step_info.F90 +++ b/src/write_step_info.F90 @@ -268,29 +268,29 @@ subroutine check_blowup(istep, mesh) write(*,*) write(*,*) 'm_ice = ',m_ice(n),', m_ice_old = ',m_ice_old(n) write(*,*) 'a_ice = ',a_ice(n),', a_ice_old = ',a_ice_old(n) - write(*,*) 'thdgr = ',thdgr(n),', thdgr_old = ',thdgr_old(n) - write(*,*) 'thdgrsn = ',thdgrsn(n) +!!PS write(*,*) 'thdgr = ',thdgr(n),', thdgr_old = ',thdgr_old(n) +!!PS write(*,*) 'thdgrsn = ',thdgrsn(n) write(*,*) - if (lcurt_stress_surf) then - write(*,*) 'curl_stress_surf = ',curl_stress_surf(n) - write(*,*) - endif - do el=1,nod_in_elem2d_num(n) - elidx = nod_in_elem2D(el,n) - write(*,*) ' elem#=',el,', elemidx=',elidx - write(*,*) ' pgf_x =',pgf_x(:,elidx) - write(*,*) ' pgf_y =',pgf_y(:,elidx) -! write(*,*) ' U =',UV(1,:,elidx) -! write(*,*) ' V =',UV(2,:,elidx) - write(*,*) - enddo - write(*,*) 'Wvel(1, n) = ',Wvel(1,n) - write(*,*) 'Wvel(:, n) = ',Wvel(:,n) +!!PS if (lcurt_stress_surf) then +!!PS write(*,*) 'curl_stress_surf = ',curl_stress_surf(n) +!!PS write(*,*) +!!PS endif +!!PS do el=1,nod_in_elem2d_num(n) +!!PS elidx = nod_in_elem2D(el,n) +!!PS write(*,*) ' elem#=',el,', elemidx=',elidx +!!PS write(*,*) ' pgf_x =',pgf_x(:,elidx) +!!PS write(*,*) ' pgf_y =',pgf_y(:,elidx) +!!PS ! write(*,*) ' U =',UV(1,:,elidx) +!!PS ! write(*,*) ' V =',UV(2,:,elidx) +!!PS write(*,*) +!!PS enddo +!!PS write(*,*) 'Wvel(1, n) = ',Wvel(,n) + write(*,*) 'Wvel(:, n) = ',Wvel(ulevels_nod2D(n):nlevels_nod2D(n),n) write(*,*) - write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) + write(*,*) 'CFL_z(:,n) = ',CFL_z(ulevels_nod2D(n):nlevels_nod2D(n),n) write(*,*) - write(*,*) 'hnode(1, n) = ',hnode(1,n) - write(*,*) 'hnode(:, n) = ',hnode(:,n) +!!PS write(*,*) 'hnode(1, n) = ',hnode(1,n) + write(*,*) 'hnode(:, n) = ',hnode(ulevels_nod2D(n):nlevels_nod2D(n),n) write(*,*) endif From d4f351ee6a8b2e3826a7eb613a036935df0b03cd Mon Sep 17 00:00:00 2001 From: Patrick Date: Sat, 13 Mar 2021 18:32:47 +0100 Subject: [PATCH 14/40] switch fluxes back on --- src/oce_ale.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 967c74dc2..360243118 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -2555,10 +2555,10 @@ subroutine oce_timestep_ale(n, mesh) t0=MPI_Wtime() - water_flux = 0.0_WP - heat_flux = 0.0_WP - stress_surf= 0.0_WP - stress_node_surf= 0.0_WP +!!PS water_flux = 0.0_WP +!!PS heat_flux = 0.0_WP +!!PS stress_surf= 0.0_WP +!!PS stress_node_surf= 0.0_WP !___________________________________________________________________________ ! calculate equation of state, density, pressure and mixed layer depths @@ -2745,7 +2745,7 @@ subroutine oce_timestep_ale(n, mesh) !___________________________________________________________________________ ! solve tracer equation if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call solve_tracers_ale'//achar(27)//'[0m' -!!PS call solve_tracers_ale(mesh) + call solve_tracers_ale(mesh) t8=MPI_Wtime() !___________________________________________________________________________ From 0cead93d603a4cf05d7b3f317105a93cd657199c Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 15 Mar 2021 12:40:25 +0100 Subject: [PATCH 15/40] add routine to compute total vertice/elemental ocean volume --- src/oce_mesh.F90 | 72 ++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 66 insertions(+), 6 deletions(-) diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 2d8324ac3..7cd5cf5c5 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -90,7 +90,14 @@ subroutine find_levels_min_e2n(mesh) end subroutine end interface end module - +module check_total_volume_interface + interface + subroutine check_total_volume(mesh) + use mod_mesh + type(t_mesh), intent(inout) , target :: mesh + end subroutine + end interface +end module ! Driving routine. The distributed mesh information and mesh proper ! are read from files. @@ -2440,9 +2447,9 @@ SUBROUTINE mesh_auxiliary_arrays(mesh) endif END SUBROUTINE mesh_auxiliary_arrays - -!=================================================================== - +! +! +!_______________________________________________________________________________ SUBROUTINE check_mesh_consistency(mesh) USE MOD_MESH USE o_PARAM @@ -2475,7 +2482,7 @@ SUBROUTINE check_mesh_consistency(mesh) do elem=1, myDim_elem2D elnodes=mesh%elem2D_nodes(:, elem) if (elnodes(1) > myDim_nod2D) CYCLE - do nz=mesh%ulevels(elem), mesh%nlevels(elem) + do nz=mesh%ulevels(elem), mesh%nlevels(elem)-1 aux(nz)=aux(nz)+mesh%elem_area(elem) end do end do @@ -2493,4 +2500,57 @@ SUBROUTINE check_mesh_consistency(mesh) !call par_ex !stop END SUBROUTINE check_mesh_consistency -!================================================================== +! +! +!_______________________________________________________________________________ +subroutine check_total_volume(mesh) + USE MOD_MESH + USE o_PARAM + USE g_PARSUP + use g_comm_auto + use o_ARRAYS + + IMPLICIT NONE + type(t_mesh), intent(inout), target :: mesh + integer :: nz, n, elem , elnodes(3) + real(kind=WP) :: vol_n, vol_e, aux + +#include "associate_mesh.h" + + !___________________________________________________________________________ + vol_n=0._WP + vol_e=0._WP + !___________________________________________________________________________ + ! total ocean volume on nodes + aux=0._WP + do n=1, myDim_nod2D + do nz=ulevels_nod2D(n), nlevels_nod2D(n)-1 + aux=aux+areasvol(nz, n)*hnode(nz,n) + end do + end do + call MPI_AllREDUCE(aux, vol_n, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + !___________________________________________________________________________ + ! total ocean volume on elements + aux=0._WP + do elem=1, myDim_elem2D + elnodes=elem2D_nodes(:, elem) + if (elnodes(1) > myDim_nod2D) cycle + do nz=ulevels(elem), nlevels(elem)-1 + aux=aux+elem_area(elem)*helem(nz,elem) + end do + end do + call MPI_AllREDUCE(aux, vol_e, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + + !___write mesh statistics___________________________________________________ + if (mype==0) then + write(*,*) '____________________________________________________________________' + write(*,*) ' --> ocean volume check:', mype + write(*,*) ' > Total ocean volume node:', vol_n, ' m^3' + write(*,*) ' > Total ocean volume elem:', vol_e, ' m^3' + + end if + +end subroutine check_total_volume +! +! +!_______________________________________________________________________________ From ade88b707453f52bd8c1161fef0112ede0bca55c Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 15 Mar 2021 12:43:09 +0100 Subject: [PATCH 16/40] check for total ocean volume in init_thickness_ale --- src/oce_ale.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 360243118..580f36e4d 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -790,6 +790,8 @@ subroutine init_thickness_ale(mesh) !___________________________________________________________________________ hnode_new=hnode ! Should be initialized, because only variable part is updated. + !!PS call check_total_volume(mesh) + end subroutine init_thickness_ale ! ! From 8144937bec32aa3088fbb95403cd7538d41e7064 Mon Sep 17 00:00:00 2001 From: Patrick Date: Sat, 20 Mar 2021 16:33:13 +0100 Subject: [PATCH 17/40] exclude ice operations on cavity vertices and elements --- src/ice_EVP.F90 | 4 ++-- src/ice_fct.F90 | 53 ++++++++++++++++++++++++------------------ src/ice_maEVP.F90 | 24 +++++++++++-------- src/ice_setup_step.F90 | 8 +++++-- src/ice_thermo_oce.F90 | 39 ++++++++++++++++++++----------- 5 files changed, 78 insertions(+), 50 deletions(-) diff --git a/src/ice_EVP.F90 b/src/ice_EVP.F90 index e9357e80d..f5267bd69 100755 --- a/src/ice_EVP.F90 +++ b/src/ice_EVP.F90 @@ -458,7 +458,7 @@ subroutine EVPdynamics(mesh) elnodes = elem2D_nodes(:,el) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle !_______________________________________________________________________ @@ -508,7 +508,7 @@ subroutine EVPdynamics(mesh) elnodes = elem2D_nodes(:,el) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle !_______________________________________________________________________ diff --git a/src/ice_fct.F90 b/src/ice_fct.F90 index e82f2e11e..501b9ac36 100755 --- a/src/ice_fct.F90 +++ b/src/ice_fct.F90 @@ -62,11 +62,13 @@ subroutine ice_TG_rhs(mesh) ! Velocities at nodes do elem=1,myDim_elem2D !assembling rhs over elements + elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if cavity element skip it if (ulevels(elem)>1) cycle + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + - elnodes=elem2D_nodes(:,elem) !derivatives dx=gradient_sca(1:3,elem) dy=gradient_sca(4:6,elem) @@ -80,6 +82,7 @@ subroutine ice_TG_rhs(mesh) diff=ice_diff*sqrt(elem_area(elem)/scale_area) DO n=1,3 row=elnodes(n) +!!PS if (ulevels_nod2D(row)>1) cycle DO q = 1,3 !entries(q)= vol*dt*((dx(n)*um+dy(n)*vm)/3.0_WP - & ! diff*(dx(n)*dx(q)+ dy(n)*dy(q))- & @@ -372,7 +375,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) !_______________________________________________________________________ ! if cavity cycle over - !!PS if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 !_______________________________________________________________________ @@ -421,6 +424,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) !========================== if (tr_array_id==1) then do row=1, myDim_nod2D + if (ulevels_nod2d(row)>1) cycle n=nn_num(row) tmax(row)=maxval(m_icel(nn_pos(1:n,row))) tmin(row)=minval(m_icel(nn_pos(1:n,row))) @@ -432,6 +436,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) if (tr_array_id==2) then do row=1, myDim_nod2D + if (ulevels_nod2d(row)>1) cycle n=nn_num(row) tmax(row)=maxval(a_icel(nn_pos(1:n,row))) tmin(row)=minval(a_icel(nn_pos(1:n,row))) @@ -443,6 +448,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) if (tr_array_id==3) then do row=1, myDim_nod2D + if (ulevels_nod2d(row)>1) cycle n=nn_num(row) tmax(row)=maxval(m_snowl(nn_pos(1:n,row))) tmin(row)=minval(m_snowl(nn_pos(1:n,row))) @@ -455,6 +461,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) #if defined (__oifs) if (tr_array_id==4) then do row=1, myDim_nod2D + if (ulevels_nod2d(row)>1) cycle n=nn_num(row) tmax(row)=maxval(m_templ(nn_pos(1:n,row))) tmin(row)=minval(m_templ(nn_pos(1:n,row))) @@ -476,7 +483,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) !_______________________________________________________________________ ! if cavity cycle over - !!PS if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 !_______________________________________________________________________ @@ -525,7 +532,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) !_______________________________________________________________________ ! if cavity cycle over - !!PS if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 !_______________________________________________________________________ @@ -548,12 +555,13 @@ subroutine ice_fem_fct(tr_array_id, mesh) m_ice(n)=m_icel(n) end do do elem=1, myDim_elem2D + elnodes=elem2D_nodes(:,elem) + !___________________________________________________________________ ! if cavity cycle over - !PS if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 - elnodes=elem2D_nodes(:,elem) do q=1,3 n=elnodes(q) m_ice(n)=m_ice(n)+icefluxes(elem,q) @@ -567,12 +575,13 @@ subroutine ice_fem_fct(tr_array_id, mesh) a_ice(n)=a_icel(n) end do do elem=1, myDim_elem2D + elnodes=elem2D_nodes(:,elem) + !___________________________________________________________________ ! if cavity cycle over - !!PS if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 - elnodes=elem2D_nodes(:,elem) do q=1,3 n=elnodes(q) a_ice(n)=a_ice(n)+icefluxes(elem,q) @@ -587,14 +596,15 @@ subroutine ice_fem_fct(tr_array_id, mesh) end do do elem=1, myDim_elem2D elnodes=elem2D_nodes(:,elem) + !___________________________________________________________________ ! if cavity cycle over - !!PS if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 do q=1,3 - n=elnodes(q) - m_snow(n)=m_snow(n)+icefluxes(elem,q) + n=elnodes(q) + m_snow(n)=m_snow(n)+icefluxes(elem,q) end do end do end if @@ -602,12 +612,14 @@ subroutine ice_fem_fct(tr_array_id, mesh) #if defined (__oifs) if(tr_array_id==4) then do n=1,myDim_nod2D + if(ulevels_nod2D(n)>1) cycle !LK89140 ice_temp(n)=m_templ(n) end do do elem=1, myDim_elem2D elnodes=elem2D_nodes(:,elem) !___________________________________________________________________ ! if cavity cycle over + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 do q=1,3 @@ -656,17 +668,11 @@ SUBROUTINE ice_mass_matrix_fill(mesh) DO elem=1,myDim_elem2D elnodes=elem2D_nodes(:,elem) - !___________________________________________________________________ - ! if cavity cycle over - if(ulevels(elem)>1) cycle - + !_______________________________________________________________________ do n=1,3 row=elnodes(n) if(row>myDim_nod2D) cycle !___________________________________________________________________ - ! if node is cavity cycle over - if(ulevels_nod2d(row)>1) cycle - ! Global-to-local neighbourhood correspondence DO q=1,nn_num(row) col_pos(nn_pos(q,row))=q @@ -674,9 +680,9 @@ SUBROUTINE ice_mass_matrix_fill(mesh) offset=ssh_stiff%rowptr(row)-ssh_stiff%rowptr(1) DO q=1,3 col=elnodes(q) - !___________________________________________________________________ - ! if node is cavity cycle over - if(ulevels_nod2d(col)>1) cycle + !_______________________________________________________________ + ! if element is cavity cycle over + if(ulevels(elem)>1) cycle ipos=offset+col_pos(col) mass_matrix(ipos)=mass_matrix(ipos)+elem_area(elem)/12.0_WP @@ -749,12 +755,12 @@ subroutine ice_TG_rhs_div(mesh) #endif /* (__oifs) */ END DO do elem=1,myDim_elem2D !assembling rhs over elements + elnodes=elem2D_nodes(:,elem) !___________________________________________________________________________ ! if cavity element skip it if (ulevels(elem)>1) cycle + if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 - !! elem=myList_elem2D(m) - elnodes=elem2D_nodes(:,elem) !derivatives dx=gradient_sca(1:3,elem) dy=gradient_sca(4:6,elem) @@ -769,6 +775,7 @@ subroutine ice_TG_rhs_div(mesh) c4=sum(dx*u_ice(elnodes)+dy*v_ice(elnodes)) DO n=1,3 row=elnodes(n) +!!PS if(ulevels_nod2D(row)>1) cycle !LK89140 DO q = 1,3 entries(q)= vol*ice_dt*((1.0_WP-0.5_WP*ice_dt*c4)*(dx(n)*(um+u_ice(elnodes(q)))+ & dy(n)*(vm+v_ice(elnodes(q))))/12.0_WP - & diff --git a/src/ice_maEVP.F90 b/src/ice_maEVP.F90 index b10d4ff11..cfe1b0bac 100644 --- a/src/ice_maEVP.F90 +++ b/src/ice_maEVP.F90 @@ -61,7 +61,7 @@ subroutine stress_tensor_m(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle msum=sum(m_ice(elnodes))*val3 @@ -147,7 +147,7 @@ subroutine ssh2rhs(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle !_______________________________________________________________________ @@ -174,7 +174,7 @@ subroutine ssh2rhs(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle vol=elem_area(elem) @@ -224,7 +224,7 @@ subroutine stress2rhs_m(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle if(sum(a_ice(elnodes)) < 0.01_WP) cycle !DS @@ -325,7 +325,7 @@ subroutine EVPdynamics_m(mesh) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle !_______________________________________________________________________ @@ -354,7 +354,7 @@ subroutine EVPdynamics_m(mesh) elnodes=elem2D_nodes(:,el) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle vol=elem_area(el) @@ -401,7 +401,7 @@ subroutine EVPdynamics_m(mesh) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle msum=sum(m_ice(elnodes))*val3 @@ -434,7 +434,7 @@ subroutine EVPdynamics_m(mesh) do el=1,myDim_elem2D !__________________________________________________________________________ if (ulevels(el)>1) cycle - + if ( any(ulevels_nod2d(elnodes)>1) ) cycle !__________________________________________________________________________ if(ice_el(el)) then @@ -580,7 +580,7 @@ subroutine find_alpha_field_a(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle msum=sum(m_ice(elnodes))*val3 @@ -647,7 +647,7 @@ subroutine stress_tensor_a(mesh) do elem=1,myDim_elem2D !__________________________________________________________________________ ! if element has any cavity node skip it - !!PS if ( any(ulevels_nod2d(elnodes)>1) ) cycle + if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle !__________________________________________________________________________ @@ -798,6 +798,10 @@ subroutine find_beta_field_a(mesh) #include "associate_mesh.h" DO n=1, myDim_nod2D + !_______________________________________________________________________ + ! if element has any cavity node skip it + if (ulevels_nod2d(n)>1) cycle + ! ============== ! FESOM1.4 and stand-alone FESIM ! beta_evp_array(n) = maxval(alpha_evp_array(nod_in_elem2D(n)%addresses(1:nod_in_elem2D(n)%nmb))) diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90 index 50a8f4cc3..ff9ea9944 100755 --- a/src/ice_setup_step.F90 +++ b/src/ice_setup_step.F90 @@ -216,7 +216,7 @@ subroutine ice_timestep(step, mesh) end do #endif /* (__oifs) */ if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call cut_off...'//achar(27)//'[0m' - call cut_off + call cut_off(mesh) if (use_cavity) call cavity_ice_clean_ma(mesh) t2=MPI_Wtime() @@ -270,7 +270,11 @@ subroutine ice_initial_state(mesh) do i=1,myDim_nod2D+eDim_nod2D !_______________________________________________________________________ - if (ulevels_nod2d(i)>1) cycle ! --> if cavity, no sea ice, no initial state + if (ulevels_nod2d(i)>1) then + !!PS m_ice(i) = 1.0e15_WP + !!PS m_snow(i) = 0.1e15_WP + cycle ! --> if cavity, no sea ice, no initial state + endif !_______________________________________________________________________ if (tr_arr(1,i,1)< 0.0_WP) then diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index a6fef3ea1..63e0f9074 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -1,13 +1,19 @@ !=================================================================== -subroutine cut_off() -use o_param -use i_arrays -implicit none +subroutine cut_off(mesh) + use o_param + use i_arrays + use MOD_MESH + use g_parsup + implicit none + type(t_mesh), intent(in) , target :: mesh -where(a_ice>1.0_WP) - a_ice=1.0_WP -end where +#include "associate_mesh.h" +!_______________________________________________________________________________ +! lower cutoff: a_ice +where(a_ice>1.0_WP) a_ice=1.0_WP + +! upper cutoff: a_ice where(a_ice<0.1e-8_WP) a_ice=0.0_WP #if defined (__oifs) @@ -17,6 +23,8 @@ subroutine cut_off() #endif /* (__oifs) */ end where +!_______________________________________________________________________________ +! lower cutoff: m_ice where(m_ice<0.1e-8_WP) m_ice=0.0_WP #if defined (__oifs) @@ -25,17 +33,22 @@ subroutine cut_off() ice_temp=273.15_WP #endif /* (__oifs) */ end where +! upper cutoff: m_ice +where(m_ice>10.0_WP .and. ulevels_nod2d==1) m_ice=10.0_WP + +!_______________________________________________________________________________ +! lower cutoff: m_snow +where(m_snow<0.1e-8_WP) m_snow=0.0_WP +! upper cutoff: m_snow +where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP +!_______________________________________________________________________________ #if defined (__oifs) -where(ice_temp>273.15_WP) - ice_temp=273.15_WP -end where +where(ice_temp>273.15_WP) ice_temp=273.15_WP #endif /* (__oifs) */ #if defined (__oifs) -where(ice_temp < 173.15_WP .and. a_ice >= 0.1e-8_WP) - ice_temp=271.35_WP -end where +where(ice_temp < 173.15_WP .and. a_ice >= 0.1e-8_WP) ice_temp=271.35_WP #endif /* (__oifs) */ end subroutine cut_off From 92cb917d85349246a0ca23d53b0297ec545d3ab5 Mon Sep 17 00:00:00 2001 From: Patrick Date: Sat, 20 Mar 2021 16:40:36 +0100 Subject: [PATCH 18/40] ommit term from ssh_rhs_old when computing ssh_rhs, after talking to sergey --> leads to rigid lid approximation under the cavity --> no moving surface possible --- src/oce_ale.F90 | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 580f36e4d..1498b2bd8 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1210,7 +1210,7 @@ subroutine init_stiff_mat_ale(mesh) if (el(i)<1) cycle ! if el(i)<1, it means its an outer boundary edge this ! has only one triangle element to which it contribute - + ! which three nodes span up triangle el(i) ! elnodes ... node indices elnodes=elem2D_nodes(:,el(i)) @@ -1555,11 +1555,13 @@ subroutine compute_ssh_rhs_ale(mesh) ! shown in eq (11) rhs of "FESOM2: from finite elements to finte volumes, S. Danilov..." eq. (11) rhs if ( .not. trim(which_ALE)=='linfs') then do n=1,myDim_nod2D + if (ulevels_nod2D(n)>1) cycle nzmin = ulevels_nod2D(n) ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n)+(1.0_WP-alpha)*ssh_rhs_old(n) end do else do n=1,myDim_nod2D + if (ulevels_nod2D(n)>1) cycle ssh_rhs(n)=ssh_rhs(n)+(1.0_WP-alpha)*ssh_rhs_old(n) end do end if @@ -1639,15 +1641,10 @@ subroutine compute_hbar_ale(mesh) !_______________________________________________________________________ ssh_rhs_old(enodes(1))=ssh_rhs_old(enodes(1))+(c1+c2) ssh_rhs_old(enodes(2))=ssh_rhs_old(enodes(2))-(c1+c2) - end do !___________________________________________________________________________ ! take into account water flux -!!PS if (.not. trim(which_ALE)=='linfs') then -!!PS ssh_rhs_old(1:myDim_nod2D)=ssh_rhs_old(1:myDim_nod2D)-water_flux(1:myDim_nod2D)*area(1,1:myDim_nod2D) -!!PS call exchange_nod(ssh_rhs_old) -!!PS end if if (.not. trim(which_ALE)=='linfs') then do n=1,myDim_nod2D ssh_rhs_old(n)=ssh_rhs_old(n)-water_flux(n)*areasvol(ulevels_nod2D(n),n) @@ -1657,9 +1654,6 @@ subroutine compute_hbar_ale(mesh) !___________________________________________________________________________ ! update the thickness -!!PS hbar_old=hbar -!!PS hbar(1:myDim_nod2D)=hbar_old(1:myDim_nod2D)+ssh_rhs_old(1:myDim_nod2D)*dt/area(1,1:myDim_nod2D) -!!PS call exchange_nod(hbar) hbar_old=hbar do n=1,myDim_nod2D hbar(n)=hbar_old(n)+ssh_rhs_old(n)*dt/areasvol(ulevels_nod2D(n),n) @@ -2669,6 +2663,11 @@ subroutine oce_timestep_ale(n, mesh) !___________________________________________________________________________ if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_vel_rhs'//achar(27)//'[0m' +!!PS if (any(UV_rhs/=UV_rhs)) write(*,*) mype,' --> found NaN UV_rhs before compute_vel_rhs' +!!PS if (any(UV/=UV)) write(*,*) mype,' --> found NaN UV before compute_vel_rhs' +!!PS if (any(ssh_rhs/=ssh_rhs))write(*,*) mype,' --> found NaN ssh_rhs before compute_ssh_rhs_ale' +!!PS if (any(ssh_rhs_old/=ssh_rhs_old))write(*,*) mype,' --> found NaN ssh_rhs_old before compute_ssh_rhs_ale' + if(mom_adv/=3) then call compute_vel_rhs(mesh) else @@ -2676,7 +2675,6 @@ subroutine oce_timestep_ale(n, mesh) end if !___________________________________________________________________________ - if (any(UV_rhs/=UV_rhs)) write(*,*) ' --> found NaN UV_rhs MARK 2' call viscosity_filter(visc_option, mesh) !___________________________________________________________________________ @@ -2717,7 +2715,7 @@ subroutine oce_timestep_ale(n, mesh) ! Current dynamic elevation alpha*hbar(n+1/2)+(1-alpha)*hbar(n-1/2) ! equation (14) Danlov et.al "the finite volume sea ice ocean model FESOM2 ! ...if we do it here we don't need to write hbar_old into a restart file... - eta_n=alpha*hbar+(1.0_WP-alpha)*hbar_old + where(ulevels_nod2D==1) eta_n=alpha*hbar+(1.0_WP-alpha)*hbar_old ! --> eta_(n) ! call zero_dynamics !DS, zeros several dynamical variables; to be used for testing new implementations! From 98ecb193df9e339bf2e5a198febb9fc1f8f43142 Mon Sep 17 00:00:00 2001 From: Patrick Date: Sat, 20 Mar 2021 16:42:39 +0100 Subject: [PATCH 19/40] add also d_eta = NaN to blowup check --- src/write_step_info.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/write_step_info.F90 b/src/write_step_info.F90 index 3be6a08cd..35a5ba69d 100644 --- a/src/write_step_info.F90 +++ b/src/write_step_info.F90 @@ -237,7 +237,8 @@ subroutine check_blowup(istep, mesh) !___________________________________________________________________ ! check ssh if ( ((eta_n(n) /= eta_n(n)) .or. & - eta_n(n)<-50.0 .or. eta_n(n)>50.0)) then + eta_n(n)<-50.0 .or. eta_n(n)>50.0 .or. & + (d_eta(n) /= d_eta(n)) ) ) then !!PS eta_n(n)<-10.0 .or. eta_n(n)>10.0)) then found_blowup_loc=1 write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep @@ -292,7 +293,8 @@ subroutine check_blowup(istep, mesh) !!PS write(*,*) 'hnode(1, n) = ',hnode(1,n) write(*,*) 'hnode(:, n) = ',hnode(ulevels_nod2D(n):nlevels_nod2D(n),n) write(*,*) - endif + + endif !___________________________________________________________________ ! check surface vertical velocity --> in case of zlevel and zstar From 72c6343bef77ec4981005dd12b6cc8b4ad45f020 Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 22 Mar 2021 22:02:02 +0100 Subject: [PATCH 20/40] clean up ice_thermo_oce.F90 --- src/ice_thermo_oce.F90 | 62 ++++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 29 deletions(-) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index 63e0f9074..a37a30a82 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -9,46 +9,50 @@ subroutine cut_off(mesh) #include "associate_mesh.h" -!_______________________________________________________________________________ -! lower cutoff: a_ice -where(a_ice>1.0_WP) a_ice=1.0_WP + !___________________________________________________________________________ + ! lower cutoff: a_ice + where(a_ice>1.0_WP) a_ice=1.0_WP -! upper cutoff: a_ice -where(a_ice<0.1e-8_WP) - a_ice=0.0_WP + ! upper cutoff: a_ice + where(a_ice<0.1e-8_WP) + a_ice=0.0_WP #if defined (__oifs) - m_ice=0.0_WP - m_snow=0.0_WP - ice_temp=273.15_WP + m_ice=0.0_WP + m_snow=0.0_WP + ice_temp=273.15_WP #endif /* (__oifs) */ -end where + end where -!_______________________________________________________________________________ -! lower cutoff: m_ice -where(m_ice<0.1e-8_WP) - m_ice=0.0_WP + !___________________________________________________________________________ + ! lower cutoff: m_ice + where(m_ice<0.1e-8_WP) + m_ice=0.0_WP #if defined (__oifs) - m_snow=0.0_WP - a_ice=0.0_WP - ice_temp=273.15_WP + m_snow=0.0_WP + a_ice=0.0_WP + ice_temp=273.15_WP #endif /* (__oifs) */ -end where -! upper cutoff: m_ice -where(m_ice>10.0_WP .and. ulevels_nod2d==1) m_ice=10.0_WP - -!_______________________________________________________________________________ -! lower cutoff: m_snow -where(m_snow<0.1e-8_WP) m_snow=0.0_WP -! upper cutoff: m_snow -where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP + end where + + ! upper cutoff SH: m_ice + where(m_ice>5.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)<0.0_WP) m_ice=5.0_WP + ! upper cutoff NH: m_ice + where(m_ice>10.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)>0.0_WP) m_ice=10.0_WP + + !___________________________________________________________________________ + ! lower cutoff: m_snow + where(m_snow<0.1e-8_WP) m_snow=0.0_WP + + ! upper cutoff: m_snow + where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP -!_______________________________________________________________________________ + !___________________________________________________________________________ #if defined (__oifs) -where(ice_temp>273.15_WP) ice_temp=273.15_WP + where(ice_temp>273.15_WP) ice_temp=273.15_WP #endif /* (__oifs) */ #if defined (__oifs) -where(ice_temp < 173.15_WP .and. a_ice >= 0.1e-8_WP) ice_temp=271.35_WP + where(ice_temp < 173.15_WP .and. a_ice >= 0.1e-8_WP) ice_temp=271.35_WP #endif /* (__oifs) */ end subroutine cut_off From 2ea94531e2b65c642cdb9c64b6d0230ae26cff55 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 23 Mar 2021 13:55:54 +0100 Subject: [PATCH 21/40] fix bug regarding cavity rigid lid approximation for zstar --- src/oce_ale.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 1498b2bd8..7fa249ff5 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1555,9 +1555,12 @@ subroutine compute_ssh_rhs_ale(mesh) ! shown in eq (11) rhs of "FESOM2: from finite elements to finte volumes, S. Danilov..." eq. (11) rhs if ( .not. trim(which_ALE)=='linfs') then do n=1,myDim_nod2D - if (ulevels_nod2D(n)>1) cycle nzmin = ulevels_nod2D(n) - ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n)+(1.0_WP-alpha)*ssh_rhs_old(n) + if (ulevels_nod2D(n)>1) then + ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n) + else + ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n)+(1.0_WP-alpha)*ssh_rhs_old(n) + end if end do else do n=1,myDim_nod2D From ced4a3e77b6a7736bf492c557cef114f7148787e Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 6 Apr 2021 12:30:53 +0200 Subject: [PATCH 22/40] improve partitioning routine, so that the cavity geometric constrains fully converge --- src/fvom_init.F90 | 574 +++++++++++++++++++++---------------- src/gen_modules_config.F90 | 1 + 2 files changed, 329 insertions(+), 246 deletions(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 2a2c73972..7069fa29e 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -149,6 +149,11 @@ subroutine read_mesh_ini(mesh) ! =================== ! Surface mesh ! =================== + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: read elem2d.out & nod2d.out '//achar(27)//'[0m' + end if + open (20,file=trim(meshpath)//'nod2d.out', status='old') open (21,file=trim(meshpath)//'elem2d.out', status='old') READ(20,*) mesh%nod2D @@ -217,8 +222,8 @@ subroutine read_mesh_cavity(mesh) !___________________________________________________________________________ if (mype==0) then - write(*,*) '____________________________________________________________' - write(*,*) ' --> read cavity depth' + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: read cavity depth '//achar(27)//'[0m' end if !___________________________________________________________________________ @@ -324,6 +329,10 @@ end subroutine elem_center ! (a) find edges. To make the procedure fast ! one needs neighbourhood arrays ! ==================== +if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute edge connectivity '//achar(27)//'[0m' +end if allocate(ne_num(nod2d), ne_pos(MAX_ADJACENT, nod2D), nn_num(nod2D)) ne_num=0 @@ -632,47 +641,59 @@ END SUBROUTINE find_edges_ini !> 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 -use g_parsup -implicit none -INTEGER :: nodes(3), elems(3), eledges(3) -integer :: elem, elem1, j, n, q, node, enum,count1,count2,exit_flag,i,nz,fileID=111 -real(kind=WP) :: x,dmean -integer :: thers_lev=5 -character*200 :: file_name -type(t_mesh), intent(inout), target :: mesh + use g_config + use mod_mesh + use g_parsup + implicit none + INTEGER :: nodes(3), elems(3), eledges(3) + integer :: elem, elem1, j, n, nneighb,q , node, i, nz + integer :: count_iter, count_neighb_open, exit_flag, fileID=111 + real(kind=WP) :: x,dmean + integer :: max_iter=1000 + character*200 :: file_name + type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: read bottom depth '//achar(27)//'[0m' + end if + + ALLOCATE(mesh%depth(nod2D)) + depth => mesh%depth !required after the allocation, otherwise the pointer remains undefined + file_name=trim(meshpath)//'aux3d.out' + open(fileID, file=file_name) + read(fileID,*) nl ! the number of levels + allocate(mesh%zbar(nl)) ! their standard depths + + zbar => mesh%zbar !required after the allocation, otherwise the pointer remains undefined + read(fileID,*) zbar + if(zbar(2)>0) zbar=-zbar ! zbar is negative + + 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 + DO n=1,nod2D + read(fileID,*) x + if (x>0) x=-x + if (x>zbar(thers_zbar_lev)) x=zbar(thers_zbar_lev) !TODO KK threshholding for depth + depth(n)=x + END DO + close(fileID) + + if(depth(2)>0) depth=-depth ! depth is negative -ALLOCATE(mesh%depth(nod2D)) -depth => mesh%depth !required after the allocation, otherwise the pointer remains undefined -file_name=trim(meshpath)//'aux3d.out' -open(fileID, file=file_name) -read(fileID,*) nl ! the number of levels -allocate(mesh%zbar(nl)) ! their standard depths -zbar => mesh%zbar !required after the allocation, otherwise the pointer remains undefined -read(fileID,*) zbar -if(zbar(2)>0) zbar=-zbar ! zbar is negative -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 -DO n=1,nod2D - read(fileID,*) x - if (x>0) x=-x - if (x>zbar(thers_lev)) x=zbar(thers_lev) !TODO KK threshholding for depth - depth(n)=x -END DO -close(fileID) - -if(depth(2)>0) depth=-depth ! depth is negative - - -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 + !___________________________________________________________________________ + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute elem, vertice bottom depth index '//achar(27)//'[0m' + end if + + 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 @@ -699,25 +720,38 @@ subroutine find_levels(mesh) end if end do if((exit_flag==0).and.(dmean<0)) nlevels(n)=nl - if(dmean>=0) nlevels(n)=thers_lev + if(dmean>=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 + if (mype==0) then + !_______________________________________________________________________ + file_name=trim(meshpath)//'elvls_raw.out' + open(fileID, file=file_name) + do n=1,elem2D + write(fileID,*) nlevels(n) + end do + close(fileID) + endif !___________________________________________________________________________ ! 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=4,nl + do nz=thers_zbar_lev+1,nl exit_flag=0 - count1=0 + count_iter=0 !_______________________________________________________________________ ! iteration loop within each layer - do while((exit_flag==0).and.(count1<1000)) + 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 - q = merge(3,4,elem2D_nodes(1,n) == elem2D_nodes(4,n)) + nneighb = merge(3,4,elem2D_nodes(1,n) == elem2D_nodes(4,n)) ! ! +---isolated bottom cell ! ._______________ | _______________________. @@ -736,16 +770,15 @@ subroutine find_levels(mesh) ! |###|###|###|###|###|###|###|###|###|###|###|###|###|###| ! if (nlevels(n)>=nz) then - count2=0 + count_neighb_open=0 elems=elem_neighbors(1:3,n) - !___________________________________________________________ ! loop over neighbouring triangles - do i=1,q + do i=1,nneighb if (elems(i)>1) then if (nlevels(elems(i))>=nz) then !count neighbours - count2=count2+1 + count_neighb_open=count_neighb_open+1 endif endif enddo @@ -754,15 +787,18 @@ subroutine find_levels(mesh) ! 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 (count2<2) then + 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-11) nlevels(elems(i)) = max(nlevels(elems(i)),nz) + 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 @@ -774,6 +810,7 @@ subroutine find_levels(mesh) end if ! --> if (nlevels(n)>=nz) then end do ! --> do n=1,elem2D end do ! --> do while((exit_flag==0).and.(count1<1000)) + write(*,*) ' -[iter ]->: nlevel',count_iter,'/',max_iter,', nz=',nz end do ! --> do nz=4,nl !___________________________________________________________________________ @@ -827,8 +864,9 @@ subroutine find_levels_cavity(mesh) use g_parsup implicit none integer :: nodes(3), elems(3) - integer :: elem, node, nz, j - integer :: exit_flag, count_iter, count_neighb_open, nneighb, cavity_maxlev + integer :: elem, node, nz, j, idx + integer :: count_neighb_open, nneighb, cavity_maxlev, count_isoelem + integer :: exit_flag1, count_iter, max_iter=1000, exit_flag2, count_iter2, max_iter2=10 real(kind=WP) :: dmean character*200 :: file_name integer, allocatable, dimension(:,:) :: numelemtonode, idxelemtonode @@ -836,13 +874,17 @@ subroutine find_levels_cavity(mesh) type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute elem,vertice cavity depth index '//achar(27)//'[0m' + end if + !___________________________________________________________________________ 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 @@ -859,217 +901,248 @@ subroutine find_levels_cavity(mesh) !_______________________________________________________________________ ! vertical elem level index of cavity-ocean boundary - exit_flag=0 ulevels(elem) = 1 + if (dmean<0.0_WP) ulevels(elem) = 2 + do nz=1,nlevels(elem)-1 - !!PS if(Z(nz) I need 4 valid full depth, 3 valid mid-depth levels ulevels(elem)=nz ! to compute shechpetkin PGF exit end if end do - if ((exit_flag==0).and.(dmean<0)) ulevels(elem)=nlevels(elem) cavity_maxlev = max(cavity_maxlev,ulevels(elem)) end do + !___________________________________________________________________________ + ! write out cavity mesh files for vertice and elemental position of + ! vertical cavity-ocean boundary before the iterative geometric adaption to + ! eliminate isolated cells + if (mype==0) then + ! write out elemental cavity-ocean boundary level + file_name=trim(meshpath)//'cavity_elvls_raw.out' + open(20, file=file_name) + do elem=1,elem2D + write(20,*) ulevels(elem) + enddo + close(20) + endif !___________________________________________________________________________ ! Eliminate cells that have two cavity boundary faces --> should not be ! possible in FESOM2.0 ! loop over all cavity levels allocate(elemreducelvl(elem2d)) - do nz=1,cavity_maxlev - exit_flag=0 - count_iter=0 + allocate(numelemtonode(nl,nod2D),idxelemtonode(nl,nod2D)) + + !___________________________________________________________________________ + ! outer iteration loop + count_iter2 = 0 + exit_flag2 = 0 + 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)) - ! - ! .___________________________.~~~~~~~~~~~~~~~~~~~~~~~~~~ - ! |###|###|###|###|###|###|###| - ! |# CAVITY |###| . |###|###| OCEAN - ! |###|###|###| /|\|###| - ! |###|###| | - ! |###| +-- Not good can lead to isolated cells - ! - if (nz >= ulevels(elem)) then - count_neighb_open=0 - elems=elem_neighbors(1:3,elem) - - !___________________________________________________________ - ! loop over neighbouring triangles - do j = 1, nneighb - if (elems(j)>0) then ! if its a valid boundary triangle, 0=missing value - ! check for isolated cell - if (ulevels(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 - -!!PS if (elem==133438) then -!!PS write(*,*) -!!PS write(*,*) 'nz =', nz -!!PS write(*,*) 'elem =', elem -!!PS write(*,*) 'ulevels(elem) =', ulevels(elem) -!!PS write(*,*) 'nlevels(elem) =', nlevels(elem) -!!PS write(*,*) 'elemreducelvl(elem)=',elemreducelvl(elem) -!!PS write(*,*) 'elems =', elems -!!PS write(*,*) 'ulevels(elems) =', ulevels(elems) -!!PS write(*,*) 'nlevels(elems) =', nlevels(elems) -!!PS write(*,*) 'nlvl-ulvl =', nlevels(elems)-ulevels(elems) -!!PS write(*,*) 'elemreducelvl(elems)=',elemreducelvl(elems) -!!PS write(*,*) 'count_neighb_open =',count_neighb_open -!!PS write(*,*) -!!PS end if - - !___________________________________________________________ - ! 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<2) then - ! if cell is isolated convert it to a deeper ocean levels - ! except when this levels would remain less than 3 valid - ! bottom levels --> in case make the levels of all sorounding - ! one level shallower - if (nlevels(elem)-(nz+1)<3) then - do j = 1, nneighb - if (elems(j)>0) then - if (ulevels(elems(j))>1 .and. ulevels(elems(j))>ulevels(elem) ) then - ulevels(elems(j)) = min(ulevels(elems(j)),nz) - elemreducelvl(elems(j))=1 - end if - end if - end do + ! iteration loop within each layer + do while((exit_flag1==0).and.(count_iter 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 + ! + if ( nz >= ulevels(elem) .and. nz0) then - if (ulevels(elems(j))>1 .and. ulevels(elems(j))>ulevels(elem) ) then - ulevels(elems(j)) = min(ulevels(elems(j)),nz) - elemreducelvl(elems(j))=1 - end if - end if - end do - end if + !_______________________________________________________ + ! loop over neighbouring triangles + do j = 1, nneighb + if (elems(j)>0) 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 - !force recheck for all current ocean cells - exit_flag=0 - endif ! --> if (count_neighb_open<2) then - - end if ! --> if (nz >= ulevels(elem)) then - end do ! --> do elem=1,elem2D - end do ! --> do while((exit_flag==0).and.(count_iter<1000)) - end do ! --> do nz=1,cavity_maxlev - deallocate(elemreducelvl) - - !___________________________________________________________________________ - ! vertical vertice level index of cavity_ocean boundary - ulevels_nod2D = nl - do elem=1,elem2D - nneighb = merge(3,4,elem2D_nodes(1,elem) == elem2D_nodes(4,elem)) + !_______________________________________________________ + ! 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<2) then + count_isoelem = count_isoelem+1 + ! if cell is isolated convert it to a deeper ocean levels + ! 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)==0 ) then + ulevels(elem)=nz+1 + else + ! --> can not increase depth anymore to eleminate isolated + ! cell, otherwise lessthan 3 valid layers + ! --> therefor reduce depth of ONE!!! of the neighbouring + ! triangles. Choose trinagle whos depth is already closest + ! to nz + idx = minloc(ulevels(elems)-nz, 1, MASK=( (elems>0) .and. ((ulevels(elems)-nz)>0) ) ) + ulevels(elems(idx)) = nz-1 + elemreducelvl(elems(idx)) = elemreducelvl(elems(idx))+1 + 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)) then + end do ! --> do elem=1,elem2D + end do ! --> do while((exit_flag==0).and.(count_iter<1000)) + write(*,*) ' -[iter ]->: ulevel',count_iter,'/',max_iter,', nz=',nz + end do ! --> do nz=1,cavity_maxlev + !_______________________________________________________________________ - ! 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 - - !___________________________________________________________________________ - ! check ulevels if ulevels=nlevels(elem)) then - if (mype==0) write(*,*) ' ERROR: found element cavity depth deeper or equal bottom depth' - call par_ex(0) - end if - if (nlevels(elem)-ulevels(elem)<3) then - write(*,*) ' ERROR: found less than three valid element ocean layers' - write(*,*) ' ulevels,nlevels = ',ulevels(elem), nlevels(elem) - write(*,*) ' ulevels(neighb) = ',ulevels(elem_neighbors(1:3,elem)) - write(*,*) ' nlevels(neighb) = ',nlevels(elem_neighbors(1:3,elem)) - call par_ex(0) - end if - end do - - !___________________________________________________________________________ - ! check ulevels_nod2d if ulevels_nod2d=nlevels_nod2D(elem)) then - if (mype==0) write(*,*) ' ERROR: found vertice cavity depth deeper or equal bottom depth' - call par_ex(0) - end if - if (nlevels_nod2D(elem)-ulevels_nod2D(elem)<3) then - if (mype==0) write(*,*) ' ERROR: found less than three valid vertice ocean layers' - end if - end do - - do elem=1,elem2D - if (ulevels(elem)< maxval(ulevels_nod2D(elem2D_nodes(:,elem))) ) then - if (mype==0) then - write(*,*) ' ERROR: found element cavity depth that is shallower than its valid maximum cavity vertice depths' - write(*,*) ' ule | uln = ',ulevels(elem),' | ',ulevels_nod2D(elem2D_nodes(:,elem)) - end if - call par_ex(0) - end if - end do - - !___________________________________________________________________________ - ! check how many triangle elements contribute to every vertice in every layer - ! every vertice in every layer should be connected to at least two triangle - ! elements ! - allocate(numelemtonode(nl,nod2D),idxelemtonode(nl,nod2D)) - numelemtonode=0 - idxelemtonode=0 - do node=1, nod2D - do j=1,nod_in_elem2D_num(node) - elem=nod_in_elem2D(j,node) - do nz=ulevels(elem),nlevels(elem)-1 - numelemtonode(nz,node) = numelemtonode(nz,node) + 1 - idxelemtonode(nz,node) = elem + ! vertical vertice level index of cavity_ocean boundary + 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 - end do - - exit_flag = 0 - do node=1, nod2D - do nz=1,nl - if (numelemtonode(nz,node)== 1) then - write(*,*) 'ERROR: found vertice with just one triangle:', mype, nz, 'node=',node, ulevels_nod2d(node), nlevels_nod2D(node), & - 'elem=', idxelemtonode(nz,node), ulevels(idxelemtonode(nz,node)), nlevels(idxelemtonode(nz,node)) - exit_flag = 1 + + !_______________________________________________________________________ + ! check ulevels if ulevels=nlevels(elem)) then + if (mype==0) 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 + !!PS write(*,*) ' elem = ',elem + !!PS write(*,*) ' elem_neighbors = ',elem_neighbors(1:3,elem) + !!PS write(*,*) ' ulevels,nlevels = ',ulevels(elem), nlevels(elem) + !!PS write(*,*) ' ulevels(neighb) = ',ulevels(elem_neighbors(1:3,elem)) + !!PS write(*,*) ' nlevels(neighb) = ',nlevels(elem_neighbors(1:3,elem)) + !!PS write(*,*) ' lon_node = ',mesh%coord_nod2D(1,elem2D_nodes(1:3,elem))/rad + !!PS write(*,*) ' lat_node = ',mesh%coord_nod2D(2,elem2D_nodes(1:3,elem))/rad + !!PS write(*,*) ' lon_elem = ',sum(mesh%coord_nod2D(1,elem2D_nodes(1:3,elem))/rad)/3.0 + !!PS write(*,*) ' lat_elem = ',sum(mesh%coord_nod2D(2,elem2D_nodes(1:3,elem))/rad)/3.0 + exit_flag2 = 0 + + end if + end do + + !_______________________________________________________________________ + ! check ulevels_nod2d if ulevels_nod2d=nlevels_nod2D(elem)) then + if (mype==0) write(*,*) ' -[check]->: vertice cavity depth deeper or equal bottom depth, node=', elem2d + exit_flag2 = 0 + + end if + if (nlevels_nod2D(elem)-ulevels_nod2D(elem)<3) then + if (mype==0) write(*,*) ' -[check]->: less than three valid vertice ocean layers, node=', elem2d + exit_flag2 = 0 + + end if + end do + + do elem=1,elem2D + if (ulevels(elem)< maxval(ulevels_nod2D(elem2D_nodes(:,elem))) ) then + if (mype==0) then + write(*,*) ' -[check]->: found elem cavity shallower than its valid maximum cavity vertice depths, elem=', elem2d + !!PS write(*,*) ' ule | uln = ',ulevels(elem),' | ',ulevels_nod2D(elem2D_nodes(:,elem)) + end if + exit_flag2 = 0 + + end if + end do + + !_______________________________________________________________________ + ! check how many triangle elements contribute to every vertice in every layer + ! every vertice in every layer should be connected to at least two triangle + ! elements ! + numelemtonode=0 + idxelemtonode=0 + do node=1, nod2D + do j=1,nod_in_elem2D_num(node) + elem=nod_in_elem2D(j,node) + do nz=ulevels(elem),nlevels(elem)-1 + numelemtonode(nz,node) = numelemtonode(nz,node) + 1 + idxelemtonode(nz,node) = elem + end do + end do + end do + + count_iter=0 + do node=1, nod2D + do nz=1,nl + if (numelemtonode(nz,node)== 1) then + exit_flag2 = 0 + count_iter = count_iter+1 + write(*,*) ' -[check]->: node has only 1 triangle: n=', node, ',nz=',nz +!!PS write(*,*) ' node=', node, ulevels_nod2d(node), nlevels_nod2D(node) +!!PS write(*,*) ' elem=', idxelemtonode(nz,node), ulevels(idxelemtonode(nz,node)), nlevels(idxelemtonode(nz,node)) +!!PS write(*,*) ' lon/lat_n = ',mesh%coord_nod2D(1,node)/rad, mesh%coord_nod2D(2,node)/rad +!!PS write(*,*) ' n_eneighb =', nod_in_elem2D(1:nod_in_elem2D_num(node),node) +!!PS write(*,*) ' ulvls(n_eneighb)=', ulevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) +!!PS write(*,*) ' nlvls(n_eneighb)=', nlevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) +!!PS write(*,*) + exit + end if + end do end do - end do + + !_______________________________________________________________________ + 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(*,*) ' iter step ', count_iter2,' out of ', max_iter2 + write(*,*) + end if + + !_______________________________________________________________________ + end do + deallocate(elemreducelvl) deallocate(numelemtonode,idxelemtonode) - if (exit_flag == 1) call par_ex(0) !___________________________________________________________________________ -!!PS ! compute nodal cavity flag: 1 yes cavity/ 0 no cavity -!!PS cavity_flag = 0 -!!PS do node=1,nod2D -!!PS if (ulevels_nod2D(node)>1) cavity_flag(node)=1 -!!PS end do + if (exit_flag2 == 0) then + write(*,*) + print *, achar(27)//'[31m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;31m'//' -[ERROR]->: Cavity geometry constrains did not converge !!!'//achar(27)//'[0m' + write(*,*) + call par_ex(0) + else + write(*,*) + print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' + print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge!!!' + write(*,*) + end if + !___________________________________________________________________________ ! write out cavity mesh files for vertice and elemental position of @@ -1086,14 +1159,10 @@ subroutine find_levels_cavity(mesh) ! write out vertice cavity-ocean boundary level + yes/no cavity flag file_name=trim(meshpath)//'cavity_nlvls.out' open(20, file=file_name) -!!PS file_name=trim(meshpath)//'cavity_flag.out' -!!PS open(21, file=file_name) do node=1,nod2D write(20,*) ulevels_nod2D(node) -!!PS write(21,*) cavity_flag(node) enddo close(20) -!!PS close(21) endif end subroutine find_levels_cavity @@ -1322,6 +1391,10 @@ subroutine communication_ini(mesh) type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" + if (mype==0) then + print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' + print *, achar(27)//'[7;1m' //' -->: compute communication arrays '//achar(27)//'[0m' + end if ! Create the distributed mesh subdirectory write(npes_string,"(I10)") npes dist_mesh_dir=trim(meshpath)//'dist_'//trim(ADJUSTL(npes_string))//'/' @@ -1381,6 +1454,11 @@ end subroutine partit end interface #include "associate_mesh_ini.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...' @@ -1446,7 +1524,11 @@ subroutine check_partitioning(mesh) integer, allocatable :: ne_part(:), ne_part_num(:), ne_part_load(:,:) type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" - + + 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) diff --git a/src/gen_modules_config.F90 b/src/gen_modules_config.F90 index 25ba9be26..4c2df68b8 100755 --- a/src/gen_modules_config.F90 +++ b/src/gen_modules_config.F90 @@ -80,6 +80,7 @@ module g_config real(kind=WP) :: gammaEuler=-90. ! then around new z. ! Set to zeros to work with ! geographical coordinates + integer :: thers_zbar_lev=5 ! minimum number of levels to be character(len=5) :: which_depth_n2e='mean' namelist /geometry/ cartesian, fplane, & cyclic_length, rotated_grid, alphaEuler, betaEuler, gammaEuler, force_rotation, which_depth_n2e From 4b9b0d13ebaec087aa24abe773c4d5e133a234e8 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 6 Apr 2021 12:31:51 +0200 Subject: [PATCH 23/40] fix bug in zlevel found by thomas rackow --- src/oce_ale.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 7fa249ff5..52fcd0a51 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1870,7 +1870,7 @@ subroutine vert_vel_ale(mesh) ! layer reached already minimum layerthickness) max_dhbar2distr = 0.0_WP !max_dhbar2distr = (zbar(1:lzstar_lev)-zbar(2:lzstar_lev+1))*min_hnode - hnode(1:lzstar_lev,n); - max_dhbar2distr = (zbar(nzmin:nzmin+lzstar_lev-1)-zbar(nzmin:nzmin+lzstar_lev-1+1))*min_hnode - hnode(nzmin:nzmin+lzstar_lev-1,n); + max_dhbar2distr = (zbar(nzmin:nzmin+lzstar_lev-1)-zbar(nzmin+1:nzmin+lzstar_lev-1+1))*min_hnode - hnode(nzmin:nzmin+lzstar_lev-1,n); where (max_dhbar2distr>=0.0_WP) max_dhbar2distr=0.0_WP !_______________________________________________________________ From b5aea4cfc68fcfcb6a56693dc99f3959f08010a6 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 8 Apr 2021 11:07:00 +0200 Subject: [PATCH 24/40] fix bug in cavity partitioning --- src/fvom_init.F90 | 135 ++++++++++++++++++++++++++-------------------- 1 file changed, 78 insertions(+), 57 deletions(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 7069fa29e..cf690ae03 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -810,7 +810,7 @@ subroutine find_levels(mesh) end if ! --> if (nlevels(n)>=nz) then end do ! --> do n=1,elem2D end do ! --> do while((exit_flag==0).and.(count1<1000)) - write(*,*) ' -[iter ]->: nlevel',count_iter,'/',max_iter,', nz=',nz + write(*,"(A, I5, A, i5, A, I3)") ' -[iter ]->: nlevel, iter/maxiter=',count_iter,'/',max_iter,', nz=',nz end do ! --> do nz=4,nl !___________________________________________________________________________ @@ -854,7 +854,8 @@ subroutine find_levels(mesh) endif end subroutine find_levels - +! +! !_______________________________________________________________________________ ! finds elemental and nodal levels of cavity-ocean boundary. ! Creates 2 files: cavity_elvls.out, cavity_nlvls.out @@ -870,10 +871,10 @@ subroutine find_levels_cavity(mesh) real(kind=WP) :: dmean character*200 :: file_name integer, allocatable, dimension(:,:) :: numelemtonode, idxelemtonode - integer, allocatable, dimension(:) :: elemreducelvl + integer, allocatable, dimension(:) :: elemreducelvl, elemfixlvl type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" - + !___________________________________________________________________________ if (mype==0) then print *, achar(27)//'[1m' //'____________________________________________________________'//achar(27)//'[0m' print *, achar(27)//'[7;1m' //' -->: compute elem,vertice cavity depth index '//achar(27)//'[0m' @@ -931,21 +932,22 @@ subroutine find_levels_cavity(mesh) ! Eliminate cells that have two cavity boundary faces --> should not be ! possible in FESOM2.0 ! loop over all cavity levels - allocate(elemreducelvl(elem2d)) + allocate(elemreducelvl(elem2d),elemfixlvl(elem2d)) allocate(numelemtonode(nl,nod2D),idxelemtonode(nl,nod2D)) !___________________________________________________________________________ ! outer iteration loop count_iter2 = 0 exit_flag2 = 0 + elemfixlvl = 0 do while((exit_flag2==0) .and. (count_iter20) then ! if its a valid boundary triangle, 0=missing value ! check for isolated cell - if ( ulevels(elems(j))<= nz .and. & + 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 @@ -993,7 +995,7 @@ subroutine find_levels_cavity(mesh) ! 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)==0 ) then + if ( (nlevels(elem)-(nz+1))>=3 .and. elemreducelvl(elem)==0 .and. elemfixlvl(elem)==0) then ulevels(elem)=nz+1 else ! --> can not increase depth anymore to eleminate isolated @@ -1009,10 +1011,13 @@ subroutine find_levels_cavity(mesh) !force recheck for all current ocean cells exit_flag1=0 end if ! --> if (count_neighb_open<2) then - end if ! --> if (nz >= ulevels(elem)) 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(*,*) ' -[iter ]->: ulevel',count_iter,'/',max_iter,', nz=',nz + write(*,"(A, I5, A, i5, A, I3)") ' -[iter ]->: ulevel, iter/maxiter=',count_iter,'/',max_iter,', nz=',nz end do ! --> do nz=1,cavity_maxlev !_______________________________________________________________________ @@ -1026,7 +1031,7 @@ subroutine find_levels_cavity(mesh) node=elem2D_nodes(j,elem) ulevels_nod2D(node)=min(ulevels_nod2D(node),ulevels(elem)) end do - end do + end do ! --> do elem=1,elem2D !_______________________________________________________________________ ! check ulevels if ulevels: less than three valid elem ocean layers, elem=',elem - !!PS write(*,*) ' elem = ',elem - !!PS write(*,*) ' elem_neighbors = ',elem_neighbors(1:3,elem) - !!PS write(*,*) ' ulevels,nlevels = ',ulevels(elem), nlevels(elem) - !!PS write(*,*) ' ulevels(neighb) = ',ulevels(elem_neighbors(1:3,elem)) - !!PS write(*,*) ' nlevels(neighb) = ',nlevels(elem_neighbors(1:3,elem)) - !!PS write(*,*) ' lon_node = ',mesh%coord_nod2D(1,elem2D_nodes(1:3,elem))/rad - !!PS write(*,*) ' lat_node = ',mesh%coord_nod2D(2,elem2D_nodes(1:3,elem))/rad - !!PS write(*,*) ' lon_elem = ',sum(mesh%coord_nod2D(1,elem2D_nodes(1:3,elem))/rad)/3.0 - !!PS write(*,*) ' lat_elem = ',sum(mesh%coord_nod2D(2,elem2D_nodes(1:3,elem))/rad)/3.0 exit_flag2 = 0 end if - end do + end do ! --> do elem=1,elem2D !_______________________________________________________________________ ! check ulevels_nod2d if ulevels_nod2d=nlevels_nod2D(elem)) then - if (mype==0) write(*,*) ' -[check]->: vertice cavity depth deeper or equal bottom depth, node=', elem2d + do node=1,nod2D + !___________________________________________________________________ + if (ulevels_nod2D(node)>=nlevels_nod2D(node)) then + if (mype==0) write(*,*) ' -[check]->: vertice cavity depth deeper or equal bottom depth, node=', node exit_flag2 = 0 - end if - if (nlevels_nod2D(elem)-ulevels_nod2D(elem)<3) then - if (mype==0) write(*,*) ' -[check]->: less than three valid vertice ocean layers, node=', elem2d + + !___________________________________________________________________ + if (nlevels_nod2D(node)-ulevels_nod2D(node)<3) then + if (mype==0) write(*,*) ' -[check]->: less than three valid vertice ocean layers, node=', node exit_flag2 = 0 - end if - end do + end do ! --> do node=1,nod2D do elem=1,elem2D if (ulevels(elem)< maxval(ulevels_nod2D(elem2D_nodes(:,elem))) ) then - if (mype==0) then - write(*,*) ' -[check]->: found elem cavity shallower than its valid maximum cavity vertice depths, elem=', elem2d - !!PS write(*,*) ' ule | uln = ',ulevels(elem),' | ',ulevels_nod2D(elem2D_nodes(:,elem)) - end if + if (mype==0) write(*,*) ' -[check]->: found elem cavity shallower than its valid maximum cavity vertice depths, elem=', elem2d exit_flag2 = 0 - end if - end do + end do ! --> do elem=1,elem2D !_______________________________________________________________________ - ! check how many triangle elements contribute to every vertice in every layer - ! every vertice in every layer should be connected to at least two triangle - ! elements ! + ! compute how many triangle elements contribute to every vertice in every layer numelemtonode=0 idxelemtonode=0 do node=1, nod2D @@ -1095,32 +1086,59 @@ subroutine find_levels_cavity(mesh) idxelemtonode(nz,node) = elem end do end do - end do + end do ! --> do node=1, nod2D + !_______________________________________________________________________ + ! check if every vertice in every layer should be connected to at least + ! two triangle elements ! count_iter=0 do node=1, nod2D - do nz=1,nl - if (numelemtonode(nz,node)== 1) then + + !___________________________________________________________________ + do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 + + !_______________________________________________________________ + ! nodes has zero neighbouring triangles and is completely isolated + ! need to adapt ulevels by hand --> inflicts another outher + ! iteration loop (exit_flag2=0) + if (numelemtonode(nz,node)==0) then exit_flag2 = 0 count_iter = count_iter+1 - write(*,*) ' -[check]->: node has only 1 triangle: n=', node, ',nz=',nz -!!PS write(*,*) ' node=', node, ulevels_nod2d(node), nlevels_nod2D(node) -!!PS write(*,*) ' elem=', idxelemtonode(nz,node), ulevels(idxelemtonode(nz,node)), nlevels(idxelemtonode(nz,node)) -!!PS write(*,*) ' lon/lat_n = ',mesh%coord_nod2D(1,node)/rad, mesh%coord_nod2D(2,node)/rad -!!PS write(*,*) ' n_eneighb =', nod_in_elem2D(1:nod_in_elem2D_num(node),node) -!!PS write(*,*) ' ulvls(n_eneighb)=', ulevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) -!!PS write(*,*) ' nlvls(n_eneighb)=', nlevels(nod_in_elem2D(1:nod_in_elem2D_num(node),node)) -!!PS write(*,*) - exit + write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz,node) ,' triangle: n=', node, ', nz=',nz + !___________________________________________________________ + ! if node has no neighboring triangle somewhere in the middle + ! of the water column at nz (can happen but seldom) than set + ! all ulevels(elem) of sorounding trinagles whos ulevel is + ! depper than nz, equal to nz and fix that value to forbit it + ! to be changed (elemfixlvl > 0) + do j=1,nod_in_elem2D_num(node) + elem=nod_in_elem2D(j,node) + if (ulevels(elem)>nz) then + ulevels(elem) = nz + elemfixlvl(elem) = elemfixlvl(elem)+1 + end if + end do + end if + + !_______________________________________________________________ + ! nodes has just one neighbouring triangle --> but needs two --> + ! inflicts another outher iteration loop (exit_flag2=0) + if (numelemtonode(nz,node)==1) then + exit_flag2 = 0 + count_iter = count_iter+1 + write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz,node) ,' triangle: n=', node, ', nz=',nz end if - end do - end do + + end do ! --> do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 + + 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(*,*) ' iter step ', count_iter2,' out of ', max_iter2 + write(*,"(A, I3, A, I3)") ' iter step ', count_iter2,'/', max_iter2 write(*,*) end if @@ -1130,6 +1148,8 @@ subroutine find_levels_cavity(mesh) deallocate(numelemtonode,idxelemtonode) !___________________________________________________________________________ + ! check if cavity geometry totaly converged or failed to converge in the later + ! case will break up model if (exit_flag2 == 0) then write(*,*) print *, achar(27)//'[31m' //'____________________________________________________________'//achar(27)//'[0m' @@ -1139,7 +1159,8 @@ subroutine find_levels_cavity(mesh) else write(*,*) print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' - print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge!!!' + print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge, Yippee-Ki-Yay, Beep!!!' + write(*,*) end if From b489384b9d255b23d561e4cebd468b3bb0f5017f Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 8 Apr 2021 12:34:47 +0200 Subject: [PATCH 25/40] delete some check marks --- src/oce_ale.F90 | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 52fcd0a51..9d2fbb2a4 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1811,6 +1811,7 @@ subroutine vert_vel_ale(mesh) if (Fer_GM) then fer_Wvel(nz,n)=fer_Wvel(nz,n)/area(nz,n) end if + end do end do ! | @@ -2400,16 +2401,7 @@ subroutine impl_vert_visc_ale(mesh) b(nz)=b(nz)-min(0._WP, wd)*zinv c(nz)=c(nz)-max(0._WP, wd)*zinv - if (a(nz)/=a(nz) .or. b(nz)/=b(nz) .or. c(nz)/=c(nz)) then - write(*,*) ' --> found a,b,c is NaN' - write(*,*) 'mype=',mype - write(*,*) 'nz=',nz - write(*,*) 'a(nz), b(nz), c(nz)=',a(nz), b(nz), c(nz) - write(*,*) 'Av(nz,elem)=',Av(nz,elem) - write(*,*) 'Av(nz+1,elem)=',Av(nz+1,elem) - write(*,*) 'Z_n(nz-1:nz+1)=',Z_n(nz-1:nz+1) - write(*,*) 'zbar_n(nz:nz+1)=',zbar_n(nz:nz+1) - endif + end do ! The last row zinv=1.0_WP*dt/(zbar_n(nzmax-1)-zbar_n(nzmax)) @@ -2666,10 +2658,12 @@ subroutine oce_timestep_ale(n, mesh) !___________________________________________________________________________ if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_vel_rhs'//achar(27)//'[0m' -!!PS if (any(UV_rhs/=UV_rhs)) write(*,*) mype,' --> found NaN UV_rhs before compute_vel_rhs' -!!PS if (any(UV/=UV)) write(*,*) mype,' --> found NaN UV before compute_vel_rhs' -!!PS if (any(ssh_rhs/=ssh_rhs))write(*,*) mype,' --> found NaN ssh_rhs before compute_ssh_rhs_ale' -!!PS if (any(ssh_rhs_old/=ssh_rhs_old))write(*,*) mype,' --> found NaN ssh_rhs_old before compute_ssh_rhs_ale' + +!!PS if (any(UV_rhs/=UV_rhs)) write(*,*) n, mype,' --> found NaN UV_rhs before compute_vel_rhs' +!!PS if (any(UV/=UV)) write(*,*) n, mype,' --> found NaN UV before compute_vel_rhs' +!!PS if (any(ssh_rhs/=ssh_rhs)) write(*,*) n, mype,' --> found NaN ssh_rhs before compute_vel_rhs' +!!PS if (any(ssh_rhs_old/=ssh_rhs_old)) write(*,*) n, mype,' --> found NaN ssh_rhs_old before compute_vel_rhs' +!!PS if (any(abs(Wvel_e)>1.0e20)) write(*,*) n, mype,' --> found Inf Wvel_e before compute_vel_rhs' if(mom_adv/=3) then call compute_vel_rhs(mesh) @@ -2684,7 +2678,7 @@ subroutine oce_timestep_ale(n, mesh) if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call impl_vert_visc_ale'//achar(27)//'[0m' if(i_vert_visc) call impl_vert_visc_ale(mesh) t2=MPI_Wtime() - + !___________________________________________________________________________ ! >->->->->->->->->->->->-> ALE-part starts <-<-<-<-<-<-<-<-<-<-<-<- !___________________________________________________________________________ @@ -2699,6 +2693,7 @@ subroutine oce_timestep_ale(n, mesh) ! Take updated ssh matrix and solve --> new ssh! t30=MPI_Wtime() call solve_ssh_ale(mesh) + if ((toy_ocean) .AND. (TRIM(which_toy)=="soufflet")) call relax_zonal_vel(mesh) t3=MPI_Wtime() From 4568b03e0b8f87105d3daeb4806ae7737e419682 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 14 Apr 2021 17:45:19 +0200 Subject: [PATCH 26/40] change some screen output when partition --- src/fvom_init.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index cf690ae03..63ad3baf8 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -815,6 +815,7 @@ subroutine find_levels(mesh) !___________________________________________________________________________ ! 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)) @@ -1022,6 +1023,7 @@ subroutine find_levels_cavity(mesh) !_______________________________________________________________________ ! 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)) @@ -1153,13 +1155,13 @@ subroutine find_levels_cavity(mesh) if (exit_flag2 == 0) then write(*,*) print *, achar(27)//'[31m' //'____________________________________________________________'//achar(27)//'[0m' - print *, achar(27)//'[7;31m'//' -[ERROR]->: Cavity geometry constrains did not converge !!!'//achar(27)//'[0m' + print *, achar(27)//'[7;31m'//' -[ERROR]->: Cavity geometry constrains did not converge !!! *\(>︿<)/*'//achar(27)//'[0m' write(*,*) call par_ex(0) else write(*,*) print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' - print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge, Yippee-Ki-Yay, Beep!!!' + print *, ' -['//achar(27)//'[7;32m'//' OK '//achar(27)//'[0m'//']->: Cavity geometry constrains did converge !!! *\(^o^)/*' write(*,*) end if From d7201e8180fa2d81a065da6360d595b344cbcfe3 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 14 Apr 2021 17:46:05 +0200 Subject: [PATCH 27/40] small change in ../src/write_step_info.F90 --- src/write_step_info.F90 | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/write_step_info.F90 b/src/write_step_info.F90 index 35a5ba69d..01de0d6d3 100644 --- a/src/write_step_info.F90 +++ b/src/write_step_info.F90 @@ -56,11 +56,12 @@ subroutine write_step_info(istep,outfreq, mesh) loc =0. !_______________________________________________________________________ do n=1, myDim_nod2D - loc_eta = loc_eta + area(1, n)*eta_n(n) - loc_hbar = loc_hbar + area(1, n)*hbar(n) - loc_deta = loc_deta + area(1, n)*d_eta(n) - loc_dhbar = loc_dhbar + area(1, n)*(hbar(n)-hbar_old(n)) - loc_wflux = loc_wflux + area(1, n)*water_flux(n) + if (ulevels_nod2D(n)>1) cycle + loc_eta = loc_eta + area(ulevels_nod2D(n), n)*eta_n(n) + loc_hbar = loc_hbar + area(ulevels_nod2D(n), n)*hbar(n) + loc_deta = loc_deta + area(ulevels_nod2D(n), n)*d_eta(n) + loc_dhbar = loc_dhbar + area(ulevels_nod2D(n), n)*(hbar(n)-hbar_old(n)) + loc_wflux = loc_wflux + area(ulevels_nod2D(n), n)*water_flux(n) !!PS loc_hflux = loc_hflux + area(1, n)*heat_flux(n) !!PS loc_temp = loc_temp + area(1, n)*sum(tr_arr(:,n,1))/(nlevels_nod2D(n)-1) !!PS loc_salt = loc_salt + area(1, n)*sum(tr_arr(:,n,2))/(nlevels_nod2D(n)-1) @@ -75,11 +76,19 @@ subroutine write_step_info(istep,outfreq, mesh) !!PS call MPI_AllREDUCE(loc_hflux, int_hflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) !!PS call MPI_AllREDUCE(loc_temp , int_temp , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) !!PS call MPI_AllREDUCE(loc_salt , int_salt , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + int_eta = int_eta /ocean_area int_hbar = int_hbar /ocean_area int_deta = int_deta /ocean_area int_dhbar= int_dhbar/ocean_area int_wflux= int_wflux/ocean_area + +!!PS int_eta = int_eta /ocean_areawithcav +!!PS int_hbar = int_hbar /ocean_areawithcav +!!PS int_deta = int_deta /ocean_areawithcav +!!PS int_dhbar= int_dhbar/ocean_areawithcav +!!PS int_wflux= int_wflux/ocean_areawithcav + !!PS int_hflux= int_hflux/ocean_area !!PS int_temp = int_temp /ocean_area !!PS int_salt = int_salt /ocean_area @@ -300,19 +309,20 @@ subroutine check_blowup(istep, mesh) ! check surface vertical velocity --> in case of zlevel and zstar ! vertical coordinate its indicator if Volume is conserved for ! Wvel(1,n)~maschine preccision - if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) .or. abs(Wvel(1,n))>1e-12 )) then +!!PS if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) .or. abs(Wvel(1,n))>1e-12 )) then + if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) )) then found_blowup_loc=1 write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep write(*,*) ' --STOP--> found surface layer vertical velocity becomes NaN or >1e-12' write(*,*) 'mype = ',mype write(*,*) 'mstep = ',istep write(*,*) 'node = ',n + write(*,*) 'uln, nln = ',ulevels_nod2D(n), nlevels_nod2D(n) + write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad write(*,*) write(*,*) 'Wvel(1, n) = ',Wvel(1,n) write(*,*) 'Wvel(:, n) = ',Wvel(:,n) write(*,*) - write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad - write(*,*) write(*,*) 'hnode(1, n) = ',hnode(1,n) write(*,*) 'hnode(:, n) = ',hnode(:,n) write(*,*) 'hflux = ',heat_flux(n) From b07e7f9bd54b21dd631ba3fdc1a58bf5aaff4ad4 Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 19 Apr 2021 12:34:08 +0200 Subject: [PATCH 28/40] switch back compuation of total surface windstress on elements --- src/ice_oce_coupling.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 15db83c06..e85fec10f 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -52,12 +52,12 @@ subroutine oce_fluxes_mom(mesh) !_______________________________________________________________________ elnodes=elem2D_nodes(:,elem) - !!PS stress_surf(1,elem)=sum(stress_iceoce_x(elnodes)*a_ice(elnodes) + & - !!PS stress_atmoce_x(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP - !!PS stress_surf(2,elem)=sum(stress_iceoce_y(elnodes)*a_ice(elnodes) + & - !!PS stress_atmoce_y(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP - stress_surf(1,elem)=sum(stress_node_surf(1,elnodes))/3.0_WP - stress_surf(2,elem)=sum(stress_node_surf(2,elnodes))/3.0_WP + stress_surf(1,elem)=sum(stress_iceoce_x(elnodes)*a_ice(elnodes) + & + stress_atmoce_x(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP + stress_surf(2,elem)=sum(stress_iceoce_y(elnodes)*a_ice(elnodes) + & + stress_atmoce_y(elnodes)*(1.0_WP-a_ice(elnodes)))/3.0_WP + !!PS stress_surf(1,elem)=sum(stress_node_surf(1,elnodes))/3.0_WP + !!PS stress_surf(2,elem)=sum(stress_node_surf(2,elnodes))/3.0_WP END DO !___________________________________________________________________________ From 7502300e24161b3426c32e996ccb2821a1e80c88 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 20 Apr 2021 14:28:17 +0200 Subject: [PATCH 29/40] comment lower cutoff for m_snow --- src/ice_thermo_oce.F90 | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index a37a30a82..ab1e609ac 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -3,6 +3,7 @@ subroutine cut_off(mesh) use o_param use i_arrays use MOD_MESH + use g_config, only: use_cavity use g_parsup implicit none type(t_mesh), intent(in) , target :: mesh @@ -34,18 +35,22 @@ subroutine cut_off(mesh) #endif /* (__oifs) */ end where - ! upper cutoff SH: m_ice - where(m_ice>5.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)<0.0_WP) m_ice=5.0_WP - ! upper cutoff NH: m_ice - where(m_ice>10.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)>0.0_WP) m_ice=10.0_WP - !___________________________________________________________________________ - ! lower cutoff: m_snow - where(m_snow<0.1e-8_WP) m_snow=0.0_WP + if (use_cavity) then + ! upper cutoff SH: m_ice + where(m_ice>5.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)<0.0_WP) m_ice=5.0_WP + + ! upper cutoff NH: m_ice + where(m_ice>10.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)>0.0_WP) m_ice=10.0_WP + + ! upper cutoff: m_snow + where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP + + !___________________________________________________________________________ + ! lower cutoff: m_snow + !!PS where(m_snow<0.1e-8_WP) m_snow=0.0_WP + end if - ! upper cutoff: m_snow - where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP - !___________________________________________________________________________ #if defined (__oifs) where(ice_temp>273.15_WP) ice_temp=273.15_WP From 99ddb3dc47def0583ba11299dcde30b4660c5763 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 20 Apr 2021 18:11:05 +0200 Subject: [PATCH 30/40] check for numeric consistency --- src/oce_ale_tracer.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 7e57cfbc8..914abb2f1 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -524,7 +524,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) Ty =Ty *isredi Ty1=Ty1*isredi ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... - a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv*area(nz ,n)/areasvol(nz,n) +!!PS numerics a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv*area(nz ,n)/areasvol(nz,n) + a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv !!PS numerics *area(nz ,n)/areasvol(nz,n) c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=-a(nz)-c(nz)+hnode_new(nz,n) From dcf37362c2fc2f56ed505873390ea957fc6ca5e6 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 20 Apr 2021 22:36:10 +0200 Subject: [PATCH 31/40] try to fix numerical deviation from github testcase --- src/oce_ale_pressure_bv.F90 | 12 +++++++- src/oce_ale_tracer.F90 | 58 ++++++++++++++++++++++++++++++++----- 2 files changed, 62 insertions(+), 8 deletions(-) diff --git a/src/oce_ale_pressure_bv.F90 b/src/oce_ale_pressure_bv.F90 index c4ce88e0c..0feb5cd0a 100644 --- a/src/oce_ale_pressure_bv.F90 +++ b/src/oce_ale_pressure_bv.F90 @@ -2761,6 +2761,7 @@ subroutine sw_alpha_beta(TF1,SF1, mesh) use o_arrays use g_parsup use o_param + use g_comm_auto implicit none ! type(t_mesh), intent(in) , target :: mesh @@ -2815,6 +2816,8 @@ subroutine sw_alpha_beta(TF1,SF1, mesh) sw_alpha(nz,n) = a_over_b*sw_beta(nz,n) end do end do +call exchange_nod(sw_alpha) +call exchange_nod(sw_beta) end subroutine sw_alpha_beta ! ! @@ -2992,6 +2995,7 @@ SUBROUTINE density_linear(t, s, bulk_0, bulk_pz, bulk_pz2, rho_out, mesh) USE o_ARRAYS USE o_PARAM use g_PARSUP !, only: par_ex,pe_status +use g_config !, only: which_toy, toy_ocean IMPLICIT NONE real(kind=WP), intent(IN) :: t,s @@ -3005,7 +3009,13 @@ SUBROUTINE density_linear(t, s, bulk_0, bulk_pz, bulk_pz2, rho_out, mesh) bulk_0 = 1 bulk_pz = 0 bulk_pz2 = 0 - rho_out = density_0 + 0.8_WP*(s - 34.0_WP) - 0.2_WP*(t - 20.0_WP) + + IF((toy_ocean) .AND. (TRIM(which_toy)=="soufflet")) THEN + rho_out = density_0 - 0.00025_WP*(t - 10.0_WP)*density_0 + ELSE + rho_out = density_0 + 0.8_WP*(s - 34.0_WP) - 0.2*(t - 20.0_WP) + END IF + end subroutine density_linear ! ! diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 914abb2f1..a75cf140c 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -499,7 +499,17 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) !!PS v_adv =zinv*area(nz+1,n)/areasvol(nz,n) !!PS b(nz) =b(nz)+Wvel_i(nz, n)*zinv-min(0._WP, Wvel_i(nz+1, n))*v_adv !!PS c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv - v_adv =zinv*area(nz ,n)/areasvol(nz,n) + + !___________________________________________________________________ + ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) + ! is not exactly 1.0_WP when no cavity is used leads to not binary + ! identical results with github testcase + if (.not. use_cavity) then + !!PS v_adv =zinv!!PS numeric *area(nz ,n)/areasvol(nz,n) + v_adv =zinv + else + v_adv =zinv*area(nz ,n)/areasvol(nz,n) + end if b(nz) =b(nz)+Wvel_i(nz, n)*v_adv v_adv =zinv*area(nz+1,n)/areasvol(nz,n) @@ -523,9 +533,18 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) (zbar_n(nz+1)-Z_n(nz+1 ))*zinv2 *slope_tapered(3,nz+1,n)**2*Ki(nz+1,n) Ty =Ty *isredi Ty1=Ty1*isredi + ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... -!!PS numerics a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv*area(nz ,n)/areasvol(nz,n) - a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv !!PS numerics *area(nz ,n)/areasvol(nz,n) + !___________________________________________________________________ + ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) + ! is not exactly 1.0_WP when no cavity is used leads to not binary + ! identical results with github testcase + if (.not. use_cavity) then + a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv + else + a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv*area(nz ,n)/areasvol(nz,n) + end if + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=-a(nz)-c(nz)+hnode_new(nz,n) @@ -534,8 +553,15 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! update from the vertical advection if (do_wimpl) then - !!PS v_adv=zinv - v_adv=zinv*area(nz ,n)/areasvol(nz,n) + !___________________________________________________________________ + ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) + ! is not exactly 1.0_WP when no cavity is used leads to not binary + ! identical results with github testcase + if (.not. use_cavity) then + v_adv=zinv + else + v_adv=zinv*area(nz ,n)/areasvol(nz,n) + end if a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv !!PS v_adv=v_adv*areasvol(nz+1,n)/areasvol(nz,n) @@ -556,14 +582,32 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) (zbar_n(nz)-Z_n(nz)) *zinv1 *slope_tapered(3,nz,n)**2 *Ki(nz,n) Ty =Ty *isredi ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... - a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv*area(nz ,n)/areasvol(nz,n) + + !___________________________________________________________________ + ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) + ! is not exactly 1.0_WP when no cavity is used leads to not binary + ! identical results with github testcase + if (.not. use_cavity) then + a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv + else + a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv*area(nz ,n)/areasvol(nz,n) + end if c(nz)=0.0_WP b(nz)=-a(nz)+hnode_new(nz,n) ! update from the vertical advection if (do_wimpl) then !!PS v_adv=zinv - v_adv=zinv*area(nz ,n)/areasvol(nz,n) + !___________________________________________________________________ + ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) + ! is not exactly 1.0_WP when no cavity is used leads to not binary + ! identical results with github testcase + if (.not. use_cavity) then + !!PS v_adv=zinv!!PS numeric *area(nz ,n)/areasvol(nz,n) + v_adv=zinv + else + v_adv=zinv*area(nz ,n)/areasvol(nz,n) + end if a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv end if From 72f79faa71bd096968d239befe01dd66f6b99487 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 21 Apr 2021 16:41:00 +0200 Subject: [PATCH 32/40] simplify solution to overcome numerical problem when computing area(nz ,n)/areasvol(nz,n) --- src/oce_ale_tracer.F90 | 66 ++++++++++++++---------------------------- 1 file changed, 21 insertions(+), 45 deletions(-) diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 61505c100..b0ce14032 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -516,15 +516,10 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) !!PS c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv !___________________________________________________________________ - ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) - ! is not exactly 1.0_WP when no cavity is used leads to not binary - ! identical results with github testcase - if (.not. use_cavity) then - !!PS v_adv =zinv!!PS numeric *area(nz ,n)/areasvol(nz,n) - v_adv =zinv - else - v_adv =zinv*area(nz ,n)/areasvol(nz,n) - end if + ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for + ! numerical reasons, to gurante that area/areasvol in case of no + ! cavity is ==1.0_WP + v_adv =zinv* ( area(nz ,n)/areasvol(nz,n) ) b(nz) =b(nz)+Wvel_i(nz, n)*v_adv v_adv =zinv*area(nz+1,n)/areasvol(nz,n) @@ -551,15 +546,10 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... !___________________________________________________________________ - ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) - ! is not exactly 1.0_WP when no cavity is used leads to not binary - ! identical results with github testcase - if (.not. use_cavity) then - a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv - else - a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv*area(nz ,n)/areasvol(nz,n) - end if - + ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for + ! numerical reasons, to gurante that area/areasvol in case of no + ! cavity is ==1.0_WP + a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv* ( area(nz ,n)/areasvol(nz,n) ) c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) b(nz)=-a(nz)-c(nz)+hnode_new(nz,n) @@ -568,15 +558,11 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! update from the vertical advection if (do_wimpl) then - !___________________________________________________________________ - ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) - ! is not exactly 1.0_WP when no cavity is used leads to not binary - ! identical results with github testcase - if (.not. use_cavity) then - v_adv=zinv - else - v_adv=zinv*area(nz ,n)/areasvol(nz,n) - end if + !_______________________________________________________________ + ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for + ! numerical reasons, to gurante that area/areasvol in case of no + ! cavity is ==1.0_WP + v_adv=zinv* ( area(nz ,n)/areasvol(nz,n) ) a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv !!PS v_adv=v_adv*areasvol(nz+1,n)/areasvol(nz,n) @@ -599,30 +585,20 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... !___________________________________________________________________ - ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) - ! is not exactly 1.0_WP when no cavity is used leads to not binary - ! identical results with github testcase - if (.not. use_cavity) then - a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv - else - a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv*area(nz ,n)/areasvol(nz,n) - end if + ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for + ! numerical reasons, to gurante that area/areasvol in case of no + ! cavity is ==1.0_WP + a(nz)=-(Kv(nz,n)+Ty)*zinv1*zinv* ( area(nz ,n)/areasvol(nz,n) ) c(nz)=0.0_WP b(nz)=-a(nz)+hnode_new(nz,n) ! update from the vertical advection if (do_wimpl) then - !!PS v_adv=zinv !___________________________________________________________________ - ! do this here for numerical reasons since area(nz ,n)/areasvol(nz,n) - ! is not exactly 1.0_WP when no cavity is used leads to not binary - ! identical results with github testcase - if (.not. use_cavity) then - !!PS v_adv=zinv!!PS numeric *area(nz ,n)/areasvol(nz,n) - v_adv=zinv - else - v_adv=zinv*area(nz ,n)/areasvol(nz,n) - end if + ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for + ! numerical reasons, to gurante that area/areasvol in case of no + ! cavity is ==1.0_WP + v_adv=zinv* ( area(nz ,n)/areasvol(nz,n) ) a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv end if From 4a95870acaaa3fd89c93c8ce3b6dd4ce81feaa8a Mon Sep 17 00:00:00 2001 From: Dmitry Sidorenko Date: Wed, 28 Apr 2021 10:17:02 +0200 Subject: [PATCH 33/40] density / buoyancy flux computation is addad (no cost, always computed) --- src/gen_modules_diag.F90 | 12 ++++++------ src/ice_oce_coupling.F90 | 7 ++++++- src/io_meandata.F90 | 4 +++- src/oce_modules.F90 | 1 + src/oce_setup_step.F90 | 6 ++++-- 5 files changed, 20 insertions(+), 10 deletions(-) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index 7dc0540b8..22498f7b3 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -18,7 +18,7 @@ module diagnostics public :: ldiag_solver, lcurt_stress_surf, ldiag_energy, ldiag_dMOC, ldiag_DVD, ldiag_forc, ldiag_salt3D, ldiag_curl_vel3, diag_list, & compute_diagnostics, rhs_diag, curl_stress_surf, curl_vel3, wrhof, rhof, & u_x_u, u_x_v, v_x_v, v_x_w, u_x_w, dudx, dudy, dvdx, dvdy, dudz, dvdz, utau_surf, utau_bott, av_dudz_sq, av_dudz, av_dvdz, stress_bott, u_surf, v_surf, u_bott, v_bott, & - std_dens_min, std_dens_max, std_dens_N, std_dens, std_dens_UVDZ, std_dens_DIV, std_dens_Z, std_dens_dVdT, std_dens_flux, dens_flux, & + std_dens_min, std_dens_max, std_dens_N, std_dens, std_dens_UVDZ, std_dens_DIV, std_dens_Z, std_dens_dVdT, std_dens_flux, dens_flux_e, & compute_diag_dvd_2ndmoment_klingbeil_etal_2014, compute_diag_dvd_2ndmoment_burchard_etal_2008, compute_diag_dvd ! Arrays used for diagnostics, some shall be accessible to the I/O ! 1. solver diagnostics: A*x=rhs? @@ -50,7 +50,7 @@ module diagnostics real(kind=WP), save, target :: std_dd(std_dens_N-1) real(kind=WP), save, target :: std_dens_min=1030., std_dens_max=1040. real(kind=WP), save, allocatable, target :: std_dens_UVDZ(:,:,:), std_dens_flux(:,:,:), std_dens_dVdT(:,:), std_dens_DIV(:,:), std_dens_Z(:,:) - real(kind=WP), save, allocatable, target :: dens_flux(:) + real(kind=WP), save, allocatable, target :: dens_flux_e(:) logical :: ldiag_solver =.false. logical :: lcurt_stress_surf=.false. @@ -408,7 +408,7 @@ subroutine diag_densMOC(mode, mesh) allocate(std_dens_VOL2( std_dens_N, myDim_elem2D)) allocate(std_dens_flux(3,std_dens_N, myDim_elem2D)) allocate(std_dens_Z ( std_dens_N, myDim_elem2D)) - allocate(dens_flux(elem2D)) + allocate(dens_flux_e(elem2D)) allocate(aux (nl-1)) allocate(dens (nl)) allocate(el_depth(nl)) @@ -438,7 +438,7 @@ subroutine diag_densMOC(mode, mesh) std_dens_UVDZ=0. std_dens_w =0.! temporat thing for wiighting (ageraging) mean fields within a bin std_dens_flux=0. - dens_flux =0. + dens_flux_e =0. std_dens_VOL2=0. std_dens_DIV =0. std_dens_Z =0. @@ -450,11 +450,11 @@ subroutine diag_densMOC(mode, mesh) ! density flux on elements (although not related to binning it might be usefull for diagnostic and to verify the consistency) do jj=1,3 - dens_flux(elem)= dens_flux(elem) + (sw_alpha(ulevels_nod2D(elnodes(jj)),elnodes(jj)) * heat_flux_in(elnodes(jj)) / vcpw + & + dens_flux_e(elem)=dens_flux_e(elem) + (sw_alpha(ulevels_nod2D(elnodes(jj)),elnodes(jj)) * heat_flux_in(elnodes(jj)) / vcpw + & sw_beta(ulevels_nod2D(elnodes(jj)),elnodes(jj)) * (relax_salt (elnodes(jj)) + water_flux(elnodes(jj)) * & tr_arr(ulevels_nod2D(elnodes(jj)),elnodes(jj),2))) end do - dens_flux(elem) = dens_flux(elem)/3.0_WP + dens_flux_e(elem) =dens_flux_e(elem)/3.0_WP ! density_dmoc is the sigma_2 density given at nodes. it is computed in oce_ale_pressure_bv do nz=nzmin, nzmax-1 aux(nz)=sum(density_dmoc(nz, elnodes))/3.-1000. diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 2e8bd57fc..510f432b1 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -250,7 +250,12 @@ subroutine oce_fluxes(mesh) end if virtual_salt=virtual_salt-net/ocean_area end if - + + where (ulevels_nod2d == 1) + dens_flux=sw_alpha(1,:) * heat_flux_in / vcpw + sw_beta(1, :) * (relax_salt + water_flux * tr_arr(1,:,2)) + elsewhere + dens_flux=0.0_WP + end where !___________________________________________________________________________ ! balance SSS restoring to climatology if (use_cavity) then diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 23a827d8b..86c9fdf91 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -207,6 +207,8 @@ subroutine ini_mean_io(mesh) call def_stream(nod2D, myDim_nod2D, 'alpha', 'thermal expansion', 'none', sw_alpha(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) CASE ('beta ') call def_stream(nod2D, myDim_nod2D, 'beta', 'saline contraction', 'none', sw_beta (1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) +CASE ('dens_flux ') + call def_stream(nod2D, myDim_nod2D , 'dflux', 'density flux', 'kg/(m3*s)', dens_flux(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) CASE ('runoff ') sel_forcvar(10)= 1 call def_stream(nod2D, myDim_nod2D, 'runoff', 'river runoff', 'none', runoff(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) @@ -358,7 +360,7 @@ subroutine ini_mean_io(mesh) call def_stream((/std_dens_N, nod2D /), (/std_dens_N, myDim_nod2D/), 'std_dens_DIV', 'm3/s', 'm3/s' ,std_dens_DIV(:,:), 1, 'y', i_real4, mesh) call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_dens_Z', 'm', 'm' ,std_dens_Z(:,:), 1, 'y', i_real4, mesh) call def_stream((/nl-1, nod2D /), (/nl-1, myDim_nod2D /), 'density_dMOC', 'density' , 'm', density_dmoc(:,:), 1, 'y', i_real4, mesh) - call def_stream(elem2D, myDim_elem2D , 'density_flux', 'density' , 'm', dens_flux(:), 1, 'y', i_real4, mesh) + call def_stream(elem2D, myDim_elem2D , 'density_flux_e', 'density flux at elems ', 'm', dens_flux_e(:), 1, 'y', i_real4, mesh) end if !___________________________________________________________________________________________________________________________________ CASE ('pgf_x ') diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index 6ac8506bc..6597bcee2 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -272,6 +272,7 @@ MODULE o_ARRAYS real(kind=WP), allocatable,dimension(:,:,:) :: slope_tapered real(kind=WP), allocatable,dimension(:,:,:) :: sigma_xy real(kind=WP), allocatable,dimension(:,:) :: sw_beta, sw_alpha +real(kind=WP), allocatable,dimension(:) :: dens_flux !real(kind=WP), allocatable,dimension(:,:,:) :: tsh, tsv, tsh_nodes !real(kind=WP), allocatable,dimension(:,:) :: hd_flux,vd_flux !Isoneutral diffusivities (or xy diffusivities if Redi=.false) diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index 0d9dc083b..02aa93f91 100755 --- a/src/oce_setup_step.F90 +++ b/src/oce_setup_step.F90 @@ -344,8 +344,10 @@ SUBROUTINE array_setup(mesh) ! alpha and beta in the EoS allocate(sw_beta(nl-1, node_size), sw_alpha(nl-1, node_size)) -sw_beta=0.0_WP -sw_alpha=0.0_WP +allocate(dens_flux(node_size)) +sw_beta =0.0_WP +sw_alpha =0.0_WP +dens_flux=0.0_WP if (Fer_GM) then allocate(fer_c(node_size),fer_scal(node_size), fer_gamma(2, nl, node_size), fer_K(nl, node_size)) From 38b437b6387271c36a7cc523551437d86792e532 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 29 Apr 2021 12:26:25 +0200 Subject: [PATCH 34/40] fix bug of sea-ice height accumulation along the cavity-ocean edge by ensuring uice,vice == 0 boundary condition along the cavity edge --- src/ice_EVP.F90 | 57 ++++++++++++++++++++++---------- src/ice_fct.F90 | 11 ------- src/ice_maEVP.F90 | 74 ++++++++++++++++++++++++++++++++---------- src/ice_setup_step.F90 | 15 +++++++++ src/ice_thermo_oce.F90 | 28 ++++++++-------- 5 files changed, 124 insertions(+), 61 deletions(-) diff --git a/src/ice_EVP.F90 b/src/ice_EVP.F90 index 2cf727698..fcfd8d224 100755 --- a/src/ice_EVP.F90 +++ b/src/ice_EVP.F90 @@ -58,7 +58,6 @@ subroutine stress_tensor(ice_strength, mesh) do el=1,myDim_elem2D !__________________________________________________________________________ ! if element contains cavity node skip it - !!PS if ( any(ulevels_nod2d(elem2D_nodes(:,el)) > 1) ) cycle if (ulevels(el) > 1) cycle ! ===== Check if there is ice on elem @@ -166,7 +165,6 @@ subroutine stress_tensor_no1(ice_strength, mesh) do el=1,myDim_elem2D !__________________________________________________________________________ ! if element contains cavity node skip it - !!PS if ( any(ulevels_nod2d(elem2D_nodes(:,el)) > 1) ) cycle if (ulevels(el) > 1) cycle ! ===== Check if there is ice on elem @@ -203,7 +201,10 @@ subroutine stress_tensor_no1(ice_strength, mesh) ! ===== if delta is too small or zero, viscosity will too large (unlimited) ! (limit delta_inv) - delta_inv = 1.0_WP/max(delta,delta_min) + delta_inv = 1.0_WP/max(delta,delta_min) + +!!PS delta_inv = delta/(delta+delta_min) + zeta = ice_strength(el)*delta_inv ! ===== Limiting pressure/Delta (zeta): it may still happen that pressure/Delta ! is too large in some regions and CFL criterion is violated. @@ -218,7 +219,7 @@ subroutine stress_tensor_no1(ice_strength, mesh) !end if zeta = zeta*Tevp_inv - + r1 = zeta*(eps11(el)+eps22(el)) - ice_strength(el)*Tevp_inv r2 = zeta*(eps11(el)-eps22(el))*vale r3 = zeta*eps12(el)*vale @@ -247,6 +248,7 @@ subroutine stress2rhs_e(mesh) USE i_therm_param USE i_arrays USE g_PARSUP +use g_config, only: use_cavity IMPLICIT NONE @@ -266,7 +268,13 @@ subroutine stress2rhs_e(mesh) DO ed=1,myDim_edge2D ednodes=edges(:,ed) el=edge_tri(:,ed) - if(myList_edge2D(ed)>edge2D_in) cycle + if(myList_edge2D(ed)>edge2D_in) cycle + + ! stress boundary condition at ocean cavity boundary edge ==0 + if (use_cavity) then + if ( (ulevels(el(1))>1) .or. ( el(2)>0 .and. ulevels(el(2))>1) ) cycle + end if + ! elements on both sides uc = - sigma12(el(1))*edge_cross_dxdy(1,ed) + sigma11(el(1))*edge_cross_dxdy(2,ed) & + sigma12(el(2))*edge_cross_dxdy(3,ed) - sigma11(el(2))*edge_cross_dxdy(4,ed) @@ -345,7 +353,6 @@ subroutine stress2rhs(inv_areamass, ice_strength, mesh) ! if (any(m_ice(elnodes)<= 0.) .or. any(a_ice(elnodes) <=0.)) CYCLE !____________________________________________________________________________ ! if element contains cavity node skip it - !!OS if ( any(ulevels_nod2d(elem2D_nodes(:,el)) > 1) ) cycle if (ulevels(el) > 1) cycle !____________________________________________________________________________ @@ -488,7 +495,6 @@ subroutine EVPdynamics(mesh) elnodes = elem2D_nodes(:,el) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle !_______________________________________________________________________ @@ -542,7 +548,6 @@ subroutine EVPdynamics(mesh) elnodes = elem2D_nodes(:,el) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle !_______________________________________________________________________ @@ -601,7 +606,7 @@ subroutine EVPdynamics(mesh) do n=1,myDim_nod2D !_________________________________________________________________________ - ! if cavity ndoe skip it + ! if cavity node skip it if ( ulevels_nod2d(n)>1 ) cycle !_________________________________________________________________________ @@ -631,15 +636,31 @@ subroutine EVPdynamics(mesh) end if end do - DO ed=1,myDim_edge2D - ! boundary conditions - if(myList_edge2D(ed) > edge2D_in) then - U_ice(edges(1:2,ed))=0.0_WP - V_ice(edges(1:2,ed))=0.0_WP - endif - end do - - call exchange_nod(U_ice,V_ice) + + !___________________________________________________________________________ + ! apply sea ice velocity boundary condition + DO ed=1,myDim_edge2D + !_______________________________________________________________________ + ! apply coastal sea ice velocity boundary conditions + if(myList_edge2D(ed) > edge2D_in) then + U_ice(edges(1:2,ed))=0.0_WP + V_ice(edges(1:2,ed))=0.0_WP + endif + + !_______________________________________________________________________ + ! apply sea ice velocity boundary conditions at cavity-ocean edge + if (use_cavity) then + if ( (ulevels(edge_tri(1,ed))>1) .or. & + ( edge_tri(2,ed)>0 .and. ulevels(edge_tri(2,ed))>1) ) then + U_ice(edges(1:2,ed))=0.0_WP + V_ice(edges(1:2,ed))=0.0_WP + end if + end if + + end do + + !___________________________________________________________________________ + call exchange_nod(U_ice,V_ice) END DO diff --git a/src/ice_fct.F90 b/src/ice_fct.F90 index 501b9ac36..fb91b370c 100755 --- a/src/ice_fct.F90 +++ b/src/ice_fct.F90 @@ -66,8 +66,6 @@ subroutine ice_TG_rhs(mesh) !_______________________________________________________________________ ! if cavity element skip it if (ulevels(elem)>1) cycle - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 - !derivatives dx=gradient_sca(1:3,elem) @@ -82,7 +80,6 @@ subroutine ice_TG_rhs(mesh) diff=ice_diff*sqrt(elem_area(elem)/scale_area) DO n=1,3 row=elnodes(n) -!!PS if (ulevels_nod2D(row)>1) cycle DO q = 1,3 !entries(q)= vol*dt*((dx(n)*um+dy(n)*vm)/3.0_WP - & ! diff*(dx(n)*dx(q)+ dy(n)*dy(q))- & @@ -375,7 +372,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) !_______________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 !_______________________________________________________________________ @@ -483,7 +479,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) !_______________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 !_______________________________________________________________________ @@ -532,7 +527,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) !_______________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 !_______________________________________________________________________ @@ -559,7 +553,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) !___________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 do q=1,3 @@ -579,7 +572,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) !___________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 do q=1,3 @@ -599,7 +591,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) !___________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 do q=1,3 @@ -619,7 +610,6 @@ subroutine ice_fem_fct(tr_array_id, mesh) elnodes=elem2D_nodes(:,elem) !___________________________________________________________________ ! if cavity cycle over - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 if(ulevels(elem)>1) cycle !LK89140 do q=1,3 @@ -759,7 +749,6 @@ subroutine ice_TG_rhs_div(mesh) !___________________________________________________________________________ ! if cavity element skip it if (ulevels(elem)>1) cycle - if(any(ulevels_nod2D(elnodes)>1)) cycle !LK89140 !derivatives dx=gradient_sca(1:3,elem) diff --git a/src/ice_maEVP.F90 b/src/ice_maEVP.F90 index 365307566..f5aedca3d 100644 --- a/src/ice_maEVP.F90 +++ b/src/ice_maEVP.F90 @@ -66,7 +66,6 @@ subroutine stress_tensor_m(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle msum=sum(m_ice(elnodes))*val3 @@ -162,7 +161,6 @@ subroutine ssh2rhs(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle !_______________________________________________________________________ @@ -189,7 +187,6 @@ subroutine ssh2rhs(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle vol=elem_area(elem) @@ -239,7 +236,6 @@ subroutine stress2rhs_m(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle if(sum(a_ice(elnodes)) < 0.01_WP) cycle !DS @@ -355,7 +351,6 @@ subroutine EVPdynamics_m(mesh) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle !_______________________________________________________________________ @@ -384,7 +379,6 @@ subroutine EVPdynamics_m(mesh) elnodes=elem2D_nodes(:,el) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle vol=elem_area(el) @@ -431,7 +425,6 @@ subroutine EVPdynamics_m(mesh) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(el) > 1) cycle msum=sum(m_ice(elnodes))*val3 @@ -468,7 +461,7 @@ subroutine EVPdynamics_m(mesh) do el=1,myDim_elem2D !__________________________________________________________________________ if (ulevels(el)>1) cycle - if ( any(ulevels_nod2d(elnodes)>1) ) cycle + !__________________________________________________________________________ if(ice_el(el)) then @@ -538,7 +531,7 @@ subroutine EVPdynamics_m(mesh) (sigma12(el)*dx(3)+sigma22(el)*dy(3) - sigma11(el)*meancos) ! metrics end if end if - end do + end do ! --> do el=1,myDim_elem2D do i=1, myDim_nod2d !__________________________________________________________________________ @@ -568,20 +561,44 @@ subroutine EVPdynamics_m(mesh) u_ice_aux(i) = det*((1.0_WP+beta_evp+drag)*rhsu +rdt*coriolis_node(i)*rhsv) v_ice_aux(i) = det*((1.0_WP+beta_evp+drag)*rhsv -rdt*coriolis_node(i)*rhsu) end if - end do + end do ! --> do i=1, myDim_nod2d - call exchange_nod_begin(u_ice_aux, v_ice_aux) + !___________________________________________________________________________ + ! apply sea ice velocity boundary condition + do ed=1,myDim_edge2D + !_______________________________________________________________________ + ! apply coastal sea ice velocity boundary conditions + if(myList_edge2D(ed) > edge2D_in) then + u_ice_aux(edges(:,ed))=0.0_WP + v_ice_aux(edges(:,ed))=0.0_WP + end if + + !_______________________________________________________________________ + ! apply sea ice velocity boundary conditions at cavity-ocean edge + if (use_cavity) then + if ( (ulevels(edge_tri(1,ed))>1) .or. & + ( edge_tri(2,ed)>0 .and. ulevels(edge_tri(2,ed))>1) ) then + u_ice_aux(edges(1:2,ed))=0.0_WP + v_ice_aux(edges(1:2,ed))=0.0_WP + end if + end if + end do ! --> do ed=1,myDim_edge2D + + !___________________________________________________________________________ + call exchange_nod_begin(u_ice_aux, v_ice_aux) - do row=1, myDim_nod2d + do row=1, myDim_nod2d u_rhs_ice(row)=0.0_WP v_rhs_ice(row)=0.0_WP - end do + end do - call exchange_nod_end - end do + call exchange_nod_end + + end do ! --> do shortstep=1, steps u_ice=u_ice_aux v_ice=v_ice_aux + end subroutine EVPdynamics_m ! ! @@ -624,7 +641,6 @@ subroutine find_alpha_field_a(mesh) elnodes=elem2D_nodes(:,elem) !_______________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle msum=sum(m_ice(elnodes))*val3 @@ -700,7 +716,6 @@ subroutine stress_tensor_a(mesh) do elem=1,myDim_elem2D !__________________________________________________________________________ ! if element has any cavity node skip it - if ( any(ulevels_nod2d(elnodes)>1) ) cycle if (ulevels(elem) > 1) cycle !__________________________________________________________________________ @@ -782,6 +797,7 @@ subroutine EVPdynamics_a(mesh) use o_PARAM use i_therm_param use g_parsup +use g_config, only: use_cavity use g_comm_auto use ice_maEVP_interfaces @@ -830,7 +846,7 @@ subroutine EVPdynamics_a(mesh) rhsv=v_ice(i)+drag*v_w(i)+rdt*(inv_thickness*stress_atmice_y(i)+v_rhs_ice(i)) rhsu=beta_evp_array(i)*u_ice_aux(i)+rhsu - rhsv=beta_evp_array(i)*v_ice_aux(i)+rhsv + rhsv=beta_evp_array(i)*v_ice_aux(i)+rhsv !solve (Coriolis and water stress are treated implicitly) fc=rdt*coriolis_node(i) det=(1.0_WP+beta_evp_array(i)+drag)**2+fc**2 @@ -838,6 +854,28 @@ subroutine EVPdynamics_a(mesh) u_ice_aux(i)=det*((1.0_WP+beta_evp_array(i)+drag)*rhsu+fc*rhsv) v_ice_aux(i)=det*((1.0_WP+beta_evp_array(i)+drag)*rhsv-fc*rhsu) end do + + !___________________________________________________________________________ + ! apply sea ice velocity boundary condition + do ed=1,myDim_edge2D + !_______________________________________________________________________ + ! apply coastal sea ice velocity boundary conditions + if(myList_edge2D(ed) > edge2D_in) then + u_ice_aux(edges(:,ed))=0.0_WP + v_ice_aux(edges(:,ed))=0.0_WP + end if + + !_______________________________________________________________________ + ! apply sea ice velocity boundary conditions at cavity-ocean edge + if (use_cavity) then + if ( (ulevels(edge_tri(1,ed))>1) .or. & + ( edge_tri(2,ed)>0 .and. ulevels(edge_tri(2,ed))>1) ) then + u_ice_aux(edges(1:2,ed))=0.0_WP + v_ice_aux(edges(1:2,ed))=0.0_WP + end if + end if + end do ! --> do ed=1,myDim_edge2D + call exchange_nod(u_ice_aux, v_ice_aux) end do diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90 index 1d56e957d..183b2ba5e 100755 --- a/src/ice_setup_step.F90 +++ b/src/ice_setup_step.F90 @@ -242,6 +242,21 @@ subroutine ice_timestep(step, mesh) if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call thermodynamics...'//achar(27)//'[0m' call thermodynamics(mesh) #endif /* (__icepack) */ + + + do i=1,myDim_nod2D+eDim_nod2D + if ( ( U_ice(i)/=0.0_WP .and. mesh%ulevels_nod2d(i)>1) .or. (V_ice(i)/=0.0_WP .and. mesh%ulevels_nod2d(i)>1) ) then + write(*,*) " --> found cavity velocity /= 0.0_WP , ", mype + write(*,*) " ulevels_nod2d(n) = ", mesh%ulevels_nod2d(i) + write(*,*) " U_ice(n) = ", U_ice(i) + write(*,*) " V_ice(n) = ", V_ice(i) + write(*,*) + end if + end do + + + + t3=MPI_Wtime() rtime_ice = rtime_ice + (t3-t0) rtime_tot = rtime_tot + (t3-t0) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index ab1e609ac..9e12224dc 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -36,20 +36,20 @@ subroutine cut_off(mesh) end where - if (use_cavity) then - ! upper cutoff SH: m_ice - where(m_ice>5.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)<0.0_WP) m_ice=5.0_WP - - ! upper cutoff NH: m_ice - where(m_ice>10.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)>0.0_WP) m_ice=10.0_WP - - ! upper cutoff: m_snow - where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP - - !___________________________________________________________________________ - ! lower cutoff: m_snow - !!PS where(m_snow<0.1e-8_WP) m_snow=0.0_WP - end if +!!PS if (use_cavity) then +!!PS ! upper cutoff SH: m_ice +!!PS where(m_ice>5.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)<0.0_WP) m_ice=5.0_WP +!!PS +!!PS ! upper cutoff NH: m_ice +!!PS where(m_ice>10.0_WP .and. ulevels_nod2d==1 .and. geo_coord_nod2D(2,:)>0.0_WP) m_ice=10.0_WP +!!PS +!!PS ! upper cutoff: m_snow +!!PS where(m_snow>2.5_WP .and. ulevels_nod2d==1) m_snow=2.5_WP +!!PS +!!PS !___________________________________________________________________________ +!!PS ! lower cutoff: m_snow +!!PS !!PS where(m_snow<0.1e-8_WP) m_snow=0.0_WP +!!PS end if !___________________________________________________________________________ #if defined (__oifs) From 2fae97e183fef007eb711e206889db4bf5ee3898 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 29 Apr 2021 13:19:24 +0200 Subject: [PATCH 35/40] trim filename in src/gen_surface_forcing.F90 since it can be very long now --- src/gen_surface_forcing.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gen_surface_forcing.F90 b/src/gen_surface_forcing.F90 index 151e2c8ca..eb9b5eb5d 100644 --- a/src/gen_surface_forcing.F90 +++ b/src/gen_surface_forcing.F90 @@ -676,7 +676,7 @@ SUBROUTINE getcoeffld(fld_idx, rdate, mesh) delta_t = 1.0_wp if (mype==0) then write(*,*) 'WARNING: no temporal extrapolation into future (nearest neighbour is used): ', trim(var_name), ' !' - write(*,*) file_name + write(*,*) trim(file_name) write(*,*) nc_time(1), nc_time(nc_Ntime), now_date end if elseif (t_indx < 1) then ! NO extrapolation back in time @@ -685,7 +685,7 @@ SUBROUTINE getcoeffld(fld_idx, rdate, mesh) delta_t = 1.0_wp if (mype==0) then write(*,*) 'WARNING: no temporal extrapolation back in time (nearest neighbour is used): ', trim(var_name), ' !' - write(*,*) file_name + write(*,*) trim(file_name) write(*,*) nc_time(1), nc_time(nc_Ntime), now_date end if end if From e9fea1b30f51e372f2e4514ba0d6003a4054562b Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 29 Apr 2021 15:43:10 +0200 Subject: [PATCH 36/40] do elemreducelvl and elemfixlvl array as logicals instead of integer, should reduce memory demand for large meshes --- src/fvom_init.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 79f49530b..1b57fa45a 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -887,7 +887,7 @@ subroutine find_levels_cavity(mesh) real(kind=WP) :: dmean character(MAX_PATH) :: file_name integer, allocatable, dimension(:,:) :: numelemtonode, idxelemtonode - integer, allocatable, dimension(:) :: elemreducelvl, elemfixlvl + logical, allocatable, dimension(:) :: elemreducelvl, elemfixlvl type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" !___________________________________________________________________________ @@ -955,10 +955,10 @@ subroutine find_levels_cavity(mesh) ! outer iteration loop count_iter2 = 0 exit_flag2 = 0 - elemfixlvl = 0 + elemfixlvl = .false. do while((exit_flag2==0) .and. (count_iter2 in case make the levels of all sorounding ! triangles shallower - if ( (nlevels(elem)-(nz+1))>=3 .and. elemreducelvl(elem)==0 .and. elemfixlvl(elem)==0) then + if ( (nlevels(elem)-(nz+1))>=3 .and. elemreducelvl(elem)==.false. .and. elemfixlvl(elem)==.false.) then ulevels(elem)=nz+1 else ! --> can not increase depth anymore to eleminate isolated @@ -1021,7 +1021,7 @@ subroutine find_levels_cavity(mesh) ! to nz idx = minloc(ulevels(elems)-nz, 1, MASK=( (elems>0) .and. ((ulevels(elems)-nz)>0) ) ) ulevels(elems(idx)) = nz-1 - elemreducelvl(elems(idx)) = elemreducelvl(elems(idx))+1 + elemreducelvl(elems(idx)) = .true. end if !force recheck for all current ocean cells @@ -1132,7 +1132,7 @@ subroutine find_levels_cavity(mesh) elem=nod_in_elem2D(j,node) if (ulevels(elem)>nz) then ulevels(elem) = nz - elemfixlvl(elem) = elemfixlvl(elem)+1 + elemfixlvl(elem) = .true. end if end do end if @@ -1161,7 +1161,7 @@ subroutine find_levels_cavity(mesh) !_______________________________________________________________________ end do - deallocate(elemreducelvl) + deallocate(elemreducelvl,elemfixlvl) deallocate(numelemtonode,idxelemtonode) !___________________________________________________________________________ From 401a225dadbad80964d87ac5f06857880e02b9ce Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 29 Apr 2021 16:36:48 +0200 Subject: [PATCH 37/40] reduce size of numelemtonode idxelemtonode in fvom_init.F90 --- src/fvom_init.F90 | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 1b57fa45a..98eea6b45 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -886,7 +886,7 @@ subroutine find_levels_cavity(mesh) integer :: exit_flag1, count_iter, max_iter=1000, exit_flag2, count_iter2, max_iter2=10 real(kind=WP) :: dmean character(MAX_PATH) :: file_name - integer, allocatable, dimension(:,:) :: numelemtonode, idxelemtonode + integer, allocatable, dimension(:) :: numelemtonode, idxelemtonode logical, allocatable, dimension(:) :: elemreducelvl, elemfixlvl type(t_mesh), intent(inout), target :: mesh #include "associate_mesh_ini.h" @@ -949,7 +949,7 @@ subroutine find_levels_cavity(mesh) ! possible in FESOM2.0 ! loop over all cavity levels allocate(elemreducelvl(elem2d),elemfixlvl(elem2d)) - allocate(numelemtonode(nl,nod2D),idxelemtonode(nl,nod2D)) + allocate(numelemtonode(nl),idxelemtonode(nl)) !___________________________________________________________________________ ! outer iteration loop @@ -1093,35 +1093,35 @@ subroutine find_levels_cavity(mesh) !_______________________________________________________________________ ! compute how many triangle elements contribute to every vertice in every layer - numelemtonode=0 - idxelemtonode=0 + count_iter=0 do node=1, nod2D + !___________________________________________________________________ + numelemtonode=0 + idxelemtonode=0 + + !___________________________________________________________________ + ! compute how many triangle elements contribute to vertice in every layer do j=1,nod_in_elem2D_num(node) elem=nod_in_elem2D(j,node) do nz=ulevels(elem),nlevels(elem)-1 - numelemtonode(nz,node) = numelemtonode(nz,node) + 1 - idxelemtonode(nz,node) = elem + numelemtonode(nz) = numelemtonode(nz) + 1 + idxelemtonode(nz) = elem end do end do - end do ! --> do node=1, nod2D - - !_______________________________________________________________________ - ! check if every vertice in every layer should be connected to at least - ! two triangle elements ! - count_iter=0 - do node=1, nod2D !___________________________________________________________________ + ! check if every vertice in every layer should be connected to at least + ! two triangle elements ! do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 !_______________________________________________________________ ! nodes has zero neighbouring triangles and is completely isolated ! need to adapt ulevels by hand --> inflicts another outher ! iteration loop (exit_flag2=0) - if (numelemtonode(nz,node)==0) then + if (numelemtonode(nz)==0) then exit_flag2 = 0 count_iter = count_iter+1 - write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz,node) ,' triangle: n=', node, ', nz=',nz + write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz) ,' triangle: n=', node, ', nz=',nz !___________________________________________________________ ! if node has no neighboring triangle somewhere in the middle ! of the water column at nz (can happen but seldom) than set @@ -1140,10 +1140,10 @@ subroutine find_levels_cavity(mesh) !_______________________________________________________________ ! nodes has just one neighbouring triangle --> but needs two --> ! inflicts another outher iteration loop (exit_flag2=0) - if (numelemtonode(nz,node)==1) then + if (numelemtonode(nz)==1) then exit_flag2 = 0 count_iter = count_iter+1 - write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz,node) ,' triangle: n=', node, ', nz=',nz + write(*,"( A, I1, A, I7, A, I3)") ' -[check]->: node has only ', numelemtonode(nz) ,' triangle: n=', node, ', nz=',nz end if end do ! --> do nz = ulevels_nod2D(node), nlevels_nod2D(node)-1 From 229b8d6fa3e67b2f52abfb365208dce6a019d891 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 29 Apr 2021 16:37:58 +0200 Subject: [PATCH 38/40] reduce size of numelemtonode in oce_mesh.F90: subroutine find_levels_cavity and ensure deallocation --- src/oce_mesh.F90 | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index fb91adc32..7d37a8372 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -916,7 +916,7 @@ subroutine find_levels_cavity(mesh) real(kind=WP) :: t0, t1 logical :: file_exist=.False. integer :: elem, elnodes(3), ule, uln(3), node, j, nz - integer, allocatable, dimension(:,:) :: numelemtonode + integer, allocatable, dimension(:) :: numelemtonode !NR Cannot include the pointers before the targets are allocated... !NR #include "associate_mesh.h" @@ -1297,27 +1297,29 @@ subroutine find_levels_cavity(mesh) !___________________________________________________________________________ - allocate(numelemtonode(mesh%nl,myDim_nod2d+eDim_nod2D)) - numelemtonode=0 + allocate(numelemtonode(mesh%nl)) do node=1, myDim_nod2D+eDim_nod2D + numelemtonode=0 + !_______________________________________________________________________ do j=1,mesh%nod_in_elem2D_num(node) elem=mesh%nod_in_elem2D(j,node) do nz=mesh%ulevels(elem),mesh%nlevels(elem)-1 - numelemtonode(nz,node) = numelemtonode(nz,node) + 1 + numelemtonode(nz) = numelemtonode(nz) + 1 end do end do - end do - - ! check how many triangle elements contribute to every vertice in every layer - ! every vertice in every layer should be connected to at least two triangle - ! elements ! - do node=1, myDim_nod2D+eDim_nod2D - do nz=1,mesh%nl - if (numelemtonode(nz,node)== 1) then + + !_______________________________________________________________________ + ! check how many triangle elements contribute to every vertice in every layer + ! every vertice in every layer should be connected to at least two triangle + ! elements ! + do nz=mesh%ulevels_nod2D(node),mesh%nlevels_nod2D(node)-1 + if (numelemtonode(nz)== 1) then write(*,*) 'ERROR A: found vertice with just one triangle:', mype, node, nz end if end do - end do + + end do + deallocate(numelemtonode) end subroutine find_levels_cavity ! From 03e7aacce2a059a3a507c7588735701b8af9448c Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 29 Apr 2021 16:44:38 +0200 Subject: [PATCH 39/40] exchange logical comparison == with .eqv. --- src/fvom_init.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/fvom_init.F90 b/src/fvom_init.F90 index 98eea6b45..14bc79742 100755 --- a/src/fvom_init.F90 +++ b/src/fvom_init.F90 @@ -1011,7 +1011,9 @@ subroutine find_levels_cavity(mesh) ! 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)==.false. .and. elemfixlvl(elem)==.false.) then + if ( (nlevels(elem)-(nz+1))>=3 .and. & + elemreducelvl(elem) .eqv. .false. .and. & + elemfixlvl(elem) .eqv. .false.) then ulevels(elem)=nz+1 else ! --> can not increase depth anymore to eleminate isolated From 745949a89d7a879f8abd951782bef31dd972023c Mon Sep 17 00:00:00 2001 From: Dmitry Sidorenko Date: Fri, 30 Apr 2021 13:04:28 +0200 Subject: [PATCH 40/40] fixed a bug in PPM. Another thing is that due to the lagrangian nature of the scheme it shall not be used in combination with Adams-Bashforth. We shall hide this option in the namelist.oce for now but it hopefully will be supported in the next release. --- src/oce_adv_tra_ver.F90 | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/src/oce_adv_tra_ver.F90 b/src/oce_adv_tra_ver.F90 index cd07d947f..06986e724 100644 --- a/src/oce_adv_tra_ver.F90 +++ b/src/oce_adv_tra_ver.F90 @@ -390,18 +390,16 @@ subroutine adv_tra_vert_ppm(ttf, w, do_Xmoment, mesh, flux, init_zero) nzmax=nlevels_nod2D(n) nzmin=ulevels_nod2D(n) - ! tracer at surface layer - tv(nzmin)=ttf(nzmin,n) - - ! tracer at surface+1 layer - ! tv(2)=-ttf(1,n)*min(sign(1.0, W(2,n)), 0._WP)+ttf(2,n)*max(sign(1.0, W(2,n)), 0._WP) - tv(nzmin+1)=0.5*(ttf(nzmin,n)+ttf(nzmin+1,n)) - - ! tacer at bottom-1 layer - !tv(nzmax-1)=-ttf(nzmax-2,n)*min(sign(1.0, W(nzmax-1,n)), 0._WP)+ttf(nzmax-1,n)*max(sign(1.0, W(nzmax-1,n)), 0._WP) - tv(nzmax-1)=0.5_WP*(ttf(nzmax-2,n)+ttf(nzmax-1,n)) - - ! tracer at bottom layer + ! tracer at surface level + tv(nzmin)=ttf(nzmin,n) + ! tracer at surface+1 level +! tv(2)=-ttf(1,n)*min(sign(1.0, W(2,n)), 0._WP)+ttf(2,n)*max(sign(1.0, W(2,n)), 0._WP) +! tv(3)=-ttf(2,n)*min(sign(1.0, W(3,n)), 0._WP)+ttf(3,n)*max(sign(1.0, W(3,n)), 0._WP) + tv(nzmin+1)=0.5*(ttf(nzmin, n)+ttf(nzmin+1,n)) + ! tacer at bottom-1 level + tv(nzmax-1)=-ttf(nzmax-2,n)*min(sign(1.0, W(nzmax-1,n)), 0._WP)+ttf(nzmax-1,n)*max(sign(1.0, W(nzmax-1,n)), 0._WP) +! tv(nzmax-1)=0.5_WP*(ttf(nzmax-2,n)+ttf(nzmax-1,n)) + ! tracer at bottom level tv(nzmax)=ttf(nzmax-1,n) !_______________________________________________________________________ @@ -409,7 +407,7 @@ subroutine adv_tra_vert_ppm(ttf, w, do_Xmoment, mesh, flux, init_zero) ! see Colella and Woodward, JCP, 1984, 174-201 --> equation (1.9) ! loop over layers (segments) !!PS do nz=3, nzmax-3 - do nz=nzmin+2, nzmax-2 + do nz=nzmin+1, nzmax-3 !___________________________________________________________________ ! for uniform spaced vertical grids --> piecewise parabolic method (ppm) ! equation (1.9) @@ -419,10 +417,10 @@ subroutine adv_tra_vert_ppm(ttf, w, do_Xmoment, mesh, flux, init_zero) ! for non-uniformity spaced vertical grids --> piecewise parabolic ! method (ppm) see see Colella and Woodward, JCP, 1984, 174-201 ! --> full equation (1.6), (1.7) and (1.8) - dzjm1 = hnode_new(nz-2,n) - dzj = hnode_new(nz-1,n) - dzjp1 = hnode_new(nz,n) - dzjp2 = hnode_new(nz+1,n) + dzjm1 = hnode_new(nz-1,n) + dzj = hnode_new(nz ,n) + dzjp1 = hnode_new(nz+1,n) + dzjp2 = hnode_new(nz+2,n) ! Be carefull here vertical operation have to be done on NEW vertical mesh !!! !___________________________________________________________________