From b7390f7e04f2e10cbadeb0c7e146c7485981905c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 27 Apr 2022 14:25:43 -0600 Subject: [PATCH 1/7] Makes set_u_at_v and set_v_at_u public --- src/parameterizations/vertical/MOM_set_viscosity.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index fb969953c4..367cf44d58 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -34,7 +34,7 @@ module MOM_set_visc #include public set_viscous_BBL, set_viscous_ML, set_visc_init, set_visc_end -public set_visc_register_restarts +public set_visc_register_restarts, set_u_at_v, set_v_at_u ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with From 9c103f1f9701796005e9a1fecee2d72e1a7daaa5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 27 Apr 2022 16:43:28 -0600 Subject: [PATCH 2/7] First draft for fpmix --- src/core/MOM_dynamics_split_RK2.F90 | 78 ++- .../vertical/MOM_vert_friction.F90 | 610 ++++++++++++++++++ 2 files changed, 687 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 06d828de96..8c3612f50b 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -68,6 +68,9 @@ module MOM_dynamics_split_RK2 use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units use MOM_wave_interface, only: wave_parameters_CS, Stokes_PGF +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member implicit none ; private @@ -131,6 +134,8 @@ module MOM_dynamics_split_RK2 real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure !! anomaly in each layer due to free surface height !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]. + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to ge + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean !! to the seafloor [R L Z T-2 ~> Pa] @@ -159,10 +164,12 @@ module MOM_dynamics_split_RK2 !! Euler (1) [nondim]. 0 is often used. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debug_OBC !< If true, do debugging calls for open boundary conditions. + logical :: fpmix !< If true, apply profiles of MTM flux magnitude and direction. logical :: module_is_initialized = .false. !< Record whether this module has been initialized. !>@{ Diagnostic IDs + integer :: id_uold = -1, id_vold = -1 integer :: id_uh = -1, id_vh = -1 integer :: id_umo = -1, id_vmo = -1 integer :: id_umo_2d = -1, id_vmo_2d = -1 @@ -320,6 +327,11 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! eta_pred is the predictor value of the free surface height or column mass, ! [H ~> m or kg m-2]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold + ! uold and vold are the velocities before vert_visc is applied. These arrays + ! are only used if fpmix is enabled [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_old_rad_OBC real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_OBC ! u_old_rad_OBC and v_old_rad_OBC are the starting velocities, which are @@ -348,8 +360,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1]. h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth from Cvmix real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. - + logical :: LU_pred ! Controls if it is predictor step or not logical :: dyn_p_surf logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the ! relative weightings of the layers in calculating @@ -629,10 +642,41 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%debug) then call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, scale=US%L_T_to_m_s) endif + + if (CS%fpmix) then + uold(:,:,:) = 0.0 + vold(:,:,:) = 0.0 + do k = 1, nz + do j = js , je + do I = Isq, Ieq + uold(I,j,k) = up(I,j,k) + enddo + enddo + do J = Jsq, Jeq + do i = is, ie + vold(i,J,k) = vp(i,J,k) + enddo + enddo + enddo + endif + call vertvisc_coef(up, vp, h, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, & CS%OBC) call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + + if (CS%fpmix) then + LU_pred = .true. + hbl(:,:) = 0.0 + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) + if (ASSOCIATED(CS%energetic_PBL_CSp)) & + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) + call vertFPmix(LU_pred, up, vp, uold, vold, hbl, h, forces, & + dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) + call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & + GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + endif + if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)") if (G%nonblocking_updates) then call cpu_clock_end(id_clock_vertvisc) @@ -847,9 +891,36 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! u <- u + dt d/dz visc d/dz u ! u_av <- u_av + dt d/dz visc d/dz u_av call cpu_clock_begin(id_clock_vertvisc) + + if (CS%fpmix) then + uold(:,:,:) = 0.0 + vold(:,:,:) = 0.0 + do k = 1, nz + do j = js , je + do I = Isq, Ieq + uold(I,j,k) = u(I,j,k) + enddo + enddo + do J = Jsq, Jeq + do i = is, ie + vold(i,J,k) = v(i,J,k) + enddo + enddo + enddo + endif + call vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves) + + if (CS%fpmix) then + LU_pred = .false. + call vertFPmix(LU_pred, u, v, uold, vold, hbl, h, forces, dt, & + G, GV, US, CS%vertvisc_CSp, CS%OBC) + call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & + CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + endif + if (G%nonblocking_updates) then call cpu_clock_end(id_clock_vertvisc) call start_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass) @@ -914,6 +985,11 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s enddo ; enddo enddo + if (CS%fpmix) then + if (CS%id_uold > 0) call post_data(CS%id_uold , uold, CS%diag) + if (CS%id_vold > 0) call post_data(CS%id_vold , vold, CS%diag) + endif + ! The time-averaged free surface height has already been set by the last call to btstep. ! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index d384500c3d..d5a7aa9804 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -3,6 +3,7 @@ module MOM_vert_friction ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domains, only : pass_var, To_All, Omit_corners +use MOM_domains, only : pass_vector, Scalar_Pair use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : post_product_u, post_product_sum_u use MOM_diag_mediator, only : post_product_v, post_product_sum_v @@ -31,6 +32,7 @@ module MOM_vert_friction public vertvisc, vertvisc_remnant, vertvisc_coef public vertvisc_limit_vel, vertvisc_init, vertvisc_end public updateCFLtruncationValue +public vertFPmix ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -126,6 +128,9 @@ module MOM_vert_friction integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 integer :: id_du_dt_str = -1, id_dv_dt_str = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 + integer :: id_FPmask_u = -1, id_FPmask_v = -1 , id_FPhbl_u = -1, id_FPhbl_v = -1 + integer :: id_tauFP_u = -1, id_tauFP_v = -1 , id_FPtau2x_u = -1, id_FPtau2x_v = -1 + integer :: id_FPtau2s_u = -1, id_FPtau2s_v = -1, id_FPtau2w_u = -1, id_FPtau2w_v = -1 integer :: id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1 @@ -142,6 +147,579 @@ module MOM_vert_friction contains +!> Add nonlocal momentum flux profile increments +!! TODO: add more description +subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OBC) ! FPmix + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: ui !< Zonal velocity after vertvisc [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: vi !< Meridional velocity after vertvisc [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: uold !< Old Zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: vold !< Old Meridional velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: hbl_h ! boundary layer depth + logical, intent(inout) :: LU_pred !w predictor step or NOT + type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces + real, intent(in) :: dt !< Time increment [T ~> s] + type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure + type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure + + ! local variables + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: mask3d_u !Test Plots @ 3-D centers + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: mask3d_v + real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !2-D + real, dimension(SZI_(G),SZJB_(G)) :: hbl_v + integer, dimension(SZIB_(G),SZJ_(G)) :: kbl_u + integer, dimension(SZI_(G),SZJB_(G)) :: kbl_v + real, dimension(SZI_(G),SZJ_(G)) :: ustar2_h !2-D surface + real, dimension(SZIB_(G),SZJ_(G)) :: ustar2_u + real, dimension(SZI_(G),SZJB_(G)) :: ustar2_v + real, dimension(SZIB_(G),SZJ_(G)) :: taux_u + real, dimension(SZI_(G),SZJB_(G)) :: tauy_v + real, dimension(SZIB_(G),SZJ_(G)) :: tauy_u + real, dimension(SZI_(G),SZJB_(G)) :: taux_v + real, dimension(SZIB_(G),SZJ_(G)) :: omega_w2x_u + real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v + + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tau_u !3-D interfaces + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tau_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauxDG_u + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauyDG_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauxDG_v + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauyDG_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2x_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2x_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2s_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2s_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_s2x_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_s2x_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_s2w_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_s2w_v + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: du_rot !3-D centers + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: dv_rot + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: vi_u + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: ui_v + + real :: tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, MAXinc, MINthick + real :: du, dv, du_v, dv_u , dup, dvp , uZero, vZero + real :: fEQband, Cemp_SS , Cemp_LS , Cemp_CG, Cemp_DG , Wgt_SS + real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG + real :: pi, tmp, cos_tmp, sin_tmp, depth, taux, tauy, tauk, tauxI , tauyI, sign_f + real :: tauxh, tauyh, tauh, omega_s2xh, omega_s2wh, omega_tau2xh, omega_tau2wh + real :: taux0, tauy0, tau0, sigma, G_sig, Wind_x, Wind_y, omega_w2s, omega_tau2s,omega_s2x + real :: omega_tau2x, omega_tau2w, omega_SS, omega_LS, omega_tmp, omega_s2xI, omega_s2w + integer :: kblmin, kbld, kp, km, kp1, L19 ,jNseam + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz + is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke + + pi = 4. * atan2(1.,1.) + L19 = 1 !w Options A = 1, B = 2, C = 3 + Cemp_CG = 3.6 !w L91 cross-gradient + Cemp_DG = 1.0 !w L91 down-gradient + MAXinc = -1.0 !w if positive + MINthick= 0.01 !w GV%H_subroundoff !w 0.5 + kblmin = 1 + jNseam = 457 !w north seam = SZJ_(G) + +if(LU_pred ) then !w predictor step only, surface forcing + ustar2_h(:,:) = 0. + do j = js,je !w ?GMM -1,+1 with forces% + do i = is,ie + ustar2_h(i,j) = forces%ustar(i,j) * forces%ustar(i,j) + !w omega_w2x_h(i,j) = forces%omega_w2x(i,j) + enddo + enddo + call pass_var(ustar2_h ,G%Domain) ! update halos ?GMM + call pass_var( hbl_h ,G%Domain) + +! SURFACE ustar2 and x-stress to u points and ustar2 and y-stress to v points + ustar2_u(:,:) = 0. + ustar2_v(:,:) = 0. + hbl_u(:,:) = 0. + hbl_v(:,:) = 0. + taux_u(:,:) = 0. + tauy_v(:,:) = 0. + do j = js,je + do I = Isq,Ieq + tmp = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) + ustar2_u(I,j)=(G%mask2dT(i,j)*ustar2_h(i,j)+G%mask2dT(i+1,j)*ustar2_h(i+1,j))/tmp + hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j)* hbl_h(i+1,j)) /tmp + taux_u(I,j) = forces%taux(I,j) / GV%H_to_RZ + enddo + enddo + do J = Jsq,Jeq + do i = is,ie + tmp = MAX ( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1) ) ) + ustar2_u(I,j)=(G%mask2dT(i,j)*ustar2_h(i,j)+G%mask2dT(i+1,j)*ustar2_h(i+1,j))/tmp + hbl_v(i,J) = (G%mask2dT(i,j)* hbl_h(i,J) + G%mask2dT(i,j+1)* hbl_h(i,j+1)) /tmp + if( j > jNseam-1 ) then + ustar2_v(i,J) = ustar2_h(i,j ) !w ( j > 456 ) j >= 457 + hbl_v(i,J) = hbl_h(i,j) + endif + tauy_v(i,J) = forces%tauy(i,J) / GV%H_to_RZ + enddo + enddo + call pass_vector(taux_u , tauy_v, G%Domain, To_All+Scalar_Pair) + if (CS%debug) then + call uvchksum("ustar2 ",ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.) + call uvchksum(" hbl ", hbl_u , hbl_v , G%HI, haloshift=0, scalar_pair=.true.) + call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=0, scalar_pair=.true.) + endif +!W endif !w predictor step + +!w surface tauy_u , taux_v and omega_w2x_[u,v] & Implicit interface stresses tauxDG_u and tauyDG_v + tauy_u(:,:) = 0.0 + taux_v(:,:) = 0.0 + kbl_u(:,:) = 0 + kbl_v(:,:) = 0 + omega_w2x_u(:,:) = 0.0 + omega_w2x_v(:,:) = 0.0 + tauxDG_u(:,:,:) = 0.0 + tauyDG_v(:,:,:) = 0.0 + do j = js,je + do I = Isq,Ieq + if( (G%mask2dCu(I,j) > 0.5) ) then + tauy = 0.0 + tmp = MAX(1.0, (G%mask2dCv(i,j) + G%mask2dCv(i,j-1) + G%mask2dCv(i+1,j) + G%mask2dCv(i+1,j-1) ) ) + if ( G%mask2dCv(i ,j ) > 0.5 ) tauy = tauy + tauy_v(i ,j ) + if ( G%mask2dCv(i ,j-1) > 0.5 ) tauy = tauy + tauy_v(i ,j-1) + if ( G%mask2dCv(i+1,j ) > 0.5 ) tauy = tauy + tauy_v(i+1,j ) + if ( G%mask2dCv(i+1,j-1) > 0.5 ) tauy = tauy + tauy_v(i+1,j-1) + tauy = tauy / tmp + tauy_u(I,j) = (tauy/(abs(tauy)+GV%H_subroundoff)) * sqrt(MAX(GV%H_subroundoff,ustar2_u(I,j)*ustar2_u(I,j)-taux_u(I,j)*taux_u(I,j) )) + omega_w2x_u(I,j) = atan2( tauy_u(I,j) , taux_u(I,j) ) + tauxDG_u(I,j,1) = taux_u(I,j) !w ustar2_u(I,j) * cos(omega_w2x_u(I,j)) + depth = 0.0 + do k = 1, nz + depth = depth + CS%h_u(I,j,k) + if( (depth .ge. hbl_u(I,j)) .and. (kbl_u(I,j) .eq. 0 ) .and. (k > (kblmin-1)) ) then + kbl_u(I,j) = k + hbl_u(I,j) = depth + endif + enddo + endif + enddo + enddo + do J = Jsq,Jeq + do i = is,ie + if( (G%mask2dCv(i,J) > 0.5) ) then + taux = 0.0 + if ( j < 457 ) then + tmp = MAX(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i,j+1) + G%mask2dCu(i-1,j) + G%mask2dCu(i-1,j+1) ) ) + if ( G%mask2dCu(i ,j ) > 0.5 ) taux = taux + taux_u(i ,j ) + if ( G%mask2dCu(i ,j+1) > 0.5 ) taux = taux + taux_u(i ,j+1) + if ( G%mask2dCu(i-1,j ) > 0.5 ) taux = taux + taux_u(i-1,j ) + if ( G%mask2dCu(i-1,j+1) > 0.5 ) taux = taux + taux_u(i-1,j+1) + else + tmp = MAX(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i-1,j) ) ) + if ( G%mask2dCu(i ,j ) > 0.5 ) taux = taux + taux_u(i ,j ) + if ( G%mask2dCu(i-1,j ) > 0.5 ) taux = taux + taux_u(i-1,j ) + endif + taux = taux / tmp + taux_v(i,J) = (taux/(abs(taux)+GV%H_subroundoff)) * sqrt(MAX(GV%H_subroundoff,ustar2_v(i,J)*ustar2_v(i,J)-tauy_v(i,J)*tauy_v(i,J) )) + omega_w2x_v(i,J) = atan2( tauy_v(i,J) , taux_v(i,J) ) + tauyDG_v(i,J,1) = tauy_v(i,J) !w ustar2_v(i,J) * cos(omega_w2x_v(i,J)) + depth = 0.0 + do k = 1, nz + depth = depth + CS%h_v(i,J,k) + if( (depth .ge. hbl_v(i,J)) .and. (kbl_v(i,J) .eq. 0 ) .and. (k > (kblmin-1)) ) then + kbl_v(i,J) = k + hbl_v(i,J) = depth + endif + enddo + endif + enddo + enddo +endif !w predictor step + +! Thickness weighted diagnostic interpolations ! Copy Implicit [uv]i to [uv]old + call pass_vector(ui,vi, G%Domain, To_All+Scalar_Pair) + vi_u(:,:,:) = 0. + ui_v(:,:,:) = 0. + tauxDG_u(:,:,:) = 0.0 + tauyDG_v(:,:,:) = 0.0 + tauxDG_v(:,:,:) = 0. + tauyDG_u(:,:,:) = 0. + do k = 1, nz + kp = MIN( k+1 , nz) + do j = js-1 ,je+1 + do I = Isq-1, Ieq+1 + tauxDG_u(I,j,k+1) = CS%a_u(I,j,kp) * (ui(I,j,k) - ui(I,j,kp)) + enddo + enddo + do J = Jsq-1, Jeq+1 + do i = is-1, ie+1 + tauyDG_v(i,J,k+1) = CS%a_v(i,J,kp) * (vi(i,J,k) - vi(i,J,kp)) + enddo + enddo + + ! v to u points + do j = js , je + do I = Isq, Ieq + vi_u(I,j,k) = set_v_at_u(vi, h, G, GV, I, j, k, G%mask2dCv, OBC) + tauyDG_u(I,j,k)= set_v_at_u(tauyDG_v, h, G, GV, I, j, k, G%mask2dCv, OBC) + enddo + enddo + ! u to v points + do J = Jsq, Jeq + do i = is, ie + ui_v(I,j,k) = set_u_at_v(ui, h, G, GV, i, J, k, G%mask2dCu, OBC) + tauxDG_v(i,J,k)= set_u_at_v(tauxDG_u, h, G, GV, i, J, k, G%mask2dCu, OBC) + enddo + enddo + enddo + if (CS%debug) then + call uvchksum(" vi_u ui_v ", vi_u , ui_v , G%HI, haloshift=0, scalar_pair=.true.) + endif + +! compute angles, tau2x_[u,v], tau2w_[u,v], tau2s_[u,v], s2x_[u,v], s2w_[u,v] and stress mag tau_[u,v] + omega_tau2x_u(:,:,:) = 0.0 + omega_tau2x_v(:,:,:) = 0.0 + omega_tau2w_u(:,:,:) = 0.0 + omega_tau2w_v(:,:,:) = 0.0 + omega_tau2s_u(:,:,:) = 0.0 + omega_tau2s_v(:,:,:) = 0.0 + omega_s2x_u(:,:,:) = 0.0 + omega_s2x_v(:,:,:) = 0.0 + omega_s2w_u(:,:,:) = 0. + omega_s2w_v(:,:,:) = 0. + tau_u(:,:,:) = 0.0 + tau_v(:,:,:) = 0.0 + +!w Default implicit (I) stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] Profiles + do j = js,je + do I = Isq,Ieq + if( (G%mask2dCu(I,j) > 0.5) ) then + tauyDG_u(I,j,1) = tauy_u(I,j) ! SURFACE + tau_u(I,j,1) = ustar2_u(I,j) !w stress magnitude + Omega_tau2w_u(I,j,1) = 0.0 + Omega_tau2x_u(I,j,1) = omega_w2x_u(I,j) + Omega_tau2s_u(I,j,1) = 0.0 + omega_s2x_u(I,j,1) = omega_w2x_u(I,j) + omega_s2w_u(I,j,1) = 0.0 + + do k=1,nz + kp1 = MIN(k+1 , nz) + tau_u(I,j,k+1) = sqrt( tauxDG_u(I,j,k+1)*tauxDG_u(I,j,k+1) + tauyDG_u(I,j,k+1)*tauyDG_u(I,j,k+1)) + Omega_tau2x_u(I,j,k+1) = atan2( tauyDG_u(I,j,k+1) , tauxDG_u(I,j,k+1) ) + + du = ui(i,J,k) - ui(i,J,kp1) + dv = vi_u(i,J,k) - vi_u(i,J,kp1) + omega_s2x_u(I,j,k+1) = atan2( dv , du) !w ~ Omega_tau2x + + omega_tmp = Omega_tau2x_u(I,j,k+1) - omega_w2x_u(I,j) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2w_u(I,j,k+1) = omega_tmp + + omega_tmp = Omega_tau2x_u(I,j,k+1) - omega_s2x_u(I,j,k+1) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2s_u(I,j,k+1) = omega_tmp !w ~ 0 + + omega_tmp = omega_s2x_u(I,j,k+1) - omega_w2x_u(I,j) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + omega_s2w_u(I,j,k+1) = omega_tmp !w ~ Omega_tau2w + + enddo + endif + enddo + enddo + do J = Jsq, Jeq + do i = is, ie + if( (G%mask2dCv(i,J) > 0.5) ) then + tauxDG_v(i,J,1) = taux_v(i,J) ! SURFACE + tau_v(i,J,1) = ustar2_v(i,J) + Omega_tau2w_v(i,J,1) = 0.0 + Omega_tau2x_v(i,J,1) = omega_w2x_v(i,J) + Omega_tau2s_v(i,J,1) = 0.0 + omega_s2x_v(i,J,1) = omega_w2x_v(i,J) + omega_s2w_v(i,J,1) = 0.0 + + do k=1,nz-1 + kp1 = MIN(k+1 , nz) + tau_v(i,J,k+1) = sqrt ( tauxDG_v(i,J,k+1)*tauxDG_v(i,J,k+1) + tauyDG_v(i,J,k+1)*tauyDG_v(i,J,k+1) ) + Omega_tau2x_v(i,J,k+1) = atan2( tauyDG_v(i,J,k+1) , tauxDG_v(i,J,k+1) ) + + du = ui_v(i,J,k) - ui_v(i,J,kp1) + dv = vi(i,J,k) - vi(i,J,kp1) + omega_s2x_v(i,J,k+1) = atan2( dv , du ) !~ Omega_tau2x + + omega_tmp = Omega_tau2x_v(i,J,k+1) - omega_w2x_v(i,J) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2w_v(i,J,k+1) = omega_tmp + + omega_tmp = Omega_tau2x_v(i,J,k+1) - omega_s2x_v(i,J,k+1) + if (omega_tmp .gt. pi ) omega_tmp = omega_tmp - 2.*pi + if (omega_tmp .le. (0.-pi) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2s_v(i,J,k+1) = omega_tmp !w ~ 0 + + omega_tmp = omega_s2x_v(i,J,k+1) - omega_w2x_v(i,J) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + omega_s2w_v(i,J,k+1) = omega_tmp !w ~ Omega_tau2w + + enddo + endif + enddo + enddo +! ********************************************************************************************** +!w Parameterized stress orientation from the wind at interfaces (tau2x) and centers (tau2x) OVERWRITE to kbl-interface above hbl + du_rot(:,:,:) = 0.0 + dv_rot(:,:,:) = 0.0 + mask3d_u(:,:,:) = 0.0 + mask3d_v(:,:,:) = 0.0 + do j = js,je !w U-points + do I = Isq,Ieq + if( (G%mask2dCu(I,j) > 0.5) ) then + kbld = MIN( (kbl_u(I,j)) , (nz-2) ) + if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1 + !w if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1 + + tauh = tau_u(I,j,kbld+1) + GV%H_subroundoff + omega_tau2wh = omega_tau2w_u(I,j,kbld+1) + + depth = 0. ! surface boundary conditions + tauNLup = 0.0 + do k=1, kbld + depth = depth + CS%h_u(I,j,k) + if ( (L19 > 0) ) then + sigma = MIN ( 1.0 , depth / hbl_u(i,j) ) + G_sig = MIN ( 0.287 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (1.74392*sigma - 2.58538) ) ) + + tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) !w linear stress mag + omega_s2x = Omega_tau2x_u(I,j,k+1) + cos_tmp = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) + sin_tmp = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) + Wind_x = ustar2_u(I,j) * cos(omega_w2x_u(I,j)) !w taux_u primary + Wind_y = ustar2_u(I,j) * sin(omega_w2x_u(I,j)) !w tauy_u interpolated + tauNL_DG = ( Wind_x *cos_tmp + Wind_y *sin_tmp ) !wind in x' + tauNL_CG = ( Wind_y *cos_tmp - Wind_x *sin_tmp ) !WCG in y' + omega_w2s = atan2( tauNL_CG , tauNL_DG ) !W wind to shear x' (limiter) + omega_s2w = 0.0-omega_w2s + tauNL_CG = Cemp_CG * G_sig * tauNL_CG +!OPTIONS + if(L19 .eq. 1) then !A L19=1 + tau_MAG = MAX( tau_MAG , tauNL_CG ) + tauNL_DG = sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - tau_u(I,j,k+1) + endif + + if(L19 .eq. 2) then !B L19=2 + tauNL_CG = MIN( tauNL_CG , tau_MAG ) + tauNL_DG = sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - tau_u(I,j,k+1) + endif + + if(L19 .eq. 3) then !C L19=3 + tauNL_DG = tau_MAG - tau_u(I,j,k+1) + tau_MAG = sqrt( tau_MAG*tau_MAG + tauNL_CG*tauNL_CG ) + endif + omega_tmp = atan2( tauNL_CG , (tau_u(I,j,k+1)+tauNL_DG) ) !W Limiters + + tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp ) !w back to x,y coordinates + tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp ) + tauNLdn = tauNL_X ! SOLUTION + du_rot(I,j,k) = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff) + tauNLup = tauNLdn + + mask3d_u(I,j,k) = tauNL_CG / (tau_MAG) !W (tauNLup - tauNLdn) + mask3d_v(i,j,k) = (tau_u(I,j,k+1)+tauNL_DG) / (tau_MAG) + ! DIAGNOSTICS + tau_u(I,j,k+1) = sqrt( (tauxDG_u(I,j,k+1) + tauNL_X)**2 + (tauyDG_u(I,j,k+1) + tauNL_Y)**2 ) + omega_tau2x = atan2((tauyDG_u(I,j,k+1) + tauNL_Y) , (tauxDG_u(I,j,k+1) + tauNL_X) ) + + omega_tau2w = omega_tau2x - omega_w2x_u(I,j) + if (omega_tau2w .gt. pi ) omega_tau2w = omega_tau2w - 2.*pi + if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi + Omega_tau2w_u(I,j,k+1) = omega_tau2w + Omega_tau2s_u(I,j,k+1) = omega_tmp !W omega_tau2x - Omega_tau2x_u(I,j,k+1) + Omega_tau2x_u(I,j,k+1) = 0.0 - omega_w2s !W omega_s2x !W 0.0 - omega_w2s !W omega_tau2x + + endif + enddo + endif + enddo + enddo +!w V-point dv increment %%%%%%%%%%%%%%%%%%%%%%%%%%%% + do J = Jsq,Jeq + do i = is,ie + if( (G%mask2dCv(i,J) > 0.5) ) then + kbld = MIN( (kbl_v(i,J)) , (nz-2) ) + if ( tau_v(i,J,kbld+2) > tau_v(i,J,kbld+1) ) kbld = kbld + 1 + tauh = tau_v(i,J,kbld+1) + omega_tau2wh = omega_tau2w_u(I,j,kbld+1) + + depth = 0. !surface boundary conditions + tauNLup = 0.0 + do k=1, kbld + depth = depth + CS%h_v(i,J,k) + if ( (L19 > 0) ) then + sigma = MIN ( 1.0 , (depth ) / hbl_v(I,J) ) + G_sig = MIN ( 0.287 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (1.74392*sigma - 2.58538) ) ) + + tau_MAG = (ustar2_v(i,J) * (1.-sigma) ) + (tauh * sigma ) !w linear stress + omega_s2x = Omega_tau2x_v(i,J,k+1) + cos_tmp = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) + sin_tmp = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) + Wind_x = ustar2_v(i,J) * cos(omega_w2x_v(i,J)) !w taux_v interpolated + Wind_y = ustar2_v(i,J) * sin(omega_w2x_v(i,J)) !w tauy_v primary + tauNL_DG = ( Wind_x *cos_tmp + Wind_y *sin_tmp ) + tauNL_CG = ( Wind_y *cos_tmp - Wind_x *sin_tmp ) !w WCG + omega_w2s = atan2( tauNL_CG , tauNL_DG ) ! tau2x' limiter + omega_s2w = 0.0 - omega_w2s + tauNL_CG = Cemp_CG * G_sig * tauNL_CG +!OPTIONS + if(L19 .eq. 1) then !A L19=1 + tau_MAG = MAX( tau_MAG , tauNL_CG ) + tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) + endif + + if(L19 .eq. 2) then !B L19=2 + tauNL_CG = MIN( tauNL_CG , tau_MAG ) + tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) + endif + + if(L19 .eq. 3) then !C L19=3 + tauNL_DG = 0.0 - tau_v(i,J,k+1) + tau_MAG + tau_MAG = sqrt( tau_MAG*tau_MAG + tauNL_CG*tauNL_CG ) + endif + + omega_tmp = atan2( tauNL_CG , tau_v(i,J,k+1) + tauNL_DG ) !W LIMITERS as (tauNL_CG / tau_MAG) + + tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp ) ! back to x,y coordinate + tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp ) + tauNLdn = tauNL_Y + dv_rot(i,J,k) = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) ) ! SOLUTION + tauNLup = tauNLdn + ! DIAGNOSTICS + tau_v(i,J,k+1) = sqrt( (tauxDG_v(i,J,k+1) + tauNL_X)**2 + (tauyDG_v(i,J,k+1) + tauNL_Y)**2 ) + omega_tau2x = atan2( (tauyDG_v(i,J,k+1) + tauNL_Y) , (tauxDG_v(i,J,k+1) + tauNL_X) ) + omega_tau2w = omega_tau2x - omega_w2x_v(i,J) + if (omega_tau2w .gt. pi ) omega_tau2w = omega_tau2w - 2.*pi + if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi + + Omega_tau2w_v(i,J,k+1) = omega_tau2w + Omega_tau2s_v(i,J,k+1) = omega_tmp !W omega_tau2x - Omega_tau2x_v(i,J,k+1) + Omega_tau2x_v(i,J,k+1) = 0.0 - omega_w2s !W omega_s2x !W 0.0 - omega_w2s !W omega_tau2x + endif + enddo + endif + enddo + enddo + if (CS%debug) then + call uvchksum("FP-tau_[uv] ", tau_u, tau_v, G%HI, haloshift=0, scalar_pair=.true.) + call uvchksum("FP-omega_s2x ",omega_s2x_u,omega_s2x_v,G%HI,haloshift=0,scalar_pair=.true.) + call uvchksum("FP-omega_s2w ",omega_s2w_u,omega_s2w_v,G%HI,haloshift=0,scalar_pair=.true.) + call uvchksum("FP-omega_t2w ",omega_tau2x_u,omega_tau2x_v,G%HI,haloshift=0,scalar_pair=.true.) + call uvchksum("FP-omega_t2x ",omega_tau2x_u ,omega_tau2x_v ,G%HI,haloshift=0,scalar_pair=.true.) + call uvchksum("FP-d[uv]_rot ",du_rot, dv_rot, G%HI, haloshift=0,scalar_pair=.true.) + call uvchksum("FP-d[uv]_out ",uold , vold , G%HI, haloshift=0,scalar_pair=.true.) + endif + +!w OUTPUT + do k=1,nz + do j = js,je + do I = Isq,Ieq + ui(I,j,k) = uold(I,j,k) + du_rot(I,j,k) + uold(I,j,k) = du_rot(I,j,k) + enddo + enddo + do J = Jsq,Jeq + do i = is,ie + vi(i,J,k) = vold(i,J,k) + dv_rot(i,J,k) + vold(i,J,k) = dv_rot(i,J,k) + enddo + enddo + enddo + +if( LU_pred .eq. .false. ) then !W CONDITION DIAGNOSTIC OUTPUT THEN POST + do j = js,je + do I = Isq,Ieq + if( (G%mask2dCu(I,j) > 0.5) ) then + kbld = kbl_u(I,j) + ustar2 = ustar2_u(I,j) + tau_u(I,j,1) = tau_u(I,j,1) / ustar2 + Omega_tau2w_u(I,j,1) = Omega_tau2w_u(I,j,1) / pi + Omega_tau2x_u(I,j,1) = Omega_tau2x_u(I,j,1) / pi + Omega_tau2s_u(I,j,1) = Omega_tau2s_u(I,j,1) / pi + do k=1,nz + !w mask3d_u(I,j,k) = + tau_u(I,j,k+1) = tau_u(I,j,k+1) / ustar2 + Omega_tau2w_u(I,j,k+1) = Omega_tau2w_u(I,j,k+1) /pi + Omega_tau2x_u(I,j,k+1) = Omega_tau2x_u(I,j,k+1) /pi + Omega_tau2s_u(I,j,k+1) = Omega_tau2s_u(I,j,k+1) /pi + if( k .eq. kbld+2) then + tau_u(I,j,k) = 0.0 - tau_u(I,j,k) + Omega_tau2w_u(I,j,k) = 1.05 + Omega_tau2x_u(I,j,k) = 1.05 + Omega_tau2s_u(I,j,k) = 1.05 + endif + enddo + Omega_tau2x_u(I,j,nz+1) = omega_w2x_u(I,j) / pi + mask3d_u(I,j,nz) = ustar2_u(I,j) + mask3d_u(I,j,nz-1) = sqrt(taux_u(I,j)*taux_u(I,j) + tauy_u(I,j)*tauy_u(I,j) ) + endif + enddo + enddo + do J = Jsq,Jeq !w v-points + do i = is,ie + if( (G%mask2dCv(i,J) > 0.5) ) then + kbld = kbl_v(i,J) + ustar2 = ustar2_v(i,J) + tau_v(i,J,1) = tau_v(i,J,1) / ustar2 + Omega_tau2w_v(i,J,1) = Omega_tau2w_v(i,J,1) / pi + Omega_tau2x_v(i,J,1) = Omega_tau2x_v(i,J,1) / pi + Omega_tau2s_v(i,J,1) = Omega_tau2s_v(i,J,1) / pi + do k=1,nz + !w mask3d_v(i,J,k) = tauxDG_v(i,J,k) !w vi(i,J,k) - v(i,J,k) !w dv_rot(i,J,k) + tau_v(i,J,k+1) = tau_v(i,J,k+1) / ustar2 + Omega_tau2w_v(i,J,k+1) = Omega_tau2w_v(i,J,k+1) /pi + Omega_tau2x_v(i,J,k+1) = Omega_tau2x_v(i,J,k+1) /pi + Omega_tau2s_v(i,J,k+1) = Omega_tau2s_v(i,J,k+1) /pi + if( k .eq. kbld+2) then + tau_v(i,J,k) = 0.0 - tau_v(i,J,k) + Omega_tau2w_v(i,J,k) = 1.05 + Omega_tau2x_v(i,J,k) = 1.05 + Omega_tau2s_v(i,J,k) = 1.05 + endif + enddo + Omega_tau2x_v(i,J,nz+1) = omega_w2x_v(i,J) / pi + mask3d_v(i,J,nz) = ustar2_v(i,J) + mask3d_v(i,J,nz-1) = sqrt(taux_v(i,J)*taux_v(i,J) + tauy_v(i,J)*tauy_v(i,J) ) + endif + enddo + enddo + + if (CS%id_tauFP_u > 0) call post_data(CS%id_tauFP_u, tau_u, CS%diag) + if (CS%id_tauFP_v > 0) call post_data(CS%id_tauFP_v, tau_v, CS%diag) + if (CS%id_FPtau2s_u > 0) call post_data(CS%id_FPtau2s_u, omega_tau2s_u, CS%diag) + if (CS%id_FPtau2s_v > 0) call post_data(CS%id_FPtau2s_v, omega_tau2s_v, CS%diag) + if (CS%id_FPtau2w_u > 0) call post_data(CS%id_FPtau2w_u, omega_tau2w_u, CS%diag) + if (CS%id_FPtau2w_v > 0) call post_data(CS%id_FPtau2w_v, omega_tau2w_v, CS%diag) + if (CS%id_FPtau2x_u > 0) call post_data(CS%id_FPtau2x_u, omega_tau2x_u, CS%diag) + if (CS%id_FPtau2x_v > 0) call post_data(CS%id_FPtau2x_v, omega_tau2x_v, CS%diag) + if (CS%id_FPmask_u > 0) call post_data(CS%id_FPmask_u, mask3d_u, CS%diag) + if (CS%id_FPmask_v > 0) call post_data(CS%id_FPmask_v, mask3d_v, CS%diag) + if (CS%id_FPhbl_u > 0) call post_data(CS%id_FPhbl_u, hbl_u, CS%diag) + if (CS%id_FPhbl_v > 0) call post_data(CS%id_FPhbl_v, hbl_v, CS%diag) + + if (cs%debug) then + call uvchksum("post viscFPmix [ui,vi]",ui,vi,G%HI,haloshift=0,scalar_pair=.true.) + endif +endif ! LU_pred = false + +end subroutine vertFPmix + !> Perform a fully implicit vertical diffusion !! of momentum. Stress top and bottom boundary conditions are used. !! @@ -1828,6 +2406,38 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', & thickness_units, conversion=GV%H_to_MKS) + !w FPmix + CS%id_FPhbl_u = register_diag_field('ocean_model', 'FPhbl_u', diag%axesCu1, Time, & + 'Boundary-Layer Depth (u-points)','m') !w , conversion=GV%H_to_MKS) + CS%id_FPhbl_v = register_diag_field('ocean_model', 'FPhbl_v', diag%axesCv1, Time, & + 'Boundary-Layer Depth (v-points)','m') + + CS%id_FPmask_u = register_diag_field('ocean_model', 'FPmask_u', diag%axesCuL, Time, & + 'FP overwrite mask (u-points)','binary') + CS%id_FPmask_v = register_diag_field('ocean_model', 'FPmask_v', diag%axesCvL, Time, & + 'FP overwrite mask (v-points)','binary') + + CS%id_tauFP_u = register_diag_field('ocean_model', 'tauFP_u', diag%axesCui, Time, & + 'Stress Mag Profile (u-points)', 'm2 s-2') !w , conversion=GV%H_to_MKS) + CS%id_tauFP_v = register_diag_field('ocean_model', 'tauFP_v', diag%axesCvi, Time, & + 'Stress Mag Profile (v-points)', 'm2 s-2') !w , conversion=GV%H_to_MKS) + + CS%id_FPtau2s_u = register_diag_field('ocean_model', 'FPtau2s_u', diag%axesCui, Time, & + 'stress from shear direction (u-points)', 'pi ') + CS%id_FPtau2s_v = register_diag_field('ocean_model', 'FPtau2s_v', diag%axesCvi, Time, & + 'stress from shear direction (v-points)', 'pi ') + + CS%id_FPtau2w_u = register_diag_field('ocean_model', 'FPtau2w_u', diag%axesCui, Time, & + 'stress from wind direction (u-points)', 'pi ') + CS%id_FPtau2w_v = register_diag_field('ocean_model', 'FPtau2w_v', diag%axesCvi, Time, & + 'stress from wind direction (v-points)', 'pi ') + + CS%id_FPtau2x_u = register_diag_field('ocean_model', 'FPs2w_u', diag%axesCui, Time, & + 'shear from wind (u-points)', 'pi ') + CS%id_FPtau2x_v = register_diag_field('ocean_model', 'FPs2w_v', diag%axesCvi, Time, & + 'shear from wind (v-points)', 'pi ' + ! w - end + CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, & 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_du_dt_visc > 0) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) From cb65bdcefd1a51acea19767b5e12502798aba7ec Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 May 2022 14:33:09 -0600 Subject: [PATCH 3/7] Change name of logical Replaces LU_pred to L_diag, since now this logical only controls if diagnostics should be posted. --- src/core/MOM_dynamics_split_RK2.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 8c3612f50b..f6cf456f98 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -362,7 +362,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth from Cvmix real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. - logical :: LU_pred ! Controls if it is predictor step or not + logical :: L_diag ! Controls if diagostics are posted in the vertFPmix logical :: dyn_p_surf logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the ! relative weightings of the layers in calculating @@ -666,12 +666,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) if (CS%fpmix) then - LU_pred = .true. + L_diag = .false. hbl(:,:) = 0.0 if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) if (ASSOCIATED(CS%energetic_PBL_CSp)) & call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) - call vertFPmix(LU_pred, up, vp, uold, vold, hbl, h, forces, & + call vertFPmix(L_diag, up, vp, uold, vold, hbl, h, forces, & dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) @@ -914,8 +914,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves) if (CS%fpmix) then - LU_pred = .false. - call vertFPmix(LU_pred, u, v, uold, vold, hbl, h, forces, dt, & + L_diag = .true. + call vertFPmix(L_diag, u, v, uold, vold, hbl, h, forces, dt, & G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) From 143d117527746736707c49444ad554cf3f3901d6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 May 2022 15:22:07 -0600 Subject: [PATCH 4/7] Updates to vertFPmix This commit adds the latest updates to the vertFPmix subroutine after Bill Large did some cleaning. We have highlight places in the code where work must be done. --- .../vertical/MOM_vert_friction.F90 | 648 ++++++------------ 1 file changed, 227 insertions(+), 421 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index d5a7aa9804..1d4f7bf646 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -128,8 +128,8 @@ module MOM_vert_friction integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 integer :: id_du_dt_str = -1, id_dv_dt_str = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 - integer :: id_FPmask_u = -1, id_FPmask_v = -1 , id_FPhbl_u = -1, id_FPhbl_v = -1 - integer :: id_tauFP_u = -1, id_tauFP_v = -1 , id_FPtau2x_u = -1, id_FPtau2x_v = -1 + integer :: id_FPdiag_u = -1, id_FPdiag_v = -1 , id_FPw2x = -1 !W id_FPhbl_u = -1, id_FPhbl_v = -1 + integer :: id_tauFP_u = -1, id_tauFP_v = -1 !W, id_FPtau2x_u = -1, id_FPtau2x_v = -1 integer :: id_FPtau2s_u = -1, id_FPtau2s_v = -1, id_FPtau2w_u = -1, id_FPtau2w_v = -1 integer :: id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 @@ -147,9 +147,8 @@ module MOM_vert_friction contains -!> Add nonlocal momentum flux profile increments -!! TODO: add more description -subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OBC) ! FPmix +!> Add nonlocal stress increments to u^n (uold) and v^n (vold) using ui and vi. +subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OBC) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -164,123 +163,77 @@ subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, U real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: hbl_h ! boundary layer depth - logical, intent(inout) :: LU_pred !w predictor step or NOT - type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces - real, intent(in) :: dt !< Time increment [T ~> s] - type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure - type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure + logical, intent(in) :: L_diag !< controls if diagnostics should be posted + type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces + real, intent(in) :: dt !< Time increment [T ~> s] + type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure + type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure ! local variables - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: mask3d_u !Test Plots @ 3-D centers - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: mask3d_v - real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !2-D + ! WGL; TODO: add description to local variables + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: FPdiag_u !< this is for ... + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: FPdiag_v + real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u real, dimension(SZI_(G),SZJB_(G)) :: hbl_v integer, dimension(SZIB_(G),SZJ_(G)) :: kbl_u integer, dimension(SZI_(G),SZJB_(G)) :: kbl_v - real, dimension(SZI_(G),SZJ_(G)) :: ustar2_h !2-D surface real, dimension(SZIB_(G),SZJ_(G)) :: ustar2_u real, dimension(SZI_(G),SZJB_(G)) :: ustar2_v real, dimension(SZIB_(G),SZJ_(G)) :: taux_u real, dimension(SZI_(G),SZJB_(G)) :: tauy_v - real, dimension(SZIB_(G),SZJ_(G)) :: tauy_u - real, dimension(SZI_(G),SZJB_(G)) :: taux_v - real, dimension(SZIB_(G),SZJ_(G)) :: omega_w2x_u - real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v + real, dimension(SZIB_(G),SZJ_(G)) :: omega_w2x_u + real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tau_u !3-D interfaces + ! GMM; TODO: make arrays allocatable if possible + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tau_u real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tau_v real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauxDG_u real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauyDG_u real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauxDG_v real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauyDG_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2x_u - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2x_v real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2s_u real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2s_v real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_s2x_u - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_s2x_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_s2w_u - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_s2w_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: du_rot !3-D centers - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: dv_rot - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: vi_u - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: ui_v - - real :: tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, MAXinc, MINthick - real :: du, dv, du_v, dv_u , dup, dvp , uZero, vZero - real :: fEQband, Cemp_SS , Cemp_LS , Cemp_CG, Cemp_DG , Wgt_SS + + real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp, omega_tmp + real :: du, dv, depth, sigma, Wind_x, Wind_y + real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG - real :: pi, tmp, cos_tmp, sin_tmp, depth, taux, tauy, tauk, tauxI , tauyI, sign_f - real :: tauxh, tauyh, tauh, omega_s2xh, omega_s2wh, omega_tau2xh, omega_tau2wh - real :: taux0, tauy0, tau0, sigma, G_sig, Wind_x, Wind_y, omega_w2s, omega_tau2s,omega_s2x - real :: omega_tau2x, omega_tau2w, omega_SS, omega_LS, omega_tmp, omega_s2xI, omega_s2w - integer :: kblmin, kbld, kp, km, kp1, L19 ,jNseam + real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w + integer :: kblmin, kbld, kp1 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke pi = 4. * atan2(1.,1.) - L19 = 1 !w Options A = 1, B = 2, C = 3 - Cemp_CG = 3.6 !w L91 cross-gradient - Cemp_DG = 1.0 !w L91 down-gradient - MAXinc = -1.0 !w if positive - MINthick= 0.01 !w GV%H_subroundoff !w 0.5 + Cemp_CG = 3.6 kblmin = 1 - jNseam = 457 !w north seam = SZJ_(G) + FPdiag_u(:,:,:) = 0.0 + FPdiag_v(:,:,:) = 0.0 + taux_u(:,:) = 0. + tauy_v(:,:) = 0. -if(LU_pred ) then !w predictor step only, surface forcing - ustar2_h(:,:) = 0. - do j = js,je !w ?GMM -1,+1 with forces% - do i = is,ie - ustar2_h(i,j) = forces%ustar(i,j) * forces%ustar(i,j) - !w omega_w2x_h(i,j) = forces%omega_w2x(i,j) - enddo - enddo - call pass_var(ustar2_h ,G%Domain) ! update halos ?GMM - call pass_var( hbl_h ,G%Domain) - -! SURFACE ustar2 and x-stress to u points and ustar2 and y-stress to v points - ustar2_u(:,:) = 0. - ustar2_v(:,:) = 0. - hbl_u(:,:) = 0. - hbl_v(:,:) = 0. - taux_u(:,:) = 0. - tauy_v(:,:) = 0. do j = js,je do I = Isq,Ieq - tmp = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) - ustar2_u(I,j)=(G%mask2dT(i,j)*ustar2_h(i,j)+G%mask2dT(i+1,j)*ustar2_h(i+1,j))/tmp - hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j)* hbl_h(i+1,j)) /tmp - taux_u(I,j) = forces%taux(I,j) / GV%H_to_RZ + taux_u(I,j) = forces%taux(I,j) / GV%H_to_RZ !W rho0=1035. enddo enddo + do J = Jsq,Jeq do i = is,ie - tmp = MAX ( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1) ) ) - ustar2_u(I,j)=(G%mask2dT(i,j)*ustar2_h(i,j)+G%mask2dT(i+1,j)*ustar2_h(i+1,j))/tmp - hbl_v(i,J) = (G%mask2dT(i,j)* hbl_h(i,J) + G%mask2dT(i,j+1)* hbl_h(i,j+1)) /tmp - if( j > jNseam-1 ) then - ustar2_v(i,J) = ustar2_h(i,j ) !w ( j > 456 ) j >= 457 - hbl_v(i,J) = hbl_h(i,j) - endif - tauy_v(i,J) = forces%tauy(i,J) / GV%H_to_RZ + tauy_v(i,J) = forces%tauy(i,J) / GV%H_to_RZ enddo enddo - call pass_vector(taux_u , tauy_v, G%Domain, To_All+Scalar_Pair) - if (CS%debug) then - call uvchksum("ustar2 ",ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.) - call uvchksum(" hbl ", hbl_u , hbl_v , G%HI, haloshift=0, scalar_pair=.true.) - call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=0, scalar_pair=.true.) - endif -!W endif !w predictor step -!w surface tauy_u , taux_v and omega_w2x_[u,v] & Implicit interface stresses tauxDG_u and tauyDG_v - tauy_u(:,:) = 0.0 - taux_v(:,:) = 0.0 - kbl_u(:,:) = 0 - kbl_v(:,:) = 0 + call pass_var( hbl_h ,G%Domain, halo=1 ) + call pass_vector(taux_u , tauy_v, G%Domain, To_All ) + ustar2_u(:,:) = 0. + ustar2_v(:,:) = 0. + hbl_u(:,:) = 0. + hbl_v(:,:) = 0. + kbl_u(:,:) = 0 + kbl_v(:,:) = 0 omega_w2x_u(:,:) = 0.0 omega_w2x_v(:,:) = 0.0 tauxDG_u(:,:,:) = 0.0 @@ -288,16 +241,14 @@ subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, U do j = js,je do I = Isq,Ieq if( (G%mask2dCu(I,j) > 0.5) ) then - tauy = 0.0 - tmp = MAX(1.0, (G%mask2dCv(i,j) + G%mask2dCv(i,j-1) + G%mask2dCv(i+1,j) + G%mask2dCv(i+1,j-1) ) ) - if ( G%mask2dCv(i ,j ) > 0.5 ) tauy = tauy + tauy_v(i ,j ) - if ( G%mask2dCv(i ,j-1) > 0.5 ) tauy = tauy + tauy_v(i ,j-1) - if ( G%mask2dCv(i+1,j ) > 0.5 ) tauy = tauy + tauy_v(i+1,j ) - if ( G%mask2dCv(i+1,j-1) > 0.5 ) tauy = tauy + tauy_v(i+1,j-1) - tauy = tauy / tmp - tauy_u(I,j) = (tauy/(abs(tauy)+GV%H_subroundoff)) * sqrt(MAX(GV%H_subroundoff,ustar2_u(I,j)*ustar2_u(I,j)-taux_u(I,j)*taux_u(I,j) )) - omega_w2x_u(I,j) = atan2( tauy_u(I,j) , taux_u(I,j) ) - tauxDG_u(I,j,1) = taux_u(I,j) !w ustar2_u(I,j) * cos(omega_w2x_u(I,j)) + tmp = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) + hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j) * hbl_h(i+1,j)) /tmp + tmp = MAX(1.0, (G%mask2dCv(i,j) + G%mask2dCv(i,j-1) + G%mask2dCv(i+1,j) + G%mask2dCv(i+1,j-1) ) ) + tauy = ( G%mask2dCv(i ,j )*tauy_v(i ,j ) + G%mask2dCv(i ,j-1)*tauy_v(i ,j-1) & + + G%mask2dCv(i+1,j )*tauy_v(i+1,j ) + G%mask2dCv(i+1,j-1)*tauy_v(i+1,j-1) ) / tmp + ustar2_u(I,j) = sqrt( taux_u(I,j)*taux_u(I,j) + tauy*tauy ) + omega_w2x_u(I,j) = atan2( tauy , taux_u(I,j) ) + tauxDG_u(I,j,1) = taux_u(I,j) depth = 0.0 do k = 1, nz depth = depth + CS%h_u(I,j,k) @@ -312,22 +263,14 @@ subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, U do J = Jsq,Jeq do i = is,ie if( (G%mask2dCv(i,J) > 0.5) ) then - taux = 0.0 - if ( j < 457 ) then - tmp = MAX(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i,j+1) + G%mask2dCu(i-1,j) + G%mask2dCu(i-1,j+1) ) ) - if ( G%mask2dCu(i ,j ) > 0.5 ) taux = taux + taux_u(i ,j ) - if ( G%mask2dCu(i ,j+1) > 0.5 ) taux = taux + taux_u(i ,j+1) - if ( G%mask2dCu(i-1,j ) > 0.5 ) taux = taux + taux_u(i-1,j ) - if ( G%mask2dCu(i-1,j+1) > 0.5 ) taux = taux + taux_u(i-1,j+1) - else - tmp = MAX(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i-1,j) ) ) - if ( G%mask2dCu(i ,j ) > 0.5 ) taux = taux + taux_u(i ,j ) - if ( G%mask2dCu(i-1,j ) > 0.5 ) taux = taux + taux_u(i-1,j ) - endif - taux = taux / tmp - taux_v(i,J) = (taux/(abs(taux)+GV%H_subroundoff)) * sqrt(MAX(GV%H_subroundoff,ustar2_v(i,J)*ustar2_v(i,J)-tauy_v(i,J)*tauy_v(i,J) )) - omega_w2x_v(i,J) = atan2( tauy_v(i,J) , taux_v(i,J) ) - tauyDG_v(i,J,1) = tauy_v(i,J) !w ustar2_v(i,J) * cos(omega_w2x_v(i,J)) + tmp = MAX ( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1) ) ) + hbl_v(i,J) = (G%mask2dT(i,j)* hbl_h(i,J) + G%mask2dT(i,j+1) * hbl_h(i,j+1)) /tmp + tmp = MAX(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i,j+1) + G%mask2dCu(i-1,j) + G%mask2dCu(i-1,j+1) ) ) + taux = ( G%mask2dCu(i ,j )*taux_u(i ,j ) + G%mask2dCu(i ,j+1)*taux_u(i ,j+1) & + + G%mask2dCu(i-1,j )*taux_u(i-1,j ) + G%mask2dCu(i-1,j+1)*taux_u(i-1,j+1) ) / tmp + ustar2_v(i,J) = sqrt( tauy_v(i,J)*tauy_v(i,J) + taux*taux ) + omega_w2x_v(i,J) = atan2( tauy_v(i,J) , taux ) + tauyDG_v(i,J,1) = tauy_v(i,J) depth = 0.0 do k = 1, nz depth = depth + CS%h_v(i,J,k) @@ -339,98 +282,80 @@ subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, U endif enddo enddo -endif !w predictor step - -! Thickness weighted diagnostic interpolations ! Copy Implicit [uv]i to [uv]old - call pass_vector(ui,vi, G%Domain, To_All+Scalar_Pair) - vi_u(:,:,:) = 0. - ui_v(:,:,:) = 0. - tauxDG_u(:,:,:) = 0.0 - tauyDG_v(:,:,:) = 0.0 - tauxDG_v(:,:,:) = 0. - tauyDG_u(:,:,:) = 0. + + if (CS%debug) then + call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=1, scalar_pair=.true.) + call uvchksum("ustar2 ",ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.) + call uvchksum(" hbl ", hbl_u , hbl_v , G%HI, haloshift=0, scalar_pair=.true.) + endif + + ! Compute downgradient stresses do k = 1, nz - kp = MIN( k+1 , nz) - do j = js-1 ,je+1 - do I = Isq-1, Ieq+1 - tauxDG_u(I,j,k+1) = CS%a_u(I,j,kp) * (ui(I,j,k) - ui(I,j,kp)) + kp1 = MIN( k+1 , nz) + do j = js ,je + do I = Isq , Ieq + tauxDG_u(I,j,k+1) = CS%a_u(I,j,kp1) * (ui(I,j,k) - ui(I,j,kp1)) enddo enddo - do J = Jsq-1, Jeq+1 - do i = is-1, ie+1 - tauyDG_v(i,J,k+1) = CS%a_v(i,J,kp) * (vi(i,J,k) - vi(i,J,kp)) + do J = Jsq , Jeq + do i = is , ie + tauyDG_v(i,J,k+1) = CS%a_v(i,J,kp1) * (vi(i,J,k) - vi(i,J,kp1)) enddo enddo + enddo + + call pass_vector(tauxDG_u, tauyDG_v , G%Domain, To_All) + call pass_vector(ui,vi, G%Domain, To_All) + tauxDG_v(:,:,:) = 0. + tauyDG_u(:,:,:) = 0. + ! Thickness weighted interpolations + do k = 1, nz ! v to u points do j = js , je do I = Isq, Ieq - vi_u(I,j,k) = set_v_at_u(vi, h, G, GV, I, j, k, G%mask2dCv, OBC) - tauyDG_u(I,j,k)= set_v_at_u(tauyDG_v, h, G, GV, I, j, k, G%mask2dCv, OBC) + tauyDG_u(I,j,k) = set_v_at_u(tauyDG_v, h, G, GV, I, j, k, G%mask2dCv, OBC) enddo enddo ! u to v points do J = Jsq, Jeq do i = is, ie - ui_v(I,j,k) = set_u_at_v(ui, h, G, GV, i, J, k, G%mask2dCu, OBC) - tauxDG_v(i,J,k)= set_u_at_v(tauxDG_u, h, G, GV, i, J, k, G%mask2dCu, OBC) + tauxDG_v(i,J,k) = set_u_at_v(tauxDG_u, h, G, GV, i, J, k, G%mask2dCu, OBC) enddo enddo enddo if (CS%debug) then - call uvchksum(" vi_u ui_v ", vi_u , ui_v , G%HI, haloshift=0, scalar_pair=.true.) + call uvchksum(" tauyDG_u tauxDG_v",tauyDG_u,tauxDG_v, G%HI, haloshift=0, scalar_pair=.true.) endif -! compute angles, tau2x_[u,v], tau2w_[u,v], tau2s_[u,v], s2x_[u,v], s2w_[u,v] and stress mag tau_[u,v] - omega_tau2x_u(:,:,:) = 0.0 - omega_tau2x_v(:,:,:) = 0.0 + ! compute angles, tau2x_[u,v], tau2w_[u,v], tau2s_[u,v], s2w_[u,v] and stress mag tau_[u,v] omega_tau2w_u(:,:,:) = 0.0 omega_tau2w_v(:,:,:) = 0.0 omega_tau2s_u(:,:,:) = 0.0 omega_tau2s_v(:,:,:) = 0.0 - omega_s2x_u(:,:,:) = 0.0 - omega_s2x_v(:,:,:) = 0.0 - omega_s2w_u(:,:,:) = 0. - omega_s2w_v(:,:,:) = 0. tau_u(:,:,:) = 0.0 tau_v(:,:,:) = 0.0 -!w Default implicit (I) stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] Profiles + !w Default implicit (I) stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] Profiles do j = js,je do I = Isq,Ieq if( (G%mask2dCu(I,j) > 0.5) ) then - tauyDG_u(I,j,1) = tauy_u(I,j) ! SURFACE - tau_u(I,j,1) = ustar2_u(I,j) !w stress magnitude + ! SURFACE + tauyDG_u(I,j,1) = ustar2_u(I,j) * cos(omega_w2x_u(I,j)) + tau_u(I,j,1) = ustar2_u(I,j) Omega_tau2w_u(I,j,1) = 0.0 - Omega_tau2x_u(I,j,1) = omega_w2x_u(I,j) Omega_tau2s_u(I,j,1) = 0.0 - omega_s2x_u(I,j,1) = omega_w2x_u(I,j) - omega_s2w_u(I,j,1) = 0.0 + ! WGL; TODO: can we use set_v_at_u to get tauyDG_u? do k=1,nz kp1 = MIN(k+1 , nz) tau_u(I,j,k+1) = sqrt( tauxDG_u(I,j,k+1)*tauxDG_u(I,j,k+1) + tauyDG_u(I,j,k+1)*tauyDG_u(I,j,k+1)) - Omega_tau2x_u(I,j,k+1) = atan2( tauyDG_u(I,j,k+1) , tauxDG_u(I,j,k+1) ) - - du = ui(i,J,k) - ui(i,J,kp1) - dv = vi_u(i,J,k) - vi_u(i,J,kp1) - omega_s2x_u(I,j,k+1) = atan2( dv , du) !w ~ Omega_tau2x - - omega_tmp = Omega_tau2x_u(I,j,k+1) - omega_w2x_u(I,j) + Omega_tau2x = atan2( tauyDG_u(I,j,k+1) , tauxDG_u(I,j,k+1) ) + omega_tmp = Omega_tau2x - omega_w2x_u(I,j) if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi Omega_tau2w_u(I,j,k+1) = omega_tmp - - omega_tmp = Omega_tau2x_u(I,j,k+1) - omega_s2x_u(I,j,k+1) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - Omega_tau2s_u(I,j,k+1) = omega_tmp !w ~ 0 - - omega_tmp = omega_s2x_u(I,j,k+1) - omega_w2x_u(I,j) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - omega_s2w_u(I,j,k+1) = omega_tmp !w ~ Omega_tau2w - + Omega_tau2s_u(I,j,k+1) = 0.0 enddo endif enddo @@ -438,49 +363,30 @@ subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, U do J = Jsq, Jeq do i = is, ie if( (G%mask2dCv(i,J) > 0.5) ) then - tauxDG_v(i,J,1) = taux_v(i,J) ! SURFACE + ! SURFACE + tauxDG_v(i,J,1) = ustar2_v(i,J) * sin(omega_w2x_v(i,J)) tau_v(i,J,1) = ustar2_v(i,J) Omega_tau2w_v(i,J,1) = 0.0 - Omega_tau2x_v(i,J,1) = omega_w2x_v(i,J) Omega_tau2s_v(i,J,1) = 0.0 - omega_s2x_v(i,J,1) = omega_w2x_v(i,J) - omega_s2w_v(i,J,1) = 0.0 + ! WGL; TODO: can we use set_u_at_v to get tauxDG_v? do k=1,nz-1 kp1 = MIN(k+1 , nz) tau_v(i,J,k+1) = sqrt ( tauxDG_v(i,J,k+1)*tauxDG_v(i,J,k+1) + tauyDG_v(i,J,k+1)*tauyDG_v(i,J,k+1) ) - Omega_tau2x_v(i,J,k+1) = atan2( tauyDG_v(i,J,k+1) , tauxDG_v(i,J,k+1) ) - - du = ui_v(i,J,k) - ui_v(i,J,kp1) - dv = vi(i,J,k) - vi(i,J,kp1) - omega_s2x_v(i,J,k+1) = atan2( dv , du ) !~ Omega_tau2x - - omega_tmp = Omega_tau2x_v(i,J,k+1) - omega_w2x_v(i,J) + omega_tau2x = atan2( tauyDG_v(i,J,k+1) , tauxDG_v(i,J,k+1) ) + omega_tmp = omega_tau2x - omega_w2x_v(i,J) if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi Omega_tau2w_v(i,J,k+1) = omega_tmp - - omega_tmp = Omega_tau2x_v(i,J,k+1) - omega_s2x_v(i,J,k+1) - if (omega_tmp .gt. pi ) omega_tmp = omega_tmp - 2.*pi - if (omega_tmp .le. (0.-pi) ) omega_tmp = omega_tmp + 2.*pi - Omega_tau2s_v(i,J,k+1) = omega_tmp !w ~ 0 - - omega_tmp = omega_s2x_v(i,J,k+1) - omega_w2x_v(i,J) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - omega_s2w_v(i,J,k+1) = omega_tmp !w ~ Omega_tau2w - + Omega_tau2s_v(i,J,k+1) = 0.0 enddo endif enddo enddo -! ********************************************************************************************** -!w Parameterized stress orientation from the wind at interfaces (tau2x) and centers (tau2x) OVERWRITE to kbl-interface above hbl - du_rot(:,:,:) = 0.0 - dv_rot(:,:,:) = 0.0 - mask3d_u(:,:,:) = 0.0 - mask3d_v(:,:,:) = 0.0 - do j = js,je !w U-points + + ! Parameterized stress orientation from the wind at interfaces (tau2x) + ! and centers (tau2x) OVERWRITE to kbl-interface above hbl + do j = js,je do I = Isq,Ieq if( (G%mask2dCu(I,j) > 0.5) ) then kbld = MIN( (kbl_u(I,j)) , (nz-2) ) @@ -488,238 +394,151 @@ subroutine vertFPmix(LU_pred, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, U !w if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1 tauh = tau_u(I,j,kbld+1) + GV%H_subroundoff - omega_tau2wh = omega_tau2w_u(I,j,kbld+1) - - depth = 0. ! surface boundary conditions + ! surface boundary conditions + depth = 0. tauNLup = 0.0 do k=1, kbld depth = depth + CS%h_u(I,j,k) - if ( (L19 > 0) ) then - sigma = MIN ( 1.0 , depth / hbl_u(i,j) ) - G_sig = MIN ( 0.287 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (1.74392*sigma - 2.58538) ) ) - - tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) !w linear stress mag - omega_s2x = Omega_tau2x_u(I,j,k+1) - cos_tmp = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) - sin_tmp = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) - Wind_x = ustar2_u(I,j) * cos(omega_w2x_u(I,j)) !w taux_u primary - Wind_y = ustar2_u(I,j) * sin(omega_w2x_u(I,j)) !w tauy_u interpolated - tauNL_DG = ( Wind_x *cos_tmp + Wind_y *sin_tmp ) !wind in x' - tauNL_CG = ( Wind_y *cos_tmp - Wind_x *sin_tmp ) !WCG in y' - omega_w2s = atan2( tauNL_CG , tauNL_DG ) !W wind to shear x' (limiter) - omega_s2w = 0.0-omega_w2s - tauNL_CG = Cemp_CG * G_sig * tauNL_CG -!OPTIONS - if(L19 .eq. 1) then !A L19=1 - tau_MAG = MAX( tau_MAG , tauNL_CG ) - tauNL_DG = sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - tau_u(I,j,k+1) - endif - - if(L19 .eq. 2) then !B L19=2 - tauNL_CG = MIN( tauNL_CG , tau_MAG ) - tauNL_DG = sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - tau_u(I,j,k+1) - endif - - if(L19 .eq. 3) then !C L19=3 - tauNL_DG = tau_MAG - tau_u(I,j,k+1) - tau_MAG = sqrt( tau_MAG*tau_MAG + tauNL_CG*tauNL_CG ) - endif - omega_tmp = atan2( tauNL_CG , (tau_u(I,j,k+1)+tauNL_DG) ) !W Limiters - - tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp ) !w back to x,y coordinates - tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp ) - tauNLdn = tauNL_X ! SOLUTION - du_rot(I,j,k) = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff) - tauNLup = tauNLdn - - mask3d_u(I,j,k) = tauNL_CG / (tau_MAG) !W (tauNLup - tauNLdn) - mask3d_v(i,j,k) = (tau_u(I,j,k+1)+tauNL_DG) / (tau_MAG) - ! DIAGNOSTICS - tau_u(I,j,k+1) = sqrt( (tauxDG_u(I,j,k+1) + tauNL_X)**2 + (tauyDG_u(I,j,k+1) + tauNL_Y)**2 ) - omega_tau2x = atan2((tauyDG_u(I,j,k+1) + tauNL_Y) , (tauxDG_u(I,j,k+1) + tauNL_X) ) - - omega_tau2w = omega_tau2x - omega_w2x_u(I,j) - if (omega_tau2w .gt. pi ) omega_tau2w = omega_tau2w - 2.*pi - if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi - Omega_tau2w_u(I,j,k+1) = omega_tau2w - Omega_tau2s_u(I,j,k+1) = omega_tmp !W omega_tau2x - Omega_tau2x_u(I,j,k+1) - Omega_tau2x_u(I,j,k+1) = 0.0 - omega_w2s !W omega_s2x !W 0.0 - omega_w2s !W omega_tau2x - - endif + sigma = MIN ( 1.0 , depth / hbl_u(i,j) ) + + ! linear stress mag + tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) + cos_tmp = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) + sin_tmp = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) + + ! rotate to wind coordinates + Wind_x = ustar2_u(I,j) * cos(omega_w2x_u(I,j)) + Wind_y = ustar2_u(I,j) * sin(omega_w2x_u(I,j)) + tauNL_DG = ( Wind_x *cos_tmp + Wind_y *sin_tmp ) + tauNL_CG = ( Wind_y *cos_tmp - Wind_x *sin_tmp ) + omega_w2s = atan2( tauNL_CG , tauNL_DG ) + omega_s2w = 0.0-omega_w2s + tauNL_CG = Cemp_CG * G_sig(sigma) * tauNL_CG + tau_MAG = MAX( tau_MAG , tauNL_CG ) + tauNL_DG = sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - tau_u(I,j,k+1) + + ! back to x,y coordinates + tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp ) + tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp ) + tauNLdn = tauNL_X + + ! nonlocal increment and update to uold + du = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff) + ui(I,j,k) = uold(I,j,k) + du + uold(I,j,k) = du + tauNLup = tauNLdn + + ! diagnostics + FPdiag_u(I,j,k+1) = tauNL_CG / (tau_MAG + GV%H_subroundoff) + Omega_tau2s_u(I,j,k+1) = atan2( tauNL_CG , (tau_u(I,j,k+1)+tauNL_DG) ) + tau_u(I,j,k+1) = sqrt( (tauxDG_u(I,j,k+1) + tauNL_X)**2 + (tauyDG_u(I,j,k+1) + tauNL_Y)**2 ) + omega_tau2x = atan2((tauyDG_u(I,j,k+1) + tauNL_Y) , (tauxDG_u(I,j,k+1) + tauNL_X) ) + omega_tau2w = omega_tau2x - omega_w2x_u(I,j) + if (omega_tau2w >= pi ) omega_tau2w = omega_tau2w - 2.*pi + if (omega_tau2w <= (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi + Omega_tau2w_u(I,j,k+1) = omega_tau2w + enddo + do k= kbld+1, nz + ui(I,j,k) = uold(I,j,k) + uold(I,j,k) = 0.0 enddo endif enddo enddo -!w V-point dv increment %%%%%%%%%%%%%%%%%%%%%%%%%%%% + + ! v-point dv increment do J = Jsq,Jeq do i = is,ie if( (G%mask2dCv(i,J) > 0.5) ) then kbld = MIN( (kbl_v(i,J)) , (nz-2) ) if ( tau_v(i,J,kbld+2) > tau_v(i,J,kbld+1) ) kbld = kbld + 1 tauh = tau_v(i,J,kbld+1) - omega_tau2wh = omega_tau2w_u(I,j,kbld+1) - depth = 0. !surface boundary conditions + !surface boundary conditions + depth = 0. tauNLup = 0.0 do k=1, kbld depth = depth + CS%h_v(i,J,k) - if ( (L19 > 0) ) then - sigma = MIN ( 1.0 , (depth ) / hbl_v(I,J) ) - G_sig = MIN ( 0.287 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (1.74392*sigma - 2.58538) ) ) - - tau_MAG = (ustar2_v(i,J) * (1.-sigma) ) + (tauh * sigma ) !w linear stress - omega_s2x = Omega_tau2x_v(i,J,k+1) - cos_tmp = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) - sin_tmp = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) - Wind_x = ustar2_v(i,J) * cos(omega_w2x_v(i,J)) !w taux_v interpolated - Wind_y = ustar2_v(i,J) * sin(omega_w2x_v(i,J)) !w tauy_v primary - tauNL_DG = ( Wind_x *cos_tmp + Wind_y *sin_tmp ) - tauNL_CG = ( Wind_y *cos_tmp - Wind_x *sin_tmp ) !w WCG - omega_w2s = atan2( tauNL_CG , tauNL_DG ) ! tau2x' limiter - omega_s2w = 0.0 - omega_w2s - tauNL_CG = Cemp_CG * G_sig * tauNL_CG -!OPTIONS - if(L19 .eq. 1) then !A L19=1 - tau_MAG = MAX( tau_MAG , tauNL_CG ) - tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - endif - - if(L19 .eq. 2) then !B L19=2 - tauNL_CG = MIN( tauNL_CG , tau_MAG ) - tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) - endif - - if(L19 .eq. 3) then !C L19=3 - tauNL_DG = 0.0 - tau_v(i,J,k+1) + tau_MAG - tau_MAG = sqrt( tau_MAG*tau_MAG + tauNL_CG*tauNL_CG ) - endif + sigma = MIN ( 1.0 , (depth ) / hbl_v(I,J) ) + + ! linear stress + tau_MAG = (ustar2_v(i,J) * (1.-sigma) ) + (tauh * sigma ) + cos_tmp = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) + sin_tmp = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) + + ! rotate into wind coordinate + Wind_x = ustar2_v(i,J) * cos(omega_w2x_v(i,J)) + Wind_y = ustar2_v(i,J) * sin(omega_w2x_v(i,J)) + tauNL_DG = ( Wind_x *cos_tmp + Wind_y *sin_tmp ) + tauNL_CG = ( Wind_y *cos_tmp - Wind_x *sin_tmp ) + omega_w2s = atan2( tauNL_CG , tauNL_DG ) + omega_s2w = 0.0 - omega_w2s + tauNL_CG = Cemp_CG * G_sig(sigma) * tauNL_CG + tau_MAG = MAX( tau_MAG , tauNL_CG ) + tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt( tau_MAG*tau_MAG - tauNL_CG*tauNL_CG ) + + ! back to x,y coordinate + tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp ) + tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp ) + tauNLdn = tauNL_Y + dv = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) ) + vi(i,J,k) = vold(i,J,k) + dv + vold(i,J,k) = dv + tauNLup = tauNLdn + + ! diagnostics + FPdiag_v(i,j,k+1) = tau_MAG / tau_v(i,J,k+1) + Omega_tau2s_v(i,J,k+1) = atan2( tauNL_CG , tau_v(i,J,k+1) + tauNL_DG ) + tau_v(i,J,k+1) = sqrt( (tauxDG_v(i,J,k+1) + tauNL_X)**2 + (tauyDG_v(i,J,k+1) + tauNL_Y)**2 ) + omega_tau2x = atan2( (tauyDG_v(i,J,k+1) + tauNL_Y) , (tauxDG_v(i,J,k+1) + tauNL_X) ) + omega_tau2w = omega_tau2x - omega_w2x_v(i,J) + if (omega_tau2w .gt. pi ) omega_tau2w = omega_tau2w - 2.*pi + if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi + Omega_tau2w_v(i,J,k+1) = omega_tau2w + enddo - omega_tmp = atan2( tauNL_CG , tau_v(i,J,k+1) + tauNL_DG ) !W LIMITERS as (tauNL_CG / tau_MAG) - - tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp ) ! back to x,y coordinate - tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp ) - tauNLdn = tauNL_Y - dv_rot(i,J,k) = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) ) ! SOLUTION - tauNLup = tauNLdn - ! DIAGNOSTICS - tau_v(i,J,k+1) = sqrt( (tauxDG_v(i,J,k+1) + tauNL_X)**2 + (tauyDG_v(i,J,k+1) + tauNL_Y)**2 ) - omega_tau2x = atan2( (tauyDG_v(i,J,k+1) + tauNL_Y) , (tauxDG_v(i,J,k+1) + tauNL_X) ) - omega_tau2w = omega_tau2x - omega_w2x_v(i,J) - if (omega_tau2w .gt. pi ) omega_tau2w = omega_tau2w - 2.*pi - if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi - - Omega_tau2w_v(i,J,k+1) = omega_tau2w - Omega_tau2s_v(i,J,k+1) = omega_tmp !W omega_tau2x - Omega_tau2x_v(i,J,k+1) - Omega_tau2x_v(i,J,k+1) = 0.0 - omega_w2s !W omega_s2x !W 0.0 - omega_w2s !W omega_tau2x - endif + do k= kbld+1, nz + vi(i,J,k) = vold(i,J,k) + vold(i,J,k) = 0.0 enddo endif enddo enddo + if (CS%debug) then call uvchksum("FP-tau_[uv] ", tau_u, tau_v, G%HI, haloshift=0, scalar_pair=.true.) - call uvchksum("FP-omega_s2x ",omega_s2x_u,omega_s2x_v,G%HI,haloshift=0,scalar_pair=.true.) - call uvchksum("FP-omega_s2w ",omega_s2w_u,omega_s2w_v,G%HI,haloshift=0,scalar_pair=.true.) - call uvchksum("FP-omega_t2w ",omega_tau2x_u,omega_tau2x_v,G%HI,haloshift=0,scalar_pair=.true.) - call uvchksum("FP-omega_t2x ",omega_tau2x_u ,omega_tau2x_v ,G%HI,haloshift=0,scalar_pair=.true.) - call uvchksum("FP-d[uv]_rot ",du_rot, dv_rot, G%HI, haloshift=0,scalar_pair=.true.) - call uvchksum("FP-d[uv]_out ",uold , vold , G%HI, haloshift=0,scalar_pair=.true.) endif -!w OUTPUT - do k=1,nz - do j = js,je - do I = Isq,Ieq - ui(I,j,k) = uold(I,j,k) + du_rot(I,j,k) - uold(I,j,k) = du_rot(I,j,k) - enddo - enddo - do J = Jsq,Jeq - do i = is,ie - vi(i,J,k) = vold(i,J,k) + dv_rot(i,J,k) - vold(i,J,k) = dv_rot(i,J,k) - enddo - enddo - enddo - -if( LU_pred .eq. .false. ) then !W CONDITION DIAGNOSTIC OUTPUT THEN POST - do j = js,je - do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - kbld = kbl_u(I,j) - ustar2 = ustar2_u(I,j) - tau_u(I,j,1) = tau_u(I,j,1) / ustar2 - Omega_tau2w_u(I,j,1) = Omega_tau2w_u(I,j,1) / pi - Omega_tau2x_u(I,j,1) = Omega_tau2x_u(I,j,1) / pi - Omega_tau2s_u(I,j,1) = Omega_tau2s_u(I,j,1) / pi - do k=1,nz - !w mask3d_u(I,j,k) = - tau_u(I,j,k+1) = tau_u(I,j,k+1) / ustar2 - Omega_tau2w_u(I,j,k+1) = Omega_tau2w_u(I,j,k+1) /pi - Omega_tau2x_u(I,j,k+1) = Omega_tau2x_u(I,j,k+1) /pi - Omega_tau2s_u(I,j,k+1) = Omega_tau2s_u(I,j,k+1) /pi - if( k .eq. kbld+2) then - tau_u(I,j,k) = 0.0 - tau_u(I,j,k) - Omega_tau2w_u(I,j,k) = 1.05 - Omega_tau2x_u(I,j,k) = 1.05 - Omega_tau2s_u(I,j,k) = 1.05 - endif - enddo - Omega_tau2x_u(I,j,nz+1) = omega_w2x_u(I,j) / pi - mask3d_u(I,j,nz) = ustar2_u(I,j) - mask3d_u(I,j,nz-1) = sqrt(taux_u(I,j)*taux_u(I,j) + tauy_u(I,j)*tauy_u(I,j) ) - endif - enddo - enddo - do J = Jsq,Jeq !w v-points - do i = is,ie - if( (G%mask2dCv(i,J) > 0.5) ) then - kbld = kbl_v(i,J) - ustar2 = ustar2_v(i,J) - tau_v(i,J,1) = tau_v(i,J,1) / ustar2 - Omega_tau2w_v(i,J,1) = Omega_tau2w_v(i,J,1) / pi - Omega_tau2x_v(i,J,1) = Omega_tau2x_v(i,J,1) / pi - Omega_tau2s_v(i,J,1) = Omega_tau2s_v(i,J,1) / pi - do k=1,nz - !w mask3d_v(i,J,k) = tauxDG_v(i,J,k) !w vi(i,J,k) - v(i,J,k) !w dv_rot(i,J,k) - tau_v(i,J,k+1) = tau_v(i,J,k+1) / ustar2 - Omega_tau2w_v(i,J,k+1) = Omega_tau2w_v(i,J,k+1) /pi - Omega_tau2x_v(i,J,k+1) = Omega_tau2x_v(i,J,k+1) /pi - Omega_tau2s_v(i,J,k+1) = Omega_tau2s_v(i,J,k+1) /pi - if( k .eq. kbld+2) then - tau_v(i,J,k) = 0.0 - tau_v(i,J,k) - Omega_tau2w_v(i,J,k) = 1.05 - Omega_tau2x_v(i,J,k) = 1.05 - Omega_tau2s_v(i,J,k) = 1.05 - endif - enddo - Omega_tau2x_v(i,J,nz+1) = omega_w2x_v(i,J) / pi - mask3d_v(i,J,nz) = ustar2_v(i,J) - mask3d_v(i,J,nz-1) = sqrt(taux_v(i,J)*taux_v(i,J) + tauy_v(i,J)*tauy_v(i,J) ) - endif - enddo - enddo - - if (CS%id_tauFP_u > 0) call post_data(CS%id_tauFP_u, tau_u, CS%diag) - if (CS%id_tauFP_v > 0) call post_data(CS%id_tauFP_v, tau_v, CS%diag) - if (CS%id_FPtau2s_u > 0) call post_data(CS%id_FPtau2s_u, omega_tau2s_u, CS%diag) - if (CS%id_FPtau2s_v > 0) call post_data(CS%id_FPtau2s_v, omega_tau2s_v, CS%diag) - if (CS%id_FPtau2w_u > 0) call post_data(CS%id_FPtau2w_u, omega_tau2w_u, CS%diag) - if (CS%id_FPtau2w_v > 0) call post_data(CS%id_FPtau2w_v, omega_tau2w_v, CS%diag) - if (CS%id_FPtau2x_u > 0) call post_data(CS%id_FPtau2x_u, omega_tau2x_u, CS%diag) - if (CS%id_FPtau2x_v > 0) call post_data(CS%id_FPtau2x_v, omega_tau2x_v, CS%diag) - if (CS%id_FPmask_u > 0) call post_data(CS%id_FPmask_u, mask3d_u, CS%diag) - if (CS%id_FPmask_v > 0) call post_data(CS%id_FPmask_v, mask3d_v, CS%diag) - if (CS%id_FPhbl_u > 0) call post_data(CS%id_FPhbl_u, hbl_u, CS%diag) - if (CS%id_FPhbl_v > 0) call post_data(CS%id_FPhbl_v, hbl_v, CS%diag) - - if (cs%debug) then - call uvchksum("post viscFPmix [ui,vi]",ui,vi,G%HI,haloshift=0,scalar_pair=.true.) + ! GMM; TODO: can you make the arrays used below allocatable? + if(L_diag) then + if (CS%id_tauFP_u > 0) call post_data(CS%id_tauFP_u, tau_u, CS%diag) + if (CS%id_tauFP_v > 0) call post_data(CS%id_tauFP_v, tau_v, CS%diag) + if (CS%id_FPtau2s_u > 0) call post_data(CS%id_FPtau2s_u, omega_tau2s_u, CS%diag) + if (CS%id_FPtau2s_v > 0) call post_data(CS%id_FPtau2s_v, omega_tau2s_v, CS%diag) + if (CS%id_FPtau2w_u > 0) call post_data(CS%id_FPtau2w_u, omega_tau2w_u, CS%diag) + if (CS%id_FPtau2w_v > 0) call post_data(CS%id_FPtau2w_v, omega_tau2w_v, CS%diag) + if (CS%id_FPdiag_u > 0) call post_data(CS%id_FPdiag_u, FPdiag_u, CS%diag) + if (CS%id_FPdiag_v > 0) call post_data(CS%id_FPdiag_v, FPdiag_v, CS%diag) + if (CS%id_FPw2x > 0) call post_data(CS%id_FPw2x, forces%omega_w2x , CS%diag) endif -endif ! LU_pred = false end subroutine vertFPmix +!> Returns the empirical shape-function given sigma. +real function G_sig(sigma) + real , intent(in) :: sigma !< non-dimensional normalized boundary layer depth [m] + + ! local variables + real :: p1, c2, c3 !< parameters used to fit and match empirycal shape-functions. + + ! parabola + p1 = 0.287 + ! cubic function + c2 = 1.74392 + c3 = 2.58538 + G_sig = MIN ( p1 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (c2*sigma - c3) ) ) +end function G_sig + !> Perform a fully implicit vertical diffusion !! of momentum. Stress top and bottom boundary conditions are used. !! @@ -2406,37 +2225,24 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', & thickness_units, conversion=GV%H_to_MKS) - !w FPmix - CS%id_FPhbl_u = register_diag_field('ocean_model', 'FPhbl_u', diag%axesCu1, Time, & - 'Boundary-Layer Depth (u-points)','m') !w , conversion=GV%H_to_MKS) - CS%id_FPhbl_v = register_diag_field('ocean_model', 'FPhbl_v', diag%axesCv1, Time, & - 'Boundary-Layer Depth (v-points)','m') - - CS%id_FPmask_u = register_diag_field('ocean_model', 'FPmask_u', diag%axesCuL, Time, & - 'FP overwrite mask (u-points)','binary') - CS%id_FPmask_v = register_diag_field('ocean_model', 'FPmask_v', diag%axesCvL, Time, & - 'FP overwrite mask (v-points)','binary') - + CS%id_FPw2x = register_diag_field('ocean_model', 'FPw2x', diag%axesT1, Time, & + 'Wind direction from x-axis','radians') + CS%id_FPdiag_u = register_diag_field('ocean_model', 'FPdiag_u', diag%axesCui, Time, & + 'FP diagmostic (u-points)','binary') + CS%id_FPdiag_v = register_diag_field('ocean_model', 'FPdiag_v', diag%axesCvi, Time, & + 'FP diagnostic (v-points)','binary') CS%id_tauFP_u = register_diag_field('ocean_model', 'tauFP_u', diag%axesCui, Time, & - 'Stress Mag Profile (u-points)', 'm2 s-2') !w , conversion=GV%H_to_MKS) + 'Stress Mag Profile (u-points)', 'm2 s-2') CS%id_tauFP_v = register_diag_field('ocean_model', 'tauFP_v', diag%axesCvi, Time, & - 'Stress Mag Profile (v-points)', 'm2 s-2') !w , conversion=GV%H_to_MKS) - + 'Stress Mag Profile (v-points)', 'm2 s-2') CS%id_FPtau2s_u = register_diag_field('ocean_model', 'FPtau2s_u', diag%axesCui, Time, & - 'stress from shear direction (u-points)', 'pi ') + 'stress from shear direction (u-points)', 'radians ') CS%id_FPtau2s_v = register_diag_field('ocean_model', 'FPtau2s_v', diag%axesCvi, Time, & - 'stress from shear direction (v-points)', 'pi ') - + 'stress from shear direction (v-points)', 'radians') CS%id_FPtau2w_u = register_diag_field('ocean_model', 'FPtau2w_u', diag%axesCui, Time, & - 'stress from wind direction (u-points)', 'pi ') + 'stress from wind direction (u-points)', 'radians') CS%id_FPtau2w_v = register_diag_field('ocean_model', 'FPtau2w_v', diag%axesCvi, Time, & - 'stress from wind direction (v-points)', 'pi ') - - CS%id_FPtau2x_u = register_diag_field('ocean_model', 'FPs2w_u', diag%axesCui, Time, & - 'shear from wind (u-points)', 'pi ') - CS%id_FPtau2x_v = register_diag_field('ocean_model', 'FPs2w_v', diag%axesCvi, Time, & - 'shear from wind (v-points)', 'pi ' - ! w - end + 'stress from wind direction (v-points)', 'radians') CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, & 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) From 7624a83b810e21616b089a772398d7d287ca7feb Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 24 May 2022 14:51:31 -0600 Subject: [PATCH 5/7] Add missing use for vertFPmix --- src/core/MOM_dynamics_split_RK2.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index f6cf456f98..b74df389b3 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -64,7 +64,7 @@ module MOM_dynamics_split_RK2 use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant use MOM_vert_friction, only : vertvisc_init, vertvisc_end, vertvisc_CS -use MOM_vert_friction, only : updateCFLtruncationValue +use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units use MOM_wave_interface, only: wave_parameters_CS, Stokes_PGF From 864506e850d5ec4f72457e01bf3dea6305c8eb8c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 24 May 2022 14:56:17 -0600 Subject: [PATCH 6/7] Add omega_w2x to fluxes and forces omega_w2x is the counter-clockwise angle of the wind stress with respect to the horizontal abscissa (x-coordinate) at tracer points [rad]. This variable is needed in the vertPFmix subroutine. --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 3 +++ src/core/MOM_forcing_type.F90 | 24 ++++++++++++++++++- .../vertical/MOM_vert_friction.F90 | 1 + 3 files changed, 27 insertions(+), 1 deletion(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 69841bf84a..41572b969e 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -311,6 +311,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., & press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, & cfc=CS%use_CFC, hevap=CS%enthalpy_cpl) + call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) @@ -721,6 +722,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed) call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed) + call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed) if (CS%rigid_sea_ice) then call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed) @@ -880,6 +882,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)) + forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j)) enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) else ! C-grid wind stresses. diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index d4afabc2de..9d95e7159f 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -67,6 +67,7 @@ module MOM_forcing_type ! surface stress components and turbulent velocity scale real, pointer, dimension(:,:) :: & + omega_w2x => NULL(), & !< the counter-clockwise angle of the wind stress with respect ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1]. ustar_gustless => NULL() !< surface friction velocity scale without any !! any augmentation for gustiness [Z T-1 ~> m s-1]. @@ -221,7 +222,9 @@ module MOM_forcing_type taux => NULL(), & !< zonal wind stress [R L Z T-2 ~> Pa] tauy => NULL(), & !< meridional wind stress [R L Z T-2 ~> Pa] ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1]. - net_mass_src => NULL() !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1] + net_mass_src => NULL(), & !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1] + omega_w2x => NULL() !< the counter-clockwise angle of the wind stress with respect + !! to the horizontal abscissa (x-coordinate) at tracer points [rad]. ! applied surface pressure from other component models (e.g., atmos, sea ice, land ice) real, pointer, dimension(:,:) :: p_surf_full => NULL() @@ -357,6 +360,7 @@ module MOM_forcing_type integer :: id_taux = -1 integer :: id_tauy = -1 integer :: id_ustar = -1 + integer :: id_omega_w2x = -1 integer :: id_psurf = -1 integer :: id_TKE_tidal = -1 @@ -1320,6 +1324,9 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', & 'm s-1', conversion=US%Z_to_m*US%s_to_T) + handles%id_omega_w2x = register_diag_field('ocean_model', 'omega_w2x', diag%axesT1, Time, & + 'Counter-clockwise angle of the wind stress from the horizontal axis.', 'rad') + if (present(use_berg_fluxes)) then if (use_berg_fluxes) then handles%id_ustar_berg = register_diag_field('ocean_model', 'ustar_berg', diag%axesT1, Time, & @@ -2164,6 +2171,11 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) enddo ; enddo endif + if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then + do j=js,je ; do i=is,ie + fluxes%omega_w2x(i,j) = forces%omega_w2x(i,j) + enddo ; enddo + endif if (do_pres) then if (associated(forces%p_surf) .and. associated(fluxes%p_surf)) then do j=js,je ; do i=is,ie @@ -2295,6 +2307,11 @@ subroutine copy_back_forcing_fields(fluxes, forces, G) enddo ; enddo endif + if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then + do j=js,je ; do i=is,ie + forces%omega_w2x(i,j) = fluxes%omega_w2x(i,j) + enddo ; enddo + endif end subroutine copy_back_forcing_fields !> Offer mechanical forcing fields for diagnostics for those @@ -2948,6 +2965,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) & call post_data(handles%id_ustar, fluxes%ustar, diag) + if ((handles%id_omega_w2x > 0) .and. associated(fluxes%omega_w2x)) & + call post_data(handles%id_omega_w2x, fluxes%omega_w2x, diag) + if ((handles%id_ustar_berg > 0) .and. associated(fluxes%ustar_berg)) & call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag) @@ -3264,6 +3284,7 @@ end subroutine myAlloc subroutine deallocate_forcing_type(fluxes) type(forcing), intent(inout) :: fluxes !< Forcing fields structure + if (associated(fluxes%omega_w2x)) deallocate(fluxes%omega_w2x) if (associated(fluxes%ustar)) deallocate(fluxes%ustar) if (associated(fluxes%ustar_gustless)) deallocate(fluxes%ustar_gustless) if (associated(fluxes%buoy)) deallocate(fluxes%buoy) @@ -3325,6 +3346,7 @@ subroutine deallocate_mech_forcing(forces) if (associated(forces%taux)) deallocate(forces%taux) if (associated(forces%tauy)) deallocate(forces%tauy) if (associated(forces%ustar)) deallocate(forces%ustar) + if (associated(forces%omega_w2x)) deallocate(forces%omega_w2x) if (associated(forces%p_surf)) deallocate(forces%p_surf) if (associated(forces%p_surf_full)) deallocate(forces%p_surf_full) if (associated(forces%net_mass_src)) deallocate(forces%net_mass_src) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 1d4f7bf646..605fda5dce 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -25,6 +25,7 @@ module MOM_vert_friction use MOM_variables, only : ocean_internal_state use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS +use MOM_set_visc, only : set_v_at_u, set_u_at_v implicit none ; private #include From 9b4bd84b5e3c7ac1cf67a19670e7197c9ea4cdf5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 21 Jun 2022 15:17:51 -0600 Subject: [PATCH 7/7] Add mssing call to get_param for FPMIX This line of code was lost during the last merge. --- src/core/MOM_dynamics_split_RK2.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index b74df389b3..288d7d9092 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -164,7 +164,7 @@ module MOM_dynamics_split_RK2 !! Euler (1) [nondim]. 0 is often used. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debug_OBC !< If true, do debugging calls for open boundary conditions. - logical :: fpmix !< If true, apply profiles of MTM flux magnitude and direction. + logical :: fpmix !< If true, applies profiles of momentum flux magnitude and direction. logical :: module_is_initialized = .false. !< Record whether this module has been initialized. @@ -327,6 +327,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! eta_pred is the predictor value of the free surface height or column mass, ! [H ~> m or kg m-2]. + ! GMM, TODO: make these allocatable? real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! uold and vold are the velocities before vert_visc is applied. These arrays @@ -1278,6 +1279,9 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param "If true, use the summed layered fluxes plus an "//& "adjustment due to the change in the barotropic velocity "//& "in the barotropic continuity equation.", default=.true.) + call get_param(param_file, mdl, "FPMIX", CS%fpmix, & + "If true, apply profiles of momentum flux magnitude and "//& + " direction", default=.false.) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.)