diff --git a/config/namelist.dyn b/config/namelist.dyn index 43ab9d0ba..545a7fe67 100644 --- a/config/namelist.dyn +++ b/config/namelist.dyn @@ -22,5 +22,15 @@ wsplit_maxcfl= 1.0 ! maximum allowed CFL criteria in vertical (0.5 < w_max_c ! in older FESOM it used to be w_exp_max=1.e-3 ldiag_KE=.false. ! activates energy diagnostics AB_order=2 + +use_ssh_se_subcycl = .false. +se_BTsteps = 50 +se_BTtheta = 0.14 ! default: 0.14, +se_bottdrag = .true. +se_bdrag_si = .true. ! bottomdrag semi-implicite/explicite +se_visc = .true. +se_visc_gamma0 = 10 +se_visc_gamma1 = 19500 !19500 (core2@32spd), 2750 (core2@72spd) +se_visc_gamma2 = 0 / diff --git a/src/MOD_DYN.F90 b/src/MOD_DYN.F90 index b77437bbf..6f2a0c329 100644 --- a/src/MOD_DYN.F90 +++ b/src/MOD_DYN.F90 @@ -42,7 +42,7 @@ MODULE MOD_DYN ! set main structure for dynamicss, contains viscosity options and parameters + ! option for momentum advection TYPE T_DYN -!___________________________________________________________________________ + !___________________________________________________________________________ ! instant zonal merdional velocity & Adams-Bashfort rhs real(kind=WP), allocatable, dimension(:,:,:) :: uv, uv_rhs, fer_uv real(kind=WP), allocatable, dimension(:,:,:,:) :: uv_rhsAB @@ -55,6 +55,24 @@ MODULE MOD_DYN ! sea surface height arrays real(kind=WP), allocatable, dimension(:) :: eta_n, d_eta, ssh_rhs, ssh_rhs_old + + !___arrays for split explicite ssh computation______________________________ + ! se_uvh...transport velocity, + real(kind=WP), allocatable, dimension(:,:,:):: se_uvh + !se_uv_rhs...vertical integral of transport velocity rhs, se_uvBT_4AB... + ! barotropic transport velocities (vertically integrated), contains actual + ! timestep (1:2) and previous timestep (3:4) for adams-bashfort interpolation + real(kind=WP), allocatable, dimension(:,:) :: se_uvBT_rhs, se_uvBT_4AB + + ! se_uvBT...barotropic trnasport velocities from barotropic time stepping + ! se_uvBT_theta...velocities for dissipative time stepping of thickness equation + ! UBTmean_mean... Mean BT velocity to trim 3D velocity in tracers + real(kind=WP), allocatable, dimension(:,:) :: se_uvBT, se_uvBT_theta, se_uvBT_mean, se_uvBT_12 + + ! array that are needed for viscosity and bottomdrag stabilization of + ! split-expl subcycling method + real(kind=WP), allocatable, dimension(:,:) :: se_uvBT_stab_hvisc + real(kind=WP), allocatable, dimension(:) :: se_uvBT_stab_bdrag ! LA: 2023-05-17 iceberg arrays real(kind=WP), allocatable, dimension(:) :: eta_n_ib ! kh 18.03.21 additional array for asynchronous iceberg computations @@ -99,7 +117,24 @@ MODULE MOD_DYN logical :: use_wsplit = .false. ! maximum allowed CFL criteria in vertical (0.5 < w_max_cfl < 1.) ! in older FESOM it used to be w_exp_max=1.e-3 - real(kind=WP) :: wsplit_maxcfl= 1.0 + real(kind=WP) :: wsplit_maxcfl = 1.0 + + ! switch between ssh computation, by solver or split explicite subcycling + ! use_ssh_se_subcycl = .false. --> solver + ! use_ssh_se_subcycl = .true. --> split explicite subcycling + logical :: use_ssh_se_subcycl = .false. + + ! barotropic subcycling time-steps and dissipation parameter + integer :: se_BTsteps = 50 + real(kind=WP) :: se_BTtheta = 0.14_WP + logical :: se_bottdrag = .true. + logical :: se_bdrag_si = .true. + logical :: se_visc = .true. + real(kind=WP) :: se_visc_gamma0 = 10 + real(kind=WP) :: se_visc_gamma1 = 2750 + real(kind=WP) :: se_visc_gamma2 = 0 + + !___________________________________________________________________________ ! energy diagnostic part: will be computed inside the model ("hard integration"): logical :: ldiag_ke = .true. ! different contributions to velocity change. will be computed inside the code. @@ -229,7 +264,8 @@ subroutine WRITE_T_DYN(dynamics, unit, iostat, iomsg) write(unit, iostat=iostat, iomsg=iomsg) dynamics%use_freeslip write(unit, iostat=iostat, iomsg=iomsg) dynamics%use_wsplit write(unit, iostat=iostat, iomsg=iomsg) dynamics%wsplit_maxcfl - + write(unit, iostat=iostat, iomsg=iomsg) dynamics%use_ssh_se_subcycl + !___________________________________________________________________________ call dynamics%solverinfo%WRITE_T_SOLVERINFO(unit) @@ -249,6 +285,17 @@ subroutine WRITE_T_DYN(dynamics, unit, iostat, iomsg) call write_bin_array(dynamics%fer_w , unit, iostat, iomsg) call write_bin_array(dynamics%fer_uv, unit, iostat, iomsg) end if + if (dynamics%use_ssh_se_subcycl) then + call write_bin_array(dynamics%se_uvh , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_rhs , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_4AB , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_theta, unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_mean , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_12 , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_stab_hvisc , unit, iostat, iomsg) + call write_bin_array(dynamics%se_uvBT_stab_bdrag , unit, iostat, iomsg) + end if end subroutine WRITE_T_DYN @@ -275,6 +322,7 @@ subroutine READ_T_DYN(dynamics, unit, iostat, iomsg) read(unit, iostat=iostat, iomsg=iomsg) dynamics%use_freeslip read(unit, iostat=iostat, iomsg=iomsg) dynamics%use_wsplit read(unit, iostat=iostat, iomsg=iomsg) dynamics%wsplit_maxcfl + read(unit, iostat=iostat, iomsg=iomsg) dynamics%use_ssh_se_subcycl !___________________________________________________________________________ call dynamics%solverinfo%READ_T_SOLVERINFO(unit) @@ -295,6 +343,17 @@ subroutine READ_T_DYN(dynamics, unit, iostat, iomsg) call read_bin_array(dynamics%fer_w , unit, iostat, iomsg) call read_bin_array(dynamics%fer_uv , unit, iostat, iomsg) end if + if (dynamics%use_ssh_se_subcycl) then + call read_bin_array(dynamics%se_uvh , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_rhs , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_4AB , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_theta, unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_mean , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_12 , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_stab_hvisc , unit, iostat, iomsg) + call read_bin_array(dynamics%se_uvBT_stab_bdrag , unit, iostat, iomsg) + end if end subroutine READ_T_DYN diff --git a/src/associate_mesh_ass.h b/src/associate_mesh_ass.h index d2fa85123..d48a84500 100644 --- a/src/associate_mesh_ass.h +++ b/src/associate_mesh_ass.h @@ -57,14 +57,14 @@ Z_3d_n(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D) => mesh%Z_3d_n(:,:) #if defined(__async_icebergs) Z_3d_n_ib(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D) => mesh%Z_3d_n_ib(:,:) #endif -helem(1:mesh%nl-1, 1:myDim_elem2D) => mesh%helem(:,:) -bottom_elem_thickness(1:myDim_elem2D) => mesh%bottom_elem_thickness(:) +helem(1:mesh%nl-1, 1:myDim_elem2D+eDim_elem2D) => mesh%helem(:,:) +bottom_elem_thickness(1:myDim_elem2D+eDim_elem2D) => mesh%bottom_elem_thickness(:) bottom_node_thickness(1:myDim_nod2D+eDim_nod2D) => mesh%bottom_node_thickness(:) dhe(1:myDim_elem2D) => mesh%dhe(:) hbar(1:myDim_nod2D+eDim_nod2D) => mesh%hbar(:) hbar_old(1:myDim_nod2D+eDim_nod2D) => mesh%hbar_old(:) -!zbar_n(1:mesh%nl) => mesh%zbar_n -!Z_n(1:mesh%nl-1) => mesh%Z_n +!zbar_n(1:mesh%nl) => mesh%zbar_n +!Z_n(1:mesh%nl-1) => mesh%Z_n zbar_n_bot(1:myDim_nod2D+eDim_nod2D) => mesh%zbar_n_bot(:) zbar_e_bot(1:myDim_elem2D+eDim_elem2D) => mesh%zbar_e_bot(:) zbar_n_srf(1:myDim_nod2D+eDim_nod2D) => mesh%zbar_n_srf(:) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index b9c6f8745..8ecdda816 100644 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -80,9 +80,10 @@ module diagnostics logical :: ldiag_extflds =.false. logical :: ldiag_ice =.false. logical :: ldiag_trflx =.false. - - namelist /diag_list/ ldiag_solver, lcurt_stress_surf, ldiag_curl_vel3, ldiag_Ri, ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, & - ldiag_salt3D, ldiag_forc, ldiag_vorticity, ldiag_extflds, ldiag_trflx + + namelist /diag_list/ ldiag_solver, lcurt_stress_surf, ldiag_curl_vel3, ldiag_Ri, & + ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, ldiag_salt3D, ldiag_forc, & + ldiag_vorticity, ldiag_extflds, ldiag_trflx, ldiag_ice contains @@ -181,69 +182,78 @@ subroutine diag_curl_vel3(mode, dynamics, partit, mesh) #include "associate_mesh_ass.h" UV => dynamics%uv(:,:,:) - !===================== + !___________________________________________________________________________ if (firstcall) then !allocate the stuff at the first call allocate(curl_vel3(nl-1, myDim_nod2D+eDim_nod2D)) firstcall=.false. if (mode==0) return end if + !___________________________________________________________________________ curl_vel3=0. - - DO ed=1,myDim_edge2D + do ed=1,myDim_edge2D enodes=edges(:,ed) el=edge_tri(:,ed) + !_______________________________________________________________________ nl1=nlevels(el(1))-1 nu1=ulevels(el(1)) deltaX1=edge_cross_dxdy(1,ed) deltaY1=edge_cross_dxdy(2,ed) - nl2=0 - nu2=0 + + !_______________________________________________________________________ if (el(2)>0) then deltaX2=edge_cross_dxdy(3,ed) deltaY2=edge_cross_dxdy(4,ed) nl2=nlevels(el(2))-1 nu2=ulevels(el(2)) - end if - - nl12 = min(nl1,nl2) - nu12 = min(nu1,nu2) - DO nz=nu1,nu12-1 - c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1)) - curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 - curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 - END DO - if (nu2>0) then - DO nz=nu2,nu12-1 + nl12 = min(nl1,nl2) + nu12 = max(nu1,nu2) + !___________________________________________________________________ + do nz=nu1,nu12-1 + c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1)) + curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 + curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 + end do + do nz=nu2,nu12-1 c1= -deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2)) curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 - END DO + end do + !___________________________________________________________________ + do nz=nu12,nl12 + c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))- & + deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2)) + curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 + curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 + end do + !___________________________________________________________________ + do nz=nl12+1,nl1 + c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1)) + curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 + curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 + end do + do nz=nl12+1,nl2 + c1= -deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2)) + curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 + curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 + end do + !_______________________________________________________________________ + else + do nz=nu1,nl1 + c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1)) + curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 + curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 + end do end if - DO nz=nu12,nl12 - c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))- & - deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2)) - curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 - curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 - END DO - DO nz=nl12+1,nl1 - c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1)) - curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 - curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 - END DO - DO nz=nl12+1,nl2 - c1= -deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2)) - curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1 - curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1 - END DO - END DO - - DO n=1, myDim_nod2D - !!PS DO nz=1, nlevels_nod2D(n)-1 - DO nz=ulevels_nod2D(n), nlevels_nod2D(n)-1 + end do + + !___________________________________________________________________________ + do n=1, myDim_nod2D + do nz=ulevels_nod2D(n), nlevels_nod2D(n)-1 curl_vel3(nz,n)=curl_vel3(nz,n)/areasvol(nz,n) - END DO - END DO + end do + end do + end subroutine diag_curl_vel3 ! ============================================================== ! diff --git a/src/io_blowup.F90 b/src/io_blowup.F90 index 87e570c3b..9c2a1f8cc 100644 --- a/src/io_blowup.F90 +++ b/src/io_blowup.F90 @@ -102,12 +102,27 @@ subroutine ini_blowup_io(year, ice, dynamics, tracers, partit, mesh) !___Define the netCDF variables for 2D fields_______________________________ !___SSH_____________________________________________________________________ call def_variable(bid, 'eta_n' , (/nod2D/) , 'sea surface elevation', 'm', dynamics%eta_n); - call def_variable(bid, 'd_eta' , (/nod2D/) , 'change in ssh from solver', 'm', dynamics%d_eta); !___ALE related fields______________________________________________________ call def_variable(bid, 'hbar' , (/nod2D/) , 'ALE surface elevation hbar_n+0.5', 'm', hbar); !!PS call def_variable(bid, 'hbar_old' , (/nod2D/) , 'ALE surface elevation hbar_n-0.5', 'm', hbar_old); - call def_variable(bid, 'ssh_rhs' , (/nod2D/) , 'RHS for the elevation', '?', dynamics%ssh_rhs); - call def_variable(bid, 'ssh_rhs_old', (/nod2D/) , 'RHS for the elevation', '?', dynamics%ssh_rhs_old); + if (.not. dynamics%use_ssh_se_subcycl) then + call def_variable(bid, 'd_eta' , (/nod2D/) , 'change in ssh from solver', 'm', dynamics%d_eta); + call def_variable(bid, 'ssh_rhs' , (/nod2D/) , 'RHS for the elevation', '?', dynamics%ssh_rhs); + call def_variable(bid, 'ssh_rhs_old', (/nod2D/) , 'RHS for the elevation', '?', dynamics%ssh_rhs_old); + + else + call def_variable(bid, 'ubt_rhs' , (/elem2D/), 'zonal RHS barotr. transp. equation' , '?' , dynamics%se_uvBT_rhs( 1,:)); + call def_variable(bid, 'vbt_rhs' , (/elem2D/), 'merid. RHS barotr. transp. equation', '?' , dynamics%se_uvBT_rhs( 2,:)); + call def_variable(bid, 'ubt' , (/elem2D/), 'zonal barotr. transp.' , '?' , dynamics%se_uvBT( 1,:)); + call def_variable(bid, 'vbt' , (/elem2D/), 'merid. barotr. transp.' , '?' , dynamics%se_uvBT( 2,:)); + call def_variable(bid, 'ubt_theta', (/elem2D/), 'zonal barotr. theta term.' , '?' , dynamics%se_uvBT_theta(1,:)); + call def_variable(bid, 'vbt_theta', (/elem2D/), 'merid. barotr. theta term' , '?' , dynamics%se_uvBT_theta(2,:)); + call def_variable(bid, 'ubt_mean' , (/elem2D/), 'zonal barotr. mean term.' , '?' , dynamics%se_uvBT_mean( 1,:)); + call def_variable(bid, 'vbt_mean' , (/elem2D/), 'merid. barotr. mean term' , '?' , dynamics%se_uvBT_mean( 2,:)); + call def_variable(bid, 'uh' , (/nl-1, elem2D/), 'zonal velocity' , 'm/s', dynamics%se_uvh(1,:,:)); + call def_variable(bid, 'vh' , (/nl-1, elem2D/), 'meridional velocity' , 'm/s', dynamics%se_uvh(2,:,:)); + end if + !___Define the netCDF variables for 3D fields_______________________________ call def_variable(bid, 'hnode' , (/nl-1, nod2D/) , 'ALE stuff', '?', hnode); call def_variable(bid, 'helem' , (/nl-1, elem2D/) , 'Element layer thickness', 'm/s', helem(:,:)); diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 4232da250..ff789ea9d 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -564,6 +564,8 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) if (Fer_GM) then call def_stream( nod2D , myDim_nod2D , 'fer_C', 'GM, depth independent speed', 'm/s' , fer_c(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) end if +CASE ('cfl_z ') + call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'cfl_z', 'vertical CFL criteria', '?', dynamics%cfl_z(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) !_______________________________________________________________________________ ! Density MOC diagnostics @@ -702,7 +704,32 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) if (sel_forcvar(14)==0) call def_stream(nod2D , myDim_nod2D , 'relaxsalt', 'relaxation salt flux' , 'm/s*psu', relax_salt(:) , io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) #endif if (sel_forcvar(15)==0) call def_stream(nod2D , myDim_nod2D , 'realsalt' , 'real salt flux from sea ice' , 'm/s*psu', real_salt_flux(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) - end if + end if +!_______________________________________________________________________________ +! Split-Explicite subcycling varaibles +CASE ('SPLIT-EXPL') + if (dynamics%use_ssh_se_subcycl) then + call def_stream(elem2D, myDim_elem2D , 'ubt' , 'zonal barotrop. transport' , '?' , dynamics%se_uvBT(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'vbt' , 'merid. barotrop. transport', '?' , dynamics%se_uvBT(2,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'ubt_rhs', 'zonal barotrop. transport RHS' , '?' , dynamics%se_uvBT_rhs(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'vbt_rhs', 'merid. barotrop. transport RHS', '?' , dynamics%se_uvBT_rhs(2,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'ubt_12' , 'zonal barotrop. transport 12' , '?' , dynamics%se_uvBT_12(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'vbt_12' , 'merid. barotrop. transport 12', '?' , dynamics%se_uvBT_12(2,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'ubt_theta', 'zonal barotrop. transport theta' , '?' , dynamics%se_uvBT_theta(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'vbt_theta', 'merid. barotrop. transport theta', '?' , dynamics%se_uvBT_theta(2,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'ubt_mean' , 'zonal barotrop. transport mean' , '?' , dynamics%se_uvBT_mean(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'vbt_mean' , 'merid. barotrop. transport mean', '?' , dynamics%se_uvBT_mean(2,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream((/nl-1,elem2D/), (/nl-1,myDim_elem2D/) , 'u_rhs' , 'zonal transport rhs' , '?' , dynamics%uv_rhs(1,:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream((/nl-1,elem2D/), (/nl-1,myDim_elem2D/) , 'v_rhs' , 'merid. transport rhs', '?' , dynamics%uv_rhs(2,:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream((/nl-1,elem2D/), (/nl-1,myDim_elem2D/) , 'uh' , 'zonal transport' , '?' , dynamics%se_uvh(1,:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream((/nl-1,elem2D/), (/nl-1,myDim_elem2D/) , 'vh' , 'merid. transport', '?' , dynamics%se_uvh(2,:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + if (dynamics%se_visc) then + call def_stream(elem2D, myDim_elem2D , 'ubt_hvisc' , 'zonal hvisc. stabil.' , '?' , dynamics%se_uvBT_stab_hvisc(1,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(elem2D, myDim_elem2D , 'vbt_hvisc' , 'merid. hvisc. stabil.' , '?' , dynamics%se_uvBT_stab_hvisc(2,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + end if + + !_______________________________________________________________________________ CASE DEFAULT if (mype==0) write(*,*) 'stream ', io_list(i)%id, ' is not defined !' @@ -781,73 +808,74 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) !___________________________________________________________________________ if (dynamics%ldiag_ke) then - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_u_xVEL', 'work of advection [u]', 'm2/s2', dynamics%ke_adv_xVEL(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_v_xVEL', 'work of advection [v]', 'm2/s2', dynamics%ke_adv_xVEL(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_u_xVEL', 'work of advection [u]', 'm2/s2', dynamics%ke_adv_xVEL(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_v_xVEL', 'work of advection [v]', 'm2/s2', dynamics%ke_adv_xVEL(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_u_xVEL', 'work Coriolis :) [u]', 'm2/s2', dynamics%ke_cor_xVEL(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_v_xVEL', 'work Coriolis :) [v]', 'm2/s2', dynamics%ke_cor_xVEL(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_u_xVEL', 'work Coriolis :) [u]', 'm2/s2', dynamics%ke_cor_xVEL(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_v_xVEL', 'work Coriolis :) [v]', 'm2/s2', dynamics%ke_cor_xVEL(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_u_xVEL', 'work of pressure gradient [x]', 'm2/s2', dynamics%ke_pre_xVEL(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_v_xVEL', 'work of pressure gradient [y]', 'm2/s2', dynamics%ke_pre_xVEL(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_u_xVEL', 'work of pressure gradient [x]', 'm2/s2', dynamics%ke_pre_xVEL(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_v_xVEL', 'work of pressure gradient [y]', 'm2/s2', dynamics%ke_pre_xVEL(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_u_xVEL', 'work of hor. visc. [x]', 'm2/s2', dynamics%ke_hvis_xVEL(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_v_xVEL', 'work of hor. visc. [y]', 'm2/s2', dynamics%ke_hvis_xVEL(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_u_xVEL', 'work of hor. visc. [x]', 'm2/s2', dynamics%ke_hvis_xVEL(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_v_xVEL', 'work of hor. visc. [y]', 'm2/s2', dynamics%ke_hvis_xVEL(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_u_xVEL', 'work of ver. visc. [x]', 'm2/s2', dynamics%ke_vvis_xVEL(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_v_xVEL', 'work of ver. visc. [y]', 'm2/s2', dynamics%ke_vvis_xVEL(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_u_xVEL', 'work of ver. visc. [x]', 'm2/s2', dynamics%ke_vvis_xVEL(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_v_xVEL', 'work of ver. visc. [y]', 'm2/s2', dynamics%ke_vvis_xVEL(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_du2', 'KE change [0.5 du2]', 'm2/s2', dynamics%ke_du2(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_dv2', 'KE change [0.5 dv2]', 'm2/s2', dynamics%ke_du2(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_du2', 'KE change [0.5 du2]', 'm2/s2', dynamics%ke_du2(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_dv2', 'KE change [0.5 dv2]', 'm2/s2', dynamics%ke_du2(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'W_x_RHO', 'buoyancy work = ke_pre x VEL', 'm2/s2', dynamics%ke_wrho(:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'W_x_RHO', 'buoyancy work = ke_pre x VEL', 'm2/s2', dynamics%ke_wrho(:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_wind_x_xVEL', 'work of wind [x]', 'm2/s2', dynamics%ke_wind_xVEL(1,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_wind_y_xVEL', 'work of wind [y]', 'm2/s2', dynamics%ke_wind_xVEL(2,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_wind_x_xVEL', 'work of wind [x]', 'm2/s2', dynamics%ke_wind_xVEL(1,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_wind_y_xVEL', 'work of wind [y]', 'm2/s2', dynamics%ke_wind_xVEL(2,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_drag_x_xVEL', 'work of drag [x]', 'm2/s2', dynamics%ke_drag_xVEL(1,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_drag_y_xVEL', 'work of drag [y]', 'm2/s2', dynamics%ke_drag_xVEL(2,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_drag_x_xVEL', 'work of drag [x]', 'm2/s2', dynamics%ke_drag_xVEL(1,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_drag_y_xVEL', 'work of drag [y]', 'm2/s2', dynamics%ke_drag_xVEL(2,:), io_list(i)%freq, 'y', 8, partit, mesh) ! same as above but without multiplying with UMEAN (for later computation of turbulence fluxes) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_u', 'advection *dT [u]', 'm/s', dynamics%ke_adv(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_v', 'advection *dT [v]', 'm/s', dynamics%ke_adv(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_u', 'advection *dT [u]', 'm/s', dynamics%ke_adv(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_adv_v', 'advection *dT [v]', 'm/s', dynamics%ke_adv(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_u', 'Coriolis *dT [X]', 'm/s', dynamics%ke_cor(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_v', 'Coriolis *dT [Y]', 'm/s', dynamics%ke_cor(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_u', 'Coriolis *dT [X]', 'm/s', dynamics%ke_cor(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_cor_v', 'Coriolis *dT [Y]', 'm/s', dynamics%ke_cor(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_u', 'pressure gradient *dT [x]', 'm/s', dynamics%ke_pre(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_v', 'pressure gradient *dT [y]', 'm/s', dynamics%ke_pre(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_u', 'pressure gradient *dT [x]', 'm/s', dynamics%ke_pre(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_pre_v', 'pressure gradient *dT [y]', 'm/s', dynamics%ke_pre(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_u', 'hor. visc. [x] *dT', 'm/s', dynamics%ke_hvis(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_v', 'hor. visc. [y] *dT', 'm/s', dynamics%ke_hvis(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_u', 'hor. visc. [x] *dT', 'm/s', dynamics%ke_hvis(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_hvis_v', 'hor. visc. [y] *dT', 'm/s', dynamics%ke_hvis(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_u', 'ver. visc. [x] *dT', 'm/s', dynamics%ke_vvis(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_v', 'ver. visc. [y] *dT', 'm/s', dynamics%ke_vvis(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_u', 'ver. visc. [x] *dT', 'm/s', dynamics%ke_vvis(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_vvis_v', 'ver. visc. [y] *dT', 'm/s', dynamics%ke_vvis(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_Umean', 'mean U', 'm/s', dynamics%ke_umean(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_Vmean', 'mean V', 'm/s', dynamics%ke_umean(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_Umean', 'mean U', 'm/s', dynamics%ke_umean(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_Vmean', 'mean V', 'm/s', dynamics%ke_umean(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_U2mean', 'U2 mean', 'm/s', dynamics%ke_u2mean(1,:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_V2mean', 'V2 mean', 'm/s', dynamics%ke_u2mean(2,:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_U2mean', 'U2 mean', 'm/s', dynamics%ke_u2mean(1,:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'ke_V2mean', 'V2 mean', 'm/s', dynamics%ke_u2mean(2,:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'ke_dW', 'd/dz (vertical VEL)', 'm/s', dynamics%ke_dW(:,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'ke_PFULL', 'full Pressure', 'm/s', dynamics%ke_Pfull(:,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'ke_dW', 'd/dz (vertical VEL)', 'm/s', dynamics%ke_dW(:,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'ke_PFULL', 'full Pressure', 'm/s', dynamics%ke_Pfull(:,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_wind_x', 'wind [x] *dT', 'm/s', dynamics%ke_wind(1,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_wind_y', 'wind [y] *dT', 'm/s', dynamics%ke_wind(2,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_wind_x', 'wind [x] *dT', 'm/s', dynamics%ke_wind(1,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_wind_y', 'wind [y] *dT', 'm/s', dynamics%ke_wind(2,:), io_list(i)%freq, 'y', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_drag_x', 'drag [x] *dT', 'm/s', dynamics%ke_drag(1,:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(elem2D, myDim_elem2D, 'ke_drag_y', 'drag [y] *dT', 'm/s', dynamics%ke_drag(2,:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_drag_x', 'drag [x] *dT', 'm/s', dynamics%ke_drag(1,:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(elem2D, myDim_elem2D, 'ke_drag_y', 'drag [y] *dT', 'm/s', dynamics%ke_drag(2,:), io_list(i)%freq, 'y', 8, partit, mesh) ! surface fields for APE input computations... - call def_stream(nod2D , myDim_nod2D , 'ke_J', 'surface temperature flux [Js]','°C/s', dynamics%ke_J(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_G', 'surface salinity flux [Gs]','PSU/s', dynamics%ke_G(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_D', 'surface density', 'kg/m^3', dynamics%ke_D(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_D2', 'surface density squared', 'kg^2/m^6', dynamics%ke_D2(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_n0', 'dRHO/dz', 'kg/m^4', dynamics%ke_n0(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_JD', 'surface temperature flux [Js] * RHO','°C*kg/s/m^3', dynamics%ke_JD(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_GD', 'surface salinity flux [Gs] * RHO','PSU*kg/s/m^3', dynamics%ke_GD(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_swA', 'Thermal expansion coeff (alpha)', '1/°C', dynamics%ke_swA(:), io_list(i)%freq, 'm', 8, partit, mesh) - call def_stream(nod2D , myDim_nod2D , 'ke_swB', 'Taline contraction coeff (beta)', '1/PSU', dynamics%ke_swB(:), io_list(i)%freq, 'm', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_J', 'surface temperature flux [Js]','°C/s', dynamics%ke_J(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_G', 'surface salinity flux [Gs]','PSU/s', dynamics%ke_G(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_D', 'surface density', 'kg/m^3', dynamics%ke_D(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_D2', 'surface density squared', 'kg^2/m^6', dynamics%ke_D2(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_n0', 'dRHO/dz', 'kg/m^4', dynamics%ke_n0(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_JD', 'surface temperature flux [Js] * RHO','°C*kg/s/m^3', dynamics%ke_JD(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_GD', 'surface salinity flux [Gs] * RHO','PSU*kg/s/m^3', dynamics%ke_GD(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_swA', 'Thermal expansion coeff (alpha)', '1/°C', dynamics%ke_swA(:), io_list(i)%freq, 'y', 8, partit, mesh) + call def_stream(nod2D , myDim_nod2D , 'ke_swB', 'Taline contraction coeff (beta)', '1/PSU', dynamics%ke_swB(:), io_list(i)%freq, 'y', 8, partit, mesh) end if + end subroutine ! ! diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index dbdb8259e..734e7f478 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1,3 +1,46 @@ + +module compute_CFLz_interface + interface + subroutine compute_CFLz(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + end interface +end module + +module compute_Wvel_split_interface + interface + subroutine compute_Wvel_split(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + end interface +end module + +module compute_vert_vel_transpv_interface + interface + subroutine compute_vert_vel_transpv(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + end interface +end module + module oce_ale_interfaces interface subroutine init_bottom_elem_thickness(partit, mesh) @@ -146,6 +189,8 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) end subroutine end interface end module + + ! CONTENT: ! ------------ ! subroutine ale_init @@ -212,7 +257,7 @@ subroutine init_ale(dynamics, partit, mesh) allocate(mesh%hbar_old(myDim_nod2D+eDim_nod2D)) ! helem: layer thickness at elements. It is interpolated from hnode. - allocate(mesh%helem(1:nl-1, myDim_elem2D)) + allocate(mesh%helem(1:nl-1, myDim_elem2D+eDim_nod2D)) ! dhe: The increment of total fluid depth on elements. It is used to update the matrix ! of the ssh operator. @@ -255,7 +300,7 @@ subroutine init_ale(dynamics, partit, mesh) end if ! bottom_elem_tickness: changed bottom layer thinkness due to partial cells - allocate(mesh%bottom_elem_thickness(myDim_elem2D)) + allocate(mesh%bottom_elem_thickness(myDim_elem2D+eDim_nod2D)) allocate(mesh%zbar_e_bot(myDim_elem2D+eDim_elem2D)) allocate(mesh%zbar_e_srf(myDim_elem2D+eDim_elem2D)) @@ -270,9 +315,9 @@ subroutine init_ale(dynamics, partit, mesh) hnode(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D) => mesh%hnode(:,:) hnode_new(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D) => mesh%hnode_new(:,:) zbar_3d_n(1:mesh%nl, 1:myDim_nod2D+eDim_nod2D) => mesh%zbar_3d_n(:,:) - !Z_3d_n(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D) => mesh%Z_3d_n(:,:) - helem(1:mesh%nl-1, 1:myDim_elem2D) => mesh%helem(:,:) - bottom_elem_thickness(1:myDim_elem2D) => mesh%bottom_elem_thickness(:) + Z_3d_n(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D) => mesh%Z_3d_n(:,:) + helem(1:mesh%nl-1, 1:myDim_elem2D+eDim_nod2D) => mesh%helem(:,:) + bottom_elem_thickness(1:myDim_elem2D+eDim_nod2D) => mesh%bottom_elem_thickness(:) bottom_node_thickness(1:myDim_nod2D+eDim_nod2D) => mesh%bottom_node_thickness(:) dhe(1:myDim_elem2D) => mesh%dhe(:) hbar(1:myDim_nod2D+eDim_nod2D) => mesh%hbar(:) @@ -470,6 +515,7 @@ subroutine init_bottom_elem_thickness(partit, mesh) !___________________________________________________________________________ call exchange_elem(zbar_e_bot, partit) + call exchange_elem(bottom_elem_thickness, partit) end subroutine init_bottom_elem_thickness ! @@ -746,6 +792,7 @@ subroutine init_thickness_ale(dynamics, partit, mesh) USE MOD_PARTIT USE MOD_PARSUP USE MOD_DYN + use g_comm_auto implicit none type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit @@ -989,7 +1036,7 @@ subroutine init_thickness_ale(dynamics, partit, mesh) !___________________________________________________________________________ hnode_new=hnode ! Should be initialized, because only variable part is updated. - + call exchange_elem(helem, partit) !!PS call check_total_volume(partit, mesh) end subroutine init_thickness_ale @@ -1004,6 +1051,7 @@ subroutine update_thickness_ale(partit, mesh) USE MOD_PARSUP use o_ARRAYS use g_config,only: which_ale,lzstar_lev,min_hnode + use g_comm_auto implicit none type(t_partit), intent(inout), target :: partit type(t_mesh) , intent(inout), target :: mesh @@ -1159,6 +1207,13 @@ subroutine update_thickness_ale(partit, mesh) end do !$OMP END PARALLEL DO endif + + !___________________________________________________________________________ +!$OMP MASTER + call exchange_elem(helem, partit) +!$OMP END MASTER +!$OMP BARRIER + end subroutine update_thickness_ale ! ! @@ -2021,6 +2076,8 @@ subroutine vert_vel_ale(dynamics, partit, mesh) use g_comm_auto use io_RESTART !!PS use g_forcing_arrays !!PS + use compute_Wvel_split_interface + use compute_CFLz_interface implicit none type(t_dyn) , intent(inout), target :: dynamics type(t_partit), intent(inout), target :: partit @@ -2049,14 +2106,20 @@ subroutine vert_vel_ale(dynamics, partit, mesh) Wvel_e => dynamics%w_e(:,:) Wvel_i => dynamics%w_i(:,:) CFL_z => dynamics%cfl_z(:,:) - ssh_rhs => dynamics%ssh_rhs(:) - ssh_rhs_old => dynamics%ssh_rhs_old(:) eta_n => dynamics%eta_n(:) - d_eta => dynamics%d_eta(:) if (Fer_GM) then fer_UV => dynamics%fer_uv(:,:,:) fer_Wvel=> dynamics%fer_w(:,:) - end if + end if + + if ( .not. dynamics%use_ssh_se_subcycl) then + d_eta => dynamics%d_eta(:) + ssh_rhs => dynamics%ssh_rhs(:) + ssh_rhs_old => dynamics%ssh_rhs_old(:) + else + + end if + !___________________________________________________________________________ ! Contributions from levels in divergence !$OMP PARALLEL DO @@ -2124,6 +2187,8 @@ subroutine vert_vel_ale(dynamics, partit, mesh) !_______________________________________________________________________ ! if ed is not a boundary edge --> calc div(u_vec*h) for every layer ! for el(2) + c1 = 0.0_WP + c2 = 0.0_WP if(el(2)>0)then deltaX2=edge_cross_dxdy(3,ed) deltaY2=edge_cross_dxdy(4,ed) @@ -2157,7 +2222,7 @@ subroutine vert_vel_ale(dynamics, partit, mesh) #else !$OMP END ORDERED #endif - end if + end if !~-> if(el(2)>0)then end do ! --> do ed=1, myDim_edge2D !$OMP END DO !$OMP END PARALLEL @@ -2180,6 +2245,14 @@ subroutine vert_vel_ale(dynamics, partit, mesh) fer_Wvel(nz,n)=fer_Wvel(nz,n)+fer_Wvel(nz+1,n) end if end do + + !_______________________________________________________________________ + if ( any(Wvel(nzmin:nzmax, n)/=Wvel(nzmin:nzmax, n)) ) then + write(*,*) ' --> subroutine vert_vel_ale --> found Nan in Wvel after cumulativ summation(...)' + write(*,*) ' mype =', mype + write(*,*) ' node =', n + write(*,*) ' Wvel(nzmin:nzmax, n)=', Wvel(nzmin:nzmax, n) + end if end do !$OMP END PARALLEL DO !___________________________________________________________________________ @@ -2197,6 +2270,14 @@ subroutine vert_vel_ale(dynamics, partit, mesh) end if end do + + !_______________________________________________________________________ + if ( any(Wvel(nzmin:nzmax, n)/=Wvel(nzmin:nzmax, n)) ) then + write(*,*) ' --> subroutine vert_vel_ale --> found Nan in Wvel after divide with area' + write(*,*) ' mype =', mype + write(*,*) ' node =', n + write(*,*) ' Wvel(nzmin:nzmax, n)=', Wvel(nzmin:nzmax, n) + end if end do !$OMP END PARALLEL DO ! | @@ -2472,6 +2553,14 @@ subroutine vert_vel_ale(dynamics, partit, mesh) ! continutity Wvel(nzmin,n)=Wvel(nzmin,n)-water_flux(n) + !_______________________________________________________________________ + if ( any(Wvel(nzmin:nzmax, n)/=Wvel(nzmin:nzmax, n)) ) then + write(*,*) ' --> subroutine vert_vel_ale --> found Nan in Wvel after ssh is distributed' + write(*,*) ' mype =', mype + write(*,*) ' node =', n + write(*,*) ' Wvel(nzmin:nzmax, n)=', Wvel(nzmin:nzmax, n) + end if + endif ! --> if (nzmin==1) then end do ! --> do n=1, myDim_nod2D @@ -2493,10 +2582,12 @@ subroutine vert_vel_ale(dynamics, partit, mesh) write(*,*) write(*,*) 'water_flux = ', water_flux(n) write(*,*) - write(*,*) "ssh_rhs = ", ssh_rhs(n) - write(*,*) "ssh_rhs_old = ", ssh_rhs_old(n) write(*,*) "eta_n = ", eta_n(n) - write(*,*) "d_eta = ", d_eta(n) + if ( .not. dynamics%use_ssh_se_subcycl) then + write(*,*) "d_eta = ", d_eta(n) + write(*,*) "ssh_rhs = ", ssh_rhs(n) + write(*,*) "ssh_rhs_old = ", ssh_rhs_old(n) + end if write(*,*) write(*,*) "zbar_3d_n(1,n) = ", zbar_3d_n(1,n) write(*,*) "dd1 = ", dd1 @@ -2525,65 +2616,375 @@ subroutine vert_vel_ale(dynamics, partit, mesh) call exchange_nod(hnode_new, partit) ! Or extend cycles above if (Fer_GM) call exchange_nod(fer_Wvel, partit) !$OMP BARRIER + + !___________________________________________________________________________ + ! compute vertical CFL_z criteria + call compute_CFLz(dynamics, partit, mesh) + + !___________________________________________________________________________ + ! compute implicite explicite splitting of vetical velocity Wvel according + ! to CFL_z criteria + call compute_Wvel_split(dynamics, partit, mesh) + +end subroutine vert_vel_ale + +! +! +!_______________________________________________________________________________ +! calculate vertical velocity from eq.3 in S. Danilov et al. : FESOM2: from +! finite elements to finite volumes. +! +! dh_k/dt + grad(u*h)_k + (w^t - w^b) + water_flux_k=1 = 0 +! +! w^t = w^b - dh_k/dt - grad(u*h)_k - water_flux=1 +! --> do cumulativ summation from bottom to top +subroutine compute_vert_vel_transpv(dynamics, partit, mesh) + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_MESH + USE MOD_DYN + use o_ARRAYS, only: water_flux + use g_config, only: dt, which_ale + use g_comm_auto + use compute_Wvel_split_interface + use compute_CFLz_interface + implicit none + !___________________________________________________________________________ + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + integer :: node, elem, nz, ed, nzmin, nzmax, ednodes(2), edelem(2) + real(kind=WP) :: hh_inv, deltaX1, deltaX2, deltaY1, deltaY2 + real(kind=WP) :: c1(mesh%nl-1), c2(mesh%nl-1) + + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:,:), pointer :: UVh, fer_UV + real(kind=WP), dimension(:,:) , pointer :: Wvel, Wvel_e, fer_Wvel +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + UVh => dynamics%se_uvh(:,:,:) + Wvel => dynamics%w(:,:) + Wvel_e => dynamics%w_e(:,:) + if (Fer_GM) then + fer_UV => dynamics%fer_uv(:,:,:) + fer_Wvel=> dynamics%fer_w(:,:) + end if + + !___________________________________________________________________________ +!$OMP PARALLEL DO + do node=1, myDim_nod2D+eDim_nod2D + Wvel(:, node)=0.0_WP + if (Fer_GM) then + fer_Wvel(:, node)=0.0_WP + end if + end do ! --> do node=1, myDim_nod2D+eDim_nod2D +!$OMP END PARALLEL DO + + !___________________________________________________________________________ +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(ed, ednodes, edelem, nz, nzmin, nzmax, & +!$OMP deltaX1, deltaY1, deltaX2, deltaY2, c1, c2) +!$OMP DO + do ed=1, myDim_edge2D + ! local indice of nodes that span up edge ed + ednodes=edges(:,ed) + + ! local index of element that contribute to edge + edelem=edge_tri(:,ed) + + ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid edelem(1) to + ! center of edge --> needed to calc flux perpedicular to edge from elem edelem(1) + deltaX1=edge_cross_dxdy(1,ed) + deltaY1=edge_cross_dxdy(2,ed) + + !_______________________________________________________________________ + ! calc div(u_vec*h) for every layer + ! do it with gauss-law: int( div(u_vec)*dV) = int( u_vec * n_vec * dS ) + nzmin = ulevels(edelem(1)) + nzmax = nlevels(edelem(1))-1 + ! we introduced c1 & c2 as arrays here to avoid deadlocks when in OpenMP mode + do nz = nzmax, nzmin, -1 + ! --> h * u_vec * n_vec + ! --> e_vec = (dx,dy), n_vec = (-dy,dx); + ! --> h * u*(-dy) + v*dx + c1(nz)=( UVh(2, nz, edelem(1))*deltaX1 - UVh(1, nz, edelem(1))*deltaY1 ) + ! inflow(outflow) "flux" to control volume of node enodes1 + ! is equal to outflow(inflow) "flux" to control volume of node enodes2 + end do ! --> do nz=nzmax,nzmin,-1 + if (Fer_GM) then + do nz = nzmax, nzmin, -1 + c2(nz)=(fer_UV(2, nz, edelem(1))*deltaX1 - fer_UV(1, nz, edelem(1))*deltaY1)*helem(nz, edelem(1)) + end do + end if +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock (partit%plock(ednodes(1))) +#else +!$OMP ORDERED +#endif + Wvel(nzmin:nzmax, ednodes(1))= Wvel(nzmin:nzmax, ednodes(1))+c1(nzmin:nzmax) + if (Fer_GM) then + fer_Wvel(nzmin:nzmax, ednodes(1))= fer_Wvel(nzmin:nzmax, ednodes(1))+c2(nzmin:nzmax) + end if +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(1))) + call omp_set_lock (partit%plock(ednodes(2))) +#endif + Wvel(nzmin:nzmax, ednodes(2))= Wvel(nzmin:nzmax, ednodes(2))-c1(nzmin:nzmax) + if (Fer_GM) then + fer_Wvel(nzmin:nzmax, ednodes(2))= fer_Wvel(nzmin:nzmax, ednodes(2))-c2(nzmin:nzmax) + end if +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(2))) +#else +!$OMP END ORDERED +#endif + + !_______________________________________________________________________ + ! if ed is not a boundary edge --> calc div(u_vec*h) for every layer + ! for edelem(2) + c1 = 0.0_WP + c2 = 0.0_WP + if(edelem(2)>0)then + deltaX2=edge_cross_dxdy(3,ed) + deltaY2=edge_cross_dxdy(4,ed) + nzmin = ulevels(edelem(2)) + nzmax = nlevels(edelem(2))-1 + do nz = nzmax, nzmin, -1 + c1(nz)=-(UVh(2, nz, edelem(2))*deltaX2 - UVh(1, nz, edelem(2))*deltaY2) + end do ! --> do nz=nzmax,nzmin,-1 + if (Fer_GM) then + do nz = nzmax, nzmin, -1 + c2(nz)=-(fer_UV(2, nz, edelem(2))*deltaX2-fer_UV(1, nz, edelem(2))*deltaY2)*helem(nz, edelem(2)) + end do ! --> do nz=nzmax,nzmin,-1 + end if +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock (partit%plock(ednodes(1))) +#else +!$OMP ORDERED +#endif + Wvel(nzmin:nzmax, ednodes(1))= Wvel(nzmin:nzmax, ednodes(1))+c1(nzmin:nzmax) + if (Fer_GM) then + fer_Wvel(nzmin:nzmax, ednodes(1))= fer_Wvel(nzmin:nzmax, ednodes(1))+c2(nzmin:nzmax) + end if +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(1))) + call omp_set_lock (partit%plock(ednodes(2))) +#endif + Wvel(nzmin:nzmax, ednodes(2))= Wvel(nzmin:nzmax, ednodes(2))-c1(nzmin:nzmax) + if (Fer_GM) then + fer_Wvel(nzmin:nzmax, ednodes(2))= fer_Wvel(nzmin:nzmax, ednodes(2))-c2(nzmin:nzmax) + end if +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(2))) +#else +!$OMP END ORDERED +#endif + end if !--> if(edelem(2)>0)then + + end do ! --> do ed=1, myDim_edge2D +!$OMP END DO +!$OMP END PARALLEL + + !___________________________________________________________________________ + ! add the contribution from -dh/dt * area +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax) + do node=1, myDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2d(node)-1 + do nz=nzmax,nzmin,-1 + Wvel(nz, node)=Wvel(nz, node)-(hnode_new(nz, node)-hnode(nz, node))*area(nz, node)/dt + end do ! --> do nz=nzmax,nzmin,-1 + end do ! --> do node=1, myDim_nod2D +!$OMP END PARALLEL DO + + !___________________________________________________________________________ + ! Sum up to get W*area + ! cumulative summation of div(u_vec*h) vertically + ! W_k = W_k+1 - div(h_k*u_k) + ! W_k ... vertical flux troughdo node=1, myDim_nod2D +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax) + do node=1, myDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2d(node)-1 + do nz=nzmax,nzmin,-1 + Wvel(nz, node)=Wvel(nz, node)+Wvel(nz+1, node) + if (Fer_GM) then + fer_Wvel(nz, node)=fer_Wvel(nz, node)+fer_Wvel(nz+1, node) + end if + end do ! --> do nz=nzmax,nzmin,-1 + end do ! --> do node=1, myDim_nod2D +!$OMP END PARALLEL DO + + !___________________________________________________________________________ + ! divide with depth dependent cell area to convert from Vertical flux to + ! physical vertical velocities in units m/s +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax) + do node=1, myDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2d(node)-1 + do nz=nzmin,nzmax + Wvel(nz, node)=Wvel(nz, node)/area(nz, node) + if (Fer_GM) then + fer_Wvel(nz, node)=fer_Wvel(nz, node)/area(nz, node) + end if + end do ! --> do nz=nzmax,nzmin,-1 + end do ! --> do node=1, myDim_nod2D +!$OMP END PARALLEL DO + + !___________________________________________________________________________ + ! Add surface fresh water flux as upper boundary condition for + ! continutity + if (.not. (trim(which_ale)=='linfs' )) then +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nzmin) + do node=1, myDim_nod2D + nzmin = ulevels_nod2D(node) + if (nzmin==1) Wvel(nzmin, node)=Wvel(nzmin, node)-water_flux(node) + end do ! --> do node=1, myDim_nod2D +!$OMP END PARALLEL DO + end if + + + !___________________________________________________________________________ +!$OMP MASTER + call exchange_nod(Wvel, partit) + if (Fer_GM) call exchange_nod(fer_Wvel, partit) +!$OMP END MASTER +!$OMP BARRIER + + !___________________________________________________________________________ + ! compute vertical CFL_z criteria + call compute_CFLz(dynamics, partit, mesh) + + !___________________________________________________________________________ + ! compute implicite explicite splitting of vetical velocity Wvel according + ! to CFL_z criteria + call compute_Wvel_split(dynamics, partit, mesh) + +end subroutine compute_vert_vel_transpv +! +! +!_______________________________________________________________________________ +! compute vertical CFL_z criteria and print out warning when critical value over +! stepped +subroutine compute_CFLz(dynamics, partit, mesh) + use g_config, only: dt, flag_warn_cflz + use MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + use o_PARAM + use g_comm_auto + implicit none + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh), intent(inout), target :: mesh + !___________________________________________________________________________ + integer :: node, nz, nzmin, nzmax + real(kind=WP) :: cflmax, c1, c2 + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:) , pointer :: Wvel, CFL_z +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + Wvel => dynamics%w(:,:) + CFL_z => dynamics%cfl_z(:,:) + !___________________________________________________________________________ ! calc vertical CFL criteria for debugging purpose and vertical Wvel splitting !$OMP PARALLEL DO - do n=1, myDim_nod2D+eDim_nod2D - CFL_z(1,n)=0._WP + do node=1, myDim_nod2D+eDim_nod2D + CFL_z(1,node)=0._WP end do !$OMP END PARALLEL DO -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n, nz, nzmin, nzmax, c1, c2) - do n=1, myDim_nod2D+eDim_nod2D - nzmin = ulevels_nod2D(n) - nzmax = nlevels_nod2D(n)-1 + !___________________________________________________________________________ +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax, c1, c2) + do node=1, myDim_nod2D+eDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2D(node)-1 do nz=nzmin,nzmax - c1(1)=abs(Wvel(nz,n) *dt/hnode_new(nz,n)) !c1->c1(1) is made for the sake of reproducibility with the master branch (rounding error) - c2(1)=abs(Wvel(nz+1,n)*dt/hnode_new(nz,n)) !otherwise just add these terms (c(1) & c(2)) to CFL_z, respectively! + c1=abs(Wvel(nz,node) *dt/hnode_new(nz,node)) !c1->c1(1) is made for the sake of reproducibility with the master branch (rounding error) + c2=abs(Wvel(nz+1,node)*dt/hnode_new(nz,node)) !otherwise just add these terms (c(1) & c(2)) to CFL_z, respectively! ! strong condition: ! total volume change induced by the vertical motion ! no matter, upwind or downwind ! - CFL_z(nz, n)=CFL_z(nz,n)+c1(1) - CFL_z(nz+1,n)= c2(1) - end do - end do + CFL_z(nz, node)=CFL_z(nz,node)+c1 + CFL_z(nz+1,node)= c2 + end do ! --> do nz=nzmin,nzmax + end do ! --> do node=1, myDim_nod2D+eDim_nod2D !$OMP END PARALLEL DO -cflmax=0. -!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(n) REDUCTION(max:cflmax) + + !___________________________________________________________________________ + cflmax=0. +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(node) REDUCTION(max:cflmax) !$OMP DO - do n=1, myDim_nod2D+eDim_nod2D - cflmax=max(cflmax, maxval(CFL_z(:, n))) - end do + do node=1, myDim_nod2D+eDim_nod2D + cflmax=max(cflmax, maxval(CFL_z(:, node))) + end do ! --> do node=1, myDim_nod2D+eDim_nod2D !$OMP END DO !$OMP END PARALLEL + !___________________________________________________________________________ if (cflmax > 1.0_WP .and. flag_warn_cflz) then -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n, nz, nzmin, nzmax) - do n=1, myDim_nod2D - nzmin = ulevels_nod2D(n) - nzmax = nlevels_nod2D(n)-1 - do nz=nzmin,nzmax - if (abs(CFL_z(nz,n)-cflmax) < 1.e-12 .and. CFL_z(nz,n) > 1.75_WP .and. CFL_z(nz,n)<=2.5_WP ) then - print '(A, A, F4.2, A, I6, A, F7.2,A,F6.2, A, I3,I3)', achar(27)//'[33m'//' --> WARNING CFLz>1.75:'//achar(27)//'[0m',& - 'CFLz_max=',cflmax,',mstep=',mstep,',glon/glat=',geo_coord_nod2D(1,n)/rad,'/',geo_coord_nod2D(2,n)/rad,& - ',nz/nzmin=',nz,nzmin - elseif (abs(CFL_z(nz,n)-cflmax) < 1.e-12 .and. CFL_z(nz,n) > 2.5_WP) then - print '(A, A, F4.2, A, I6, A, F7.2,A,F6.2, A, I3,I3)', achar(27)//'[31m'//' --> WARNING CFLz>2.5:'//achar(27)//'[0m',& - 'CFLz_max=',cflmax,',mstep=',mstep,',glon/glat=',geo_coord_nod2D(1,n)/rad,'/',geo_coord_nod2D(2,n)/rad,& - ',nz/nzmin=',nz,nzmin +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax) + do node=1, myDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2D(node)-1 + do nz=nzmin, nzmax + if (abs(CFL_z(nz,node)-cflmax) < 1.e-12 .and. CFL_z(nz,node) > 1.75_WP .and. CFL_z(nz,node)<=2.5_WP ) then + print '(A, A, F4.2, A, I6, A, F7.2,A,F6.2, A, I3,I3,I3)', achar(27)//'[33m'//' --> WARNING CFLz>1.75:'//achar(27)//'[0m',& + 'CFLz_max=',cflmax,',mstep=',mstep,',glon/glat=',geo_coord_nod2D(1,node)/rad,'/',geo_coord_nod2D(2,node)/rad,& + ',nz/nzmin/nzmax=',nz,nzmin,nzmax + elseif (abs(CFL_z(nz,node)-cflmax) < 1.e-12 .and. CFL_z(nz,node) > 2.5_WP) then + print '(A, A, F4.2, A, I6, A, F7.2,A,F6.2, A, I3,I3,I3)', achar(27)//'[31m'//' --> WARNING CFLz>2.5:'//achar(27)//'[0m',& + 'CFLz_max=',cflmax,',mstep=',mstep,',glon/glat=',geo_coord_nod2D(1,node)/rad,'/',geo_coord_nod2D(2,node)/rad,& + ',nz/nzmin/nzmax=',nz,nzmin,nzmax !!PS write(*,*) '***********************************************************' !!PS write(*,*) 'max. CFL_z = ', cflmax, ' mype = ', mype !!PS write(*,*) 'mstep = ', mstep - !!PS write(*,*) 'glon, glat = ', geo_coord_nod2D(:,n)/rad - !!PS write(*,*) '2D node = ', myList_nod2D(n) + !!PS write(*,*) 'glon, glat = ', geo_coord_nod2D(:,node)/rad + !!PS write(*,*) '2D node = ', myList_nod2D(node) !!PS write(*,*) 'nz = ', nz !!PS write(*,*) '***********************************************************' end if - end do - end do + end do ! --> do nz=nzmin,nzmax + end do ! --> do node=1, myDim_nod2D !$OMP END PARALLEL DO - end if - + end if ! --> if (cflmax > 1.0_WP .and. flag_warn_cflz) then +end subroutine compute_CFLz +! +! +!_______________________________________________________________________________ +subroutine compute_Wvel_split(dynamics, partit, mesh) + use MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + use o_PARAM + use g_comm_auto + implicit none + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh), intent(inout), target :: mesh + !___________________________________________________________________________ + integer :: node, nz, nzmin, nzmax + real(kind=WP) :: dd + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:) , pointer :: Wvel, Wvel_e, Wvel_i, CFL_z +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + Wvel => dynamics%w( :,:) + Wvel_e => dynamics%w_e( :,:) + Wvel_i => dynamics%w_i( :,:) + CFL_z => dynamics%cfl_z(:,:) + !___________________________________________________________________________ ! Split implicit vertical velocity onto implicit and explicit components using CFL criteria: ! wsplit_maxcfl constrains the allowed explicit w according to the CFL at this place @@ -2591,23 +2992,23 @@ subroutine vert_vel_ale(dynamics, partit, mesh) ! wsplit_maxcfl=0 means w_exp is zero (everything computed implicitly) ! wsplit_maxcfl=inf menas w_impl is zero (everything computed explicitly) ! a guess for optimal choice of wsplit_maxcfl would be 0.95 -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n, nz, nzmin, nzmax, dd) - do n=1, myDim_nod2D+eDim_nod2D - nzmin = ulevels_nod2D(n) - nzmax = nlevels_nod2D(n) - do nz=nzmin,nzmax - Wvel_e(nz,n)=Wvel(nz,n) - Wvel_i(nz,n)=0.0_WP - if (dynamics%use_wsplit .and. (CFL_z(nz, n) > dynamics%wsplit_maxcfl)) then - dd=max((CFL_z(nz, n)-dynamics%wsplit_maxcfl), 0.0_WP)/max(dynamics%wsplit_maxcfl, 1.e-12) - Wvel_e(nz,n)=(1.0_WP/(1.0_WP+dd))*Wvel(nz,n) !explicit part =1. if dd=0. - Wvel_i(nz,n)=(dd /(1.0_WP+dd))*Wvel(nz,n) !implicit part =1. if dd=inf +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax, dd) + do node=1, myDim_nod2D+eDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2D(node) + do nz=nzmin, nzmax + Wvel_e(nz, node)=Wvel(nz, node) + Wvel_i(nz, node)=0.0_WP + if (dynamics%use_wsplit .and. (CFL_z(nz, node) > dynamics%wsplit_maxcfl)) then + dd=max((CFL_z(nz, node)-dynamics%wsplit_maxcfl), 0.0_WP)/max(dynamics%wsplit_maxcfl, 1.e-12) + Wvel_e(nz, node)=(1.0_WP/(1.0_WP+dd))*Wvel(nz, node) !explicit part =1. if dd=0. + Wvel_i(nz, node)=(dd /(1.0_WP+dd))*Wvel(nz, node) !implicit part =1. if dd=inf end if - end do - end do + end do ! --> do nz=nzmin,nzmax + end do ! --> do node=1, myDim_nod2D+eDim_nod2D !$OMP END PARALLEL DO -end subroutine vert_vel_ale - +end subroutine compute_Wvel_split +! ! ! !=============================================================================== @@ -2878,6 +3279,8 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) use g_cvmix_tidal use Toy_Channel_Soufflet use oce_ale_interfaces + use compute_vert_vel_transpv_interface + use compute_ssh_split_explicit_interface use pressure_bv_interface use pressure_force_4_linfs_interface use pressure_force_4_zxxxx_interface @@ -2886,6 +3289,8 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) use write_step_info_interface use check_blowup_interface use fer_solve_interface + use impl_vert_visc_ale_vtransp_interface + IMPLICIT NONE integer , intent(in) :: n type(t_dyn) , intent(inout), target :: dynamics @@ -2908,10 +3313,10 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, 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 @@ -2976,9 +3381,9 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call oce_mixing_KPP'//achar(27)//'[0m' call oce_mixing_KPP(Av, Kv_double, dynamics, tracers, partit, mesh) !$OMP PARALLEL DO - DO node=1, myDim_nod2D+eDim_nod2D + do node=1, myDim_nod2D+eDim_nod2D Kv(:, node)=Kv_double(:, node, 1) - END DO + end do !$OMP END PARALLEL DO call mo_convect(ice, partit, mesh) @@ -3029,115 +3434,209 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) t1=MPI_Wtime() !___________________________________________________________________________ + ! add contribution from momentum advection, coriolis and pressure gradient | + ! force to UV_rhs + ! UV_rhs = dt*[ (R_advec + R_coriolis)_AB2^n + R_pressure^n ] if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_vel_rhs'//achar(27)//'[0m' call compute_vel_rhs(ice, dynamics, partit, mesh) !___________________________________________________________________________ - if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call viscosity_filter'//achar(27)//'[0m' + ! Energy diagnostic contribution if (dynamics%ldiag_ke) then + ! if use solver + if (.not. dynamics%use_ssh_se_subcycl) then !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) - do elem=1, myDim_elem2D - nzmax = nlevels(elem) - nzmin = ulevels(elem) - do nz=nzmin,nzmax-1 - dynamics%ke_rhs_bak(:,nz,elem)=dynamics%UV_rhs(:,nz,elem) - end do - end do + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_rhs_bak(:,nz,elem)=dynamics%UV_rhs(:,nz,elem) + end do + end do !$OMP END PARALLEL DO + + ! if use splitexpl subcycl. UV_rhs in units of transport --> therefor *1/helem + else +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_rhs_bak(:,nz,elem)=dynamics%UV_rhs(:,nz,elem)/helem(nz,elem) + end do + end do +!$OMP END PARALLEL DO + end if end if + + !___________________________________________________________________________ + ! add contribution from horizontal viscosity to UV_rhs + ! UV_rhs = dt*[ (R_advec + R_coriolis)^n + R_pressure + R_hviscos] + ! UV_rhs = UV_rhs + dt*R_hviscos + if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call viscosity_filter'//achar(27)//'[0m' call viscosity_filter(dynamics%opt_visc, dynamics, partit, mesh) + + !___________________________________________________________________________ + ! Energy diagnostic contribution if (dynamics%ldiag_ke) then + ! if use solver + if (.not. dynamics%use_ssh_se_subcycl) then !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) - do elem=1, myDim_elem2D - nzmax = nlevels(elem) - nzmin = ulevels(elem) - do nz=nzmin,nzmax-1 - dynamics%ke_hvis(:,nz,elem)=dynamics%UV_rhs(:,nz,elem)-dynamics%ke_rhs_bak(:,nz,elem) - end do - end do + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_hvis(:,nz,elem)=dynamics%UV_rhs(:,nz,elem) - dynamics%ke_rhs_bak(:,nz,elem) + end do + end do !$OMP END PARALLEL DO +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_rhs_bak(:,nz,elem)=dynamics%UV_rhs(:,nz,elem) + end do + end do +!$OMP END PARALLEL DO + + ! if use splitexpl subcycl. UV_rhs in units of transport --> therefor *1/helem + else +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_hvis(:,nz,elem)=dynamics%UV_rhs(:,nz,elem)/helem(nz,elem) - dynamics%ke_rhs_bak(:,nz,elem) + end do + end do +!$OMP END PARALLEL DO +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_rhs_bak(:,nz,elem)=dynamics%UV_rhs(:,nz,elem)/helem(nz,elem) + end do + end do +!$OMP END PARALLEL DO + + end if end if !___________________________________________________________________________ if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call impl_vert_visc_ale'//achar(27)//'[0m' - if (dynamics%ldiag_ke) then -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) - do elem=1, myDim_elem2D - nzmax = nlevels(elem) - nzmin = ulevels(elem) - do nz=nzmin,nzmax-1 - dynamics%ke_rhs_bak(:,nz,elem)=dynamics%UV_rhs(:,nz,elem) - end do - end do -!$OMP END PARALLEL DO + if(dynamics%use_ivertvisc) then + if ( .not. dynamics%use_ssh_se_subcycl ) then + call impl_vert_visc_ale(dynamics,partit, mesh) + else + call impl_vert_visc_ale_vtransp(dynamics, partit, mesh) + end if end if - if(dynamics%use_ivertvisc) call impl_vert_visc_ale(dynamics,partit, mesh) + + !___________________________________________________________________________ if (dynamics%ldiag_ke) then + ! if use solver + if (.not. dynamics%use_ssh_se_subcycl) then !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) - do elem=1, myDim_elem2D - nzmax = nlevels(elem) - nzmin = ulevels(elem) - do nz=nzmin,nzmax-1 - dynamics%ke_vvis(:,nz,elem)=dynamics%UV_rhs(:,nz,elem)-dynamics%ke_rhs_bak(:,nz,elem) - end do - end do + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_vvis(:,nz,elem)=dynamics%UV_rhs(:,nz,elem) - dynamics%ke_rhs_bak(:,nz,elem) + end do + end do !$OMP END PARALLEL DO + + ! if use splitexpl subcycl. UV_rhs in units of transport --> therefor *1/helem + else + do elem=1, myDim_elem2D + nzmax = nlevels(elem) + nzmin = ulevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_vvis(:,nz,elem)=dynamics%UV_rhs(:,nz,elem)/helem(nz,elem) - dynamics%ke_rhs_bak(:,nz,elem) + end do + end do + end if end if t2=MPI_Wtime() !___________________________________________________________________________ ! >->->->->->->->->->->->-> ALE-part starts <-<-<-<-<-<-<-<-<-<-<-<- !___________________________________________________________________________ + + !___________________________________________________________________________ + ! Compute SSH via solver ! Update stiffness matrix by dhe=hbar(n+1/2)-hbar(n-1/2) on elements, only ! needed for zlevel and zstar - if (.not. trim(which_ale)=='linfs') call update_stiff_mat_ale(partit, mesh) - if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_ssh_rhs_ale'//achar(27)//'[0m' - ! ssh_rhs=-alpha*\nabla\int(U_n+U_rhs)dz-(1-alpha)*... - ! see "FESOM2: from finite elements to finte volumes, S. Danilov..." eq. (18) rhs - call compute_ssh_rhs_ale(dynamics, partit, mesh) - - ! Take updated ssh matrix and solve --> new ssh! - t30=MPI_Wtime() - call solve_ssh_ale(dynamics, partit, mesh) - - if ((toy_ocean) .AND. (TRIM(which_toy)=="soufflet")) then - if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call relax_zonal_vel'//achar(27)//'[0m' - call relax_zonal_vel(dynamics, partit, mesh) - end if - t3=MPI_Wtime() - - ! estimate new horizontal velocity u^(n+1) - ! u^(n+1) = u* + [-g * tau * theta * grad(eta^(n+1)-eta^(n)) ] - if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call update_vel'//achar(27)//'[0m' - ! ke will be computed inside there if dynamics%ldiag_ke is .TRUE. - call update_vel(dynamics, partit, mesh) - - ! --> eta_(n) --> eta_(n+1) = eta_(n) + deta = eta_(n) + (eta_(n+1) + eta_(n)) - t4=MPI_Wtime() - - ! Update to hbar(n+3/2) and compute dhe to be used on the next step - if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_hbar_ale'//achar(27)//'[0m' - call compute_hbar_ale(dynamics, partit, 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... - ! - where(ulevels_nod2D==1) is used here because of the rigid lid - ! approximation under the cavity - ! - at points in the cavity the time derivative term in ssh matrix will be - ! omitted; and (14) will not be applied at cavity points. Additionally, - ! since there is no real elevation, but only surface pressure, there is - ! no layer motion under the cavity. In this case the ice sheet acts as a - ! rigid lid. + if (.not. dynamics%use_ssh_se_subcycl) then + if (.not. trim(which_ale)=='linfs') call update_stiff_mat_ale(partit, mesh) + if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_ssh_rhs_ale'//achar(27)//'[0m' + ! ssh_rhs=-alpha*\nabla\int(U_n+U_rhs)dz-(1-alpha)*... + ! see "FESOM2: from finite elements to finte volumes, S. Danilov..." eq. (18) rhs + call compute_ssh_rhs_ale(dynamics, partit, mesh) + + ! Take updated ssh matrix and solve --> new ssh! + t30=MPI_Wtime() + call solve_ssh_ale(dynamics, partit, mesh) + + if ((toy_ocean) .AND. (TRIM(which_toy)=="soufflet")) then + if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call relax_zonal_vel'//achar(27)//'[0m' + call relax_zonal_vel(dynamics, partit, mesh) + end if + t3=MPI_Wtime() + + ! estimate new horizontal velocity u^(n+1) + ! u^(n+1) = u* + [-g * tau * theta * grad(eta^(n+1)-eta^(n)) ] + if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call update_vel'//achar(27)//'[0m' + ! ke will be computed inside there if dynamics%ldiag_ke is .TRUE. + call update_vel(dynamics, partit, mesh) + + ! --> eta_(n) --> eta_(n+1) = eta_(n) + deta = eta_(n) + (eta_(n+1) + eta_(n)) + t4=MPI_Wtime() + + ! Update to hbar(n+3/2) and compute dhe to be used on the next step + if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call compute_hbar_ale'//achar(27)//'[0m' + call compute_hbar_ale(dynamics, partit, 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... + ! - where(ulevels_nod2D==1) is used here because of the rigid lid + ! approximation under the cavity + ! - at points in the cavity the time derivative term in ssh matrix will be + ! omitted; and (14) will not be applied at cavity points. Additionally, + ! since there is no real elevation, but only surface pressure, there is + ! no layer motion under the cavity. In this case the ice sheet acts as a + ! rigid lid. !$OMP PARALLEL DO - do node=1, myDim_nod2D+eDim_nod2D - eta_n(node)=alpha*hbar(node)+(1.0_WP-alpha)*hbar_old(node) - end do + do node=1, myDim_nod2D+eDim_nod2D + eta_n(node)=alpha*hbar(node)+(1.0_WP-alpha)*hbar_old(node) + end do !$OMP END PARALLEL DO - ! --> eta_(n) - ! call zero_dynamics !DS, zeros several dynamical variables; to be used for testing new implementations! - t5=MPI_Wtime() + ! --> eta_(n) + ! call zero_dynamics !DS, zeros several dynamical variables; to be used for testing new implementations! + t5=MPI_Wtime() + + !___________________________________________________________________________ + ! Compute SSH via split-explicite subcycling + else + ! Compute vertical integral of transport velocity rhs omitting the contributions from + ! the elevation and Coriolis. + t30=MPI_Wtime() + call compute_BT_rhs_SE_vtransp(dynamics, partit, mesh) + + ! Do barotropic step, get eta_{n+1} and BT transport + call compute_BT_step_SE_ale(dynamics, partit, mesh) + t3=MPI_Wtime() + + ! Trim U to be consistent with BT transport + call update_trim_vel_ale_vtransp(1, dynamics, partit, mesh) + t4=MPI_Wtime() + t5=t4 + end if ! --> if (.not. dynamics%use_ssh_se_subcycl) then + !___________________________________________________________________________ ! Do horizontal and vertical scaling of GM/Redi diffusivity if (Fer_GM .or. Redi) then @@ -3154,23 +3653,38 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) t6=MPI_Wtime() !___________________________________________________________________________ - ! The main step of ALE procedure --> this is were the magic happens --> here - ! is decided how change in hbar is distributed over the vertical layers - if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call vert_vel_ale'//achar(27)//'[0m' ! keep the old vertical velocity for computation of the mean between the timesteps (is used in compute_ke_wrho) if (dynamics%ldiag_ke) then !$OMP PARALLEL DO - do node=1, myDim_nod2D+eDim_nod2D - dynamics%w_old(:, node)=dynamics%w(:, node) - end do + do node=1, myDim_nod2D+eDim_nod2D + dynamics%w_old(:, node)=dynamics%w(:, node) + end do !$OMP END PARALLEL DO end if - call vert_vel_ale(dynamics, partit, mesh) + + !___________________________________________________________________________ + ! The main step of ALE procedure --> this is were the magic happens --> here + ! is decided how change in hbar is distributed over the vertical layers + if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call vert_vel_ale'//achar(27)//'[0m' + if ( .not. dynamics%use_ssh_se_subcycl) then + call vert_vel_ale(dynamics, partit, mesh) + else + if (trim(which_ale)=='zstar' ) then + call compute_thickness_zstar(dynamics, partit, mesh) + else + hnode_new = hnode + end if + call compute_vert_vel_transpv(dynamics, partit, mesh) + end if t7=MPI_Wtime() + + !___________________________________________________________________________ + ! energy diagnostic computation if (dynamics%ldiag_ke) then call compute_ke_wrho(dynamics, partit, mesh) call compute_apegen (dynamics, tracers, partit, mesh) end if + !___________________________________________________________________________ ! solve tracer equation if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call solve_tracers_ale'//achar(27)//'[0m' @@ -3182,13 +3696,22 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call update_thickness_ale'//achar(27)//'[0m' call update_thickness_ale(partit, mesh) t9=MPI_Wtime() + +!PS !___________________________________________________________________________ +!PS ! Trim to make velocity consistent with BT velocity at n+1/2 +!PS ! Velocity will be used to advance momentum +!PS if (dynamics%use_ssh_se_subcycl) then +!PS call update_trim_vel_ale_vtransp(2, dynamics, partit, mesh) +!PS end if + !___________________________________________________________________________ ! write out global fields for debugging if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call write_step_info'//achar(27)//'[0m' call write_step_info(n,logfile_outfreq, ice, dynamics, tracers, partit, mesh) + !___________________________________________________________________________ ! write energy diagnostic info (dynamics%ldiag_ke = .true.) - if (dynamics%ldiag_ke) then + if ( (dynamics%ldiag_ke) .and. (mod(n,logfile_outfreq)==0) ) then call write_enegry_info(dynamics, partit, mesh) end if diff --git a/src/oce_ale_ssh_splitexpl_subcycl.F90 b/src/oce_ale_ssh_splitexpl_subcycl.F90 new file mode 100644 index 000000000..f2d795865 --- /dev/null +++ b/src/oce_ale_ssh_splitexpl_subcycl.F90 @@ -0,0 +1,1514 @@ +! +! +!_______________________________________________________________________________ +module momentum_adv_scalar_transpv_interface + interface + subroutine momentum_adv_scalar_transpv(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(in) , target :: mesh + end subroutine + end interface +end module +! +! +!_______________________________________________________________________________ +module impl_vert_visc_ale_vtransp_interface + interface + subroutine impl_vert_visc_ale_vtransp(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + end interface +end module +! +! +!_______________________________________________________________________________ +module compute_ssh_split_explicit_interface + interface + subroutine compute_BT_rhs_SE_vtransp(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + + subroutine compute_BT_step_SE_ale(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + + subroutine update_trim_vel_ale_vtransp(mode, dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + integer , intent(in) :: mode + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + end subroutine + + end interface +end module +! +! +!_______________________________________________________________________________ +! Transports are used instead of velocities, Urhs, Vrhs are also for transports. +subroutine momentum_adv_scalar_transpv(dynamics, partit, mesh) + USE MOD_MESH + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + USE o_PARAM + USE g_comm_auto + IMPLICIT NONE + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(in) , target :: mesh + !___________________________________________________________________________ + integer :: node, elem, ed, nz + integer :: nl1, ul1, nl2, ul2, nl12, ul12 + real(kind=WP) :: uv12, uv1, uv2, qc, qu, qd, num_ord=0.95_WP + integer :: ednodes(2), edelem(2) + real(kind=WP) :: wu(mesh%nl), wv(mesh%nl), un1(mesh%nl), un2(mesh%nl) + + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:,:,:), pointer :: UV_rhsAB + real(kind=WP), dimension(:,:,:) , pointer :: UV, UVnode_rhs, UVnode, UVh + real(kind=WP), dimension(:,:) , pointer :: Wvel_e +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + UV =>dynamics%uv(:,:,:) + UV_rhsAB =>dynamics%uv_rhsAB(:,:,:,:) + UVnode_rhs=>dynamics%work%uvnode_rhs(:,:,:) + UVnode =>dynamics%uvnode(:,:,:) + Wvel_e =>dynamics%w_e(:,:) + UVh =>dynamics%se_uvh(:,:,:) + !___________________________________________________________________________ + ! 1st. compute vertical momentum advection component: w * du/dz, w*dv/dz +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(node, elem, ed, nz, nl1, ul1, nl2, ul2, nl12, ul12, & +!$OMP uv1, uv2, uv12, qc, qu, qd, wu, wv, & +!$OMP ednodes, edelem, un1, un2) +!$OMP DO + do node=1, myDim_nod2D + ul1 = ulevels_nod2D(node) + nl1 = nlevels_nod2D(node) + wu(1:nl1) = 0.0_WP + wv(1:nl1) = 0.0_WP + UVnode_rhs(:,:,node) = 0.0_WP + !_______________________________________________________________________ + ! surface + nz=ul1 + wu(nz)=UVnode(1, nz, node)*Wvel_e(nz, node)*area(nz, node) + wv(nz)=UVnode(2, nz, node)*Wvel_e(nz, node)*area(nz, node) + + !_______________________________________________________________________ + ! subsurface --> centered 2nd order + nz=ul1+1 ! Central differences at the second and last but one levels + wu(nz)=0.5_WP*(UVnode(1, nz, node)+UVnode(1, nz-1, node))*Wvel_e(nz, node)*area(nz, node) + wv(nz)=0.5_WP*(UVnode(2, nz, node)+UVnode(2, nz-1, node))*Wvel_e(nz, node)*area(nz, node) + + !_______________________________________________________________________ + ! bulk --> centered 4th order (num_ord=1), upwind 3rd order (num_ord=0) + do nz=ul1+2, nl1-2 + qc = (UVnode(1, nz-1, node)-UVnode(1, nz , node))/(hnode(nz-1, node)+hnode(nz , node)) ! factor 2 is accounted for later + qu = (UVnode(1, nz , node)-UVnode(1, nz+1, node))/(hnode(nz , node)+hnode(nz+1, node)) + qd = (UVnode(1, nz-2, node)-UVnode(1, nz-1, node))/(hnode(nz-2, node)+hnode(nz-1, node)) + + uv1 = UVnode(1, nz , node)+(2*qc+qu)*hnode(nz , node)/6.0_WP ! Gradient reconstruction 2(2qc+qu)(h/2)(1/6) + uv2 = UVnode(1, nz-1, node)-(2*qc+qd)*hnode(nz-1, node)/6.0_WP + uv12 = (Wvel_e(nz, node)+abs(Wvel_e(nz, node)))*uv1+ & + (Wvel_e(nz, node)-abs(Wvel_e(nz, node)))*uv2 + wu(nz)=0.5_WP*(num_ord*(uv1+uv2)*Wvel_e(nz, node)+(1.0_WP-num_ord)*uv12)*area(nz, node) + + qc = (UVnode(2, nz-1, node)-UVnode(2, nz , node))/(hnode(nz-1, node)+hnode(nz , node)) + qu = (UVnode(2, nz , node)-UVnode(2, nz+1, node))/(hnode(nz , node)+hnode(nz+1, node)) + qd = (UVnode(2, nz-2, node)-UVnode(2, nz-1, node))/(hnode(nz-2, node)+hnode(nz-1, node)) + + uv1 = UVnode(2, nz , node)+(2*qc+qu)*hnode(nz , node)/6.0_WP + uv2 = UVnode(2, nz-1, node)-(2*qc+qd)*hnode(nz-1, node)/6.0_WP + uv12 = (Wvel_e(nz, node)+abs(Wvel_e(nz, node)))*uv1+ & + (Wvel_e(nz, node)-abs(Wvel_e(nz, node)))*uv2 + wv(nz)=0.5_WP*(num_ord*(uv1+uv2)*Wvel_e(nz, node)+(1.0_WP-num_ord)*uv12)*area(nz, node) + end do ! --> do nz=ul1+2, nl1-2 + + !_______________________________________________________________________ + ! one layer above bottom --> centered 2nd order + nz=nl1-1 + wu(nz)=0.5_WP*(UVnode(1, nz, node)+UVnode(1, nz-1, node))*Wvel_e(nz, node)*area(nz, node) + wv(nz)=0.5_WP*(UVnode(2, nz, node)+UVnode(2, nz-1, node))*Wvel_e(nz, node)*area(nz, node) + + !_______________________________________________________________________ + ! bottom layer --> boundary condition + nz=nl1 + wu(nz)=0.0_WP + wv(nz)=0.0_WP + + !_______________________________________________________________________ + ! set to the rhs for transports, not velocities!!! --> No division by h + do nz=ul1, nl1-1 + UVnode_rhs(1, nz, node)= UVnode_rhs(1, nz, node) - (wu(nz)-wu(nz+1)) + UVnode_rhs(2, nz, node)= UVnode_rhs(2, nz, node) - (wv(nz)-wv(nz+1)) + end do + + end do ! --> do node=1, myDim_nod2D +!$OMP END DO + + + !___________________________________________________________________________ + ! 2nd. compute horizontal advection component: u*du/dx, u*dv/dx & v*du/dy, v*dv/dy + ! loop over triangle edges +!$OMP DO + do ed=1, myDim_edge2D + ! local indice of nodes that span up edge ed + ednodes = edges(:,ed) + + ! local index of element that contribute to edge + edelem = edge_tri(1:2,ed) + + !_______________________________________________________________________ + ! index off surface layer in case of cavity !=1 and index of mid depth + ! bottom layer + ul1 = ulevels(edelem(1)) + nl1 = nlevels(edelem(1))-1 + + !_______________________________________________________________________ + !NR --> Natalja Style + un1 = 0.0_WP + un1(ul1:nl1) = ( UVh(2, ul1:nl1, edelem(1))*edge_cross_dxdy(1,ed) & + - UVh(1, ul1:nl1, edelem(1))*edge_cross_dxdy(2,ed)) + + !_______________________________________________________________________ + ! if edelem(2)==0 than edge is boundary edge + if(edelem(2)>0) then + ul2 = ulevels(edelem(2)) + nl2 = nlevels(edelem(2))-1 + + !___________________________________________________________________ + !NR --> Natalja Style + un2 = 0.0_WP + un2(ul2:nl2) = -( UVh(2, ul2:nl2, edelem(2))*edge_cross_dxdy(3,ed) & + - UVh(1, ul2:nl2, edelem(2))*edge_cross_dxdy(4,ed)) + + !___________________________________________________________________ + ! nl12 ... minimum number of layers -1 between element edelem(1) & edelem(2) that + ! contribute to edge ed + ! nu12 ... upper index of layers between element edelem(1) & edelem(2) that + ! contribute to edge ed + ! be carefull !!! --> if ed is a boundary edge than edelem(1)~=0 and edelem(2)==0 + ! that means nl1>0, nl2==0, nl12=min(nl1,nl2)=0 !!! + ul12 = max(ul1, ul2) + nl12 = min(nl1, nl2) + + !___________________________________________________________________ + ! ensure openmp numerical reproducability +#if defined(__openmp_reproducible) +!$OMP ORDERED +#endif + !___________________________________________________________________ + !NR add contribution to first edge node --> ednodes(1) + !NR Do not calculate on Halo nodes, as the result will not be used. + !NR The "if" is cheaper than the avoided computiations. + if (ednodes(1) <= myDim_nod2d) then +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(ednodes(1))) +#endif + ! cavity domain where only edelem(1) exist + do nz=ul1 , ul12-1 + UVnode_rhs(1, nz, ednodes(1)) = UVnode_rhs(1, nz, ednodes(1)) + un1(nz)*UV(1, nz, edelem(1)) + UVnode_rhs(2, nz, ednodes(1)) = UVnode_rhs(2, nz, ednodes(1)) + un1(nz)*UV(2, nz, edelem(1)) + end do + ! cavity domain where only edelem(2) exist + do nz=ul2 , ul12-1 + UVnode_rhs(1, nz, ednodes(1)) = UVnode_rhs(1, nz, ednodes(1)) + un2(nz)*UV(1, nz, edelem(2)) + UVnode_rhs(2, nz, ednodes(1)) = UVnode_rhs(2, nz, ednodes(1)) + un2(nz)*UV(2, nz, edelem(2)) + end do + ! bulk domain where edelem(1) and edelem(2) exist + do nz=ul12, nl12 + UVnode_rhs(1, nz, ednodes(1)) = UVnode_rhs(1, nz, ednodes(1)) + (un1(nz) + un2(nz))*(UV(1, nz, edelem(1))+UV(1, nz, edelem(2)))*0.5_WP + UVnode_rhs(2, nz, ednodes(1)) = UVnode_rhs(2, nz, ednodes(1)) + (un1(nz) + un2(nz))*(UV(2, nz, edelem(1))+UV(2, nz, edelem(2)))*0.5_WP + end do + ! bottom domain where only edelem(1) exist + do nz=nl12+1, nl1 + UVnode_rhs(1, nz, ednodes(1)) = UVnode_rhs(1, nz, ednodes(1)) + un1(nz)*UV(1, nz, edelem(1)) + UVnode_rhs(2, nz, ednodes(1)) = UVnode_rhs(2, nz, ednodes(1)) + un1(nz)*UV(2, nz, edelem(1)) + end do + ! bottom domain where only edelem(2) exist + do nz=nl12+1, nl2 + UVnode_rhs(1, nz, ednodes(1)) = UVnode_rhs(1, nz, ednodes(1)) + un2(nz)*UV(1, nz, edelem(2)) + UVnode_rhs(2, nz, ednodes(1)) = UVnode_rhs(2, nz, ednodes(1)) + un2(nz)*UV(2, nz, edelem(2)) + end do +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(1))) +#endif + end if + !___________________________________________________________________ + !NR add contribution to second edge node --> ednodes(2) + !NR Do not calculate on Halo nodes, as the result will not be used. + !NR The "if" is cheaper than the avoided computiations. + if (ednodes(2) <= myDim_nod2d) then +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(ednodes(2))) +#endif + ! cavity domain where only edelem(1) exist + do nz=ul1 , ul12-1 + UVnode_rhs(1, nz, ednodes(2)) = UVnode_rhs(1, nz, ednodes(2)) - un1(nz)*UV(1, nz, edelem(1)) + UVnode_rhs(2, nz, ednodes(2)) = UVnode_rhs(2, nz, ednodes(2)) - un1(nz)*UV(2, nz, edelem(1)) + end do + ! cavity domain where only edelem(2) exist + do nz=ul2 , ul12-1 + UVnode_rhs(1, nz, ednodes(2)) = UVnode_rhs(1, nz, ednodes(2)) - un2(nz)*UV(1, nz, edelem(2)) + UVnode_rhs(2, nz, ednodes(2)) = UVnode_rhs(2, nz, ednodes(2)) - un2(nz)*UV(2, nz, edelem(2)) + end do + ! bulk domain where edelem(1) and edelem(2) exist + do nz=ul12, nl12 + UVnode_rhs(1, nz, ednodes(2)) = UVnode_rhs(1, nz, ednodes(2)) - (un1(nz) + un2(nz))*(UV(1, nz, edelem(1))+UV(1, nz, edelem(2)))*0.5_WP + UVnode_rhs(2, nz, ednodes(2)) = UVnode_rhs(2, nz, ednodes(2)) - (un1(nz) + un2(nz))*(UV(2, nz, edelem(1))+UV(2, nz, edelem(2)))*0.5_WP + end do + ! bottom domain where only edelem(1) exist + do nz=nl12+1, nl1 + UVnode_rhs(1, nz, ednodes(2)) = UVnode_rhs(1, nz, ednodes(2)) - un1(nz)*UV(1, nz, edelem(1)) + UVnode_rhs(2, nz, ednodes(2)) = UVnode_rhs(2, nz, ednodes(2)) - un1(nz)*UV(2, nz, edelem(1)) + end do + ! bottom domain where only edelem(2) exist + do nz=nl12+1, nl2 + UVnode_rhs(1, nz, ednodes(2)) = UVnode_rhs(1, nz, ednodes(2)) - un2(nz)*UV(1, nz, edelem(2)) + UVnode_rhs(2, nz, ednodes(2)) = UVnode_rhs(2, nz, ednodes(2)) - un2(nz)*UV(2, nz, edelem(2)) + end do +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(2))) +#endif + end if + + !_______________________________________________________________________ + ! if edelem(2)==0 than edge is boundary edge + else ! --> if(edelem(2)>0) then + !___________________________________________________________________ + !NR add contribution to first edge node --> ednodes(1) + !NR Do not calculate on Halo nodes, as the result will not be used. + !NR The "if" is cheaper than the avoided computiations. + if (ednodes(1) <= myDim_nod2d) then +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(ednodes(1))) +#endif + ! bulk domain where only edelem(1) exist + do nz=ul1 , nl1 + UVnode_rhs(1, nz, ednodes(1)) = UVnode_rhs(1, nz, ednodes(1)) + un1(nz)*UV(1, nz, edelem(1)) + UVnode_rhs(2, nz, ednodes(1)) = UVnode_rhs(2, nz, ednodes(1)) + un1(nz)*UV(2, nz, edelem(1)) + end do +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(1))) +#endif + end if + !___________________________________________________________________ + !NR add contribution to second edge node --> ednodes(2) + !NR Do not calculate on Halo nodes, as the result will not be used. + !NR The "if" is cheaper than the avoided computiations. + if (ednodes(2) <= myDim_nod2d) then +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(ednodes(2))) +#endif + ! bulk domain where only edelem(1) exist + do nz=ul1 , nl1 + UVnode_rhs(1, nz, ednodes(2)) = UVnode_rhs(1, nz, ednodes(2)) - un1(nz)*UV(1, nz, edelem(1)) + UVnode_rhs(2, nz, ednodes(2)) = UVnode_rhs(2, nz, ednodes(2)) - un1(nz)*UV(2, nz, edelem(1)) + end do +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(2))) +#endif + end if + end if ! --> if(edelem(2)>0) then + +#if defined(__openmp_reproducible) +!$OMP END ORDERED +#endif + + end do ! --> do ed=1, myDim_edge2D +!$OMP END DO + + !___________________________________________________________________________ + ! divide total nodal momentum advection by scalar area +!$OMP DO + do node=1,myDim_nod2d + ul1 = ulevels_nod2D(node) + nl1 = nlevels_nod2D(node)-1 + UVnode_rhs(1, ul1:nl1, node) = UVnode_rhs(1, ul1:nl1, node) * areasvol_inv(ul1:nl1, node) + UVnode_rhs(2, ul1:nl1, node) = UVnode_rhs(2, ul1:nl1, node) * areasvol_inv(ul1:nl1, node) + end do ! --> do node=1,myDim_nod2d +!$OMP END DO + + !___________________________________________________________________________ +!$OMP MASTER + call exchange_nod(UVnode_rhs, partit) +!$OMP END MASTER +!$OMP BARRIER + + !___________________________________________________________________________ + ! convert total nodal advection from vertice --> elements +!$OMP DO + do elem=1, myDim_elem2D + ul1 = ulevels(elem) + nl1 = nlevels(elem)-1 + UV_rhsAB(1, 1:2, ul1:nl1, elem) = UV_rhsAB(1, 1:2, ul1:nl1, elem) + elem_area(elem)* & + ( UVnode_rhs(1:2, ul1:nl1, elem2D_nodes(1, elem)) & + + UVnode_rhs(1:2, ul1:nl1, elem2D_nodes(2, elem)) & + + UVnode_rhs(1:2, ul1:nl1, elem2D_nodes(3, elem))) / 3.0_WP + end do ! --> do edelem=1, myDim_elem2D +!$OMP END DO + + !___________________________________________________________________________ + ! for energz diagnostic + if (dynamics%ldiag_ke) then !we repeat the computation here and there are multiple ways to speed it up +!$OMP DO + do elem=1, myDim_elem2D + ul1 = ulevels(elem) + nl1 = nlevels(elem)-1 + dynamics%ke_adv_AB(1, 1, ul1:nl1, elem) = dynamics%ke_adv_AB(1, 1, ul1:nl1, elem) + elem_area(elem)* & + ( UVnode_rhs(1, ul1:nl1, elem2D_nodes(1, elem)) & + + UVnode_rhs(1, ul1:nl1, elem2D_nodes(2, elem)) & + + UVnode_rhs(1, ul1:nl1, elem2D_nodes(3, elem))) / 3.0_WP / helem(ul1:nl1, elem) + dynamics%ke_adv_AB(1, 2, ul1:nl1, elem) = dynamics%ke_adv_AB(1, 2, ul1:nl1, elem) + elem_area(elem)* & + ( UVnode_rhs(2, ul1:nl1, elem2D_nodes(1, elem)) & + + UVnode_rhs(2, ul1:nl1, elem2D_nodes(2, elem)) & + + UVnode_rhs(2, ul1:nl1, elem2D_nodes(3, elem))) / 3.0_WP / helem(ul1:nl1, elem) + end do +!$OMP END DO + end if +!$OMP END PARALLEL +end subroutine momentum_adv_scalar_transpv + +! +! +!_______________________________________________________________________________ +!SD Transport velocity version (T in the name) +!SD Solve U**-U*=dt*D_vert u**, where U*=Uh+rhs (rhs contains dt/elem_area) +!SD Express through Delta U=U**-U* and then introduce +!SD Delta u= Delta U/helem; after Delta u is found, rhs=rhs +h*Delta u. +!SD The rhs accumulated previously is not affected in this step, this is different +!SD from FESOM, and I would recomment this variant. +! +! +! solve equation vertical viscosity implicite part: +! --> U^(n+0.5,**) - U^(n+0.5,*) = dt*( Av * d/dz * u^(n+0.5,**) )|^t_b +! | +! +-> du = u^(n+0.5,**) - u^(n+0.5,*) = dU/h^(*) +! | +! +-> du * h^(*) = dU +! | u^(n+0.5,**) = du + u^(n+0.5,*) +! V +! du * h^(*) = dt*( Av*d/dz*du )|^t_b + dt*( Av*d/dz*u^(n+0.5,*) )|^t_b +! du - dt/h^(*)*( Av*d/dz*du ) = dt/h^(*)*( Av*d/dz*u^(n+0.5,*) ) +! --> solve for du +! ^ +! /|\ nvec_up (+1) +! | +! ----------- zbar_1, A_1 *----|----* +! Z_1 o T_1, Ac_1 |\ | ./| +! ----------- zbar_2, A_2 | \ ./ | Gaus Theorem: +! Z_2 o T_2, Ac_2 | \ / | --> Flux form +! ----------- zbar_3, A_3 | | | --> normal vec outwards facing +! Z_3 o T_3, Ac_3 *---|-----* +! ----------- zbar_4 \ | | ./ +! : \ | |/ +! \|/| +! * | +! V nvec_dwn (-1) +! +! --> 1st. solve homogenouse part: +! f(du) = du - dt/h^(*)*( Av*d/dz*du ) = 0 +! +! --> 2nd. compute difference quotient at du_i using Gauss-Theorem --> flux form +! | +! +-> Gauss THeorem: int(V', div(F_vec))dV = intcircle(A', F_vec*n_vec)dA +! | +! +-> du = dt/h^(*)*( Av*d/dz*du ) | *div()_z=d/dz, *int(V',)dV +! | int(V', d/dz *du)dV = int(V', d/dz *dt/h^(*)*( Av*d/dz*du ) )dV +! | ... = intcircle(A', dt/h^(*)*( Av*d/dz*du )*n_vec)dA +! | int(V', d/dz *du)dz*dA = ... +! | int(A', du)dA = intcircle(A', dt/h^(*)*( Av*d/dz*du )*n_vec)dA +! | +! +-> Ac_i area of vector cell = are of triangle, A_i, A_i+1 area of top +! | and bottom face +! V +! +! du_i * Ac_i = dt * [ Av_i /h_i * (du_i-1 - du_i )/(Z_i-1 - Z_i ) * A_i * nvec_up(+1) +! +Av_i+1/h_i * (du_i - du_i+1)/(Z_i - Z_i+1) * A_i+1 * nvec_dwn(-1) ] +! | +! +-> since we are on triangles Ac_i = A_i = A_i+1 --> can kick out A_i/Ac_i +! | and A_i+1/Ac_i +! +-> take into account normal vector direction +! V +! f(du_i) = du_i - dt*Av_i /h_i * (du_i-1 - du_i )/(Z_i-1 - Z_i ) +! + dt*Av_i+1/h_i * (du_i - du_i+1)/(Z_i - Z_i+1) +! = 0 +! +! --> 3rd. solve for coefficents a, b, c (homogenous part): +! f(du_i) = [ a*du_i-1 + b*du_i + c*du_i+1 ] +! | +! +-> estimate a, b, c by derivation of f(du_i) +! | +! +-> a = d[f(du_i)]/d[du_i-1] = - dt*Av_i/h_i / (Z_i-1 - Z_i) +! | +! +-> c = d[f(du_i)]/d[du_i+1] = - dt*Av_i+1/h_i / (Z_i - Z_i+1) +! | +! +-> b = d[f(du_i)]/d[du_i] = 1 + dt*Av_i/h_i / (Z_i-1 - Z_i) + dt*Av_i+1/h_i / (Z_i - Z_i+1) +! 1 - a - c +! +! --> 4th. solve inhomogenous part: +! [ a*du_i-1 + b*du_i + c*du_i+1 ] = RHS/A_i +! +! RHS/A_i = dt* [ Av_i /h_i * (u^(n+0.5,*)_i-1 - u^(n+0.5,*)_i )/(Z_i-1 - Z_i ) * A_i /A_i * nvec_up(+1) +! | +Av_i+1/h_i * (u^(n+0.5,*)_i - u^(n+0.5,*)_i+1)/(Z_i - Z_i+1) * A_i+1/A_i * nvec_dwn(-1) ] +! | +! +-> since we are on triangles A_i = A_i+1 --> can kick out A_i +! +-> take into account normal vector direction +! V +! = -a*u^(n+0.5,*)_i-1 + (a+c)*u^(n+0.5,*)_i - c*u^(n+0.5,*)_i+1 +! +! --> 5th. solve for du_i --> forward sweep algorithm --> see lower +! | b_1 c_1 ... | |du_1| +! | a_2 b_2 c_2 ... | |du_2| +! | a_3 b_3 c_3 ... | * |du_3| = RHS/A_i +! | a_4 b_4 c_4 ...| |du_4| +! | : | | : | +! +subroutine impl_vert_visc_ale_vtransp(dynamics, partit, mesh) + USE MOD_MESH + USE o_PARAM + USE o_ARRAYS, only: Av, stress_surf + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_DYN + USE g_CONFIG, only: dt + IMPLICIT NONE + !___________________________________________________________________________ + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + !___________________________________________________________________________ + real(kind=WP) :: a(mesh%nl-1), b(mesh%nl-1), c(mesh%nl-1), ur(mesh%nl-1), vr(mesh%nl-1) + real(kind=WP) :: cp(mesh%nl-1), up(mesh%nl-1), vp(mesh%nl-1), uu(mesh%nl-1), vv(mesh%nl-1) + integer :: elem, nz, nzmax, nzmin, elnodes(3) + real(kind=WP) :: zinv, m, friction, wu, wd + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:,:), pointer :: UV, UV_rhs, UVh + real(kind=WP), dimension(:,:) , pointer :: Wvel_i +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + UV =>dynamics%uv(:,:,:) + UV_rhs =>dynamics%uv_rhs(:,:,:) + Wvel_i =>dynamics%w_i(:,:) + UVh =>dynamics%se_uvh(:,:,:) + + !___________________________________________________________________________ +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(nz, nzmin, nzmax, elem, elnodes, & +!$OMP a, b, c, m, ur, vr, cp, up, vp, & +!$OMP zinv, friction, wu, wd, uu, vv) + + !___________________________________________________________________________ +!$OMP DO + do elem=1,myDim_elem2D + elnodes= elem2D_nodes(:,elem) + nzmin = ulevels(elem) + nzmax = nlevels(elem) + uu = 0.0_WP + vv = 0.0_WP + + !_______________________________________________________________________ + ! New velocities: compute u_k^(n+1/2, *) = u_k^(n-1/2) + UV_rhs/helem + ! UV_rhs is here in terms of transport velocity, we need real velocity + ! therefor divide with helem + uu(nzmin:nzmax-1)=(UVh(1, nzmin:nzmax-1, elem)+UV_rhs(1, nzmin:nzmax-1, elem))/helem(nzmin:nzmax-1, elem) !! u*=U*/h + vv(nzmin:nzmax-1)=(UVh(2, nzmin:nzmax-1, elem)+UV_rhs(2, nzmin:nzmax-1, elem))/helem(nzmin:nzmax-1, elem) + + !_______________________________________________________________________ + ! Operator + rhs + ! Regular part of coefficients: + + !_____1st layer (surface)_____ + nz = nzmin + zinv = 2.0_WP*dt/helem(nz, elem) + c( nz)= -Av(nz+1,elem)/(helem(nz,elem)+helem(nz+1,elem))*zinv + a( nz)= 0.0_WP + b( nz)= -c(nz)+1.0_WP + + ! Update from the vertical advection + wu = sum(Wvel_i(nz , elnodes))/3._WP + wd = sum(Wvel_i(nz+1, elnodes))/3._WP + b( nz)= b(nz)+wu*zinv + b( nz)= b(nz)-min(0._WP, wd)*zinv + c( nz)= c(nz)-max(0._WP, wd)*zinv + + ur(nz)= +c(nz)*uu(nz)-c(nz)*uu(nz+1) + vr(nz)= +c(nz)*vv(nz)-c(nz)*vv(nz+1) + + !_____bulk layers_____ + do nz=nzmin+1, nzmax-2 + zinv = 2.0_WP*dt/helem(nz, elem) + a( nz)= -Av(nz , elem)/(helem(nz-1, elem)+helem(nz , elem))*zinv + c( nz)= -Av(nz+1, elem)/(helem(nz , elem)+helem(nz+1, elem))*zinv + b( nz)= -a(nz)-c(nz)+1.0_WP + + ! Update from the vertical advection + wu=sum(Wvel_i(nz , elnodes))/3._WP + wd=sum(Wvel_i(nz+1, elnodes))/3._WP + a( nz)= a(nz)+min(0._WP, wu)*zinv + b( nz)= b(nz)+max(0._WP, wu)*zinv + + ur(nz)= -a(nz)*uu(nz-1)+(a(nz)+c(nz))*uu(nz)-c(nz)*uu(nz+1) !! add dt*D_vert u* + vr(nz)= -a(nz)*vv(nz-1)+(a(nz)+c(nz))*vv(nz)-c(nz)*vv(nz+1) !! to the rhs (because of Delta u) + end do + + !_____bottom layer_____ + nz = nzmax-1 + zinv = 2.0_WP*dt/helem(nz,elem) + a( nz)= -Av(nz, elem)/(helem(nz-1, elem)+helem(nz,elem))*zinv + b( nz)= -a(nz)+1.0_WP + c( nz)= 0.0_WP + + ! Update from the vertical advection + wu = sum(Wvel_i(nz-1, elnodes))/3._WP + a( nz)= a(nz)+min(0._WP, wu)*zinv + b( nz)= b(nz)+max(0._WP, wu)*zinv + + ur(nz)= -a(nz)*uu(nz-1)+a(nz)*uu(nz) + vr(nz)= -a(nz)*vv(nz-1)+a(nz)*vv(nz) + + !_______________________________________________________________________ + ! Forcing to the rhs: + ! The first row contains surface forcing + nz = nzmin + zinv = 1.0_WP*dt/helem(nz, elem) + ur(nz)= ur(nz) + zinv*stress_surf(1, elem)/density_0 + vr(nz)= vr(nz) + zinv*stress_surf(2, elem)/density_0 + if (dynamics%ldiag_ke) then + dynamics%ke_wind(1,elem)=stress_surf(1, elem)/density_0*dt + dynamics%ke_wind(2,elem)=stress_surf(2, elem)/density_0*dt + end if + + ! The last row contains bottom friction + nz = nzmax-1 + zinv = 1.0_WP*dt/helem(nz, elem) + friction=-C_d*sqrt(UV(1, nz, elem)**2 + UV(2, nz, elem)**2) + ur(nz)= ur(nz) + zinv*friction*UV(1, nz, elem) + vr(nz)= vr(nz) + zinv*friction*UV(2, nz, elem) + if (dynamics%ldiag_ke) then + dynamics%ke_drag(1, elem)=friction*UV(1, nz, elem)*dt + dynamics%ke_drag(2, elem)=friction*UV(2, nz, elem)*dt + end if + + !_______________________________________________________________________ + ! The sweep algorithm --> solve for du_i + ! initialize c-prime and s,t-prime + cp(nzmin) = c(nzmin)/b(nzmin) + up(nzmin) = ur(nzmin)/b(nzmin) + vp(nzmin) = vr(nzmin)/b(nzmin) + + ! sove for vectors c-prime and t, s-prime + do nz=nzmin+1, nzmax-1 + m = b(nz)-cp(nz-1)*a(nz) + cp(nz) = c(nz)/m + up(nz) = (ur(nz)-up(nz-1)*a(nz))/m + vp(nz) = (vr(nz)-vp(nz-1)*a(nz))/m + end do + + ! initialize x + ur(nzmax-1) = up(nzmax-1) + vr(nzmax-1) = vp(nzmax-1) + + ! solve for x from the vectors c-prime and d-prime + do nz=nzmax-2, nzmin, -1 + ur(nz) = up(nz)-cp(nz)*ur(nz+1) + vr(nz) = vp(nz)-cp(nz)*vr(nz+1) + end do + + !_______________________________________________________________________ + ! RHS update + ! ur and vr are in terms of du velocity difference, UV_rhs are in terms + ! of transport velocity + do nz=nzmin, nzmax-1 + UV_rhs(1, nz, elem)=UV_rhs(1, nz, elem)+ur(nz)*helem(nz, elem) !! The rhs is for transport + UV_rhs(2, nz, elem)=UV_rhs(2, nz, elem)+vr(nz)*helem(nz, elem) + end do + + !_______________________________________________________________ + + end do ! --> do elem=1,myDim_elem2D +!$OMP END DO +!$OMP END PARALLEL +end subroutine impl_vert_visc_ale_vtransp +! +! +!_______________________________________________________________________________ +!SD Compute vertical integral of transport velocity rhs omitting the contributions from +!SD the elevation and Coriolis. The elevation and Coriolis are accounted for +!SD explicitly in BT equations, and should therefore be removed from the vertically integrated rhs. +subroutine compute_BT_rhs_SE_vtransp(dynamics, partit, mesh) + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_MESH + USE MOD_DYN + USE g_config, only: dt, r_restart + USE g_comm_auto + IMPLICIT NONE + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout) , target :: mesh + !___________________________________________________________________________ + real(kind=WP) :: vert_sum_u, vert_sum_v, Fx, Fy, ab1, ab2, hh + integer :: elem, nz, nzmin, nzmax, elnodes(3), ed, el(2) + logical, save :: sfirst + real(kind=WP) :: update_ubt, update_vbt, vi, len + + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:,:), pointer :: UV_rhs + real(kind=WP), dimension(:) , pointer :: eta_n + real(kind=WP), dimension(:,:) , pointer :: UVBT_4AB, UVBT_rhs +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + eta_n =>dynamics%eta_n(:) + UV_rhs =>dynamics%uv_rhs(:,:,:) + UVBT_rhs =>dynamics%se_uvBT_rhs(:,:) + UVBT_4AB =>dynamics%se_uvBT_4AB(:,:) + + !___________________________________________________________________________ + ab1=(1.5_WP+epsilon) + if(sfirst .and. (.not. r_restart)) then + ab1=1.0_WP + sfirst=.false. + end if + ab2=1.0_WP-ab1 + + !___________________________________________________________________________ +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, elnodes, nz, nzmin, nzmax, & +!$OMP Fx, Fy, hh, vert_sum_u, vert_sum_v) +!$OMP DO + do elem=1, myDim_elem2D + elnodes= elem2D_nodes(:,elem) + nzmin = ulevels(elem) + nzmax = nlevels(elem)-1 + + !_______________________________________________________________________ + ! vertically integrate UV_rhs --> for barotropic equatiobn + vert_sum_u=0.0_WP + vert_sum_v=0.0_WP + do nz=nzmin, nzmax + vert_sum_u=vert_sum_u + UV_rhs(1, nz, elem) + vert_sum_v=vert_sum_v + UV_rhs(2, nz, elem) + end do + + !_______________________________________________________________________ + ! Remove the contribution from the elevation will be accounted explicitely + ! for in the barotropic equation + Fx = g*dt*sum(gradient_sca(1:3,elem)*eta_n(elnodes)) + Fy = g*dt*sum(gradient_sca(4:6,elem)*eta_n(elnodes)) + + ! total ocean depth H + !PS hh = sum(helem(nzmin:nzmax, elem)) + !PS hh = -zbar_e_bot(elem) + sum(eta_n(elnodes))/3.0_WP + hh = -zbar_e_bot(elem) + sum(eta_n(elnodes))/3.0_WP + + vert_sum_u=vert_sum_u + Fx*hh + vert_sum_v=vert_sum_v + Fy*hh + + !_______________________________________________________________________ + ! Remove the contribution from the Coriolis will be accounted explicitely + ! for in the barotropic equation + ! UVBT_rhs ... baroclinic forcing term in barotropic equation R_b + ! --> d/dt*U_bt + f*e_z x U_bt + g*H* grad(eta) = R_bt + ! --> from AB2 in oce_ale_vel_rhs.F90 --> ab2=-0.5 (from previouse time + ! step) --> ab1=1.5 or 0.0 (first tstep) + UVBT_rhs(1, elem)=vert_sum_u - dt*mesh%coriolis(elem)* & + (ab1*UVBT_4AB(2,elem)+ab2*UVBT_4AB(4,elem)) ! for AB-interpolated + UVBT_rhs(2, elem)=vert_sum_v + dt*mesh%coriolis(elem)* & + (ab1*UVBT_4AB(1,elem)+ab2*UVBT_4AB(3,elem)) + + !_______________________________________________________________________ + ! save actual of vertical integrated transport velocity UVBT_4AB(1:2,elem) + ! timestep for next time loop step + ! UVBT_4AB ... is U_bt in barotropic equation + ! --> d/dt*U_bt + f*e_z x U_bt + g*H*grad_H(eta) = R_bt + UVBT_4AB(3:4,elem)=UVBT_4AB(1:2,elem) + + end do +!$OMP END DO +!$OMP END PARALLEL +end subroutine compute_BT_rhs_SE_vtransp +! +! +!_______________________________________________________________________________ +! Barotropic time stepping with Forward-Backward dissipative method +! (Demange et al. 2019) is used where eta and U_bt = sum_k(U_k) = sum_k(u_k*h_k) +! are estimated from the equations ... +! --> d/dt*U_bt + f*e_z x U_bt + g*H*grad_H(eta) = R_bt +! --> d/dt*eta + div_h(U_bt) + FW = 0 +! +! +! --> forward-backward dissipative time stepping by Demagne et al. 2019 +! --> equation (6) in T. Banerjee et al.,Split-Explicite external +! mode solver in FESOM2, +! +! Ubt^(n+(m+1)/M) = Ubt^(n+(m)/M) - dt/M*[ +! + 0.5*f*e_z x (Ubt^(n+(m+1)/M) + Ubt^(n+(m)/M)) +! - h*H^m*grad_H*eta^((n+m)/M) +! - Rbt-->UVBT_rhs ] +! +! eta^(n+(m+1)/M) = eta^(n+(m)/M) - dt/M * div_H * [(1+theta)*Ubt^(n+(m+1)/M) - theta*Ubt^(n+(m)/M)] +! +subroutine compute_BT_step_SE_ale(dynamics, partit, mesh) + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_MESH + USE MOD_DYN + USE g_comm_auto + USE g_config, only: dt, which_ALE, use_cavity_fw2press + USE g_support, only: integrate_nod + use o_ARRAYS, only: water_flux + IMPLICIT NONE + !___________________________________________________________________________ + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout) , target :: mesh + !___________________________________________________________________________ + real(kind=WP) :: dtBT, BT_inv, hh, ff, rx, ry, a, b, d, c1, c2, ax, ay + real(kind=WP) :: deltaX1, deltaY1, deltaX2, deltaY2, thetaBT + integer :: step, elem, edge, node, elnodes(3), ednodes(2), edelem(2), nzmax + real(kind=WP) :: update_ubt, update_vbt, vi, len + + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:) , pointer :: eta_n, bottomdrag + real(kind=WP), dimension(:,:) , pointer :: UVBT_rhs, UVBT, UVBT_theta, UVBT_mean, UVBT_12, UVBT_harmvisc + real(kind=WP), dimension(:,:,:), pointer :: UV +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + eta_n =>dynamics%eta_n(:) + UVBT_rhs =>dynamics%se_uvBT_rhs(:,:) + UVBT =>dynamics%se_uvBT(:,:) + UVBT_theta=>dynamics%se_uvBT_theta(:,:) + UVBT_mean =>dynamics%se_uvBT_mean(:,:) + UVBT_12 =>dynamics%se_uvBT_12(:,:) + if (dynamics%se_bottdrag) then + UV =>dynamics%uv(:,:,:) + bottomdrag =>dynamics%se_uvBT_stab_bdrag(:) + end if + if (dynamics%se_visc) then + UVBT_harmvisc=>dynamics%se_uvBT_stab_hvisc(:,:) + end if + + !___________________________________________________________________________ + ! Dissipation parameter of FB dissipative method 0.14 is the default value + ! from Demange et al. + thetaBT= dynamics%se_BTtheta + + ! BTsteps should be 30 or 40. + dtBT = dt/dynamics%se_BTsteps + BT_inv = 1.0_WP/(1.0_WP*dynamics%se_BTsteps) + + !___SPLIT-EXPLICITE STABILIZATION___________________________________________ + ! trim R (UVBT_rhs): + ! UVBT_rhs --> UVBT_rhs - div_h Ah H^n div_h(Ubt^n/H^n) + Cd*|Ubot|* Ubt^n/H^n + ! The intention here is to approximately remove the term thatwill be + ! added later to the barotropic momentum equation (see subroutine + ! compute_BT_step_SE_ale) + ! --> use only harmonmic viscosity operator applied to the barotropic + ! velocity + if (dynamics%se_visc) then + !_______________________________________________________________________ + ! remove viscosity +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(edge, edelem, ednodes, hh, len, & +!$OMP vi, update_ubt, update_vbt) +!$OMP DO + do edge=1, myDim_edge2D+eDim_edge2D + + ! if ed is an outer boundary edge, skip it + if(myList_edge2D(edge)>edge2D_in) cycle + + ! elem indices that participate in edge + edelem = edge_tri(:,edge) + ednodes = edges(:,edge) + + ! total ocean depth H + !PS nzmax = minval(nlevels(edelem)) + !PS hh = -zbar(nzmax) + !PS hh = minval(-zbar_e_bot(edelem)) + hh = -sum(zbar_e_bot(edelem))*0.5_WP + sum(hbar(ednodes))*0.5_WP + + len = sqrt(sum(elem_area(edelem))) + update_ubt=(UVBT(1, edelem(1))-UVBT(1, edelem(2)))/hh + update_vbt=(UVBT(2, edelem(1))-UVBT(2, edelem(2)))/hh + vi=update_ubt*update_ubt + update_vbt*update_vbt + vi=-dt*sqrt(max(dynamics%se_visc_gamma0, & + max(dynamics%se_visc_gamma1*sqrt(vi), & + dynamics%se_visc_gamma2*vi) & + )*len) + update_ubt=update_ubt*vi + update_vbt=update_vbt*vi + + !___________________________________________________________________ +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(edelem(1))) +#else +!$OMP ORDERED +#endif + UVBT_rhs(1, edelem(1))=UVBT_rhs(1, edelem(1))-update_ubt/elem_area(edelem(1))*hh + UVBT_rhs(2, edelem(1))=UVBT_rhs(2, edelem(1))-update_vbt/elem_area(edelem(1))*hh +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(edelem(1))) + call omp_set_lock (partit%plock(edelem(2))) +#endif + UVBT_rhs(1, edelem(2))=UVBT_rhs(1, edelem(2))+update_ubt/elem_area(edelem(2))*hh + UVBT_rhs(2, edelem(2))=UVBT_rhs(2, edelem(2))+update_vbt/elem_area(edelem(2))*hh +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(edelem(2))) +#else +!$OMP END ORDERED +#endif + end do ! --> do edge=1, myDim_edge2D+eDim_edge2D +!$OMP END DO +!$OMP END PARALLEL + end if ! --> if (dynamics%se_visc) then + + !___________________________________________________________________________ + ! remove bottom drag + if (dynamics%se_bottdrag) then +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, elnodes, nzmax, hh) +!$OMP DO + do elem=1, myDim_elem2D + elnodes= elem2D_nodes(:,elem) + nzmax = nlevels(elem) + + ! total ocean depth H + !PS hh = -zbar(nzmax)+sum(eta_n(elnodes))/3.0_WP + !PS hh = -zbar(nzmax) + !PS hh = -zbar_e_bot(elem) + hh = -zbar_e_bot(elem) + sum(hbar(elnodes))/3.0_WP + + bottomdrag(elem) = dt*C_d*sqrt(UV(1, nzmax-1, elem)**2 + UV(2, nzmax-1, elem)**2) + UVBT_rhs(1, elem)=UVBT_rhs(1, elem) + bottomdrag(elem)*UVBT(1, elem)/hh + UVBT_rhs(2, elem)=UVBT_rhs(2, elem) + bottomdrag(elem)*UVBT(2, elem)/hh + end do +!$OMP END DO +!$OMP END PARALLEL + end if ! --> if (dynamics%se_bottdrag) then + + !___________________________________________________________________________ + ! initialise UVBT_mean with zeros --> OMP style +!$OMP PARALLEL DO + do elem=1, myDim_elem2D+eDim_elem2D + UVBT_mean(:, elem) = 0.0_WP + end do +!$OMP END PARALLEL DO + + !___________________________________________________________________________ + ! eta_n elevation used in BT stepping, it is just a copy of eta_n + ! UBT and VBT are transport velocities + do step=1, dynamics%se_BTsteps + !####################################################################### + !########## Dissipative forward--backward time stepping ########## + !####################################################################### + + !_______________________________________________________________________ + ! compute harmonic viscosity for stability + if (dynamics%se_visc) then + + !___________________________________________________________________ + ! initialise UVBT_harmvisc with zeros --> OMP style +!$OMP PARALLEL DO + do elem=1, myDim_elem2D+eDim_elem2D + UVBT_harmvisc(:, elem) = 0.0_WP + end do +!$OMP END PARALLEL DO + + !___________________________________________________________________ +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(edge, edelem, ednodes, hh, len, & +!$OMP vi, update_ubt, update_vbt) +!$OMP DO + do edge=1, myDim_edge2D+eDim_edge2D + + ! if ed is an outer boundary edge, skip it + if(myList_edge2D(edge)>edge2D_in) cycle + + ! elem indices that participate in edge + edelem = edge_tri(:, edge) + ednodes = edges(:, edge) + + ! total ocean depth H + !PS nzmax = minval(nlevels(edelem)) + !PS hh = -zbar(nzmax) + !PS hh = minval(-zbar_e_bot(edelem)) + hh = -sum(zbar_e_bot(edelem))*0.5_WP + sum(hbar(ednodes))*0.5_WP + + len = sqrt(sum(elem_area(edelem))) + update_ubt=(UVBT(1, edelem(1))-UVBT(1, edelem(2)))/hh + update_vbt=(UVBT(2, edelem(1))-UVBT(2, edelem(2)))/hh + vi=update_ubt*update_ubt + update_vbt*update_vbt + vi=dt*sqrt(max(dynamics%se_visc_gamma0, & + max(dynamics%se_visc_gamma1*sqrt(vi), & + dynamics%se_visc_gamma2*vi) & + )*len) + update_ubt=update_ubt*vi + update_vbt=update_vbt*vi + + !_______________________________________________________________ +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(edelem(1))) +#else +!$OMP ORDERED +#endif + UVBT_harmvisc(1, edelem(1))=UVBT_harmvisc(1, edelem(1))-update_ubt/elem_area(edelem(1))*hh + UVBT_harmvisc(2, edelem(1))=UVBT_harmvisc(2, edelem(1))-update_vbt/elem_area(edelem(1))*hh +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(edelem(1))) + call omp_set_lock (partit%plock(edelem(2))) +#endif + UVBT_harmvisc(1, edelem(2))=UVBT_harmvisc(1, edelem(2))+update_ubt/elem_area(edelem(2))*hh + UVBT_harmvisc(2, edelem(2))=UVBT_harmvisc(2, edelem(2))+update_vbt/elem_area(edelem(2))*hh +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(edelem(2))) +#else +!$OMP END ORDERED +#endif + end do ! --> do edge=1, myDim_edge2D+eDim_edge2D +!$OMP END DO +!$OMP END PARALLEL + end if ! -> if (dynamics%se_visc) then + + !_______________________________________________________________________ + ! Advance velocities. I use SI stepping for the Coriolis +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, elnodes, hh, ff, rx, ry, a, b, & +!$OMP d, ax, ay) +!$OMP DO + do elem=1, myDim_elem2D + elnodes=elem2D_nodes(:,elem) + !___________________________________________________________________ + ! compute term: + ! AAA = - dt/M*[ + 0.5*f*e_z x (Ubt^(n+(m+1)/M) + Ubt^(n+(m)/M)) + ! - h*H^m*grad_H*eta^((n+m)/M) + ! - Rbt-->UVBT_rhs ] + ! total ocean depth H + !PS hh = -zbar(nlevels(elem))+sum(eta_n(elnodes))/3.0_WP ! Total fluid depth + !PS hh = -zbar(nlevels(elem)) ! Total fluid depth + !PS hh = -zbar_e_bot(elem) + hh = -zbar_e_bot(elem) + sum(hbar(elnodes))/3.0_WP + + ff = mesh%coriolis(elem) + + rx = dtBT*(-g*hh*sum(gradient_sca(1:3,elem)*eta_n(elnodes)) + ff*UVBT(2, elem)) & + + BT_inv*UVBT_rhs(1, elem) !& +!PS + BT_inv*(UVBT_harmvisc(1, elem) - bottomdrag(elem)*UVBT(1, elem)/hh) ! <-- stabilization terms + + ry = dtBT*(-g*hh*sum(gradient_sca(4:6,elem)*eta_n(elnodes)) - ff*UVBT(1, elem)) & + + BT_inv*UVBT_rhs(2, elem) !& +!PS + BT_inv*(UVBT_harmvisc(2, elem) - bottomdrag(elem)*UVBT(2, elem)/hh) ! <-- stabilization terms + if (dynamics%se_visc) then + rx = rx + BT_inv*UVBT_harmvisc(1, elem) + ry = ry + BT_inv*UVBT_harmvisc(2, elem) + end if + if (dynamics%se_bottdrag) then + rx = rx - BT_inv*bottomdrag(elem)*UVBT(1, elem)/hh + ry = ry - BT_inv*bottomdrag(elem)*UVBT(2, elem)/hh + end if + + ! compute new velocity Ubt^(n+(m+1)/M), Vbt^(n+(m+1)/M) considering + ! in terms of Increments (deltaU) and semi-Implicit Coriolis + ! (We do it here based on increments since it saves us some significant + ! digits for the accuracy) + ! Increments: + ! deltaU = Ubt^(n+(m+1)/M)-Ubt^(n+m/M) + ! Ubt^(n+(m+1)/M)+Ubt^(n+m/M) = Ubt^(n+(m+1)/M)-Ubt^(n+m/M)+2*Ubt^(n+m/M) + ! = deltaU + 2*Ubt^(n+m/M) + ! + ! Ubt^(n+(m+1)/M)-Ubt^(n+m/M) = - dt/M*[ + 0.5*f*e_z x (Ubt^(n+(m+1)/M) + Ubt^(n+m/M)) + ! - h*H^m*grad_H*eta^(n+m/M) + ! - Rbt-->UVBT_rhs ] + ! + dt/M*[grad_h*A_h*H^n*grad_h(Ubt^(n+m/M)/H^n)] + ! - dt/M*[cd*|Ubot|*Ubt^(n+(m+1)/M)/H^n] + ! + ! --> a = dt/(2*M)*ff + ! --> b = dt/M*cd*|Ubot|/H^n + ! + ! deltaU + b*deltaU - a*deltaV = dt/M*f*Vbt^(n+m/M) - h*H^m*gradx_H*eta^(n+m/M) + Rbtx + ! + dt/M*[grad_h*A_h*H^n*grad_h(Ubt^(n+m/M)/H^n)] + ! + dt/M*[cd*|Ubot|*Ubt^(n+m/M)/H^n] + ! deltaV + b*deltaV + a*deltaU = -dt/M*f*Ubt^(n+m/M) - h*H^m*grady_H*eta^(n+m/M) + Rbty + ! + dt/M*[grad_h*A_h*H^n*grad_h(Vbt^(n+m/M)/H^n)] + ! + dt/M*[cd*|Ubot|VUbt^(n+m/M)/H^n] + ! \________________________v___________________________/ + ! Rx, Ry + ! + ! | 1+b -a | * | deltaU | = MAT* | deltaU | = | Rx | + ! | a 1+b | | deltaV | = | deltaV | = | Ry | + ! --> b' = b+1 + ! --> d = 1/(b'^2 + a^2) ; inv(MAT) = d* | b' a | + ! | -a b' | + ! + ! | deltaU | = inv(MAT) * | Rx | + ! | deltaV | = | Ry | + ! + ! Ubt^(n+(m+1)/M) = Ubt^(n+m/M) + d*( b'*Rx + a *Ry) + ! Vbt^(n+(m+1)/M) = Vbt^(n+m/M) + d*(-a *Rx + b'*Ry) + ! + ! Semi-Implicit Coriolis + a = dtBT*ff*0.5_WP + if ( (dynamics%se_bdrag_si) .and. (dynamics%se_bottdrag) ) then + b = 1.0_WP+BT_inv*bottomdrag(elem)/hh + else + b = 1.0_WP + end if + d = 1.0_WP/(b*b + a*a) + ax = d*( b*rx + a*ry ) + ay = d*( -a*rx + b*ry ) + + !___________________________________________________________________ + ! compute new velocities Ubt^(n+(m+1)/M) at barotropic time step (n+(m+1)/M) ... + ! Ubt^(n+(m+1)/M) = Ubt^(n+(m)/M) + AAA + ! equation (6) in T. Banerjee et al.,Split-Explicite external + ! mode solver in FESOM2, + ! compute barotropic velocity at time step (n+(m+1)/M) + UVBT( 1, elem) = UVBT(1, elem) + ax + UVBT( 2, elem) = UVBT(2, elem) + ay + + !___________________________________________________________________ + ! velocities for dissipative time stepping of thickness equation + ! compute: [(1+theta)*Ubt^(n+(m+1)/M) - theta*Ubt^(n+(m)/M)] by ... + ! + ! --> AAA = Ubt^(n+(m+1)/M) - Ubt^(n+(m)/M) | *theta + ! AAA * theta = (Ubt^(n+(m+1)/M) - Ubt^(n+(m)/M) )*theta + ! --> Ubt^(n+(m+1)/M) + AAA*theta = Ubt^(n+(m+1)/M) + (Ubt^(n+(m+1)/M) - Ubt^(n+(m)/M) )*theta + ! (1+theta)*Ubt^(n+(m+1)/M) + theta*Ubt^(n+(m)/M) + UVBT_theta(1, elem) = UVBT(1, elem) + thetaBT*ax ! New vel*(1+thetaBT)-old vel *thetaBT + UVBT_theta(2, elem) = UVBT(2, elem) + thetaBT*ay + + !___________________________________________________________________ + ! Mean BT velocity to trim 3D velocity in tracers, equation (10) in + ! T. Banerjee et al.,Split-Explicite external mode solver in FESOM2, + UVBT_mean( 1, elem) = UVBT_mean( 1, elem) + UVBT_theta(1, elem)*BT_inv + UVBT_mean( 2, elem) = UVBT_mean( 2, elem) + UVBT_theta(2, elem)*BT_inv + + end do +!$OMP END DO +!$OMP END PARALLEL + + !_______________________________________________________________________ +!$OMP MASTER + call exchange_elem_begin(UVBT, partit) + call exchange_elem_end(partit) +!$OMP END MASTER +!$OMP BARRIER + + !_______________________________________________________________________ + ! Store mid-step velocity (to trim 3D velocities in momentum) + if(step==dynamics%se_BTsteps/2) then + UVBT_12=UVBT + end if + + !_______________________________________________________________________ + ! Advance thickness + ! compute: dt/M * div_H * [(1+theta)*Ubt^(n+(m+1)/M) - theta*Ubt^(n+(m)/M)] + ! and advance ssh --> eta^(n+(m+1)/M) = eta^(n+(m)/M) - dt/M * div_H * [...] +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(edge, ednodes, edelem, c1, c2, & +!$OMP deltaX1, deltaX2, deltaY1, deltaY2) +!$OMP DO + do edge=1, myDim_edge2D + ednodes = edges(:,edge) + edelem = edge_tri(:,edge) + + !___________________________________________________________________ + ! compute divergence div_H * [(1+theta)*Ubt^(n+(m+1)/M) - theta*Ubt^(n+(m)/M)] + deltaX1 = edge_cross_dxdy(1,edge) + deltaY1 = edge_cross_dxdy(2,edge) + c1 = UVBT_theta(2, edelem(1))*deltaX1 - UVBT_theta(1, edelem(1))*deltaY1 + c2 = 0.0_WP + if(edelem(2)>0) then + deltaX2=edge_cross_dxdy(3,edge) + deltaY2=edge_cross_dxdy(4,edge) + c2=-(UVBT_theta(2, edelem(2))*deltaX2 - UVBT_theta(1, edelem(2))*deltaY2) + end if + + !___________________________________________________________________ + ! advance ssh --> eta^(n+(m+1)/M) = eta^(n+(m)/M) - dt/M * div_H * [...] + ! equation (6) in T. Banerjee et al.,Split-Explicite external + ! mode solver in FESOM2, +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_set_lock(partit%plock(ednodes(1))) +#else +!$OMP ORDERED +#endif + eta_n(ednodes(1))=eta_n(ednodes(1)) + (c1+c2)*dtBT/areasvol(1,ednodes(1)) +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(1))) + call omp_set_lock (partit%plock(ednodes(2))) +#endif + eta_n(ednodes(2))=eta_n(ednodes(2)) - (c1+c2)*dtBT/areasvol(1,ednodes(2)) +#if defined(_OPENMP) && !defined(__openmp_reproducible) + call omp_unset_lock(partit%plock(ednodes(2))) +#else +!$OMP END ORDERED +#endif + end do +!$OMP END DO +!$OMP END PARALLEL + + !_______________________________________________________________________ + ! Apply freshwater boundary condition + if ( .not. trim(which_ALE)=='linfs') then +!$OMP PARALLEL DO + do node=1,myDim_nod2D + eta_n(node)=eta_n(node) - dtBT*water_flux(node) + end do +!$OMP END PARALLEL DO + end if + + !_______________________________________________________________________ +!$OMP MASTER + call exchange_nod(eta_n, partit) +!$OMP END MASTER +!$OMP BARRIER + + end do ! --> do step=1, dynamics%se_BTsteps + + !___________________________________________________________________________ + hbar_old = hbar + hbar = eta_n +end subroutine compute_BT_step_SE_ale +! +! +!_______________________________________________________________________________ +! Trim U and Uh to be consistent with BT transport +subroutine update_trim_vel_ale_vtransp(mode, dynamics, partit, mesh) + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_MESH + USE MOD_DYN + use g_comm_auto + IMPLICIT NONE + !___________________________________________________________________________ + integer , intent(in) :: mode + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout) , target :: mesh + integer :: elem, nz, nzmin, nzmax + real(kind=WP) :: ubar, vbar, hh_inv + real(kind=WP) :: usum(2), udiff(2) + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:,:,:), pointer :: UVh, UV, UV_rhs + real(kind=WP), dimension(:,:) , pointer :: UVBT_mean, UVBT_12 +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + UV =>dynamics%uv(:,:,:) + UVh =>dynamics%se_uvh(:,:,:) + UV_rhs =>dynamics%uv_rhs(:,:,:) + UVBT_mean =>dynamics%se_uvBT_mean(:,:) + UVBT_12 =>dynamics%se_uvBT_12(:,:) + + !___________________________________________________________________________ + ! + ! Transport version + ! + !___________________________________________________________________________ + ! Trim full velocity to ensure that its vertical sum equals to the mean BT velocity + ! The trimmed Uh,Vh are consistent with new total height defined by eta_n + if (mode==1) then +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax, ubar, vbar, hh_inv) +!$OMP DO + do elem=1, myDim_elem2D + nzmin = ulevels(elem) + nzmax = nlevels(elem)-1 + + !___________________________________________________________________ + ubar = 0.0_WP + vbar = 0.0_WP + hh_inv = 0.0_WP + do nz=nzmin, nzmax-1 + !PS !___________________________________________________________ + !PS ! Update transport velocities: U^(n+1/2,**) = U^(n+1/2,*) + U_rhs + !PS UVh(1, nz, elem)=UVh(1, nz, elem)+UV_rhs(1, nz, elem) + !PS UVh(2, nz, elem)=UVh(2, nz, elem)+UV_rhs(2, nz, elem) + !PS + !PS !___________________________________________________________ + !PS ! vertically integrate updated transport velocity: sum(k, U_k^(n+1/2,**) ) + !PS ubar = ubar+UVh(1, nz, elem) + !PS vbar = vbar+UVh(2, nz, elem) + !PS hh_inv= hh_inv+helem(nz,elem) + !_______________________________________________________________ + ! vertically integrate updated transport velocity: sum(k, U_k^(n+1/2,**) ) + ! --> the actual update of the transport velocity is done after the + ! if (dynamics%ldiag_ke) block + ubar = ubar+UVh(1, nz, elem) + UV_rhs(1, nz, elem) + vbar = vbar+UVh(2, nz, elem) + UV_rhs(2, nz, elem) + hh_inv= hh_inv+helem(nz,elem) + end do ! --> do nz=nzmin, nzmax + ! inverse total heigh over which the trimming part is distributed --> + ! here exclude bottom cell + hh_inv=1.0_WP/hh_inv + + ! sum transport over entire water column including bottom cell need + ! the trimming with the barotropic transport + nz = nzmax + ubar = ubar+UVh(1, nz, elem) + UV_rhs(1, nz, elem) + vbar = vbar+UVh(2, nz, elem) + UV_rhs(2, nz, elem) + + !___________________________________________________________________ + if (dynamics%ldiag_ke) then + do nz=nzmin, nzmax + ! U_(n+1) - U_n = Urhs_n |* (U_(n+1)+U_n) + ! U_(n+1)^2 - U_n^2 = Urhs_n * (U_(n+1)+U_n) + ! | + ! +-> U_(n+1) = U_n+Urhs_n + ! U_(n+1)^2 - U_n^2 = Urhs_n * (2*U_n + Urhs) + ! | | + ! v v + ! udiff usum + if (nz==nzmax) then + usum(1) = 2.0_WP*UVh(1,nz,elem) + UV_rhs(1, nz, elem) + usum(2) = 2.0_WP*UVh(2,nz,elem) + UV_rhs(2, nz, elem) + else + usum(1) = 2.0_WP*UVh(1,nz,elem) + UV_rhs(1, nz, elem) + (UVBT_mean(1, elem)-ubar)*helem(nz,elem)*hh_inv + usum(2) = 2.0_WP*UVh(2,nz,elem) + UV_rhs(2, nz, elem) + (UVBT_mean(2, elem)-vbar)*helem(nz,elem)*hh_inv + end if + ! transform: transport vel --> velocity + usum(1) = usum(1)/helem(nz,elem) + usum(2) = usum(2)/helem(nz,elem) + + udiff(1) = UV_rhs(1, nz, elem) + (UVBT_mean(1, elem)-ubar)*helem(nz,elem)*hh_inv + udiff(2) = UV_rhs(2, nz, elem) + (UVBT_mean(2, elem)-vbar)*helem(nz,elem)*hh_inv + ! transform: transport vel --> velocity + udiff(1) = udiff(1)/helem(nz,elem) + udiff(2) = udiff(2)/helem(nz,elem) + + ! (U_(n+1)^2 - U_n^2)/2 = usum*udiff/2 + dynamics%ke_du2( :,nz,elem) = usum*udiff/2.0_WP + + dynamics%ke_pre_xVEL( :,nz,elem) = usum*dynamics%ke_pre (:,nz,elem)/2.0_WP + dynamics%ke_adv_xVEL( :,nz,elem) = usum*dynamics%ke_adv (:,nz,elem)/2.0_WP + dynamics%ke_cor_xVEL( :,nz,elem) = usum*dynamics%ke_cor (:,nz,elem)/2.0_WP + dynamics%ke_hvis_xVEL(:,nz,elem) = usum*dynamics%ke_hvis(:,nz,elem)/2.0_WP + dynamics%ke_vvis_xVEL(:,nz,elem) = usum*dynamics%ke_vvis(:,nz,elem)/2.0_WP + + ! U_(n+0.5) = U_n + 0.5*Urhs + dynamics%ke_umean( :,nz,elem) = usum/2.0_WP + ! U_(n+0.5)^2 + dynamics%ke_u2mean( :,nz,elem) = (usum*usum)/4.0_WP + + if (nz==nzmin) then + dynamics%ke_wind_xVEL(:,elem) = usum*dynamics%ke_wind(:,elem)/2.0_WP + end if + + if (nz==nzmax) then + dynamics%ke_drag_xVEL(:,elem) = usum*dynamics%ke_drag(:,elem)/2.0_WP + end if + end do ! --> do nz=nzmin, nzmax + end if ! --> if (dynamics%ldiag_ke) then + + !___________________________________________________________________ + ! finalize horizontal transport by making vertically integrated + ! transport equal to the value obtained from the barotropic solution + ! + ! --> equation (11) in T. Banerjee et al.,Split-Explicite external + ! mode solver in FESOM2, + ! U_k^(n+1/2) = U^(n+1/2,**) + ! + [<>^(n+1/2) - sum(k, U_k^(n+1/2,**)] * h_k^(n+1/2)/sum(k,h_k^(n+1/2)) + do nz=nzmin, nzmax-1 + !PS UVh(1, nz, elem)= UVh(1, nz, elem)+(UVBT_mean(1, elem)-ubar)*helem(nz,elem)*hh_inv + !PS UVh(2, nz, elem)= UVh(2, nz, elem)+(UVBT_mean(2, elem)-vbar)*helem(nz,elem)*hh_inv + UVh(1, nz, elem)= UVh(1, nz, elem) + UV_rhs(1, nz, elem) + (UVBT_mean(1, elem)-ubar)*helem(nz,elem)*hh_inv + UVh(2, nz, elem)= UVh(2, nz, elem) + UV_rhs(2, nz, elem) + (UVBT_mean(2, elem)-vbar)*helem(nz,elem)*hh_inv + UV( 1, nz, elem)= UVh(1, nz, elem)/helem(nz,elem) ! velocities are still needed + UV( 2, nz, elem)= UVh(2, nz, elem)/helem(nz,elem) + end do ! --> do nz=nzmin, nzmax + + ! apply no trimming in the bottom cell + nz = nzmax + UVh(1, nz, elem)= UVh(1, nz, elem) + UV_rhs(1, nz, elem) + UVh(2, nz, elem)= UVh(2, nz, elem) + UV_rhs(2, nz, elem) + UV( 1, nz, elem)= UVh(1, nz, elem)/helem(nz,elem) ! velocities are still needed + UV( 2, nz, elem)= UVh(2, nz, elem)/helem(nz,elem) + + end do ! --> do elem=1, myDim_elem2D +!$OMP END DO +!$OMP END PARALLEL + + !_______________________________________________________________________ +!$OMP MASTER + call exchange_elem(UVh, partit) ! This exchange can be avoided, but test first. + call exchange_elem(UV, partit) +!$OMP END MASTER +!$OMP BARRIER + + !___________________________________________________________________________ + ! Trim to make velocity consistent with BT velocity at n+1/2 + ! Velocity will be used to advance momentum + elseif (mode==2) then +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax, ubar, vbar, hh_inv) +!$OMP DO + do elem=1, myDim_elem2D + nzmin = ulevels(elem) + nzmax = nlevels(elem)-1 + + !___________________________________________________________________ + ubar = 0.0_WP + vbar = 0.0_WP + hh_inv = 0.0_WP + do nz=nzmin, nzmax-1 + ubar=ubar+UVh(1, nz, elem) + vbar=vbar+UVh(2, nz, elem) + hh_inv=hh_inv+helem(nz,elem) + end do + nz=nzmax + ubar=ubar+UVh(1, nz, elem) + vbar=vbar+UVh(2, nz, elem) + + !___________________________________________________________________ + hh_inv=1.0_WP/hh_inv ! Postpone 2nd order, just use available thickness + do nz=nzmin, nzmax-1 + UVh(1, nz, elem)= UVh(1, nz, elem)+(UVBT_12(1, elem)-ubar)*helem(nz,elem)*hh_inv + UVh(2, nz, elem)= UVh(2, nz, elem)+(UVBT_12(2, elem)-vbar)*helem(nz,elem)*hh_inv + UV( 1, nz, elem)= UVh(1, nz, elem)/helem(nz,elem) ! velocities are still needed + UV( 2, nz, elem)= UVh(2, nz, elem)/helem(nz,elem) ! to compute momentum advection + end do + nz=nzmax + UVh(1, nz, elem)= UVh(1, nz, elem) + UVh(2, nz, elem)= UVh(2, nz, elem) + UV( 1, nz, elem)= UVh(1, nz, elem)/helem(nz,elem) ! velocities are still needed + UV( 2, nz, elem)= UVh(2, nz, elem)/helem(nz,elem) ! to compute momentum advection + end do +!$OMP END DO +!$OMP END PARALLEL + + !_______________________________________________________________________ +!$OMP MASTER + call exchange_elem(UVh, partit) ! + call exchange_elem(UV , partit) ! Check if this is needed +!$OMP END MASTER +!$OMP BARRIER + + !___________________________________________________________________________ + ! skip trimming only adnvance velocity + ! Velocity will be used to advance momentum + elseif (mode==3) then +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax, ubar, vbar, hh_inv) +!$OMP DO + do elem=1, myDim_elem2D + nzmin = ulevels(elem) + nzmax = nlevels(elem)-1 + + !___________________________________________________________________ + do nz=nzmin, nzmax + UVh(1, nz, elem)= UVh(1, nz, elem) + UV_rhs(1, nz, elem) + UVh(2, nz, elem)= UVh(2, nz, elem) + UV_rhs(2, nz, elem) + UV( 1, nz, elem)= UVh(1, nz, elem)/helem(nz,elem) ! velocities are still needed + UV( 2, nz, elem)= UVh(2, nz, elem)/helem(nz,elem) ! to compute momentum advection + end do + end do +!$OMP END DO +!$OMP END PARALLEL + + !_______________________________________________________________________ +!$OMP MASTER + call exchange_elem(UVh, partit) ! + call exchange_elem(UV , partit) ! Check if this is needed +!$OMP END MASTER +!$OMP BARRIER + + end if ! --> if (mode==1) then +end subroutine update_trim_vel_ale_vtransp +! +! +!_______________________________________________________________________________ +! Trim U and Uh to be consistent with BT transport +subroutine compute_thickness_zstar(dynamics, partit, mesh) + USE MOD_PARTIT + USE MOD_PARSUP + USE MOD_MESH + USE MOD_DYN + use g_comm_auto + implicit none + !___________________________________________________________________________ + type(t_dyn) , intent(inout), target :: dynamics + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout) , target :: mesh + integer :: node, elem, nz, nzmin, nzmax, elnodes(3) + real(kind=WP) :: hh_inv + + !___________________________________________________________________________ + ! pointer on necessary derived types + real(kind=WP), dimension(:), pointer :: eta_n +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + eta_n => dynamics%eta_n(:) + + !___________________________________________________________________________ + ! leave bottom layer again untouched +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(node, nz, nzmin, nzmax, hh_inv) +!$OMP DO + do node=1, myDim_nod2D+eDim_nod2D + nzmin = ulevels_nod2D(node) + nzmax = nlevels_nod2D_min(node)-1 + !PS hh_inv=-1.0_WP/zbar(nzmax) + hh_inv=-1.0_WP/zbar_3d_n(nzmax, node) + do nz=nzmin, nzmax-1 + hnode_new(nz,node)=(zbar(nz)-zbar(nz+1))*(1.0_WP+hh_inv*eta_n(node)) + end do + end do ! --> do node=1, myDim_nod2D+eDim_nod2D +!$OMP END DO +!$OMP END PARALLEL + +!PS !___________________________________________________________________________ +!PS ! --> update mean layer thinkness at element +!PS !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, elnodes, nz, nzmin, nzmax) +!PS !$OMP DO +!PS do elem=1, myDim_elem2D +!PS nzmin = ulevels(elem) +!PS nzmax = nlevels(elem)-1 +!PS elnodes=elem2D_nodes(:,elem) +!PS do nz=nzmin, nzmax-1 +!PS helem(nz,elem)=sum(hnode_new(nz,elnodes))/3.0_WP +!PS end do +!PS end do +!PS !$OMP END DO +!PS !$OMP END PARALLEL +!PS +!PS !___________________________________________________________________________ +!PS !$OMP MASTER +!PS call exchange_elem(helem, partit) +!PS !$OMP END MASTER +!PS !$OMP BARRIER + +end subroutine compute_thickness_zstar + + diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 8080338e9..b6daeca4a 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -529,69 +529,90 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) if (Redi) isredi=1._WP Ty =0.0_WP Ty1 =0.0_WP - + ! solve equation diffusion equation implicite part: ! --> h^(n+0.5)* (T^(n+0.5)-Tstar) = dt*( K_33*d/dz*(T^(n+0.5)-Tstar) + K_33*d/dz*Tstar ) ! --> Tnew = T^(n+0.5)-Tstar ! --> h^(n+0.5)* (Tnew) = dt*(K_33*d/dz*Tnew) + K_33*dt*d/dz*Tstar ! --> h^(n+0.5)* (Tnew) = dt*(K_33*d/dz*Tnew) + RHS ! --> solve for T_new - ! --> V_1 (Skalar Volume), A_1 (Area of edge), . - ! no Cavity A1==V1, yes Cavity A1 !=V1 /I\ nvec_up (+1) - ! I - ! ----------- zbar_1, A_1 *----I----* - ! Z_1 o T_1, V1 |\ I ./| - ! ----------- zbar_2, A_2 | \ ./ | Gaus Theorem: - ! Z_2 o T_2, V2 | \ / | --> Flux form - ! ----------- zbar_3, A_3 | | | --> normal vec outwards facing - ! Z_3 o T_3, V3 *---|-----* - ! ----------- zbar_4 \ | I ./ - ! : \ | I/ - ! \|/I - ! * I - ! \I/ - ! * nvec_dwn (-1) + ! --> As_1 (Skalar Area), A_1 (Area of edge), + ! no Cavity A_1==As_1, yes Cavity A1 !=As_1 + ! + ! ----------- zbar_1, A_1 nvec_up (+1) + ! Z_1 o T_1, As_1 ^ + ! ----------- zbar_2, A_2 _____|_____ Gaus Theorem: + ! Z_2 o T_2, As_2 /| | |\ --> Flux form + ! ----------- zbar_3, A_3 * | | | * --> normal vec outwards facing + ! Z_3 o T_3, As_3 |\___________/| + ! ----------- zbar_4 | | | | + ! : | |_________| | + ! |/ \| + ! * | * + ! \_____|_____/ + ! | + ! V nvec_dwn (-1) + ! ! --> 1st. solve homogenouse part: - ! f(Tnew) = h^(n+0.5)* (Tnew) - dt*(K_33*dTnew/dz) = 0 + ! f(Tnew) = h^(n+0.5)* (Tnew) - dt*(K_33*dTnew/dz) = 0 ! ! --> 2nd. Difference Quotient at Tnew_i in flux form (Gaus Theorem, dont forget normal vectors!!!): - ! V_i*Tnew_i *h_i = -dt * [ K_33 * (Tnew_i-1 - Tnew_i)/(Z_i-1 - Z_i) * A_i * nvec_up - ! +K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1 * nvec_dwn ] - ! Tnew_i *h_i = -dt * [ K_33 * (Tnew_i-1 - Tnew_i)/(Z_i-1 - Z_i) * A_i /V_i * nvec_up - ! +K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1/V_i * nvec_dwn ] - ! - ! --> 3rd. solve for coefficents a, b, c: - ! f(Tnew) = [ a*dTnew_i-1 + b*dTnew_i + c*dTnew_i+1 ] + ! | + ! +-> Gauss THeorem: int(V', div(F_vec))dV = intcircle(A', F_vec*n_vec)dA + ! | + ! +-> du = dt*( K_33*d/dz*dTnew ) | *div()_z=d/dz, *int(V',)dV + ! | int(V', d/dz *du)dV = int(V', d/dz *dt*( K_33*d/dz*dTnew ) )dV + ! | ... = intcircle(A', dt*( K_33*d/dz*dTnew )*n_vec)dA + ! | int(V', d/dz *du)dz*dA = ... + ! | int(A', du)dA = intcircle(A', dt*( K_33*d/dz*dTnew )*n_vec)dA + ! | + ! +-> As_i area of scalar cell = A_i, A_i+1 area of top + ! | and bottom face of scalar cell + ! V ! - ! df(Tnew)/dTnew_i-1 = a = -dt*K_33/(Z_i-1 - Z_i) * A_i/V_i * (nvec_up =1) + ! As_i*Tnew_i *h_i = dt * [ K_33 * (Tnew_i-1 - Tnew_i )/(Z_i-1 - Z_i ) * A_i * nvec_up(+1) + ! +K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1 * nvec_dwn(-1) ] + ! Tnew_i *h_i = dt * [ K_33 * (Tnew_i-1 - Tnew_i )/(Z_i-1 - Z_i ) * A_i /As_i * nvec_up(+1) + ! +K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1/As_i * nvec_dwn(-1) ] + ! | + ! +-> since we are on scalar cell As_i != A_i != A_i+1 + ! +-> take into account normal vector direction + ! V + ! f(Tnew) = Tnew_i*h_i - dt*K_33 * (Tnew_i-1 - Tnew_i )/(Z_i-1 - Z_i ) * A_i /As_i + ! + dt*K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1/As_i + ! = 0 ! - ! df(Tnew)/dTnew_i+1 = c = dt * K_33 * 1/(Z_i - Z_i+1) * A_i+1/V_i * (nvec_dwn=-1) - ! = -dt * K_33 * 1/(Z_i - Z_i+1) * A_i+1/V_i - ! - ! df(Tnew)/dTnew_i = b = h_i + dt*K_33/(Z_i-1 - Z_i) * A_i/V_i * (nvec_up=+1) - ! - dt*K_33/(Z_i - Z_i+1) * A_i+1/V_i * (nvec_dwn=-1) - ! = h_i + dt*K_33/(Z_i-1 - Z_i) * A_i/V_i - ! + dt*K_33/(Z_i - Z_i+1) * A_i+1/V_i - ! = h_i -(a+c) + ! --> 3rd. solve for coefficents a, b, c: + ! f(Tnew) = [ a*dTnew_i-1 + b*dTnew_i + c*dTnew_i+1 ] + ! | + ! +-> estimate a, b, c by derivation of f(Tnew_i) + ! | + ! +-> a = df(Tnew)/dTnew_i-1 = -dt*K_33/(Z_i-1 - Z_i) * A_i/As_i + ! | + ! +-> c = df(Tnew)/dTnew_i+1 = -dt * K_33 * 1/(Z_i - Z_i+1) * A_i+1/As_i + ! | + ! +-> b = df(Tnew)/dTnew_i = h_i + dt*K_33/(Z_i-1 - Z_i) * A_i /As_i + ! + dt*K_33/(Z_i - Z_i+1) * A_i+1/As_i + ! = h_i -(a+c) ! ! --> 4th. solve inhomogenous part: - ! [ a*dTnew_i-1 + b*dTnew_i + c*dTnew_i+1 ] = RHS/V_i + ! [ a*dTnew_i-1 + b*dTnew_i + c*dTnew_i+1 ] = RHS/As_i ! - ! RHS = K_33*dt*d/dz*Tstar + ! RHS = K_33*dt*d/dz*Tstar ! ! --> write as Difference Quotient in flux form - ! RHS/V_i = K_33 * dt * (Tstar_i-1 - Tstar_i)/(Z_i-1 - Z_i) * A_i/V_i * (nvec_up=1) - ! + K_33 * dt * (Tstar_i - Tstar_i+1)/(Z_i - Z_i+1) * A_i+1/V_i * (nvec_dwn=-1) + ! RHS/As_i = K_33 * dt * (Tstar_i-1 - Tstar_i)/(Z_i-1 - Z_i) * A_i/As_i * (nvec_up=1) + ! + K_33 * dt * (Tstar_i - Tstar_i+1)/(Z_i - Z_i+1) * A_i+1/As_i * (nvec_dwn=-1) ! - ! = K_33*dt/(Z_i-1 - Z_i) * A_i/V_i * Tstar_i-1 - ! - K_33*dt/(Z_i-1 - Z_i) * A_i/V_i * Tstar_i - ! - K_33*dt/(Z_i - Z_i+1) * A_i+1/V_i * Tstar_i - ! + K_33*dt/(Z_i - Z_i+1) * A_i+1/V_i * Tstar_i+1 + ! = K_33*dt/(Z_i-1 - Z_i) * A_i/As_i * Tstar_i-1 + ! - K_33*dt/(Z_i-1 - Z_i) * A_i/As_i * Tstar_i + ! - K_33*dt/(Z_i - Z_i+1) * A_i+1/As_i * Tstar_i + ! + K_33*dt/(Z_i - Z_i+1) * A_i+1/As_i * Tstar_i+1 ! - ! = -a*Tstar_i-1 + (a+c)*Tstar_i - c * Tstar_i+1 - ! |-> b = h_i - (a+c), a+c = h_i-b + ! = -a*Tstar_i-1 + (a+c)*Tstar_i - c * Tstar_i+1 + ! |-> b = h_i - (a+c), a+c = h_i-b ! - ! = -a*Tstar_i-1 - (b-h_i)*Tstar_i - c * Tstar_i+1 + ! = -a*Tstar_i-1 - (b-h_i)*Tstar_i - c * Tstar_i+1 ! ! --> 5th. solve for Tnew_i --> forward sweep algorithm --> see lower ! | b_1 c_1 ... | |dTnew_1| @@ -600,12 +621,14 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) ! | a_4 b_4 c_4 ...| |dTnew_4| ! | : | | : | ! - ! --> a = -dt*K_33 / (Z_i-1 - Z_i) * A_i/V_i + ! --> a = -dt*K_33 / (Z_i-1 - Z_i) * A_i/As_i ! - ! --> c = -dt*K_33 / (Z_i - Z_i+1) * A_i+1/V_i + ! --> c = -dt*K_33 / (Z_i - Z_i+1) * A_i+1/As_i + ! + ! --> b = h^(n+0.5) -[ dt*K_33/(Z_i-1 - Z_i ) * A_i /As_i + ! +dt*K_33/(Z_i - Z_i+1) * A_i+1/As_i ] + ! = -(a+c) + h^(n+0.5) ! - ! --> b = h^(n+0.5) -[ dt*K_33/(Z_i-1 - Z_i)*A_i/V_i + dt*K_33/(Z_i - Z_i+1) * A_i+1/V_i ] = -(a+c) + h^(n+0.5) - !___________________________________________________________________________ ! loop over local nodes diff --git a/src/oce_ale_vel_rhs.F90 b/src/oce_ale_vel_rhs.F90 index 164a7df4a..6d5d76a35 100644 --- a/src/oce_ale_vel_rhs.F90 +++ b/src/oce_ale_vel_rhs.F90 @@ -46,6 +46,7 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) use g_comm_auto use g_sbf, only: l_mslp use momentum_adv_scalar_interface + use momentum_adv_scalar_transpv_interface implicit none type(t_ice) , intent(inout), target :: ice type(t_dyn) , intent(inout), target :: dynamics @@ -66,6 +67,8 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) real(kind=WP), dimension(:) , pointer :: m_ice, m_snow, a_ice real(kind=WP) , pointer :: rhoice, rhosno, inv_rhowat real(kind=WP) :: ab1, ab2, ab3 !Adams-Bashforth coefficients + real(kind=WP), dimension(:,:,:), pointer :: UVh + real(kind=WP), dimension(:,:) , pointer :: UVBT_4AB #include "associate_part_def.h" #include "associate_mesh_def.h" #include "associate_part_ass.h" @@ -80,6 +83,12 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) rhosno => ice%thermo%rhosno inv_rhowat=> ice%thermo%inv_rhowat + ! if split-explicite ssh subcycling is used + if (dynamics%use_ssh_se_subcycl) then + UVh => dynamics%se_uvh + UVBT_4AB => dynamics%se_uvBT_4AB + end if + !___________________________________________________________________________ use_pice=0 if (use_floatice .and. .not. trim(which_ale)=='linfs') use_pice=1 @@ -103,34 +112,41 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) do elem=1, myDim_elem2D nzmax = nlevels(elem) nzmin = ulevels(elem) - !___________________________________________________________________________ + !_______________________________________________________________________ ! Take care of the AB part - !!PS do nz=1,nl-1 - IF (dynamics%AB_order==2) THEN - do nz=nzmin,nzmax-1 - UV_rhs(1,nz,elem)=ab1*UV_rhsAB(1,1,nz,elem) - UV_rhs(2,nz,elem)=ab1*UV_rhsAB(1,2,nz,elem) - end do - if (dynamics%ldiag_ke) then - do nz=nzmin,nzmax-1 - dynamics%ke_adv(:,nz,elem)=ab1*dynamics%ke_adv_AB(1, :,nz,elem) - dynamics%ke_cor(:,nz,elem)=ab1*dynamics%ke_cor_AB(1, :,nz,elem) - end do - end if - ELSEIF (dynamics%AB_order==3) THEN - do nz=nzmin,nzmax-1 - UV_rhs(1,nz,elem)=ab1*UV_rhsAB(2,1,nz,elem)+ab2*UV_rhsAB(1,1,nz,elem) - UV_rhs(2,nz,elem)=ab1*UV_rhsAB(2,2,nz,elem)+ab2*UV_rhsAB(1,2,nz,elem) - end do - if (dynamics%ldiag_ke) then - do nz=nzmin,nzmax-1 - dynamics%ke_adv(:,nz,elem)=ab1*dynamics%ke_adv_AB(2,:,nz,elem)+ab2*dynamics%ke_adv_AB(1,:,nz,elem) - dynamics%ke_cor(:,nz,elem)=ab1*dynamics%ke_cor_AB(2,:,nz,elem)+ab2*dynamics%ke_cor_AB(1,:,nz,elem) - end do - end if - END IF + ! 2nd order Adams-Bashfort + if (dynamics%AB_order==2) then + do nz=nzmin,nzmax-1 + UV_rhs(1,nz,elem)=ab1*UV_rhsAB(1,1,nz,elem) + UV_rhs(2,nz,elem)=ab1*UV_rhsAB(1,2,nz,elem) + ! | + ! V + ! Here Adams-Bashfort from previouse time n-0.5 + ! Thats why (0.5_WP+epsilon)*... to achieve second + ! order AB2 --> fAB2 = 3/2*fab_n - 1/2*fab_n-1 + end do + if (dynamics%ldiag_ke) then + do nz=nzmin,nzmax-1 + dynamics%ke_adv(:,nz,elem)=ab1*dynamics%ke_adv_AB(1, :,nz,elem) + dynamics%ke_cor(:,nz,elem)=ab1*dynamics%ke_cor_AB(1, :,nz,elem) + end do + end if + + ! 3rd order Adams-Bashfort + elseif (dynamics%AB_order==3) then + do nz=nzmin,nzmax-1 + UV_rhs(1,nz,elem)=ab1*UV_rhsAB(2,1,nz,elem)+ab2*UV_rhsAB(1,1,nz,elem) + UV_rhs(2,nz,elem)=ab1*UV_rhsAB(2,2,nz,elem)+ab2*UV_rhsAB(1,2,nz,elem) + end do + if (dynamics%ldiag_ke) then + do nz=nzmin,nzmax-1 + dynamics%ke_adv(:,nz,elem)=ab1*dynamics%ke_adv_AB(2,:,nz,elem)+ab2*dynamics%ke_adv_AB(1,:,nz,elem) + dynamics%ke_cor(:,nz,elem)=ab1*dynamics%ke_cor_AB(2,:,nz,elem)+ab2*dynamics%ke_cor_AB(1,:,nz,elem) + end do + end if + end if - !___________________________________________________________________________ + !_______________________________________________________________________ ! Sea level and pressure contribution -\nabla(\eta +hpressure/rho_0) ! and the Coriolis force + metric terms elnodes=elem2D_nodes(:,elem) @@ -138,10 +154,10 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) ! p_eta=g*eta_n(elnodes)*(1-theta) !! this place needs update (1-theta)!!! p_eta = g*eta_n(elnodes) - ff = mesh%coriolis(elem)*elem_area(elem) - !mm=metric_factor(elem)*gg + ff = mesh%coriolis(elem)*elem_area(elem) + !mm=metric_factor(elem)*elem_area(elem) - !___________________________________________________________________________ + !_______________________________________________________________________ ! contribution from sea level pressure if (l_mslp) then p_air = press_air(elnodes)/1000 @@ -149,7 +165,7 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) p_air = 0.0_WP end if - !___________________________________________________________________________ + !_______________________________________________________________________ ! in case of ALE zlevel and zstar add pressure from ice to atmospheric pressure ! to account for floating ice if (use_pice > 0) then @@ -161,69 +177,105 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) p_ice = 0.0_WP endif - !____________________________________________________________________________ + !_______________________________________________________________________ ! apply pressure gradient force, as well as contributions from gradient of ! the sea surface height as well as ice pressure in case of floating sea ice ! to velocity rhs - pre = -(p_eta+p_ice+p_air) - if (use_global_tides) then pre=pre-ssh_gp(elnodes) end if - - Fx = sum(gradient_sca(1:3,elem)*pre) - Fy = sum(gradient_sca(4:6,elem)*pre) - - do nz=nzmin,nzmax-1 - ! add pressure gradient terms - UV_rhs(1,nz,elem) = UV_rhs(1,nz,elem) + (Fx-pgf_x(nz,elem))*elem_area(elem) - UV_rhs(2,nz,elem) = UV_rhs(2,nz,elem) + (Fy-pgf_y(nz,elem))*elem_area(elem) - - IF (dynamics%AB_order==2) THEN - ! add coriolis force - UV_rhsAB(1,1,nz,elem) = UV(2,nz,elem)*ff! + mm*UV(1,nz,elem)*UV(2,nz,elem) - UV_rhsAB(1,2,nz,elem) =-UV(1,nz,elem)*ff! - mm*UV(1,nz,elem)*UV(2,nz,elem) - ELSEIF (dynamics%AB_order==3) THEN - UV_rhsAB(2,1,nz,elem) = UV_rhsAB(1,1,nz,elem) - UV_rhsAB(2,2,nz,elem) = UV_rhsAB(1,2,nz,elem) - UV_rhsAB(1,1,nz,elem) = UV(2,nz,elem)*ff - UV_rhsAB(1,2,nz,elem) =-UV(1,nz,elem)*ff - END IF - end do - + Fx = sum(gradient_sca(1:3, elem)*pre) + Fy = sum(gradient_sca(4:6, elem)*pre) + + !_______________________________________________________________________ + ! when ssh split-explicite subcycling method is setted use transport velocities + ! u*h, v*h instead of u,v + if (.not. dynamics%use_ssh_se_subcycl) then + do nz=nzmin, nzmax-1 + ! add pressure gradient terms + UV_rhs( 1, nz, elem) = UV_rhs(1, nz, elem) + (Fx-pgf_x(nz, elem))*elem_area(elem) + UV_rhs( 2, nz, elem) = UV_rhs(2, nz, elem) + (Fy-pgf_y(nz, elem))*elem_area(elem) + + ! add coriolis force, initialise AB2 array of actual timestep + ! with coriolis term + if (dynamics%AB_order==2) then + UV_rhsAB(1,1,nz,elem) = UV(2,nz,elem)*ff! + mm*UV(1,nz,elem)*UV(2,nz,elem) + UV_rhsAB(1,2,nz,elem) =-UV(1,nz,elem)*ff! - mm*UV(1,nz,elem)*UV(2,nz,elem) + elseif (dynamics%AB_order==3) then + UV_rhsAB(2,1,nz,elem) = UV_rhsAB(1,1,nz,elem) + UV_rhsAB(2,2,nz,elem) = UV_rhsAB(1,2,nz,elem) + UV_rhsAB(1,1,nz,elem) = UV(2,nz,elem)*ff + UV_rhsAB(1,2,nz,elem) =-UV(1,nz,elem)*ff + end if + end do + else + UVBT_4AB(1:2, elem) = 0.0_WP + do nz=nzmin, nzmax-1 + ! add pressure gradient terms + UV_rhs( 1, nz, elem) = UV_rhs(1, nz, elem) + (Fx-pgf_x(nz, elem))*elem_area(elem)*helem(nz,elem) + UV_rhs( 2, nz, elem) = UV_rhs(2, nz, elem) + (Fy-pgf_y(nz, elem))*elem_area(elem)*helem(nz,elem) + + ! add coriolis force, initialise AB2 array of actual timestep + ! with coriolis term + if (dynamics%AB_order==2) then + UV_rhsAB(1,1,nz,elem) = UVh(2,nz,elem)*ff! + mm*UV(1,nz,elem)*UV(2,nz,elem) + UV_rhsAB(1,2,nz,elem) =-UVh(1,nz,elem)*ff! - mm*UV(1,nz,elem)*UV(2,nz,elem) + elseif (dynamics%AB_order==3) then + UV_rhsAB(2,1,nz,elem) = UV_rhsAB(1,1,nz,elem) + UV_rhsAB(2,2,nz,elem) = UV_rhsAB(1,2,nz,elem) + UV_rhsAB(1,1,nz,elem) = UVh(2,nz,elem)*ff + UV_rhsAB(1,2,nz,elem) =-UVh(1,nz,elem)*ff + end if + + ! compute barotropic velocity for adams-bashfort time stepping + ! UVBT_4AB(1:2, elem)--> actual timestep, + ! UVBT_4AB(3:4, elem)--> previous timestep (is setted in + ! call compute_BC_BT_SE_vtransp) + UVBT_4AB(1, elem) = UVBT_4AB(1, elem) + UVh(1, nz, elem) ! + UVBT_4AB(2, elem) = UVBT_4AB(2, elem) + UVh(2, nz, elem) ! + end do + end if + + !_______________________________________________________________________ if (dynamics%ldiag_ke) then - do nz=nzmin,nzmax-1 - dynamics%ke_pre(1,nz,elem)= (Fx-pgf_x(nz,elem))*dt!*elem_area(elem) !not to divide it aterwards (at the end of this subroutine) - dynamics%ke_pre(2,nz,elem)= (Fy-pgf_y(nz,elem))*dt!*elem_area(elem) !but account for DT here - - IF (dynamics%AB_order==3) THEN - dynamics%ke_cor_AB(2,1,nz,elem) = dynamics%ke_cor_AB(1,1,nz,elem) - dynamics%ke_cor_AB(2,2,nz,elem) = dynamics%ke_cor_AB(1,2,nz,elem) - - dynamics%ke_adv_AB(2,1,nz,elem)= dynamics%ke_adv_AB(1,1,nz,elem) - dynamics%ke_adv_AB(2,2,nz,elem)= dynamics%ke_adv_AB(1,2,nz,elem) - END IF - - dynamics%ke_cor_AB(1,1,nz,elem)= UV(2,nz,elem)*ff - dynamics%ke_cor_AB(1,2,nz,elem)=-UV(1,nz,elem)*ff - - dynamics%ke_adv_AB(1,1,nz,elem)= 0.0_WP - dynamics%ke_adv_AB(1,2,nz,elem)= 0.0_WP - end do + do nz=nzmin,nzmax-1 + dynamics%ke_pre(1,nz,elem)= (Fx-pgf_x(nz,elem))*dt!*elem_area(elem) !not to divide it aterwards (at the end of this subroutine) + dynamics%ke_pre(2,nz,elem)= (Fy-pgf_y(nz,elem))*dt!*elem_area(elem) !but account for DT here + + if (dynamics%AB_order==3) then + dynamics%ke_cor_AB(2,1,nz,elem) = dynamics%ke_cor_AB(1,1,nz,elem) + dynamics%ke_cor_AB(2,2,nz,elem) = dynamics%ke_cor_AB(1,2,nz,elem) + + dynamics%ke_adv_AB(2,1,nz,elem)= dynamics%ke_adv_AB(1,1,nz,elem) + dynamics%ke_adv_AB(2,2,nz,elem)= dynamics%ke_adv_AB(1,2,nz,elem) + end if + + dynamics%ke_cor_AB(1,1,nz,elem)= UV(2,nz,elem)*ff + dynamics%ke_cor_AB(1,2,nz,elem)=-UV(1,nz,elem)*ff + + dynamics%ke_adv_AB(1,1,nz,elem)= 0.0_WP + dynamics%ke_adv_AB(1,2,nz,elem)= 0.0_WP + end do end if end do !$OMP END PARALLEL DO !___________________________________________________________________________ - ! advection + ! advection --> add momentum advection to actual timerstep adams-bashfort + ! array UV_rhsAB if (dynamics%momadv_opt==1) then if (mype==0) write(*,*) 'in moment not adapted mom_adv advection typ for ALE, check your namelist' call par_ex(partit%MPI_COMM_FESOM, partit%mype, 1) elseif (dynamics%momadv_opt==2) then - call momentum_adv_scalar(dynamics, partit, mesh) + if (.not. dynamics%use_ssh_se_subcycl) then + call momentum_adv_scalar(dynamics, partit, mesh) + else + call momentum_adv_scalar_transpv(dynamics, partit, mesh) + end if end if + !___________________________________________________________________________ ! Update the rhs IF (dynamics%AB_order==2) THEN @@ -238,25 +290,35 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) end if !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) do elem=1, myDim_elem2D - nzmax = nlevels(elem) nzmin = ulevels(elem) + nzmax = nlevels(elem) do nz=nzmin,nzmax-1 UV_rhs(1,nz,elem)=dt*(UV_rhs(1,nz,elem)+UV_rhsAB(1,1,nz,elem)*ff)/elem_area(elem) UV_rhs(2,nz,elem)=dt*(UV_rhs(2,nz,elem)+UV_rhsAB(1,2,nz,elem)*ff)/elem_area(elem) + ! | | + ! V V + ! fAB = (f_pgf - 1/2*fab_n-1) +3/2*fab_n + ! + ! until here: UV_rhs = dt*[ (R_advec + R_coriolis)^n + R_pressure] + ! --> horizontal viscosity contribution still missing is added in + ! call viscosity_filter + ! --> vertical viscosity contribution still missing is added in + ! call impl_vert_visc_ale end do end do !$OMP END PARALLEL DO + !___________________________________________________________________________ if (dynamics%ldiag_ke) then !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(elem, nz, nzmin, nzmax) - do elem=1, myDim_elem2D - nzmax = nlevels(elem) - nzmin = ulevels(elem) - do nz=nzmin,nzmax-1 - dynamics%ke_adv(:,nz,elem)=dt*(dynamics%ke_adv(:,nz,elem)+dynamics%ke_adv_AB(1,:,nz,elem)*ff)/elem_area(elem) - dynamics%ke_cor(:,nz,elem)=dt*(dynamics%ke_cor(:,nz,elem)+dynamics%ke_cor_AB(1,:,nz,elem)*ff)/elem_area(elem) + do elem=1, myDim_elem2D + nzmin = ulevels(elem) + nzmax = nlevels(elem) + do nz=nzmin,nzmax-1 + dynamics%ke_adv(:,nz,elem)=dt*(dynamics%ke_adv(:,nz,elem)+dynamics%ke_adv_AB(1,:,nz,elem)*ff)/elem_area(elem) + dynamics%ke_cor(:,nz,elem)=dt*(dynamics%ke_cor(:,nz,elem)+dynamics%ke_cor_AB(1,:,nz,elem)*ff)/elem_area(elem) + end do end do - end do !$OMP END PARALLEL DO end if @@ -266,9 +328,10 @@ subroutine compute_vel_rhs(ice, dynamics, partit, mesh) END SUBROUTINE compute_vel_rhs ! +! +!_______________________________________________________________________________ ! Momentum advection on scalar control volumes with ALE adaption--> exchange zinv(nz) ! against hnode(nz,node) -!_______________________________________________________________________________ subroutine momentum_adv_scalar(dynamics, partit, mesh) USE MOD_MESH USE MOD_PARTIT @@ -396,7 +459,7 @@ subroutine momentum_adv_scalar(dynamics, partit, mesh) ul2 = ulevels(el2) un2(ul2:nl2) = - UV(2,ul2:nl2,el2)*edge_cross_dxdy(3,ed) & - + UV(1,ul2:nl2,el2)*edge_cross_dxdy(4,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, diff --git a/src/oce_dyn.F90 b/src/oce_dyn.F90 index 01cfad89b..d2a27243d 100755 --- a/src/oce_dyn.F90 +++ b/src/oce_dyn.F90 @@ -111,32 +111,46 @@ SUBROUTINE update_vel(dynamics, partit, mesh) nzmax = nlevels(elem) if (dynamics%ldiag_ke) then - DO nz=nzmin, nzmax-1 - dynamics%ke_pre(1,nz,elem) =dynamics%ke_pre(1,nz,elem) + Fx - dynamics%ke_pre(2,nz,elem) =dynamics%ke_pre(2,nz,elem) + Fy - - usum(1)=2.0_WP*UV(1,nz,elem)+(UV_rhs(1,nz,elem) + Fx) - usum(2)=2.0_WP*UV(2,nz,elem)+(UV_rhs(2,nz,elem) + Fy) - - udiff(1)=UV_rhs(1,nz,elem) + Fx - udiff(2)=UV_rhs(2,nz,elem) + Fy - dynamics%ke_du2 (:,nz,elem) = usum*udiff/2.0_WP - dynamics%ke_pre_xVEL (:,nz,elem) = usum*dynamics%ke_pre (:,nz,elem)/2.0_WP - dynamics%ke_adv_xVEL (:,nz,elem) = usum*dynamics%ke_adv (:,nz,elem)/2.0_WP - dynamics%ke_cor_xVEL (:,nz,elem) = usum*dynamics%ke_cor (:,nz,elem)/2.0_WP - dynamics%ke_hvis_xVEL(:,nz,elem) = usum*dynamics%ke_hvis(:,nz,elem)/2.0_WP - dynamics%ke_vvis_xVEL(:,nz,elem) = usum*dynamics%ke_vvis(:,nz,elem)/2.0_WP - dynamics%ke_umean(:,nz,elem) = usum/2.0_WP - dynamics%ke_u2mean(:,nz,elem) = (usum*usum)/4.0_WP - - if (nz==nzmin) then - dynamics%ke_wind_xVEL(:,elem)=usum*dynamics%ke_wind(:,elem)/2.0_WP - end if - - if (nz==nzmax-1) then - dynamics%ke_drag_xVEL(:,elem)=usum*dynamics%ke_drag(:,elem)/2.0_WP - end if - END DO + DO nz=nzmin, nzmax-1 + dynamics%ke_pre(1,nz,elem) =dynamics%ke_pre(1,nz,elem) + Fx + dynamics%ke_pre(2,nz,elem) =dynamics%ke_pre(2,nz,elem) + Fy + + + ! U_(n+1) - U_n = Urhs_n |* (U_(n+1)+U_n) + ! U_(n+1)^2 - U_n^2 = Urhs_n * (U_(n+1)+U_n) + ! | + ! +-> U_(n+1) = U_n+Urhs_n + ! U_(n+1)^2 - U_n^2 = Urhs_n * (2*U_n + Urhs) + ! | | + ! v v + ! udiff usum + usum(1) = 2.0_WP*UV(1,nz,elem)+(UV_rhs(1,nz,elem) + Fx) + usum(2) = 2.0_WP*UV(2,nz,elem)+(UV_rhs(2,nz,elem) + Fy) + udiff(1) = UV_rhs(1,nz,elem) + Fx + udiff(2) = UV_rhs(2,nz,elem) + Fy + + ! (U_(n+1)^2 - U_n^2)/2 = usum*udiff/2 + dynamics%ke_du2 (:,nz,elem) = usum*udiff/2.0_WP + + dynamics%ke_pre_xVEL (:,nz,elem) = usum*dynamics%ke_pre (:,nz,elem)/2.0_WP + dynamics%ke_adv_xVEL (:,nz,elem) = usum*dynamics%ke_adv (:,nz,elem)/2.0_WP + dynamics%ke_cor_xVEL (:,nz,elem) = usum*dynamics%ke_cor (:,nz,elem)/2.0_WP + dynamics%ke_hvis_xVEL(:,nz,elem) = usum*dynamics%ke_hvis(:,nz,elem)/2.0_WP + dynamics%ke_vvis_xVEL(:,nz,elem) = usum*dynamics%ke_vvis(:,nz,elem)/2.0_WP + + ! U_(n+0.5) = U_n + 0.5*Urhs + dynamics%ke_umean( :,nz,elem) = usum/2.0_WP + ! U_(n+0.5)^2 + dynamics%ke_u2mean( :,nz,elem) = (usum*usum)/4.0_WP + + if (nz==nzmin) then + dynamics%ke_wind_xVEL(:,elem) = usum*dynamics%ke_wind(:,elem)/2.0_WP + end if + + if (nz==nzmax-1) then + dynamics%ke_drag_xVEL(:,elem) = usum*dynamics%ke_drag(:,elem)/2.0_WP + end if + END DO end if DO nz=nzmin, nzmax-1 @@ -308,6 +322,7 @@ SUBROUTINE visc_filt_bcksct(dynamics, partit, mesh) END DO !$OMP END PARALLEL DO + !___________________________________________________________________________ !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(u1, v1, len, vi, nz, ed, el, nelem, k, elem, nzmin, nzmax, update_u, update_v) !$OMP DO DO ed=1, myDim_edge2D+eDim_edge2D @@ -346,14 +361,16 @@ SUBROUTINE visc_filt_bcksct(dynamics, partit, mesh) #endif END DO !$OMP END DO + + !___________________________________________________________________________ !$OMP MASTER call exchange_elem(U_b, partit) call exchange_elem(V_b, partit) !$OMP END MASTER !$OMP BARRIER - ! =========== - ! Compute smoothed viscous term: - ! =========== + + !___________________________________________________________________________ + ! Compute smoothed viscous term, smoth from elements to nodes !$OMP DO DO ed=1, myDim_nod2D nzmin = ulevels_nod2D(ed) @@ -373,19 +390,37 @@ SUBROUTINE visc_filt_bcksct(dynamics, partit, mesh) END DO END DO !$OMP END DO + + !___________________________________________________________________________ !$OMP MASTER call exchange_nod(U_c, V_c, partit) !$OMP END MASTER !$OMP BARRIER + + !___________________________________________________________________________ + ! smooth from nodes back to elements, take unsmooth viscos term minus + ! smoothed viscos term multiplied backscatter parameter !$OMP DO do ed=1, myDim_elem2D nelem=elem2D_nodes(:,ed) nzmin = ulevels(ed) nzmax = nlevels(ed) - Do nz=nzmin, nzmax-1 - UV_rhs(1,nz,ed)=UV_rhs(1,nz,ed)+U_b(nz,ed) -dynamics%visc_easybsreturn*sum(U_c(nz,nelem))/3.0_WP - UV_rhs(2,nz,ed)=UV_rhs(2,nz,ed)+V_b(nz,ed) -dynamics%visc_easybsreturn*sum(V_c(nz,nelem))/3.0_WP - END DO + !_______________________________________________________________________ + if (dynamics%use_ssh_se_subcycl) then + !SD Approximate update for transports. We do not care about accuracy + !SD here. --> of course, helem will be better. + Do nz=nzmin, nzmax-1 + !PS UV_rhs(1,nz,ed)=UV_rhs(1,nz,ed)+(U_b(nz,ed) -dynamics%visc_easybsreturn*sum(U_c(nz,nelem))/3.0_WP)*(zbar(nz)-zbar(nz+1)) + !PS UV_rhs(2,nz,ed)=UV_rhs(2,nz,ed)+(V_b(nz,ed) -dynamics%visc_easybsreturn*sum(V_c(nz,nelem))/3.0_WP)*(zbar(nz)-zbar(nz+1)) + UV_rhs(1,nz,ed)=UV_rhs(1,nz,ed)+(U_b(nz,ed) -dynamics%visc_easybsreturn*sum(U_c(nz,nelem))/3.0_WP)*helem(nz, ed) + UV_rhs(2,nz,ed)=UV_rhs(2,nz,ed)+(V_b(nz,ed) -dynamics%visc_easybsreturn*sum(V_c(nz,nelem))/3.0_WP)*helem(nz, ed) + END DO + else + Do nz=nzmin, nzmax-1 + UV_rhs(1,nz,ed)=UV_rhs(1,nz,ed)+U_b(nz,ed) -dynamics%visc_easybsreturn*sum(U_c(nz,nelem))/3.0_WP + UV_rhs(2,nz,ed)=UV_rhs(2,nz,ed)+V_b(nz,ed) -dynamics%visc_easybsreturn*sum(V_c(nz,nelem))/3.0_WP + END DO + end if end do !$OMP END DO !$OMP END PARALLEL @@ -436,44 +471,48 @@ SUBROUTINE visc_filt_bilapl(dynamics, partit, mesh) END DO !$OMP END PARALLEL DO + !___________________________________________________________________________ + ! Sum up velocity differences over edge with respect to elemtnal index !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(u1, v1, len, vi, ed, el, nz, nzmin, nzmax) !$OMP DO DO ed=1, myDim_edge2D+eDim_edge2D if(myList_edge2D(ed)>edge2D_in) cycle - el=edge_tri(:,ed) + el = edge_tri(:,ed) nzmin = maxval(ulevels(el)) nzmax = minval(nlevels(el)) - DO nz=nzmin,nzmax-1 + DO nz=nzmin, nzmax-1 update_u(nz)=(UV(1,nz,el(1))-UV(1,nz,el(2))) update_v(nz)=(UV(2,nz,el(1))-UV(2,nz,el(2))) END DO #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_set_lock(partit%plock(el(1))) + call omp_set_lock(partit%plock(el(1))) #else !$OMP ORDERED #endif - U_c(nzmin:nzmax-1, el(1))=U_c(nzmin:nzmax-1, el(1))-update_u(nzmin:nzmax-1) - V_c(nzmin:nzmax-1, el(1))=V_c(nzmin:nzmax-1, el(1))-update_v(nzmin:nzmax-1) + U_c(nzmin:nzmax-1, el(1))=U_c(nzmin:nzmax-1, el(1))-update_u(nzmin:nzmax-1) + V_c(nzmin:nzmax-1, el(1))=V_c(nzmin:nzmax-1, el(1))-update_v(nzmin:nzmax-1) #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_unset_lock(partit%plock(el(1))) - call omp_set_lock (partit%plock(el(2))) + call omp_unset_lock(partit%plock(el(1))) + call omp_set_lock (partit%plock(el(2))) #endif - U_c(nzmin:nzmax-1, el(2))=U_c(nzmin:nzmax-1, el(2))+update_u(nzmin:nzmax-1) - V_c(nzmin:nzmax-1, el(2))=V_c(nzmin:nzmax-1, el(2))+update_v(nzmin:nzmax-1) + U_c(nzmin:nzmax-1, el(2))=U_c(nzmin:nzmax-1, el(2))+update_u(nzmin:nzmax-1) + V_c(nzmin:nzmax-1, el(2))=V_c(nzmin:nzmax-1, el(2))+update_v(nzmin:nzmax-1) #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_unset_lock(partit%plock(el(2))) + call omp_unset_lock(partit%plock(el(2))) #else !$OMP END ORDERED #endif END DO !$OMP END DO + !___________________________________________________________________________ + ! compute viscosity on element !$OMP DO DO ed=1,myDim_elem2D - len=sqrt(elem_area(ed)) + len = sqrt(elem_area(ed)) nzmin = ulevels(ed) nzmax = nlevels(ed) - Do nz=nzmin,nzmax-1 + Do nz=nzmin, nzmax-1 ! vi has the sense of harmonic viscosity coef. because of ! division by area in the end u1=U_c(nz,ed)**2+V_c(nz,ed)**2 @@ -486,21 +525,39 @@ SUBROUTINE visc_filt_bilapl(dynamics, partit, mesh) END DO END DO !$OMP END DO + + !___________________________________________________________________________ !$OMP MASTER call exchange_elem(U_c, partit) call exchange_elem(V_c, partit) !$OMP END MASTER !$OMP BARRIER + + !___________________________________________________________________________ !$OMP DO DO ed=1, myDim_edge2D+eDim_edge2D if(myList_edge2D(ed)>edge2D_in) cycle el=edge_tri(:,ed) nzmin = maxval(ulevels(el)) nzmax = minval(nlevels(el)) - DO nz=nzmin,nzmax-1 - update_u(nz)=(U_c(nz,el(1))-U_c(nz,el(2))) - update_v(nz)=(V_c(nz,el(1))-V_c(nz,el(2))) - END DO + !_______________________________________________________________________ + if (dynamics%use_ssh_se_subcycl) then + !SD Approximate update for transports. We do not care about accuracy + !SD here. --> of course, helem will be better. + do nz=nzmin,nzmax-1 + !PS update_u(nz)=(U_c(nz,el(1))-U_c(nz,el(2)))*(zbar(nz)-zbar(nz+1)) + !PS update_v(nz)=(V_c(nz,el(1))-V_c(nz,el(2)))*(zbar(nz)-zbar(nz+1)) + update_u(nz)=(U_c(nz,el(1))-U_c(nz,el(2))) * (helem(nz, el(1))+helem(nz, el(2)))*0.5_WP + update_v(nz)=(V_c(nz,el(1))-V_c(nz,el(2))) * (helem(nz, el(1))+helem(nz, el(2)))*0.5_WP + end do + else + do nz=nzmin,nzmax-1 + update_u(nz)=(U_c(nz,el(1))-U_c(nz,el(2))) + update_v(nz)=(V_c(nz,el(1))-V_c(nz,el(2))) + end do + end if + + !_______________________________________________________________________ #if defined(_OPENMP) && ! defined(__openmp_reproducible) call omp_set_lock(partit%plock(el(1))) #else @@ -568,6 +625,8 @@ SUBROUTINE visc_filt_bidiff(dynamics, partit, mesh) V_c(:, elem) = 0.0_WP END DO !$OMP END PARALLEL DO + + !___________________________________________________________________________ !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(u1, v1, len, vi, ed, el, nz, nzmin, nzmax, update_u, update_v) !$OMP DO DO ed=1, myDim_edge2D+eDim_edge2D @@ -606,14 +665,17 @@ SUBROUTINE visc_filt_bidiff(dynamics, partit, mesh) #else !$OMP END ORDERED #endif - END DO !$OMP END DO + + !___________________________________________________________________________ !$OMP MASTER call exchange_elem(U_c, partit) call exchange_elem(V_c, partit) !$OMP END MASTER !$OMP BARRIER + + !___________________________________________________________________________ !$OMP DO DO ed=1, myDim_edge2D+eDim_edge2D if(myList_edge2D(ed)>edge2D_in) cycle @@ -621,18 +683,40 @@ SUBROUTINE visc_filt_bidiff(dynamics, partit, mesh) len=sqrt(sum(elem_area(el))) nzmin = maxval(ulevels(el)) nzmax = minval(nlevels(el)) - DO nz=nzmin,nzmax-1 - u1=(UV(1,nz,el(1))-UV(1,nz,el(2))) - v1=(UV(2,nz,el(1))-UV(2,nz,el(2))) - vi=u1*u1+v1*v1 - vi=-dt*sqrt(max(dynamics%visc_gamma0, & - max(dynamics%visc_gamma1*sqrt(vi), & - dynamics%visc_gamma2*vi) & - )*len) - ! vi=-dt*sqrt(max(dynamics%visc_gamma0, dynamics%visc_gamma1*max(sqrt(vi), dynamics%visc_gamma2*vi))*len) - update_u(nz)=vi*(U_c(nz,el(1))-U_c(nz,el(2))) - update_v(nz)=vi*(V_c(nz,el(1))-V_c(nz,el(2))) - END DO + !_______________________________________________________________________ + if (dynamics%use_ssh_se_subcycl) then + !SD Approximate update for transports. We do not care about accuracy + !SD here. --> of course, helem will be better. + do nz=nzmin,nzmax-1 + u1=(UV(1,nz,el(1))-UV(1,nz,el(2))) + v1=(UV(2,nz,el(1))-UV(2,nz,el(2))) + vi=u1*u1+v1*v1 + vi=-dt*sqrt(max(dynamics%visc_gamma0, & + max(dynamics%visc_gamma1*sqrt(vi), & + dynamics%visc_gamma2*vi) & + )*len) + ! vi=-dt*sqrt(max(dynamics%visc_gamma0, dynamics%visc_gamma1*max(sqrt(vi), dynamics%visc_gamma2*vi))*len) + !PS update_u(nz)=vi*(U_c(nz,el(1))-U_c(nz,el(2)))*(zbar(nz)-zbar(nz+1)) + !PS update_v(nz)=vi*(V_c(nz,el(1))-V_c(nz,el(2)))*(zbar(nz)-zbar(nz+1)) + update_u(nz)=vi*(U_c(nz,el(1))-U_c(nz,el(2)))*(helem(nz, el(1))+helem(nz, el(2)))*0.5_WP + update_v(nz)=vi*(V_c(nz,el(1))-V_c(nz,el(2)))*(helem(nz, el(1))+helem(nz, el(2)))*0.5_WP + end do + else + do nz=nzmin,nzmax-1 + u1=(UV(1,nz,el(1))-UV(1,nz,el(2))) + v1=(UV(2,nz,el(1))-UV(2,nz,el(2))) + vi=u1*u1+v1*v1 + vi=-dt*sqrt(max(dynamics%visc_gamma0, & + max(dynamics%visc_gamma1*sqrt(vi), & + dynamics%visc_gamma2*vi) & + )*len) + ! vi=-dt*sqrt(max(dynamics%visc_gamma0, dynamics%visc_gamma1*max(sqrt(vi), dynamics%visc_gamma2*vi))*len) + update_u(nz)=vi*(U_c(nz,el(1))-U_c(nz,el(2))) + update_v(nz)=vi*(V_c(nz,el(1))-V_c(nz,el(2))) + end do + end if + + !_______________________________________________________________________ #if defined(_OPENMP) && !defined(__openmp_reproducible) call omp_set_lock(partit%plock(el(1))) #else @@ -640,15 +724,12 @@ SUBROUTINE visc_filt_bidiff(dynamics, partit, mesh) #endif UV_rhs(1, nzmin:nzmax-1, el(1))=UV_rhs(1, nzmin:nzmax-1, el(1))-update_u(nzmin:nzmax-1)/elem_area(el(1)) UV_rhs(2, nzmin:nzmax-1, el(1))=UV_rhs(2, nzmin:nzmax-1, el(1))-update_v(nzmin:nzmax-1)/elem_area(el(1)) - #if defined(_OPENMP) && !defined(__openmp_reproducible) call omp_unset_lock(partit%plock(el(1))) call omp_set_lock (partit%plock(el(2))) #endif - UV_rhs(1, nzmin:nzmax-1, el(2))=UV_rhs(1, nzmin:nzmax-1, el(2))+update_u(nzmin:nzmax-1)/elem_area(el(2)) UV_rhs(2, nzmin:nzmax-1, el(2))=UV_rhs(2, nzmin:nzmax-1, el(2))+update_v(nzmin:nzmax-1)/elem_area(el(2)) - #if defined(_OPENMP) && !defined(__openmp_reproducible) call omp_unset_lock(partit%plock(el(2))) #else @@ -658,7 +739,9 @@ SUBROUTINE visc_filt_bidiff(dynamics, partit, mesh) !$OMP END DO !$OMP END PARALLEL end subroutine visc_filt_bidiff - +! +! +!_______________________________________________________________________________ SUBROUTINE compute_ke_wrho(dynamics, partit, mesh) USE MOD_MESH USE MOD_PARTIT @@ -701,6 +784,9 @@ SUBROUTINE compute_ke_wrho(dynamics, partit, mesh) END DO call exchange_nod(dynamics%ke_wrho, partit) END SUBROUTINE compute_ke_wrho +! +! +!_______________________________________________________________________________ ! APE generation stuff SUBROUTINE compute_apegen(dynamics, tracers, partit, mesh) USE MOD_MESH diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index c7dd8a6c9..8e4911d3d 100755 --- a/src/oce_setup_step.F90 +++ b/src/oce_setup_step.F90 @@ -477,9 +477,21 @@ SUBROUTINE dynamics_init(dynamics, partit, mesh) integer :: AB_order = 2 logical :: check_opt_visc=.true. real(kind=WP) :: wsplit_maxcfl + logical :: use_ssh_se_subcycl=.false. + integer :: se_BTsteps + real(kind=WP) :: se_BTtheta + logical :: se_visc, se_bottdrag, se_bdrag_si + real(kind=WP) :: se_visc_gamma0, se_visc_gamma1, se_visc_gamma2 + namelist /dynamics_visc / opt_visc, check_opt_visc, visc_gamma0, visc_gamma1, visc_gamma2, & use_ivertvisc, visc_easybsreturn - namelist /dynamics_general/ momadv_opt, use_freeslip, use_wsplit, wsplit_maxcfl, ldiag_KE, AB_order + + namelist /dynamics_general/ momadv_opt, use_freeslip, use_wsplit, wsplit_maxcfl, & + ldiag_KE, AB_order, & + use_ssh_se_subcycl, se_BTsteps, se_BTtheta, & + se_bottdrag, se_bdrag_si, se_visc, se_visc_gamma0, & + se_visc_gamma1, se_visc_gamma2 + !___________________________________________________________________________ ! pointer on necessary derived types #include "associate_part_def.h" @@ -503,19 +515,42 @@ SUBROUTINE dynamics_init(dynamics, partit, mesh) !___________________________________________________________________________ ! set parameters in derived type - dynamics%opt_visc = opt_visc - dynamics%check_opt_visc = check_opt_visc - dynamics%visc_gamma0 = visc_gamma0 - dynamics%visc_gamma1 = visc_gamma1 - dynamics%visc_gamma2 = visc_gamma2 - dynamics%visc_easybsreturn = visc_easybsreturn - dynamics%use_ivertvisc = use_ivertvisc - dynamics%momadv_opt = momadv_opt - dynamics%use_freeslip = use_freeslip - dynamics%use_wsplit = use_wsplit - dynamics%wsplit_maxcfl = wsplit_maxcfl - dynamics%ldiag_KE = ldiag_KE - dynamics%AB_order = AB_order + dynamics%opt_visc = opt_visc + dynamics%check_opt_visc = check_opt_visc + dynamics%visc_gamma0 = visc_gamma0 + dynamics%visc_gamma1 = visc_gamma1 + dynamics%visc_gamma2 = visc_gamma2 + dynamics%visc_easybsreturn = visc_easybsreturn + dynamics%use_ivertvisc = use_ivertvisc + dynamics%momadv_opt = momadv_opt + dynamics%use_freeslip = use_freeslip + dynamics%use_wsplit = use_wsplit + dynamics%wsplit_maxcfl = wsplit_maxcfl + dynamics%ldiag_KE = ldiag_KE + dynamics%AB_order = AB_order + dynamics%use_ssh_se_subcycl = use_ssh_se_subcycl + if (dynamics%use_ssh_se_subcycl) then + dynamics%se_BTsteps = se_BTsteps + dynamics%se_BTtheta = se_BTtheta + dynamics%se_bottdrag = se_bottdrag + dynamics%se_bdrag_si = se_bdrag_si + dynamics%se_visc = se_visc + dynamics%se_visc_gamma0 = se_visc_gamma0 + dynamics%se_visc_gamma1 = se_visc_gamma1 + dynamics%se_visc_gamma2 = se_visc_gamma2 + if (mype==0) then + write(*,*) " ___Split-Explicit barotropic subcycling_________", dynamics%se_BTsteps + write(*,*) " se_BTsteps = ", dynamics%se_BTsteps + write(*,*) " se_BTtheta = ", dynamics%se_BTtheta + write(*,*) " se_bottdrag = ", dynamics%se_bottdrag + write(*,*) " se_bdrag_si = ", dynamics%se_bdrag_si + write(*,*) " se_visc = ", dynamics%se_visc + write(*,*) " se_visc_gamma0 = ", dynamics%se_visc_gamma0 + write(*,*) " se_visc_gamma1 = ", dynamics%se_visc_gamma1 + write(*,*) " se_visc_gamma2 = ", dynamics%se_visc_gamma2 + end if + end if + !___________________________________________________________________________ ! define local vertice & elem array size elem_size=myDim_elem2D+eDim_elem2D @@ -557,14 +592,39 @@ SUBROUTINE dynamics_init(dynamics, partit, mesh) !___________________________________________________________________________ ! allocate/initialise ssh arrays in derived type allocate(dynamics%eta_n( node_size)) - allocate(dynamics%d_eta( node_size)) - allocate(dynamics%ssh_rhs( node_size)) dynamics%eta_n = 0.0_WP - dynamics%d_eta = 0.0_WP - dynamics%ssh_rhs = 0.0_WP - !!PS allocate(dynamics%ssh_rhs_old(node_size)) - !!PS dynamics%ssh_rhs_old= 0.0_WP - + if (dynamics%use_ssh_se_subcycl) then + allocate(dynamics%se_uvh( 2, nl-1, elem_size)) + allocate(dynamics%se_uvBT_rhs( 2, elem_size)) + allocate(dynamics%se_uvBT_4AB( 4, elem_size)) + allocate(dynamics%se_uvBT( 2, elem_size)) + allocate(dynamics%se_uvBT_theta(2, elem_size)) + allocate(dynamics%se_uvBT_mean( 2, elem_size)) + allocate(dynamics%se_uvBT_12( 2, elem_size)) + dynamics%se_uvh = 0.0_WP + dynamics%se_uvBT_rhs = 0.0_WP + dynamics%se_uvBT_4AB = 0.0_WP + dynamics%se_uvBT = 0.0_WP + dynamics%se_uvBT_theta = 0.0_WP + dynamics%se_uvBT_mean = 0.0_WP + dynamics%se_uvBT_12 = 0.0_WP + if (dynamics%se_visc) then + allocate(dynamics%se_uvBT_stab_hvisc(2, elem_size)) + dynamics%se_uvBT_stab_hvisc = 0.0_WP + end if + if (dynamics%se_bottdrag) then + allocate(dynamics%se_uvBT_stab_bdrag(elem_size)) + dynamics%se_uvBT_stab_bdrag = 0.0_WP + end if + else + allocate(dynamics%d_eta( node_size)) + allocate(dynamics%ssh_rhs( node_size)) + dynamics%d_eta = 0.0_WP + dynamics%ssh_rhs = 0.0_WP + !!PS allocate(dynamics%ssh_rhs_old(node_size)) + !!PS dynamics%ssh_rhs_old= 0.0_WP + end if + !___________________________________________________________________________ ! inititalise working arrays allocate(dynamics%work%uvnode_rhs(2, nl-1, node_size)) diff --git a/src/write_step_info.F90 b/src/write_step_info.F90 index cae90ae33..483983ab6 100644 --- a/src/write_step_info.F90 +++ b/src/write_step_info.F90 @@ -91,7 +91,7 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) Wvel => dynamics%w(:,:) CFL_z => dynamics%cfl_z(:,:) eta_n => dynamics%eta_n(:) - d_eta => dynamics%d_eta(:) + if ( .not. dynamics%use_ssh_se_subcycl) d_eta => dynamics%d_eta(:) m_ice => ice%data(2)%values(:) if (mod(istep,outfreq)==0) then @@ -109,7 +109,8 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) loc_deta =0. loc_dhbar =0. loc_wflux =0. - loc =0. + loc =0. + !_______________________________________________________________________ #if !defined(__openmp_reproducible) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n) REDUCTION(+:loc_eta, loc_hbar, loc_deta, loc_dhbar, loc_wflux) @@ -117,23 +118,29 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) do n=1, myDim_nod2D loc_eta = loc_eta + areasvol(ulevels_nod2D(n), n)*eta_n(n) loc_hbar = loc_hbar + areasvol(ulevels_nod2D(n), n)*hbar(n) - loc_deta = loc_deta + areasvol(ulevels_nod2D(n), n)*d_eta(n) loc_dhbar = loc_dhbar + areasvol(ulevels_nod2D(n), n)*(hbar(n)-hbar_old(n)) + if ( .not. dynamics%use_ssh_se_subcycl) then + loc_deta = loc_deta + areasvol(ulevels_nod2D(n), n)*d_eta(n) + end if loc_wflux = loc_wflux + areasvol(ulevels_nod2D(n), n)*water_flux(n) end do #if !defined(__openmp_reproducible) !$OMP END PARALLEL DO #endif + if (dynamics%use_ssh_se_subcycl) then + loc_deta=loc_dhbar + end if + !_______________________________________________________________________ call MPI_AllREDUCE(loc_eta , int_eta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) call MPI_AllREDUCE(loc_hbar , int_hbar , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_deta , int_deta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) +!PS call MPI_AllREDUCE(loc_deta , int_deta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) call MPI_AllREDUCE(loc_dhbar, int_dhbar, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) call MPI_AllREDUCE(loc_wflux, int_wflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) int_eta = int_eta /ocean_areawithcav int_hbar = int_hbar /ocean_areawithcav - int_deta = int_deta /ocean_areawithcav +!PS int_deta = int_deta /ocean_areawithcav int_dhbar= int_dhbar/ocean_areawithcav int_wflux= int_wflux/ocean_areawithcav !_______________________________________________________________________ @@ -161,7 +168,11 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) call MPI_AllREDUCE(loc , min_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(2,2,:), 1, myDim_nod2D, 'min', partit) call MPI_AllREDUCE(loc , min_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc=omp_min_max_sum1(d_eta, 1, myDim_nod2D, 'min', partit) + if ( .not. dynamics%use_ssh_se_subcycl) then + loc=omp_min_max_sum1(d_eta, 1, myDim_nod2D, 'min', partit) + else + loc=omp_min_max_sum1(hbar-hbar_old, 1, myDim_nod2D, 'min', partit) + end if call MPI_AllREDUCE(loc , min_deta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hnode(1,:), 1, myDim_nod2D, 'min', partit) call MPI_AllREDUCE(loc , min_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) @@ -193,7 +204,11 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) call MPI_AllREDUCE(loc , max_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(2,2,:), 1, myDim_nod2D, 'max', partit) call MPI_AllREDUCE(loc , max_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc=omp_min_max_sum1(d_eta, 1, myDim_nod2D, 'max', partit) + if ( .not. dynamics%use_ssh_se_subcycl) then + loc=omp_min_max_sum1(d_eta, 1, myDim_nod2D, 'max', partit) + else + loc=omp_min_max_sum1(hbar-hbar_old, 1, myDim_nod2D, 'max', partit) + end if call MPI_AllREDUCE(loc , max_deta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hnode(1, :), 1, myDim_nod2D, 'max', partit) call MPI_AllREDUCE(loc , max_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) @@ -296,6 +311,7 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) real(kind=WP), dimension(:) , pointer :: u_ice, v_ice real(kind=WP), dimension(:) , pointer :: a_ice, m_ice, m_snow real(kind=WP), dimension(:) , pointer :: a_ice_old, m_ice_old, m_snow_old + real(kind=WP), dimension(:), allocatable, target :: dhbar #include "associate_part_def.h" #include "associate_mesh_def.h" #include "associate_part_ass.h" @@ -303,10 +319,8 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) UV => dynamics%uv(:,:,:) Wvel => dynamics%w(:,:) CFL_z => dynamics%cfl_z(:,:) - ssh_rhs => dynamics%ssh_rhs(:) - ssh_rhs_old => dynamics%ssh_rhs_old(:) + eta_n => dynamics%eta_n(:) - d_eta => dynamics%d_eta(:) u_ice => ice%uice(:) v_ice => ice%vice(:) a_ice => ice%data(1)%values(:) @@ -315,13 +329,22 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) a_ice_old => ice%data(1)%values_old(:) m_ice_old => ice%data(2)%values_old(:) m_snow_old => ice%data(3)%values_old(:) + if ( .not. dynamics%use_ssh_se_subcycl) then + d_eta => dynamics%d_eta(:) + ssh_rhs => dynamics%ssh_rhs(:) + ssh_rhs_old => dynamics%ssh_rhs_old(:) + else + allocate(dhbar(myDim_nod2D+eDim_nod2D)) + dhbar = hbar-hbar_old + d_eta => dhbar + end if !___________________________________________________________________________ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n, nz) do n=1, myDim_nod2d !___________________________________________________________________ ! check ssh - if ( ((eta_n(n) /= eta_n(n)) .or. eta_n(n)<-50.0 .or. eta_n(n)>50.0 .or. (d_eta(n) /= d_eta(n)) ) ) then + if ( ((eta_n(n) /= eta_n(n)) .or. eta_n(n)<-10.0 .or. eta_n(n)>10.0 .or. (d_eta(n) /= d_eta(n)) ) ) then !$OMP CRITICAL found_blowup_loc=1 write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep @@ -338,7 +361,9 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) write(*,*) 'zbar_3d_n = ',zbar_3d_n(:,n) write(*,*) 'Z_3d_n = ',Z_3d_n(:,n) write(*,*) - write(*,*) 'ssh_rhs = ',ssh_rhs(n),', ssh_rhs_old = ',ssh_rhs_old(n) + if ( .not. dynamics%use_ssh_se_subcycl) then + write(*,*) 'ssh_rhs = ',ssh_rhs(n),', ssh_rhs_old = ',ssh_rhs_old(n) + end if write(*,*) write(*,*) 'hbar = ',hbar(n),', hbar_old = ',hbar_old(n) write(*,*) @@ -393,8 +418,10 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) write(*,*) 'd_eta(n) = ',d_eta(n) write(*,*) 'hbar = ',hbar(n) write(*,*) 'hbar_old = ',hbar_old(n) - write(*,*) 'ssh_rhs = ',ssh_rhs(n) - write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + if ( .not. dynamics%use_ssh_se_subcycl) then + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + end if write(*,*) write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) write(*,*) @@ -415,6 +442,14 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) write(*,*) 'hnode(1, n) = ',hnode(1,n) write(*,*) 'hnode(:, n) = ',hnode(:,n) write(*,*) + write(*,*) 'eta_n = ',eta_n(n) + write(*,*) 'd_eta(n) = ',d_eta(n) + write(*,*) 'hbar = ',hbar(n) + write(*,*) 'hbar_old = ',hbar_old(n) + if ( .not. dynamics%use_ssh_se_subcycl) then + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + end if write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad write(*,*) !$OMP END CRITICAL @@ -449,8 +484,10 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) write(*,*) 'd_eta(n) = ',d_eta(n) write(*,*) 'hbar = ',hbar(n) write(*,*) 'hbar_old = ',hbar_old(n) - write(*,*) 'ssh_rhs = ',ssh_rhs(n) - write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + if ( .not. dynamics%use_ssh_se_subcycl) then + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + end if write(*,*) if (use_ice) then write(*,*) 'm_ice = ',m_ice(n) @@ -499,8 +536,10 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh) write(*,*) 'd_eta(n) = ',d_eta(n) write(*,*) 'hbar = ',hbar(n) write(*,*) 'hbar_old = ',hbar_old(n) - write(*,*) 'ssh_rhs = ',ssh_rhs(n) - write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + if ( .not. dynamics%use_ssh_se_subcycl) then + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + end if write(*,*) write(*,*) 'hnode = ',hnode(:,n) write(*,*) 'hnode_new = ',hnode_new(:,n) diff --git a/work/job_albedo_chain b/work/job_albedo_chain index 1bcf0eefc..81f55c259 100755 --- a/work/job_albedo_chain +++ b/work/job_albedo_chain @@ -27,11 +27,11 @@ chain_s=1 # starting chain id # time frame of model simulation # ___COREv2___ -year_s=1948 -year_e=2009 +# year_s=1948 +# year_e=2009 # ___JRA55____ -#year_s=1958 -#year_e=2018 +year_s=1958 +year_e=2018 prescribe_rlen=0 # run length in namelist.config --> if 0 value from namelist.config is taken fedit=1 @@ -226,7 +226,7 @@ if [ $is_newsimul -eq 1 ] ; then #___BACKUP NAMELIST.* FILES INTO RESULT DIRECTORY_______________________ cp namelist.config namelist.oce namelist.ice namelist.forcing namelist.io \ - namelist.cvmix ${dname_result}/. + namelist.dyn namelist.tra namelist.cvmix ${dname_result}/. cp fesom.x ${dname_result}/. #___BACKUP SRC FILES INTO RESULT DIRECTORY______________________________