From 53c2ed3ee702bd1afbd16debb07cce5fd75ceecc Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Fri, 28 Aug 2020 18:22:02 +0200 Subject: [PATCH 1/9] Add Self-Attraction and Loading effect to surface pressure gradient (only for ocean tide scenario) --- src/oce_ale_vel_rhs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oce_ale_vel_rhs.F90 b/src/oce_ale_vel_rhs.F90 index 39585f2cd..271fc5de6 100644 --- a/src/oce_ale_vel_rhs.F90 +++ b/src/oce_ale_vel_rhs.F90 @@ -87,7 +87,7 @@ subroutine compute_vel_rhs(mesh) pre = -(p_eta+p_ice+p_air) if (use_global_tides) then - pre=pre-ssh_gp(elnodes) + pre=pre-ssh_gp(elnodes)+0.1*p_eta end if Fx = sum(gradient_sca(1:3,elem)*pre) From a2603f1d51b2a2ef2c42d3bca6df248da7807d84 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Fri, 28 Aug 2020 18:24:05 +0200 Subject: [PATCH 2/9] Add energy diagnosis for global tide work at each elem (in form of layer-integrated) --- src/gen_modules_diag.F90 | 13 ++++++++++--- src/io_meandata.F90 | 1 + 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index 7edb88a98..af750de40 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -16,7 +16,7 @@ module diagnostics private !!PS public :: ldiag_solver, lcurt_stress_surf, ldiag_energy, ldiag_dMOC, ldiag_DVD, ldiag_forc, ldiag_salt3D, ldiag_curl_vel3, diag_list, & - compute_diagnostics, rhs_diag, curl_stress_surf, curl_vel3, wrhof, rhof, & + compute_diagnostics, rhs_diag, curl_stress_surf, curl_vel3, wrhof, rhof, UV_d_Omega, & u_x_u, u_x_v, v_x_v, v_x_w, u_x_w, dudx, dudy, dvdx, dvdy, dudz, dvdz, utau_surf, utau_bott, av_dudz_sq, av_dudz, av_dvdz, stress_bott, u_surf, v_surf, u_bott, v_bott, & std_dens_min, std_dens_max, std_dens_N, std_dens, std_dens_UVDZ, std_dens_RHOZ, & compute_diag_dvd_2ndmoment_klingbeil_etal_2014, compute_diag_dvd_2ndmoment_burchard_etal_2008, compute_diag_dvd @@ -26,7 +26,7 @@ module diagnostics real(kind=WP), save, allocatable, target :: rhs_diag(:) real(kind=WP), save, allocatable, target :: curl_stress_surf(:) real(kind=WP), save, allocatable, target :: curl_vel3(:,:) - real(kind=WP), save, allocatable, target :: wrhof(:,:), rhof(:,:) + real(kind=WP), save, allocatable, target :: wrhof(:,:), rhof(:,:), UV_d_Omega(:,:) real(kind=WP), save, allocatable, target :: u_x_u(:,:), u_x_v(:,:), v_x_v(:,:), v_x_w(:,:), u_x_w(:,:) real(kind=WP), save, allocatable, target :: dudx(:,:), dudy(:,:), dvdx(:,:), dvdy(:,:), dudz(:,:), dvdz(:,:), av_dudz(:,:), av_dvdz(:,:), av_dudz_sq(:,:) real(kind=WP), save, allocatable, target :: utau_surf(:), utau_bott(:) @@ -210,13 +210,14 @@ subroutine diag_energy(mode, mesh) #include "associate_mesh.h" !===================== if (firstcall) then !allocate the stuff at the first call - allocate(wrhof(nl, myDim_nod2D), rhof(nl, myDim_nod2D)) + allocate(wrhof(nl, myDim_nod2D), rhof(nl, myDim_nod2D), UV_d_Omega(nl-1, myDim_elem2D)) allocate(u_x_u(nl-1, myDim_nod2D), u_x_v(nl-1, myDim_nod2D), v_x_v(nl-1, myDim_nod2D), v_x_w(nl-1, myDim_elem2D), u_x_w(nl-1, myDim_elem2D)) allocate(dudx(nl-1, myDim_nod2D), dudy(nl-1, myDim_nod2D), dvdx(nl-1, myDim_nod2D), dvdy(nl-1, myDim_nod2D), dudz(nl, myDim_elem2D), dvdz(nl, myDim_elem2D)) allocate(utau_surf(myDim_elem2D), utau_bott(myDim_elem2D), av_dudz_sq(nl, myDim_elem2D), av_dudz(nl, myDim_elem2D), av_dvdz(nl, myDim_elem2D)) allocate(u_surf(myDim_elem2D), v_surf(myDim_elem2D), u_bott(myDim_elem2D), v_bott(myDim_elem2D), stress_bott(2, myDim_elem2D)) rhof =0. wrhof=0. + UV_d_Omega=0. u_x_u=0. u_x_v=0. v_x_v=0. @@ -299,6 +300,12 @@ subroutine diag_energy(mode, mesh) u_x_w(nz,n)=sum(Wvel(nz, elnodes))/3.*(UV(1, iup, n)*helem(iup ,n)+UV(1, ilo, n)*helem(ilo,n))/(helem(iup ,n)+helem(ilo ,n)) v_x_w(nz,n)=sum(Wvel(nz, elnodes))/3.*(UV(2, iup, n)*helem(iup ,n)+UV(2, ilo, n)*helem(ilo,n))/(helem(iup ,n)+helem(ilo ,n)) END DO + DO nz=1, nzmax-1 + if (use_global_tides) then + UV_d_Omega(nz,n)=helem(nz,n)*(UV(1,nz,n)*sum(gradient_sca(1:3,n)*(-ssh_gp(elnodes))) & + +UV(2,nz,n)*sum(gradient_sca(4:6,n)*(-ssh_gp(elnodes))) ) + end if + END DO END DO ! this loop might be very expensive DO n=1, myDim_nod2D diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index d8999c612..47d8b3dc4 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -340,6 +340,7 @@ subroutine ini_mean_io(mesh) if (ldiag_energy) then call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'rhof', 'in-situ density at faces', 'kg/m3', rhof(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'wrhof', 'vertical velocity x density', 'kg/(s*m2)', wrhof(:,:), 1, 'm', i_real8, mesh) + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'UV_d_Omega', 'layer integrated tide work', 'm^3/s^3', UV_d_Omega(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'uu', 'u times u', 'm2/s2', u_x_u(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'uv', 'u times v', 'm2/s2', u_x_v(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'vv', 'v times v', 'm2/s2', v_x_v(:,:), 1, 'm', i_real8, mesh) From c60b4c2507ed78b582d45d02dd61088861e6e196 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Tue, 8 Sep 2020 08:54:09 +0200 Subject: [PATCH 3/9] put a switch for opening/closing global tides in namelist.oce (default: .false.) --- config/namelist.oce | 1 + src/oce_modules.F90 | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/config/namelist.oce b/config/namelist.oce index e6f6cb836..0c0ba5101 100644 --- a/config/namelist.oce +++ b/config/namelist.oce @@ -39,6 +39,7 @@ visc_sh_limit=5.0e-3 ! for KPP, max visc due to shear instability mix_scheme='KPP' ! vertical mixing scheme: KPP, PP Ricr = 0.3 ! critical bulk Richardson Number concv = 1.6 ! constant for pure convection (eqn. 23) (Large 1.5-1.6; MOM default 1.8) +use_global_tides=.false. ! add global tidal potential / &oce_tra diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index 3d2b9b393..f95e8285e 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -155,7 +155,7 @@ MODULE o_PARAM biharmonic, Abh0, scale_area, mom_adv, free_slip, i_vert_visc, w_split, w_exp_max, SPP,& Fer_GM, K_GM_max, K_GM_min, K_GM_bvref, K_GM_resscalorder, K_GM_rampmax, K_GM_rampmin, & scaling_Ferreira, scaling_Rossby, scaling_resolution, scaling_FESOM14, & - Redi, visc_sh_limit, mix_scheme, Ricr, concv, which_pgf, easy_bs_scale, easy_bs_return, visc_option + Redi, visc_sh_limit, mix_scheme, Ricr, concv, which_pgf, easy_bs_scale, easy_bs_return, visc_option, use_global_tides NAMELIST /oce_tra/ diff_sh_limit, Kv0_const, double_diffusion, K_ver, K_hor, surf_relax_T, surf_relax_S, balance_salt_water, clim_relax, & ref_sss_local, ref_sss, i_vert_diff, tracer_adv, num_tracers, tracer_ID, & From 676bab49b1ef4d2d6276abf162e3d6a200e37b92 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Mon, 28 Sep 2020 11:16:57 +0200 Subject: [PATCH 4/9] small bugs fixed --- view_pscholz/sub_fesom_data.py | 2 +- view_pscholz/sub_fesom_selectline.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/view_pscholz/sub_fesom_data.py b/view_pscholz/sub_fesom_data.py index 1f38859df..de108d2c3 100644 --- a/view_pscholz/sub_fesom_data.py +++ b/view_pscholz/sub_fesom_data.py @@ -240,7 +240,7 @@ def fesom_load_data_horiz(mesh,data,do_output=True): if data.var.find('norm')!=-1 or data.var.find('vec')!=-1: ncval2 = np.copy(ncid2.variables[aux_datavar2][:]) - if data.var.find('tuv')!=-1 or data.var.find('tuv')!=-1: + if data.var.find('tuv')!=-1 or data.var.find('suv')!=-1: ncval3 = np.copy(ncid3.variables[aux_datavar3][:]) if data.var.find('norm')!=-1: ncval3 = ncval3[:,mesh.elem0_2d_i,:] diff --git a/view_pscholz/sub_fesom_selectline.py b/view_pscholz/sub_fesom_selectline.py index de4cfb8cb..17a68135c 100644 --- a/view_pscholz/sub_fesom_selectline.py +++ b/view_pscholz/sub_fesom_selectline.py @@ -489,7 +489,7 @@ def interp_lines(self,mesh,usemidpts=True): if len(bckp_value3)!=0: data_di = np.copy(bckp_value3[:,di]) if bckp_value3.shape[0]==mesh.n2dna: - data_di[di>=mesh.nodes_2d_iz-1]=np.nan + data_di[di>=mesh.nodes_2d_izg-1]=np.nan elif bckp_value3.shape[0]==mesh.n2dea: data_di[di>=np.concatenate((mesh.elem0_2d_iz,mesh.elem0_2d_iz[mesh.pbndtri_2d_i]))-1]=np.nan if usemidpts==True: @@ -891,7 +891,7 @@ def calc_flux(self, Tref=[], Sref=35.0): # loop over sample points dz = self.zlev[1:]-self.zlev[:-1] for jj in range(0,self.value[ii].shape[0]): - self.valueflx[ii][jj,:] = self.valueflx[ii][jj,:] * -dz * rho /1000 + self.valueflx[ii][jj,:] = self.valueflx[ii][jj,:] * -dz * rho0 /1000 #_______________________________________________________________ self.unit='[10^-3 Sv]' From 29e1318b1e1d67a4314b548a922a9be80422ec10 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Tue, 17 Nov 2020 11:39:35 +0100 Subject: [PATCH 5/9] mesh focus bug --- view_pscholz/sub_fesom_mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/view_pscholz/sub_fesom_mesh.py b/view_pscholz/sub_fesom_mesh.py index d3c9b14c1..50bdde395 100644 --- a/view_pscholz/sub_fesom_mesh.py +++ b/view_pscholz/sub_fesom_mesh.py @@ -282,8 +282,8 @@ def fesom_grid_rot_r2g(self,str_mode='r2g'): self.nodes_2d_zg = self.nodes_2d_z self.nodes_2d_izg = self.nodes_2d_iz self.nodes_2d_ig = self.nodes_2d_i - if (str_mode == 'focus'): - self.fesom_remove_pbnd() + #if (str_mode == 'focus'): + #self.fesom_remove_pbnd() elif (str_mode == 'g2r'): self.nodes_2d_yr=np.degrees(np.arcsin(zg)); self.nodes_2d_xr=np.degrees(np.arctan2(yg,xg)); From 82d8f1e6c7cfca5749df51f3bf9ff74b685ffacb Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Tue, 1 Dec 2020 10:15:47 +0100 Subject: [PATCH 6/9] need repmat here --- view_pscholz/sub_fesom_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/view_pscholz/sub_fesom_mesh.py b/view_pscholz/sub_fesom_mesh.py index 50bdde395..49c6b3704 100644 --- a/view_pscholz/sub_fesom_mesh.py +++ b/view_pscholz/sub_fesom_mesh.py @@ -616,7 +616,7 @@ def fesom_interp_e2n(self,data_e): elif data_e.shape[0]==self.n2dea: area_di = np.concatenate((self.elem0_2d_area,self.elem0_2d_area[self.pbndtri_2d_i])) data_e=data_e*np.matlib.repmat(area_di,nd2,1).transpose() - data_e_area = area_di.transpose()*np.invert(data_e==0) + data_e_area = np.matlib.repmat(area_di,nd2,1).transpose()*np.invert(data_e==0) del area_di data_n=np.zeros((self.n2dna,nd2)) From c3440cbea1b8066c199b10b6dddcab00b302db2c Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Mon, 25 Jan 2021 23:31:03 +0100 Subject: [PATCH 7/9] prepare for pull --- src/gen_modules_diag.F90 | 4 ++-- src/io_meandata.F90 | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index d1eebe49f..eb1217e24 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -336,8 +336,8 @@ subroutine diag_energy(mode, mesh) DO nz=1, nzmax-1 if (use_global_tides) then - UV_d_Omega(nz,n)=helem(nz,n)*(UV(1,nz,n)*sum(gradient_sca(1:3,n)*(-ssh_gp(elnodes))) & - +UV(2,nz,n)*sum(gradient_sca(4:6,n)*(-ssh_gp(elnodes))) ) + UV_d_Omega(nz,n) = UV(1,nz,n)*sum(gradient_sca(1:3,n)*(-ssh_gp(elnodes)+0.1*9.81_WP*eta_n(elnodes))) & + + UV(2,nz,n)*sum(gradient_sca(4:6,n)*(-ssh_gp(elnodes)+0.1*9.81_WP*eta_n(elnodes))) end if END DO END DO ! --> DO n=1, myDim_elem2D diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index e05e6b401..f22b8eb51 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -367,7 +367,6 @@ subroutine ini_mean_io(mesh) if (ldiag_energy) then call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'rhof', 'in-situ density at faces', 'kg/m3', rhof(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'wrhof', 'vertical velocity x density', 'kg/(s*m2)', wrhof(:,:), 1, 'm', i_real8, mesh) - call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'UV_d_Omega', 'layer integrated tide work', 'm^3/s^3', UV_d_Omega(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'uu', 'u times u', 'm2/s2', u_x_u(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'uv', 'u times v', 'm2/s2', u_x_v(:,:), 1, 'm', i_real8, mesh) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'vv', 'v times v', 'm2/s2', v_x_v(:,:), 1, 'm', i_real8, mesh) @@ -400,6 +399,8 @@ subroutine ini_mean_io(mesh) call def_stream(elem2D, myDim_elem2D, 'ty_bot', 'bottom stress y', 'N/m2', stress_bott(2, 1:myDim_elem2D),1, 'm', i_real4, mesh) if (sel_forcvar(11)==0) call def_stream(elem2D, myDim_elem2D, 'tx_sur', 'zonal wind stress to ocean', 'm/s2', stress_surf(1, 1:myDim_elem2D),1, 'm', i_real4, mesh) ; sel_forcvar(11)=1 if (sel_forcvar(12)==0) call def_stream(elem2D, myDim_elem2D, 'ty_sur', 'meridional wind stress to ocean','m/s2', stress_surf(2, 1:myDim_elem2D),1, 'm', i_real4, mesh) ; sel_forcvar(12)=1 + + call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'UV_d_Omega', 'tide power', 'm^2/s^3', UV_d_Omega(:,:), 1, 'm', i_real8, mesh) end if From 8e1fcadaf8187fdc8ed16effd6e653a4d9bb1952 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Wed, 3 Mar 2021 20:01:05 +0100 Subject: [PATCH 8/9] delete the contribution of SAL --- src/gen_modules_diag.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index eb1217e24..1285b211e 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -336,8 +336,8 @@ subroutine diag_energy(mode, mesh) DO nz=1, nzmax-1 if (use_global_tides) then - UV_d_Omega(nz,n) = UV(1,nz,n)*sum(gradient_sca(1:3,n)*(-ssh_gp(elnodes)+0.1*9.81_WP*eta_n(elnodes))) & - + UV(2,nz,n)*sum(gradient_sca(4:6,n)*(-ssh_gp(elnodes)+0.1*9.81_WP*eta_n(elnodes))) + UV_d_Omega(nz,n) = UV(1,nz,n)*sum(gradient_sca(1:3,n)*(-ssh_gp(elnodes))) & + + UV(2,nz,n)*sum(gradient_sca(4:6,n)*(-ssh_gp(elnodes))) end if END DO END DO ! --> DO n=1, myDim_elem2D From 36e5f03d5d2e5d10bbcdc42bd9b25db5b29ff7b8 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Sun, 7 Sep 2025 19:05:36 +0200 Subject: [PATCH 9/9] Fixes #404: tide module fixed, move foreph_ini from fesom_runloop to fesom_init. --- src/fesom_module.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 734560cd5..2ed97d36e 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -285,6 +285,11 @@ subroutine fesom_init(fesom_total_nsteps) call fesom_profiler_end("ocean_setup") #endif + ! global tides + if (use_global_tides) then + call foreph_ini(yearnew, month, f%partit) + end if + ! recom setup #if defined (__recom) if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call recom_init'//achar(27)//'[0m' @@ -555,10 +560,6 @@ subroutine fesom_runloop(current_nsteps) call fesom_profiler_start("fesom_runloop_total") #endif !___MODEL TIME STEPPING LOOP________________________________________________ - if (use_global_tides) then - call foreph_ini(yearnew, month, f%partit) - end if - nstart=f%from_nstep ntotal=f%from_nstep-1+current_nsteps