From b7390f7e04f2e10cbadeb0c7e146c7485981905c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 27 Apr 2022 14:25:43 -0600 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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 06/11] 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 07/11] 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.) From 4b03d23128c364c3d4b5e966cdb85d134d325588 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 30 Aug 2023 14:25:13 -0600 Subject: [PATCH 08/11] Pass wavebands from coupler to wave_parameters_CS This commit passes the waveband information recieved from the coupler to wave_parameters_CS. This information is set to public so that it can be used elsewhere. To exercise this code the following must be set: SURFBAND = COUPLER WAVE_METHOD = SURFACE_BANDS No answer changes. --- src/user/MOM_wave_interface.F90 | 37 +++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index a548436329..02da5a0007 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -98,6 +98,21 @@ module MOM_wave_interface !! Vertical -> Mid-points real, allocatable, dimension(:,:,:), public :: & KvS !< Viscosity for Stokes Drift shear [Z2 T-1 ~> m2 s-1] + real, allocatable, dimension(:), public :: & + WaveNum_Cen !< Wavenumber bands for read/coupled [Z-1 ~> m-1] + real, allocatable, dimension(:,:,:), public :: & + UStk_Hb !< Surface Stokes Drift spectrum (zonal) [L T-1 ~> m s-1] + !! Horizontal -> H-points + !! 3rd dimension -> Freq/Wavenumber + real, allocatable, dimension(:,:,:), public :: & + VStk_Hb !< Surface Stokes Drift spectrum (meridional) [L T-1 ~> m s-1] + !! Horizontal -> H-points + !! 3rd dimension -> Freq/Wavenumber + real, allocatable, dimension(:,:), public :: & + Omega_w2x !< wind direction ccw from model x- axis [nondim radians] + integer, public :: NumBands = 0 !< Number of wavenumber/frequency partitions + !! Must match the number of bands provided + !! via either coupling or file. ! The remainder of this control structure is private integer :: WaveMethod = -99 !< Options for including wave information @@ -149,18 +164,12 @@ module MOM_wave_interface real :: LA_FracHBL !< Fraction of OSBL for averaging Langmuir number [nondim] real :: LA_HBL_min !< Minimum boundary layer depth for averaging Langmuir number [Z ~> m] logical :: LA_Misalignment = .false. !< Flag to use misalignment in Langmuir number - - integer :: NumBands = 0 !< Number of wavenumber/frequency partitions to receive - !! This needs to match the number of bands provided - !! via either coupling or file. real :: g_Earth !< The gravitational acceleration, equivalent to GV%g_Earth but with !! different dimensional rescaling appropriate for deep-water gravity !! waves [Z T-2 ~> m s-2] real :: I_g_Earth !< The inverse of the gravitational acceleration, with dimensional rescaling !! appropriate for deep-water gravity waves [T2 Z-1 ~> s2 m-1] ! Surface Wave Dependent 1d/2d/3d vars - real, allocatable, dimension(:) :: & - WaveNum_Cen !< Wavenumber bands for read/coupled [Z-1 ~> m-1] real, allocatable, dimension(:) :: & Freq_Cen !< Central frequency for wave bands, including a factor of 2*pi [T-1 ~> s-1] real, allocatable, dimension(:) :: & @@ -448,6 +457,9 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag, restar allocate( CS%WaveNum_Cen(CS%NumBands), source=0.0 ) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands), source=0.0 ) allocate( CS%STKy0(G%isd:G%ied,G%jsdB:G%jedB,CS%NumBands), source=0.0 ) + allocate( CS%UStk_Hb(G%isc:G%iec,G%jsc:G%jec,CS%NumBands), source=0.0 ) + allocate( CS%VStk_Hb(G%isc:G%iec,G%jsc:G%jec,CS%NumBands), source=0.0 ) + allocate( CS%Omega_w2x(G%isc:G%iec,G%jsc:G%jec) , source=0.0 ) CS%PartitionMode = 0 call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & "Central wavenumbers for surface Stokes drift bands.", & @@ -463,6 +475,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag, restar allocate( CS%PrescribedSurfStkY(1:CS%NumBands), source=0.0 ) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,1:CS%NumBands), source=0.0 ) allocate( CS%STKy0(G%isd:G%ied,G%jsdB:G%jedB,1:CS%NumBands), source=0.0 ) + CS%PartitionMode = 0 call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & "Central wavenumbers for surface Stokes drift bands.", & @@ -692,6 +705,15 @@ subroutine Update_Surface_Waves(G, GV, US, Time_present, dt, CS, forces) enddo call pass_vector(CS%STKx0(:,:,b),CS%STKy0(:,:,b), G%Domain) enddo + do jj=G%jsc,G%jec + do ii=G%isc,G%iec + CS%Omega_w2x(ii,jj) = forces%omega_w2x(ii,jj) + do b=1,CS%NumBands + CS%UStk_Hb(ii,jj,b) = US%m_s_to_L_T*forces%UStkb(ii,jj,b) + CS%VStk_Hb(ii,jj,b) = US%m_s_to_L_T*forces%VStkb(ii,jj,b) + enddo + enddo + enddo elseif (CS%DataSource == INPUT) then do b=1,CS%NumBands do jj=G%jsd,G%jed @@ -2009,6 +2031,9 @@ subroutine Waves_end(CS) if (allocated(CS%La_turb)) deallocate( CS%La_turb ) if (allocated(CS%STKx0)) deallocate( CS%STKx0 ) if (allocated(CS%STKy0)) deallocate( CS%STKy0 ) + if (allocated(CS%UStk_Hb)) deallocate( CS%UStk_Hb ) + if (allocated(CS%VStk_Hb)) deallocate( CS%VStk_Hb ) + if (allocated(CS%Omega_w2x)) deallocate( CS%Omega_w2x ) if (allocated(CS%KvS)) deallocate( CS%KvS ) if (allocated(CS%Us0_y)) deallocate( CS%Us0_y ) if (allocated(CS%Us0_x)) deallocate( CS%Us0_x ) From 7b7052e9b691468e650686d9ba38dade5c5b4ed5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 6 Sep 2023 13:55:23 -0600 Subject: [PATCH 09/11] Describe local variables and make code consistent --- .../vertical/MOM_vert_friction.F90 | 169 ++++++++---------- 1 file changed, 75 insertions(+), 94 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 5a62e835f8..eab2f2d29d 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -166,7 +166,7 @@ module MOM_vert_friction integer :: id_au_vv = -1, id_av_vv = -1, id_au_gl90_vv = -1, id_av_gl90_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_FPdiag_u = -1, id_FPdiag_v = -1 , id_FPw2x = -1 !W id_FPhbl_u = -1, id_FPhbl_v = -1 + integer :: 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 @@ -210,47 +210,40 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure ! local variables - ! 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(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)) :: omega_w2x_u - real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v - - ! 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_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 :: 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 :: 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 + real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth at u-pts [H ~> m] + real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth at v-pts [H ~> m] + integer, dimension(SZIB_(G),SZJ_(G)) :: kbl_u !< index of the BLD at u-pts [nondim] + integer, dimension(SZI_(G),SZJB_(G)) :: kbl_v !< index of the BLD at v-pts [nondim] + real, dimension(SZIB_(G),SZJ_(G)) :: ustar2_u !< ustar squared at u-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: ustar2_v !< ustar squared at v-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< zonal wind stress at u-pts [R L Z T-2 ~> Pa] + real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< meridional wind stress at v-pts [R L Z T-2 ~> Pa] + real, dimension(SZIB_(G),SZJ_(G)) :: omega_w2x_u !< angle between wind and x-axis at u-pts [rad] + real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v !< angle between wind and y-axis at v-pts [rad] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tau_u !< kinematic zonal mtm flux at u-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tau_v !< kinematic mer. mtm flux at v-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauxDG_u !< downgradient zonal mtm flux at u-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauyDG_u !< downgradient meri mtm flux at u-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauxDG_v !< downgradient zonal mtm flux at v-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauyDG_v !< downgradient meri mtm flux at v-pts [L2 T-2 ~> m2 s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2s_u !< angle between mtm flux and vert shear at u-pts [rad] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2s_v !< angle between mtm flux and vert shear at v-pts [rad] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u !< angle between mtm flux and wind at u-pts [rad] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v !< angle between mtm flux and wind at v-pts [rad] + + real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp, omega_tmp !< constants and dummy variables + real :: du, dv, depth, sigma, Wind_x, Wind_y !< intermediate variables + real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables + real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables + real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles + integer :: kblmin, kbld, kp1, k, nz !< vertical indices + integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq ! horizontal indices 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.) Cemp_CG = 3.6 kblmin = 1 - FPdiag_u(:,:,:) = 0.0 - FPdiag_v(:,:,:) = 0.0 taux_u(:,:) = 0. tauy_v(:,:) = 0. @@ -292,7 +285,7 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US 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 + if( (depth >= hbl_u(I,j)) .and. (kbl_u(I,j) == 0 ) .and. (k > (kblmin-1)) ) then kbl_u(I,j) = k hbl_u(I,j) = depth endif @@ -303,18 +296,18 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US do J = Jsq,Jeq do i = is,ie if( (G%mask2dCv(i,J) > 0.5) ) then - 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 ) + 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) - if( (depth .ge. hbl_v(i,J)) .and. (kbl_v(i,J) .eq. 0 ) .and. (k > (kblmin-1)) ) then + if( (depth >= hbl_v(i,J)) .and. (kbl_v(i,J) == 0) .and. (k > (kblmin-1))) then kbl_v(i,J) = k hbl_v(i,J) = depth endif @@ -331,7 +324,7 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US ! Compute downgradient stresses do k = 1, nz - kp1 = MIN( k+1 , nz) + 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)) @@ -376,7 +369,7 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US tau_u(:,:,:) = 0.0 tau_v(:,:,:) = 0.0 - !w Default implicit (I) stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] Profiles + ! stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] do j = js,je do I = Isq,Ieq if( (G%mask2dCu(I,j) > 0.5) ) then @@ -386,7 +379,6 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US Omega_tau2w_u(I,j,1) = 0.0 Omega_tau2s_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)) @@ -409,7 +401,6 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US Omega_tau2w_v(i,J,1) = 0.0 Omega_tau2s_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) ) @@ -429,9 +420,8 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US do j = js,je do I = Isq,Ieq if( (G%mask2dCu(I,j) > 0.5) ) then - kbld = MIN( (kbl_u(I,j)) , (nz-2) ) + 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 ! surface boundary conditions @@ -439,7 +429,7 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US tauNLup = 0.0 do k=1, kbld depth = depth + CS%h_u(I,j,k) - sigma = MIN ( 1.0 , depth / hbl_u(i,j) ) + sigma = min( 1.0 , depth / hbl_u(i,j) ) ! linear stress mag tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) @@ -449,32 +439,31 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US ! 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 ) + 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) + 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 + 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) + 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 + 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_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 >= 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 @@ -490,8 +479,8 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US 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 + 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) !surface boundary conditions @@ -499,27 +488,27 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US tauNLup = 0.0 do k=1, kbld depth = depth + CS%h_v(i,J,k) - sigma = MIN ( 1.0 , (depth ) / hbl_v(I,J) ) + sigma = min(1.0, depth/ hbl_v(I,J)) ! linear stress - tau_MAG = (ustar2_v(i,J) * (1.-sigma) ) + (tauh * sigma ) + 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 ) + 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 ) + 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 ) + 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 @@ -527,12 +516,11 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US 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_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 > 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 @@ -549,17 +537,14 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US call uvchksum("FP-tau_[uv] ", tau_u, tau_v, G%HI, haloshift=0, scalar_pair=.true.) endif - ! 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_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) + if (CS%id_FPw2x > 0) call post_data(CS%id_FPw2x, forces%omega_w2x , CS%diag) endif end subroutine vertFPmix @@ -576,7 +561,7 @@ real function G_sig(sigma) ! cubic function c2 = 1.74392 c3 = 2.58538 - G_sig = MIN ( p1 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (c2*sigma - c3) ) ) + G_sig = min( p1 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (c2*sigma - c3) ) ) end function G_sig !> Compute coupling coefficient associated with vertical viscosity parameterization as in Greatbatch and Lamb @@ -2875,10 +2860,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 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') CS%id_tauFP_v = register_diag_field('ocean_model', 'tauFP_v', diag%axesCvi, Time, & From 66fd876af9f702a1fc4f956ca07474613956e447 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 6 Sep 2023 15:29:57 -0600 Subject: [PATCH 10/11] Removed L_diag and moved variables in vertFPmix --- src/core/MOM_dynamics_split_RK2.F90 | 7 +--- .../vertical/MOM_vert_friction.F90 | 41 +++++++++---------- 2 files changed, 22 insertions(+), 26 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 84c84efe39..df28dc0338 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -384,7 +384,6 @@ 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]. real :: Idt_bc ! Inverse of the baroclinic timestep [T-1 ~> s-1] - 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 @@ -696,12 +695,11 @@ 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 - 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(L_diag, up, vp, uold, vold, hbl, h, forces, & + call vertFPmix(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) @@ -947,8 +945,7 @@ 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 - L_diag = .true. - call vertFPmix(L_diag, u, v, uold, vold, hbl, h, forces, dt, & + call vertFPmix(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) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index aa21f8ab89..1169126c1c 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -191,26 +191,26 @@ module MOM_vert_friction contains !> 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 +subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OBC) + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: ui !< Zonal velocity after vertvisc [L T-1 ~> m s-1] + 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)), intent(inout) :: hbl_h !< boundary layer depth [H ~> m] 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(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(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces + real, intent(in) :: dt !< Time increment [T ~> s] + 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 type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure ! local variables real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth at u-pts [H ~> m] @@ -241,6 +241,7 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles integer :: kblmin, kbld, kp1, k, nz !< vertical indices integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq ! horizontal indices + 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 @@ -321,8 +322,8 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US 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.) + 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 @@ -540,15 +541,13 @@ subroutine vertFPmix(L_diag, ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US call uvchksum("FP-tau_[uv] ", tau_u, tau_v, G%HI, haloshift=0, scalar_pair=.true.) endif - 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_FPw2x > 0) call post_data(CS%id_FPw2x, forces%omega_w2x , CS%diag) - endif + 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_FPw2x > 0) call post_data(CS%id_FPw2x, forces%omega_w2x , CS%diag) end subroutine vertFPmix From d9aa751a46b67c0d496b9baab28549b2fc679c8f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 6 Sep 2023 16:12:36 -0600 Subject: [PATCH 11/11] Revert order of variables in vertFPmix --- src/parameterizations/vertical/MOM_vert_friction.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 1169126c1c..f513f50158 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -192,7 +192,8 @@ module MOM_vert_friction !> Add nonlocal stress increments to u^n (uold) and v^n (vold) using ui and vi. subroutine vertFPmix(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 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)), & @@ -206,11 +207,9 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces real, intent(in) :: dt !< Time increment [T ~> s] - 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 - type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure - type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure + type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure ! local variables real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth at u-pts [H ~> m]