From 3f57d75f39ea82f33270157e75618c8f671e9524 Mon Sep 17 00:00:00 2001 From: Nora Loose Date: Mon, 9 Jan 2023 07:16:09 -0700 Subject: [PATCH 1/9] Add GL90 diagnostics (#293) * Add GL90 parameterization in stacked shallow water This adds a new vertical viscosity parameterization as in Greatbatch and Lamb (1990), Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization, but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient nu is computed from kappa_GM via thermal wind balance, and the following relation: nu = kappa_GM * f^2 / N^2. The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free boundary conditions at the top and bottom. In the current implementation, kappa_GM is assumed either (a) constant or as (b) having an EBT structure. A third possible formulation of nu is depth-independent: nu = f^2 * alpha The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth. Currently, the GL90 parameterization is only implemented in stacked shallow water (SSW) mode, in which case we have 1/N^2 = h/g'. More specifically, this commit adds a new subroutine that computes the couping coefficient associated with GL90 via a_cpl_gl90 = nu / h = kappa_GM * f^2 / g' or a_cpl_gl90 = nu / h = f^2 * alpha / h. Further, a_cpl_gl90 is multiplied by a function (botfn), which is 0 within the GL90 bottom boundary layer, whose depth is set by Hbbl_gl90, and 1 otherwise. This modification is necessary to avlid fluxing momentum into vanished layers that ride over steep topography. Finally, a_cpl_gl90 is added to a_cpl, where the latter is the coupling coefficient associated with the remaining vertical stresses, used in the vertical viscosity solver. More information can be found in Loose et al. (https://www.essoar.org/doi/abs/10.1002/essoar.10512867.1), Appendix B. New diagnostics: * au_gl90_visc: zonal viscous coupling coefficient associated with GL90, is contained in au_visc * av_gl90_visc: meridional viscous coupling coefficient associated with GL90, is contained in av_visc * Kv_gl90_u: GL90 vertical viscosity at u-points, is contained in Kv_u * Kv_gl90_v: GL90 vertical viscosity at v-points, is contained in Kv_v * du_dt_visc_gl90: zonal acceleration due to GL90 vertical viscosity, included in du_dt_visc * dv_dt_visc_gl90: meridional acceleration due to GL90 vertical viscosity, included in dv_dt_visc * GLwork: Kinetic Energy Source from GL90 Vertical Viscosity The energetics of the GL90 parameterization (named "GLwork") are intentionally computed in MOM_vert_friction, rather than in MOM_diagnostics, where the reamining kinetic energy budget terms are computed. We have to do the computation in MOM_vert_friction to ensure sign- definiteness when GLwork is summed in the vertical. Indeed, MOM_diagnostics does not have access to the velocities and thicknesses used in the vertical solver, but rather uses a time-mean barotropic transport [uv]h to compute the energy budget diagnostics. A detailed discussion and exploration of this issue can be found in https://github.com/ocean-eddy-cpt/MOM6/issues/25. As a result of not computing the energetics in MOM_diagnostics, GLwork is not exactly contained in KE_visc. KE_visc represents the energetics of all vertical viscosity contributions, including the GL90 vertical viscosity. We could implement a term "KE_visc_gl90" that can be 1-to-1 compared to KE_visc; that is, KE_visc - KE_visc_gl90 would represent exactly the energetics of all viscosity contributions EXCEPT the GL90 viscosity. If we implemented KE_visc_gl90, this term would in practice be very similar as GLwork, but sign-definiteness is not ensured, see above. --- src/core/MOM_variables.F90 | 4 + .../vertical/MOM_vert_friction.F90 | 184 +++++++++++++++++- 2 files changed, 182 insertions(+), 6 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 8279afa954..35cdf3038a 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -168,6 +168,10 @@ module MOM_variables PFv => NULL(), & !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2] du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2] dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2] + du_dt_visc_gl90 => NULL(), &!< Zonal acceleration due to GL90 vertical viscosity + ! (is included in du_dt_visc) [L T-2 ~> m s-2] + dv_dt_visc_gl90 => NULL(), &!< Meridional acceleration due to GL90 vertical viscosity + ! (is included in dv_dt_visc) [L T-2 ~> m s-2] du_dt_str => NULL(), & !< Zonal acceleration due to the surface stress (included !! in du_dt_visc) [L T-2 ~> m s-2] dv_dt_str => NULL(), & !< Meridional acceleration due to the surface stress (included diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index bc8ef7e893..88be824885 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -7,6 +7,8 @@ module MOM_vert_friction use MOM_diag_mediator, only : post_product_u, post_product_sum_u use MOM_diag_mediator, only : post_product_v, post_product_sum_v use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled +use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type +use MOM_domains, only : To_North, To_East use MOM_debugging, only : uvchksum, hchksum use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -156,11 +158,14 @@ module MOM_vert_friction real, allocatable, dimension(:,:) :: kappa_gl90_2d !< 2D kappa_gl90 at h-points [L2 T-1 ~> m2 s-1] !>@{ Diagnostic identifiers - integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 + integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_du_dt_visc_gl90 = -1, id_dv_dt_visc_gl90 = -1 + integer :: id_GLwork = -1 + 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_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 + integer :: id_Kv_gl90_u = -1, id_Kv_gl90_v = -1 ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1 integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1 integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1 @@ -171,6 +176,7 @@ module MOM_vert_friction type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure !! for recording accelerations leading to velocity truncations + type(group_pass_type) :: pass_KE_uv !< A handle used for group halo passes end type vertvisc_CS contains @@ -359,6 +365,12 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real :: zDS, hfr, h_a ! Temporary variables used with direct_stress. real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1]. + real, allocatable, dimension(:,:,:) :: KE_term ! A term in the kinetic energy budget + ! [H L2 T-3 ~> m3 s-3 or W m-2] + real, allocatable, dimension(:,:,:) :: KE_u ! The area integral of a KE term in a layer at u-points + ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] + real, allocatable, dimension(:,:,:) :: KE_v ! The area integral of a KE term in a layer at v-points + ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] logical :: do_i(SZIB_(G)) logical :: DoStokesMixing @@ -373,6 +385,14 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(visc): "// & "Module must be initialized before it is used.") + if (CS%id_GLwork > 0) then + allocate(KE_u(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0) + allocate(KE_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0) + allocate(KE_term(G%isd:G%ied,G%jsd:G%jed,GV%ke), source=0.0) + if (.not.G%symmetric) & + call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East) + endif + if (CS%direct_stress) then Hmix = CS%Hmix_stress I_Hmix = 1.0 / Hmix @@ -412,7 +432,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = u(I,j,k) enddo ; enddo ; endif - + if (associated(ADp%du_dt_visc_gl90)) then ; do k=1,nz ; do I=Isq,Ieq + ADp%du_dt_visc_gl90(I,j,k) = u(I,j,k) + enddo ; enddo ; endif if (associated(ADp%du_dt_str)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_str(I,j,k) = 0.0 enddo ; enddo ; endif @@ -501,6 +523,46 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & endif ; enddo ; enddo endif + ! compute vertical velocity tendency that arises from GL90 viscosity; + ! follow tridiagonal solve method as above; to avoid corrupting u, + ! use ADp%du_dt_visc_gl90 as a placeholder for updated u (due to GL90) until last do loop + if ((CS%id_du_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then + if (associated(ADp%du_dt_visc_gl90)) then + do I=Isq,Ieq ; if (do_i(I)) then + b_denom_1 = CS%h_u(I,j,1) ! CS%a_u_gl90(I,j,1) is zero + b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_u_gl90(I,j,2)) + d1(I) = b_denom_1 * b1(I) + ADp%du_dt_visc_gl90(I,j,1) = b1(I) * (CS%h_u(I,j,1) * ADp%du_dt_visc_gl90(I,j,1)) + endif ; enddo + do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then + c1(I,k) = dt_Z_to_H * CS%a_u_gl90(I,j,K) * b1(I) + b_denom_1 = CS%h_u(I,j,k) + dt_Z_to_H * (CS%a_u_gl90(I,j,K)*d1(I)) + b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_u_gl90(I,j,K+1)) + d1(I) = b_denom_1 * b1(I) + ADp%du_dt_visc_gl90(I,j,k) = (CS%h_u(I,j,k) * ADp%du_dt_visc_gl90(I,j,k) + & + dt_Z_to_H * CS%a_u_gl90(I,j,K) * ADp%du_dt_visc_gl90(I,j,k-1)) * b1(I) + endif ; enddo ; enddo + ! back substitute to solve for new velocities, held by ADp%du_dt_visc_gl90 + do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then + ADp%du_dt_visc_gl90(I,j,k) = ADp%du_dt_visc_gl90(I,j,k) + c1(I,k+1) * ADp%du_dt_visc_gl90(I,j,k+1) + endif ; enddo ; enddo ! i and k loops + do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) then + ! now fill ADp%du_dt_visc_gl90(I,j,k) with actual velocity tendency due to GL90; + ! note that on RHS: ADp%du_dt_visc(I,j,k) holds the original velocity value u(I,j,k) + ! and ADp%du_dt_visc_gl90(I,j,k) the updated velocity due to GL90 + ADp%du_dt_visc_gl90(I,j,k) = (ADp%du_dt_visc_gl90(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt + if (abs(ADp%du_dt_visc_gl90(I,j,k)) < accel_underflow) ADp%du_dt_visc_gl90(I,j,k) = 0.0 + endif ; enddo ; enddo ; + ! to compute energetics, we need to multiply by u*h, where u is original velocity before + ! velocity update; note that ADp%du_dt_visc(I,j,k) holds the original velocity value u(I,j,k) + if (CS%id_GLwork > 0) then + do k=1,nz; do I=Isq,Ieq ; if (do_i(I)) then + KE_u(I,j,k) = ADp%du_dt_visc(I,j,k) * CS%h_u(I,j,k) * G%areaCu(I,j) * ADp%du_dt_visc_gl90(I,j,k) + endif ; enddo ; enddo + endif + endif + endif + if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt if (abs(ADp%du_dt_visc(I,j,k)) < accel_underflow) ADp%du_dt_visc(I,j,k) = 0.0 @@ -542,7 +604,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = v(i,J,k) enddo ; enddo ; endif - + if (associated(ADp%dv_dt_visc_gl90)) then ; do k=1,nz ; do i=is,ie + ADp%dv_dt_visc_gl90(i,J,k) = v(i,J,k) + enddo ; enddo ; endif if (associated(ADp%dv_dt_str)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_str(i,J,k) = 0.0 enddo ; enddo ; endif @@ -601,6 +665,47 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & endif ; enddo ; enddo endif + ! compute vertical velocity tendency that arises from GL90 viscosity; + ! follow tridiagonal solve method as above; to avoid corrupting v, + ! use ADp%dv_dt_visc_gl90 as a placeholder for updated u (due to GL90) until last do loop + if ((CS%id_dv_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then + if (associated(ADp%dv_dt_visc_gl90)) then + do i=is,ie ; if (do_i(i)) then + b_denom_1 = CS%h_v(i,J,1) ! CS%a_v_gl90(i,J,1) is zero + b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_v_gl90(i,J,2)) + d1(i) = b_denom_1 * b1(i) + ADp%dv_dt_visc_gl90(I,J,1) = b1(i) * (CS%h_v(i,J,1) * ADp%dv_dt_visc_gl90(i,J,1)) + endif ; enddo + do k=2,nz ; do i=is,ie ; if (do_i(i)) then + c1(i,k) = dt_Z_to_H * CS%a_v_gl90(i,J,K) * b1(i) + b_denom_1 = CS%h_v(i,J,k) + dt_Z_to_H * (CS%a_v_gl90(i,J,K)*d1(i)) + b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_v_gl90(i,J,K+1)) + d1(i) = b_denom_1 * b1(i) + ADp%dv_dt_visc_gl90(i,J,k) = (CS%h_v(i,J,k) * ADp%dv_dt_visc_gl90(i,J,k) + & + dt_Z_to_H * CS%a_v_gl90(i,J,K) * ADp%dv_dt_visc_gl90(i,J,k-1)) * b1(i) + endif ; enddo ; enddo + ! back substitute to solve for new velocities, held by ADp%dv_dt_visc_gl90 + do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then + ADp%dv_dt_visc_gl90(i,J,k) = ADp%dv_dt_visc_gl90(i,J,k) + c1(i,k+1) * ADp%dv_dt_visc_gl90(i,J,k+1) + endif ; enddo ; enddo ! i and k loops + do k=1,nz ; do i=is,ie ; if (do_i(i)) then + ! now fill ADp%dv_dt_visc_gl90(i,J,k) with actual velocity tendency due to GL90; + ! note that on RHS: ADp%dv_dt_visc(i,J,k) holds the original velocity value v(i,J,k) + ! and ADp%dv_dt_visc_gl90(i,J,k) the updated velocity due to GL90 + ADp%dv_dt_visc_gl90(i,J,k) = (ADp%dv_dt_visc_gl90(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt + if (abs(ADp%dv_dt_visc_gl90(i,J,k)) < accel_underflow) ADp%dv_dt_visc_gl90(i,J,k) = 0.0 + endif ; enddo ; enddo ; + ! to compute energetics, we need to multiply by v*h, where u is original velocity before + ! velocity update; note that ADp%dv_dt_visc(I,j,k) holds the original velocity value v(i,J,k) + if (CS%id_GLwork > 0) then + do k=1,nz ; do i=is,ie ; if (do_i(i)) then + ! note that on RHS: ADp%dv_dt_visc(I,j,k) holds the original velocity value v(I,j,k) + KE_v(I,j,k) = ADp%dv_dt_visc(i,J,k) * CS%h_v(i,J,k) * G%areaCv(i,J) * ADp%dv_dt_visc_gl90(i,J,k) + endif ; enddo ; enddo + endif + endif + endif + if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt if (abs(ADp%dv_dt_visc(i,J,k)) < accel_underflow) ADp%dv_dt_visc(i,J,k) = 0.0 @@ -626,6 +731,23 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & enddo ! end of v-component J loop + ! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3]. + ! We do the KE-rate calculation here (rather than in MOM_diagnostics) to ensure + ! a sign-definite term. MOM_diagnostics does not have access to the velocities + ! and thicknesses used in the vertical solver, but rather uses a time-mean + ! barotropic transport [uv]h. + if (CS%id_GLwork > 0) then + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do k=1,nz + do j=js,je ; do i=is,ie + KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j,k) + KE_u(I-1,j,k) + KE_v(i,J,k) + KE_v(i,J-1,k)) + enddo ; enddo + enddo + call post_data(CS%id_GLwork, KE_term, CS%diag) + endif + call vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS) ! Here the velocities associated with open boundary conditions are applied. @@ -651,8 +773,12 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (query_averaging_enabled(CS%diag)) then if (CS%id_du_dt_visc > 0) & call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag) + if (CS%id_du_dt_visc_gl90 > 0) & + call post_data(CS%id_du_dt_visc_gl90, ADp%du_dt_visc_gl90, CS%diag) if (CS%id_dv_dt_visc > 0) & call post_data(CS%id_dv_dt_visc, ADp%dv_dt_visc, CS%diag) + if (CS%id_dv_dt_visc_gl90 > 0) & + call post_data(CS%id_dv_dt_visc_gl90, ADp%dv_dt_visc_gl90, CS%diag) if (present(taux_bot) .and. (CS%id_taux_bot > 0)) & call post_data(CS%id_taux_bot, taux_bot, CS%diag) if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) & @@ -868,6 +994,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC, VarMix) real, allocatable, dimension(:,:) :: hML_v ! Diagnostic of the mixed layer depth at v points [H ~> m or kg m-2]. real, allocatable, dimension(:,:,:) :: Kv_u !< Total vertical viscosity at u-points [Z2 T-1 ~> m2 s-1]. real, allocatable, dimension(:,:,:) :: Kv_v !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1]. + real, allocatable, dimension(:,:,:) :: Kv_gl90_u !< GL90 vertical viscosity at u-points [Z2 T-1 ~> m2 s-1]. + real, allocatable, dimension(:,:,:) :: Kv_gl90_v !< GL90 vertical viscosity at v-points [Z2 T-1 ~> m2 s-1]. real :: zcol(SZI_(G)) ! The height of an interface at h-points [H ~> m or kg m-2]. real :: botfn ! A function which goes from 1 at the bottom to 0 much more ! than Hbbl into the interior [nondim]. @@ -911,6 +1039,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC, VarMix) if (CS%id_Kv_v > 0) allocate(Kv_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0) + if (CS%id_Kv_gl90_u > 0) allocate(Kv_gl90_u(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0) + + if (CS%id_Kv_gl90_v > 0) allocate(Kv_gl90_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0) + if (CS%debug .or. (CS%id_hML_u > 0)) allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0) if (CS%debug .or. (CS%id_hML_v > 0)) allocate(hML_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0) @@ -1106,7 +1238,12 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC, VarMix) if (do_i(I)) Kv_u(I,j,k) = 0.5 * GV%H_to_Z*(CS%a_u(I,j,K)+CS%a_u(I,j,K+1)) * CS%h_u(I,j,k) enddo ; enddo endif - + ! Diagnose GL90 Kv at u-points + if (CS%id_Kv_gl90_u > 0) then + do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) Kv_gl90_u(I,j,k) = 0.5 * GV%H_to_Z*(CS%a_u_gl90(I,j,K)+CS%a_u_gl90(I,j,K+1)) * CS%h_u(I,j,k) + enddo ; enddo + endif enddo @@ -1297,7 +1434,12 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC, VarMix) if (do_i(I)) Kv_v(i,J,k) = 0.5 * GV%H_to_Z*(CS%a_v(i,J,K)+CS%a_v(i,J,K+1)) * CS%h_v(i,J,k) enddo ; enddo endif - + ! Diagnose GL90 Kv at v-points + if (CS%id_Kv_gl90_v > 0) then + do k=1,nz ; do i=is,ie + if (do_i(I)) Kv_gl90_v(i,J,k) = 0.5 * GV%H_to_Z*(CS%a_v_gl90(i,J,K)+CS%a_v_gl90(i,J,K+1)) * CS%h_v(i,J,k) + enddo ; enddo + endif enddo ! end of v-point j loop if (CS%debug) then @@ -1316,8 +1458,12 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC, VarMix) call post_data(CS%id_Kv_slow, visc%Kv_slow, CS%diag) if (CS%id_Kv_u > 0) call post_data(CS%id_Kv_u, Kv_u, CS%diag) if (CS%id_Kv_v > 0) call post_data(CS%id_Kv_v, Kv_v, CS%diag) + if (CS%id_Kv_gl90_u > 0) call post_data(CS%id_Kv_gl90_u, Kv_gl90_u, CS%diag) + if (CS%id_Kv_gl90_v > 0) call post_data(CS%id_Kv_gl90_v, Kv_gl90_v, CS%diag) if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag) if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag) + if (CS%id_au_gl90_vv > 0) call post_data(CS%id_au_gl90_vv, CS%a_u_gl90, CS%diag) + if (CS%id_av_gl90_vv > 0) call post_data(CS%id_av_gl90_vv, CS%a_v_gl90, CS%diag) if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag) if (CS%id_h_v > 0) call post_data(CS%id_h_v, CS%h_v, CS%diag) if (CS%id_hML_u > 0) call post_data(CS%id_hML_u, hML_u, CS%diag) @@ -2292,12 +2438,24 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & CS%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, Time, & 'Total vertical viscosity at v-points', 'm2 s-1', conversion=US%Z2_T_to_m2_s) + CS%id_Kv_gl90_u = register_diag_field('ocean_model', 'Kv_gl90_u', diag%axesCuL, Time, & + 'GL90 vertical viscosity at u-points', 'm2 s-1', conversion=US%Z2_T_to_m2_s) + + CS%id_Kv_gl90_v = register_diag_field('ocean_model', 'Kv_gl90_v', diag%axesCvL, Time, & + 'GL90 vertical viscosity at v-points', 'm2 s-1', conversion=US%Z2_T_to_m2_s) + CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, Time, & 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, Time, & 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=US%Z_to_m*US%s_to_T) + CS%id_au_gl90_vv = register_diag_field('ocean_model', 'au_gl90_visc', diag%axesCui, Time, & + 'Zonal Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=US%Z_to_m*US%s_to_T) + + CS%id_av_gl90_vv = register_diag_field('ocean_model', 'av_gl90_visc', diag%axesCvi, Time, & + 'Meridional Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=US%Z_to_m*US%s_to_T) + CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, Time, & 'Thickness at Zonal Velocity Points for Viscosity', & thickness_units, conversion=GV%H_to_MKS) @@ -2322,7 +2480,21 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, Time, & 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_dv_dt_visc > 0) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) - + CS%id_GLwork = register_diag_field('ocean_model', 'GLwork', diag%axesTL, Time, & + 'Kinetic Energy Source from GL90 Vertical Viscosity', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_du_dt_visc_gl90 = register_diag_field('ocean_model', 'du_dt_visc_gl90', diag%axesCuL, Time, & + 'Zonal Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + if ((CS%id_du_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then + call safe_alloc_ptr(ADp%du_dt_visc_gl90,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) + endif + CS%id_dv_dt_visc_gl90 = register_diag_field('ocean_model', 'dv_dt_visc_gl90', diag%axesCvL, Time, & + 'Meridional Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + if ((CS%id_dv_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then + call safe_alloc_ptr(ADp%dv_dt_visc_gl90,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + endif CS%id_du_dt_str = register_diag_field('ocean_model', 'du_dt_str', diag%axesCuL, Time, & 'Zonal Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_du_dt_str > 0) call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz) From 1d918b651ab2c42ba7eab540c7f1ed95d7605b44 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Thu, 12 Jan 2023 14:29:36 -0500 Subject: [PATCH 2/9] Updates to ODA driver suggested by @Hallberg-NOAA Issue#277 (#297) These changes follow suggestions from @Hallberg-NOAA for clarity in documentation and code readability. These modifications should not impact existing applications (e.g. SPEAR) which are reliant on this code. --- src/core/MOM.F90 | 2 +- src/ocean_data_assim/MOM_oda_driver.F90 | 135 +++++++++++++----------- 2 files changed, 75 insertions(+), 62 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c0ff15e858..07538cdec2 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1456,7 +1456,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & if (CS%debug) then call MOM_thermo_chksum("Pre-oda ", tv, G, US, haloshift=0) endif - call apply_oda_tracer_increments(US%T_to_s*dtdia, Time_end_thermo, G, GV, tv, h, CS%odaCS) + call apply_oda_tracer_increments(dtdia, Time_end_thermo, G, GV, tv, h, CS%odaCS) if (CS%debug) then call MOM_thermo_chksum("Post-oda ", tv, G, US, haloshift=0) endif diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 4ad11592f9..8a1aab3328 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -103,8 +103,12 @@ module MOM_oda_driver_mod type(domain2d), pointer :: mpp_domain => NULL() !< Pointer to a mpp domain object for DA type(grid_type), pointer :: oda_grid !< local tracer grid real, pointer, dimension(:,:,:) :: h => NULL() ! m or kg m-2] for DA - type(thermo_var_ptrs), pointer :: tv => NULL() !< pointer to thermodynamic variables - type(thermo_var_ptrs), pointer :: tv_bc => NULL() !< pointer to thermodynamic bias correction + real, pointer, dimension(:,:,:) :: T_tend => NULL() ! degC s-1] + real, pointer, dimension(:,:,:) :: S_tend => NULL() ! ppt s-1] + real, pointer, dimension(:,:,:) :: T_bc_tend => NULL() !< The layer temperature tendency due + !! to bias adjustment [C T-1 ~> degC s-1] + real, pointer, dimension(:,:,:) :: S_bc_tend => NULL() !< The layer salinity tendency due + !! to bias adjustment [S T-1 ~> ppt s-1] integer :: ni !< global i-direction grid size integer :: nj !< global j-direction grid size logical :: reentrant_x !< grid is reentrant in the x direction @@ -120,7 +124,7 @@ module MOM_oda_driver_mod integer :: ensemble_id = 0 !< id of the current ensemble member integer, pointer, dimension(:,:) :: ensemble_pelist !< PE list for ensemble members integer, pointer, dimension(:) :: filter_pelist !< PE list for ensemble members - integer :: assim_frequency !< analysis interval in hours + real :: assim_interval !< analysis interval [ T ~> s] ! Profiles local to the analysis domain type(ocean_profile_type), pointer :: Profiles => NULL() !< pointer to linked list of all available profiles type(ocean_profile_type), pointer :: CProfiles => NULL()!< pointer to linked list of current profiles @@ -198,8 +202,15 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) call get_param(PF, mdl, "ASSIM_METHOD", assim_method, & "String which determines the data assimilation method "//& "Valid methods are: \'EAKF\',\'OI\', and \'NO_ASSIM\'", default='NO_ASSIM') - call get_param(PF, mdl, "ASSIM_FREQUENCY", CS%assim_frequency, & - "data assimilation frequency in hours") + call get_param(PF, mdl, "ASSIM_INTERVAL", CS%assim_interval, & + "data assimilation update interval in hours",default=-1.0,units="hours",scale=3600.*US%s_to_T) + if (CS%assim_interval < 0.) then + call get_param(PF, mdl, "ASSIM_FREQUENCY", CS%assim_interval, & + "data assimilation update in hours. This parameter name will \n"//& + "be deprecated in the future. ASSIM_INTERVAL should be used instead.",default=-1.0, & + units="hours",scale=3600.*US%s_to_T) + endif + call get_param(PF, mdl, "USE_REGRIDDING", CS%use_ALE_algorithm , & "If True, use the ALE algorithm (regridding/remapping).\n"//& "If False, use the layered isopycnal algorithm.", default=.false. ) @@ -338,9 +349,9 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) ! assign thicknesses call ALE_initThicknessToCoord(CS%ALE_CS, G, CS%GV, CS%h) endif - allocate(CS%tv) - allocate(CS%tv%T(isd:ied,jsd:jed,CS%GV%ke), source=0.0) - allocate(CS%tv%S(isd:ied,jsd:jed,CS%GV%ke), source=0.0) + + allocate(CS%T_tend(isd:ied,jsd:jed,CS%GV%ke), source=0.0) + allocate(CS%S_tend(isd:ied,jsd:jed,CS%GV%ke), source=0.0) ! call set_axes_info(CS%Grid, CS%GV, CS%US, PF, CS%diag_cs, set_vertical=.true.) ! missing in Feiyu's fork allocate(CS%oda_grid) CS%oda_grid%x => CS%Grid%geolonT @@ -387,9 +398,9 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) call get_external_field_info(CS%INC_CS%T_id,size=fld_sz) CS%INC_CS%fldno = 2 if (CS%nk /= fld_sz(3)) call MOM_error(FATAL,'Increment levels /= ODA levels') - allocate(CS%tv_bc) ! storage for increment - allocate(CS%tv_bc%T(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke), source=0.0) - allocate(CS%tv_bc%S(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke), source=0.0) + + allocate(CS%T_bc_tend(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke), source=0.0) + allocate(CS%S_bc_tend(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke), source=0.0) endif call cpu_clock_end(id_clock_oda_init) @@ -468,17 +479,15 @@ end subroutine set_prior_tracer !> Returns posterior adjustments or full state !!Note that only those PEs associated with an ensemble member receive data -subroutine get_posterior_tracer(Time, CS, h, tv, increment) +subroutine get_posterior_tracer(Time, CS, increment) type(time_type), intent(in) :: Time !< the current model time type(ODA_CS), pointer :: CS !< ocean DA control structure - real, dimension(:,:,:), pointer, optional :: h !< Layer thicknesses [H ~> m or kg m-2] - type(thermo_var_ptrs), pointer, optional :: tv !< A structure pointing to various thermodynamic variables logical, optional, intent(in) :: increment !< True if returning increment only type(ocean_control_struct), pointer :: Ocean_increment=>NULL() integer :: m logical :: get_inc - integer :: seconds_per_hour = 3600. + ! return if not analysis time (retain pointers for h and tv) if (Time < CS%Time .or. CS%assim_method == NO_ASSIM) return @@ -487,7 +496,7 @@ subroutine get_posterior_tracer(Time, CS, h, tv, increment) !! switch to global pelist call set_PElist(CS%filter_pelist) call MOM_mesg('Getting posterior') - if (present(h)) h => CS%h ! get analysis thickness + !! Calculate and redistribute increments to CS%tv right after assimilation !! Retain CS%tv to calculate increments for IAU updates CS%tv_inc otherwise get_inc = .true. @@ -503,30 +512,27 @@ subroutine get_posterior_tracer(Time, CS, h, tv, increment) do m=1,CS%ensemble_size if (get_inc) then call redistribute_array(CS%mpp_domain, Ocean_increment%T(:,:,:,m),& - CS%domains(m)%mpp_domain, CS%tv%T, complete=.true.) + CS%domains(m)%mpp_domain, CS%T_tend, complete=.true.) call redistribute_array(CS%mpp_domain, Ocean_increment%S(:,:,:,m),& - CS%domains(m)%mpp_domain, CS%tv%S, complete=.true.) + CS%domains(m)%mpp_domain, CS%S_tend, complete=.true.) else call redistribute_array(CS%mpp_domain, CS%Ocean_posterior%T(:,:,:,m),& - CS%domains(m)%mpp_domain, CS%tv%T, complete=.true.) + CS%domains(m)%mpp_domain, CS%T_tend, complete=.true.) call redistribute_array(CS%mpp_domain, CS%Ocean_posterior%S(:,:,:,m),& - CS%domains(m)%mpp_domain, CS%tv%S, complete=.true.) + CS%domains(m)%mpp_domain, CS%S_tend, complete=.true.) endif enddo - if (present(tv)) tv => CS%tv - if (present(h)) h => CS%h - !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) - call pass_var(CS%tv%T,CS%domains(CS%ensemble_id)) - call pass_var(CS%tv%S,CS%domains(CS%ensemble_id)) + call pass_var(CS%T_tend,CS%domains(CS%ensemble_id)) + call pass_var(CS%S_tend,CS%domains(CS%ensemble_id)) !convert to a tendency (degC or PSU per second) - CS%tv%T = CS%tv%T / (CS%assim_frequency * seconds_per_hour) - CS%tv%S = CS%tv%S / (CS%assim_frequency * seconds_per_hour) + CS%T_tend = CS%T_tend / (CS%assim_interval) + CS%S_tend = CS%S_tend / (CS%assim_interval) end subroutine get_posterior_tracer @@ -560,21 +566,22 @@ subroutine get_bias_correction_tracer(Time, US, CS) type(ODA_CS), pointer :: CS !< ocean DA control structure ! Local variables - real, allocatable, dimension(:,:,:) :: T_bias ! Temperature biases [C ~> degC] - real, allocatable, dimension(:,:,:) :: S_bias ! Salinity biases [C ~> degC] - real, allocatable, dimension(:,:,:) :: mask_z ! Missing value mask on the horizontal model grid - ! and input-file vertical levels [nondim] + real, allocatable, dimension(:,:,:) :: T_bias ! Estimated temperature tendency bias [C T-1 ~> degC s-1] + real, allocatable, dimension(:,:,:) :: S_bias ! Estimated salinity tendency bias [S T-1 ~> ppt s-1] + real, allocatable, dimension(:,:,:) :: valid_flag ! Valid value flag on the horizontal model grid + ! and input-file vertical levels [nondim] real, allocatable, dimension(:), target :: z_in ! Cell center depths for input data [Z ~> m] real, allocatable, dimension(:), target :: z_edges_in ! Cell edge depths for input data [Z ~> m] real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc] integer, dimension(3) :: fld_sz integer :: i,j,k + call cpu_clock_begin(id_clock_bias_adjustment) call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id, Time, CS%G, T_bias, & - mask_z, z_in, z_edges_in, missing_value, scale=US%degC_to_C, spongeOngrid=.true.) + valid_flag, z_in, z_edges_in, missing_value, scale=US%degC_to_C*US%s_to_T, spongeOngrid=.true.) call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id, Time, CS%G, S_bias, & - mask_z, z_in, z_edges_in, missing_value, scale=US%ppt_to_S, spongeOngrid=.true.) + valid_flag, z_in, z_edges_in, missing_value, scale=US%ppt_to_S*US%s_to_T, spongeOngrid=.true.) ! This should be replaced to use mask_z instead of the following lines ! which are intended to zero land values using an arbitrary limit. @@ -582,17 +589,21 @@ subroutine get_bias_correction_tracer(Time, US, CS) do i=1,fld_sz(1) do j=1,fld_sz(2) do k=1,fld_sz(3) - if (T_bias(i,j,k) > 1.0E-3*US%degC_to_C) T_bias(i,j,k) = 0.0 - if (S_bias(i,j,k) > 1.0E-3*US%ppt_to_S) S_bias(i,j,k) = 0.0 +! if (T_bias(i,j,k) > 1.0E-3*US%degC_to_C) T_bias(i,j,k) = 0.0 +! if (S_bias(i,j,k) > 1.0E-3*US%ppt_to_S) S_bias(i,j,k) = 0.0 + if (valid_flag(i,j,k)==0.) then + T_bias(i,j,k)=0.0 + S_bias(i,j,k)=0.0 + endif enddo enddo enddo - CS%tv_bc%T = T_bias * CS%bias_adjustment_multiplier - CS%tv_bc%S = S_bias * CS%bias_adjustment_multiplier + CS%T_bc_tend = T_bias * CS%bias_adjustment_multiplier + CS%S_bc_tend = S_bias * CS%bias_adjustment_multiplier - call pass_var(CS%tv_bc%T, CS%domains(CS%ensemble_id)) - call pass_var(CS%tv_bc%S, CS%domains(CS%ensemble_id)) + call pass_var(CS%T_bc_tend, CS%domains(CS%ensemble_id)) + call pass_var(CS%S_bc_tend, CS%domains(CS%ensemble_id)) call cpu_clock_end(id_clock_bias_adjustment) @@ -640,8 +651,8 @@ subroutine set_analysis_time(Time,CS) integer :: yr, mon, day, hr, min, sec if (Time >= CS%Time) then - ! increment the analysis time to the next step converting to seconds - CS%Time = CS%Time + real_to_time(CS%US%T_to_s*(CS%assim_frequency*3600.)) + ! increment the analysis time to the next step + CS%Time = CS%Time + real_to_time(CS%US%T_to_s*(CS%assim_interval)) call get_date(Time, yr, mon, day, hr, min, sec) write(mesg,*) 'Model Time: ', yr, mon, day, hr, min, sec @@ -662,7 +673,7 @@ end subroutine set_analysis_time !> Apply increments to tracers subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) - real, intent(in) :: dt !< The tracer timestep [s] + real, intent(in) :: dt !< The tracer timestep [T ~> s] type(time_type), intent(in) :: Time_end !< Time at the end of the interval type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure @@ -674,12 +685,14 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) !! local variables integer :: i, j integer :: isc, iec, jsc, jec - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_inc !< an adjustment to the temperature + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_tend_inc !< an adjustment to the temperature !! tendency [C T-1 -> degC s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_inc !< an adjustment to the salinity + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_tend_inc !< an adjustment to the salinity !! tendency [S T-1 -> ppt s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: T !< The updated temperature [C ~> degC] - real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: S !< The updated salinity [S ~> ppt] + real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: T_tend !< The temperature tendency adjustment from + !! DA [C T-1 ~> degC s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: S_tend !< The salinity tendency adjustment from DA + !! [S T-1 ~> ppt s-1] real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] if (.not. associated(CS)) return @@ -687,14 +700,14 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) call cpu_clock_begin(id_clock_apply_increments) - T_inc(:,:,:) = 0.0; S_inc(:,:,:) = 0.0; T(:,:,:) = 0.0; S(:,:,:) = 0.0 + T_tend_inc(:,:,:) = 0.0; S_tend_inc(:,:,:) = 0.0; T_tend(:,:,:) = 0.0; S_tend(:,:,:) = 0.0 if (CS%assim_method > 0 ) then - T = T + CS%tv%T - S = S + CS%tv%S + T_tend = T_tend + CS%T_tend + S_tend = S_tend + CS%S_tend endif if (CS%do_bias_adjustment ) then - T = T + CS%tv_bc%T - S = S + CS%tv_bc%S + T_tend = T_tend + CS%T_bc_tend + S_tend = S_tend + CS%S_bc_tend endif if (CS%answer_date >= 20190101) then @@ -707,25 +720,25 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) isc=G%isc; iec=G%iec; jsc=G%jsc; jec=G%jec do j=jsc,jec; do i=isc,iec - call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), T(i,j,:), & - G%ke, h(i,j,:), T_inc(i,j,:), h_neglect, h_neglect_edge) - call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), S(i,j,:), & - G%ke, h(i,j,:), S_inc(i,j,:), h_neglect, h_neglect_edge) + call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), T_tend(i,j,:), & + G%ke, h(i,j,:), T_tend_inc(i,j,:), h_neglect, h_neglect_edge) + call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), S_tend(i,j,:), & + G%ke, h(i,j,:), S_tend_inc(i,j,:), h_neglect, h_neglect_edge) enddo; enddo - call pass_var(T_inc, G%Domain) - call pass_var(S_inc, G%Domain) + call pass_var(T_tend_inc, G%Domain) + call pass_var(S_tend_inc, G%Domain) - tv%T(isc:iec,jsc:jec,:) = tv%T(isc:iec,jsc:jec,:) + T_inc(isc:iec,jsc:jec,:)*dt - tv%S(isc:iec,jsc:jec,:) = tv%S(isc:iec,jsc:jec,:) + S_inc(isc:iec,jsc:jec,:)*dt + tv%T(isc:iec,jsc:jec,:) = tv%T(isc:iec,jsc:jec,:) + T_tend_inc(isc:iec,jsc:jec,:)*dt + tv%S(isc:iec,jsc:jec,:) = tv%S(isc:iec,jsc:jec,:) + S_tend_inc(isc:iec,jsc:jec,:)*dt call pass_var(tv%T, G%Domain) call pass_var(tv%S, G%Domain) call enable_averaging(dt, Time_end, CS%diag_CS) - if (CS%id_inc_t > 0) call post_data(CS%id_inc_t, T_inc, CS%diag_CS) - if (CS%id_inc_s > 0) call post_data(CS%id_inc_s, S_inc, CS%diag_CS) + if (CS%id_inc_t > 0) call post_data(CS%id_inc_t, T_tend_inc, CS%diag_CS) + if (CS%id_inc_s > 0) call post_data(CS%id_inc_s, S_tend_inc, CS%diag_CS) call disable_averaging(CS%diag_CS) call diag_update_remap_grids(CS%diag_CS) From 237194dfdeb4d20e907509645e83fee8f5c92472 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 4 Jan 2023 10:23:52 -0500 Subject: [PATCH 3/9] +Add DOME_tracer and tracer_example runtime params Added the new runtime parameters DOME_TRACER_STRIPE_WIDTH, DOME_TRACER_STRIPE_LAT and DOME_TRACER_SHEET_SPACING to specify the previously hard-coded dimensional parameters in the DOME_tracer module, and added the runtime parameters TRACER_EXAMPLE_STRIPE_WIDTH and TRACER_EXAMPLE_STRIPE_LAT to specify parameters used by tracer_example. This change requires that an ocean grid type be passed to register_DOME_tracer and USER_register_tracer_example instead of a hor_index type. Descriptions of the units of a number of internal variables were added to both modules. In addition, the confusing (and dimensionally heterogeneous) trdc array in tracer_column_physics was replaced with 3 internal variables with more suggestive names. By default all answers are bitwise identical, but there are new entries in the MOM_parameter_doc.all files for cases that have USE_DOME_TRACER=True or USE_USER_TRACER_EXAMPLE=True. --- src/tracer/DOME_tracer.F90 | 65 ++++++++++------ src/tracer/MOM_tracer_flow_control.F90 | 7 +- src/tracer/tracer_example.F90 | 104 +++++++++++++------------ 3 files changed, 100 insertions(+), 76 deletions(-) diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index d1c6ebd7bf..98788843e3 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -42,10 +42,18 @@ module DOME_tracer character(len=200) :: tracer_IC_file !< The full path to the IC file, or " " to initialize internally. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry - real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this package, in g m-3? - real :: land_val(NTR) = -1.0 !< The value of tr used where land is masked out. + real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this package, perhaps in [g kg-1] + real :: land_val(NTR) = -1.0 !< The value of tr used where land is masked out, perhaps in [g kg-1] logical :: use_sponge !< If true, sponges may be applied somewhere in the domain. + real :: stripe_width !< The meridional width of the vertical stripes in the initial condition + !! for some of the DOME tracers, in [km] or [degrees_N] or [m]. + real :: stripe_s_lat !< The southern latitude of the first vertical stripe in the initial condition + !! for some of the DOME tracers, in [km] or [degrees_N] or [m]. + real :: sheet_spacing !< The vertical spacing between successive horizontal sheets of tracer in the initial + !! conditions for some of the DOME tracers [Z ~> m], and twice the thickness of + !! these horizontal tracer sheets + integer, dimension(NTR) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and the !! surface tracer concentrations are to be provided to the coupler. @@ -58,14 +66,15 @@ module DOME_tracer contains !> Register tracer fields and subroutines to be used with MOM. -function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) - type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(DOME_tracer_CS), pointer :: CS !< A pointer that is set to point to the +function register_DOME_tracer(G, GV, US, param_file, CS, tr_Reg, restart_CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(DOME_tracer_CS), pointer :: CS !< A pointer that is set to point to the !! control structure for this module type(tracer_registry_type), pointer :: tr_Reg !< A pointer to the tracer registry. - type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct + type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct ! Local variables character(len=80) :: name, longname @@ -75,10 +84,10 @@ function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) character(len=48) :: flux_units ! The units for tracer fluxes, usually ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. character(len=200) :: inputdir - real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to the tracer field + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to one of the tracers, perhaps in [g kg-1] logical :: register_DOME_tracer integer :: isd, ied, jsd, jed, nz, m - isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke if (associated(CS)) then call MOM_error(FATAL, "DOME_register_tracer called with an "// & @@ -99,6 +108,16 @@ function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) call log_param(param_file, mdl, "INPUTDIR/DOME_TRACER_IC_FILE", & CS%tracer_IC_file) endif + call get_param(param_file, mdl, "DOME_TRACER_STRIPE_WIDTH", CS%stripe_width, & + "The meridional width of the vertical stripes in the initial condition "//& + "for the DOME tracers.", units=G%y_ax_unit_short, default=50.0) + call get_param(param_file, mdl, "DOME_TRACER_STRIPE_LAT", CS%stripe_s_lat, & + "The southern latitude of the first vertical stripe in the initial condition "//& + "for the DOME tracers.", units=G%y_ax_unit_short, default=350.0) + call get_param(param_file, mdl, "DOME_TRACER_SHEET_SPACING", CS%sheet_spacing, & + "The vertical spacing between successive horizontal sheets of tracer in the initial "//& + "conditions for the DOME tracers, and twice the thickness of these tracer sheets.", & + units="m", default=600.0, scale=US%m_to_Z) call get_param(param_file, mdl, "SPONGE", CS%use_sponge, & "If true, sponges may be applied anywhere in the domain. "//& "The exact location and properties of those sponges are "//& @@ -118,7 +137,7 @@ function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! calls. Curses on the designers and implementers of Fortran90. tr_ptr => CS%tr(:,:,:,m) ! Register the tracer for horizontal advection, diffusion, and restarts. - call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + call register_tracer(tr_ptr, tr_Reg, param_file, G%HI, GV, & name=name, longname=longname, units="kg kg-1", & registry_diags=.true., restart_CS=restart_CS, & flux_units=trim(flux_units), flux_scale=GV%H_to_MKS) @@ -154,16 +173,16 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters ! Local variables - real, allocatable :: temp(:,:,:) + real, allocatable :: temp(:,:,:) ! Target values for the tracers in the sponges, perhaps in [g kg-1] character(len=16) :: name ! A variable's name in a NetCDF file. - real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to the tracer field - real :: tr_y ! Initial zonally uniform tracer concentrations. + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to one of the tracers, perhaps in [g kg-1] + real :: tr_y ! Initial zonally uniform tracer concentrations, perhaps in [g kg-1] real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: e(SZK_(GV)+1) ! Interface heights relative to the sea surface (negative down) [Z ~> m] real :: e_top ! Height of the top of the tracer band relative to the sea surface [Z ~> m] real :: e_bot ! Height of the bottom of the tracer band relative to the sea surface [Z ~> m] - real :: d_tr ! A change in tracer concentrations, in tracer units. + real :: d_tr ! A change in tracer concentrations, in tracer units, perhaps [g kg-1] integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB @@ -194,24 +213,25 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & enddo ! This sets a stripe of tracer across the basin. - do m=2,NTR ; do j=js,je ; do i=is,ie + do m=2,min(6,NTR) ; do j=js,je ; do i=is,ie tr_y = 0.0 - if ((m <= 6) .and. (G%geoLatT(i,j) > (300.0+50.0*real(m-1))) .and. & - (G%geoLatT(i,j) < (350.0+50.0*real(m-1)))) tr_y = 1.0 + if ((G%geoLatT(i,j) > (CS%stripe_s_lat + CS%stripe_width*real(m-2))) .and. & + (G%geoLatT(i,j) < (CS%stripe_s_lat + CS%stripe_width*real(m-1)))) & + tr_y = 1.0 do k=1,nz ! This adds the stripes of tracer to every layer. CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + tr_y enddo enddo ; enddo ; enddo - if (NTR > 7) then + if (NTR >= 7) then do j=js,je ; do i=is,ie e(1) = 0.0 do k=1,nz e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z do m=7,NTR - e_top = (-600.0*real(m-1) + 3000.0) * US%m_to_Z - e_bot = (-600.0*real(m-1) + 2700.0) * US%m_to_Z + e_top = -CS%sheet_spacing * (real(m-6)) + e_bot = -CS%sheet_spacing * (real(m-6) + 0.5) if (e_top < e(K)) then if (e_top < e(K+1)) then ; d_tr = 0.0 elseif (e_bot < e(K+1)) then @@ -255,8 +275,7 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & ! do m=1,NTR do m=1,1 - ! This is needed to force the compiler not to do a copy in the sponge - ! calls. Curses on the designers and implementers of Fortran90. + ! This pointer is needed to force the compiler not to do a copy in the sponge calls. tr_ptr => CS%tr(:,:,:,m) call set_up_sponge_field(temp, tr_ptr, G, GV, nz, sponge_CSp) enddo diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index bf7076d4bb..aa3e359355 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -178,8 +178,7 @@ subroutine call_tracer_register(G, GV, US, param_file, CS, tr_Reg, restart_CS) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") - call get_param(param_file, mdl, "USE_USER_TRACER_EXAMPLE", & - CS%use_USER_tracer_example, & + call get_param(param_file, mdl, "USE_USER_TRACER_EXAMPLE", CS%use_USER_tracer_example, & "If true, use the USER_tracer_example tracer package.", & default=.false.) call get_param(param_file, mdl, "USE_DOME_TRACER", CS%use_DOME_tracer, & @@ -230,10 +229,10 @@ subroutine call_tracer_register(G, GV, US, param_file, CS, tr_Reg, restart_CS) ! tracer package registration call returns a logical false if it cannot be run ! for some reason. This then overrides the run-time selection from above. if (CS%use_USER_tracer_example) CS%use_USER_tracer_example = & - USER_register_tracer_example(G%HI, GV, param_file, CS%USER_tracer_example_CSp, & + USER_register_tracer_example(G, GV, US, param_file, CS%USER_tracer_example_CSp, & tr_Reg, restart_CS) if (CS%use_DOME_tracer) CS%use_DOME_tracer = & - register_DOME_tracer(G%HI, GV, param_file, CS%DOME_tracer_CSp, & + register_DOME_tracer(G, GV, US, param_file, CS%DOME_tracer_CSp, & tr_Reg, restart_CS) if (CS%use_ISOMIP_tracer) CS%use_ISOMIP_tracer = & register_ISOMIP_tracer(G%HI, GV, param_file, CS%ISOMIP_tracer_CSp, & diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90 index 335f82a59b..fa9b978f9c 100644 --- a/src/tracer/tracer_example.F90 +++ b/src/tracer/tracer_example.F90 @@ -39,8 +39,13 @@ module USER_tracer_example !! to initialize internally. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry - real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this subroutine, in g m-3? - real :: land_val(NTR) = -1.0 !< The value of tr that is used where land is masked out. + real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this subroutine, perhaps in [g kg-1]? + real :: land_val(NTR) = -1.0 !< The value of tr that is used where land is masked out, perhaps in [g kg-1]? + + real :: stripe_width !< The Gaussian width of the stripe in the initial condition + !! for the tracer_example tracers [L ~> m] + real :: stripe_lat !< The central latitude of the stripe in the initial condition + !! for the tracer_example tracers, in [degrees_N] or [km] or [m]. logical :: use_sponge !< If true, sponges may be applied somewhere in the domain. integer, dimension(NTR) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and the @@ -54,16 +59,17 @@ module USER_tracer_example contains !> This subroutine is used to register tracer fields and subroutines to be used with MOM. -function USER_register_tracer_example(HI, GV, param_file, CS, tr_Reg, restart_CS) - type(hor_index_type), intent(in) :: HI !< A horizontal index type structure +function USER_register_tracer_example(G, GV, US, param_file, CS, tr_Reg, restart_CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(USER_tracer_example_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module type(tracer_registry_type), pointer :: tr_Reg !< A pointer that is set to point to the control !! structure for the tracer advection and !! diffusion module - type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct + type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct ! Local variables character(len=80) :: name, longname @@ -73,10 +79,10 @@ function USER_register_tracer_example(HI, GV, param_file, CS, tr_Reg, restart_CS character(len=200) :: inputdir character(len=48) :: flux_units ! The units for tracer fluxes, usually ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. - real, pointer :: tr_ptr(:,:,:) => NULL() + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to one of the tracers, perhaps in [g kg-1] logical :: USER_register_tracer_example integer :: isd, ied, jsd, jed, nz, m - isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke if (associated(CS)) then call MOM_error(FATAL, "USER_register_tracer_example called with an "// & @@ -87,9 +93,9 @@ function USER_register_tracer_example(HI, GV, param_file, CS, tr_Reg, restart_CS ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "TRACER_EXAMPLE_IC_FILE", CS%tracer_IC_file, & - "The name of a file from which to read the initial "//& - "conditions for the DOME tracers, or blank to initialize "//& - "them internally.", default=" ") + "The name of a file from which to read the initial conditions for "//& + "the tracer_example tracers, or blank to initialize them internally.", & + default=" ") if (len_trim(CS%tracer_IC_file) >= 1) then call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") CS%tracer_IC_file = trim(slasher(inputdir))//trim(CS%tracer_IC_file) @@ -100,6 +106,12 @@ function USER_register_tracer_example(HI, GV, param_file, CS, tr_Reg, restart_CS "If true, sponges may be applied anywhere in the domain. "//& "The exact location and properties of those sponges are "//& "specified from MOM_initialization.F90.", default=.false.) + call get_param(param_file, mdl, "TRACER_EXAMPLE_STRIPE_WIDTH", CS%stripe_width, & + "The Gaussian width of the stripe in the initial condition for the "//& + "tracer_example tracers.", units="m", default=1.0e5, scale=US%m_to_L) + call get_param(param_file, mdl, "TRACER_EXAMPLE_STRIPE_LAT", CS%stripe_lat, & + "The central latitude of the stripe in the initial condition for the "//& + "tracer_example tracers.", units=G%y_ax_unit_short, default=40.0) allocate(CS%tr(isd:ied,jsd:jed,nz,NTR), source=0.0) @@ -113,11 +125,10 @@ function USER_register_tracer_example(HI, GV, param_file, CS, tr_Reg, restart_CS if (GV%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1" else ; flux_units = "kg s-1" ; endif - ! This is needed to force the compiler not to do a copy in the registration - ! calls. Curses on the designers and implementers of Fortran90. + ! This pointer is needed to force the compiler not to do a copy in the registration calls. tr_ptr => CS%tr(:,:,:,m) ! Register the tracer for horizontal advection, diffusion, and restarts. - call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + call register_tracer(tr_ptr, tr_Reg, param_file, G%HI, GV, & name=name, longname=longname, units="kg kg-1", & registry_diags=.true., flux_units=flux_units, & restart_CS=restart_CS) @@ -157,11 +168,11 @@ subroutine USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & !! for the sponges, if they are in use. ! Local variables - real, allocatable :: temp(:,:,:) + real, allocatable :: temp(:,:,:) ! Target values for the tracers in the sponges, perhaps in [g kg-1] character(len=32) :: name ! A variable's name in a NetCDF file. - real, pointer :: tr_ptr(:,:,:) => NULL() - real :: PI ! 3.1415926... calculated as 4*atan(1) - real :: tr_y ! Initial zonally uniform tracer concentrations. + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to one of the tracers, perhaps in [g kg-1] + real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] + real :: tr_y ! Initial zonally uniform tracer concentrations, perhaps in [g kg-1] real :: dist2 ! The distance squared from a line [L2 ~> m2]. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB, lntr @@ -195,9 +206,8 @@ subroutine USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & ! This sets a stripe of tracer across the basin. PI = 4.0*atan(1.0) do j=js,je - dist2 = (G%Rad_Earth_L * PI / 180.0)**2 * & - (G%geoLatT(i,j) - 40.0) * (G%geoLatT(i,j) - 40.0) - tr_y = 0.5 * exp( -dist2 / (1.0e5*US%m_to_L)**2 ) + dist2 = (G%Rad_Earth_L * PI / 180.0)**2 * (G%geoLatT(i,j) - CS%stripe_lat)**2 + tr_y = 0.5 * exp( -dist2 / CS%stripe_width**2 ) do k=1,nz ; do i=is,ie ! This adds the stripes of tracer to every layer. @@ -218,7 +228,7 @@ subroutine USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & allocate(temp(G%isd:G%ied,G%jsd:G%jed,nz)) do k=1,nz ; do j=js,je ; do i=is,ie - if (G%geoLatT(i,j) > 700.0 .and. (k > nz/2)) then + if ((G%geoLatT(i,j) > 0.5*G%len_lat + G%south_lat) .and. (k > nz/2)) then temp(i,j,k) = 1.0 else temp(i,j,k) = 0.0 @@ -227,8 +237,7 @@ subroutine USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & ! do m=1,NTR do m=1,1 - ! This is needed to force the compiler not to do a copy in the sponge - ! calls. Curses on the designers and implementers of Fortran90. + ! This pointer is needed to force the compiler not to do a copy in the sponge calls. tr_ptr => CS%tr(:,:,:,m) call set_up_sponge_field(temp, tr_ptr, G, GV, nz, sponge_CSp) enddo @@ -288,28 +297,25 @@ subroutine tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, real :: d1(SZI_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. - real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]. + real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]. + real :: diapyc_filt ! A multiplicative filter that can be set to 0 to disable diapycnal + ! advection of the tracer [nondim] + real :: dye_up ! The tracer concentration of upwelled water, perhaps in [g kg-1]? + real :: dye_down ! The tracer concentration of downwelled water, perhaps in [g kg-1]? integer :: i, j, k, is, ie, js, je, nz, m -! The following array (trdc) determines the behavior of the tracer -! diapycnal advection. The first element is 1 if tracers are -! passively advected. The second and third are the concentrations -! to which downwelling and upwelling water are set, respectively. -! For most (normal) tracers, the appropriate vales are {1,0,0}. - - real :: trdc(3) -! Uncomment the following line to dye both upwelling and downwelling. -! data trdc / 0.0,1.0,1.0 / -! Uncomment the following line to dye downwelling. -! data trdc / 0.0,1.0,0.0 / -! Uncomment the following line to dye upwelling. -! data trdc / 0.0,0.0,1.0 / -! Uncomment the following line for tracer concentrations to be set -! to zero in any diapycnal motions. -! data trdc / 0.0,0.0,0.0 / -! Uncomment the following line for most "physical" tracers, which -! are advected diapycnally in the usual manner. - data trdc / 1.0,0.0,0.0 / + ! These are the settings for most "physical" tracers, which + ! are advected diapycnally in the usual manner. + diapyc_filt = 1.0 ; dye_down = 0.0 ; dye_down = 0.0 + + ! Uncomment the following line to dye downwelling. +! diapyc_filt = 0.0 ; dye_down = 1.0 + ! Uncomment the following line to dye upwelling. +! diapyc_filt = 0.0 ; dye_up = 1.0 + ! Uncomment the following line for tracer concentrations to be set + ! to zero in any diapycnal motions. +! diapyc_filt = 0.0 ; dye_down = 0.0 ; dye_down = 0.0 + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (.not.associated(CS)) return @@ -330,21 +336,21 @@ subroutine tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, b_denom_1 = h_old(i,j,1) + ea(i,j,1) + h_neglect b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) ! d1(i) = b_denom_1 * b1(i) - d1(i) = trdc(1) * (b_denom_1 * b1(i)) + (1.0 - trdc(1)) + d1(i) = diapyc_filt * (b_denom_1 * b1(i)) + (1.0 - diapyc_filt) do m=1,NTR - CS%tr(i,j,1,m) = b1(i)*(hold0(i)*CS%tr(i,j,1,m) + trdc(3)*eb(i,j,1)) + CS%tr(i,j,1,m) = b1(i)*(hold0(i)*CS%tr(i,j,1,m) + dye_up*eb(i,j,1)) ! Add any surface tracer fluxes to the preceding line. enddo enddo do k=2,nz ; do i=is,ie - c1(i,k) = trdc(1) * eb(i,j,k-1) * b1(i) + c1(i,k) = diapyc_filt * eb(i,j,k-1) * b1(i) b_denom_1 = h_old(i,j,k) + d1(i)*ea(i,j,k) + h_neglect b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) - d1(i) = trdc(1) * (b_denom_1 * b1(i)) + (1.0 - trdc(1)) + d1(i) = diapyc_filt * (b_denom_1 * b1(i)) + (1.0 - diapyc_filt) do m=1,NTR CS%tr(i,j,k,m) = b1(i) * (h_old(i,j,k)*CS%tr(i,j,k,m) + & - ea(i,j,k)*(trdc(1)*CS%tr(i,j,k-1,m)+trdc(2)) + & - eb(i,j,k)*trdc(3)) + ea(i,j,k)*(diapyc_filt*CS%tr(i,j,k-1,m) + dye_down) + & + eb(i,j,k)*dye_up) enddo enddo ; enddo do m=1,NTR ; do k=nz-1,1,-1 ; do i=is,ie From 54b07015396d3e508db6eca774829b53d02e1dfe Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 4 Jan 2023 17:33:54 -0500 Subject: [PATCH 4/9] Named argument assignment white-space clean-up Removed extra white space around the named argument assignments for default, units or conversion arguments to align with the MOM6 style guide. Only white space is modified, and all answers are bitwise identical. --- src/core/MOM_dynamics_split_RK2.F90 | 4 ++-- src/core/MOM_dynamics_unsplit.F90 | 4 ++-- src/core/MOM_dynamics_unsplit_RK2.F90 | 4 ++-- src/core/MOM_verticalGrid.F90 | 2 +- src/framework/MOM_domains.F90 | 2 +- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 2 +- src/parameterizations/lateral/MOM_thickness_diffuse.F90 | 2 +- src/parameterizations/lateral/MOM_tidal_forcing.F90 | 2 +- src/parameterizations/vertical/MOM_CVMix_shear.F90 | 2 +- src/tracer/MOM_neutral_diffusion.F90 | 6 +++--- src/user/BFB_surface_forcing.F90 | 2 +- src/user/circle_obcs_initialization.F90 | 2 +- src/user/dumbbell_surface_forcing.F90 | 4 ++-- 14 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 748748f77f..e8909e24f9 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1545,10 +1545,10 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & - 'Effective U-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & + 'Effective U-Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, & y_cell_method='sum', v_extensive=.true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & - 'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & + 'Effective V-Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, & x_cell_method='sum', v_extensive=.true.) if (GV%Boussinesq) then CS%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, Time, & diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index a0a6633811..e6f99cc9d8 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -692,10 +692,10 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & - 'Effective U Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & + 'Effective U Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, & y_cell_method='sum', v_extensive=.true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & - 'Effective V Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & + 'Effective V Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, & x_cell_method='sum', v_extensive=.true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 85f5d6c546..fbf416d13d 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -655,10 +655,10 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & - 'Effective U-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & + 'Effective U-Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, & y_cell_method='sum', v_extensive=.true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & - 'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & + 'Effective V-Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, & x_cell_method='sum', v_extensive=.true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index 41d29488cd..f20c7bbd26 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -105,7 +105,7 @@ subroutine verticalGridInit( param_file, GV, US ) log_to_all=.true., debugging=.true.) call get_param(param_file, mdl, "G_EARTH", GV%g_Earth, & "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80, scale=US%Z_to_m*US%m_s_to_L_T**2) + units="m s-2", default=9.80, scale=US%Z_to_m*US%m_s_to_L_T**2) call get_param(param_file, mdl, "RHO_0", GV%Rho0, & "The mean ocean density used with BOUSSINESQ true to "//& "calculate accelerations and the mass for conservation "//& diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 47ac43df06..a0f3855d19 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -180,7 +180,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & !$ if (.not.MOM_thread_affinity_set()) then !$ call get_param(param_file, mdl, "OCEAN_OMP_THREADS", ocean_nthreads, & !$ "The number of OpenMP threads that MOM6 will use.", & - !$ default = 1, layoutParam=.true.) + !$ default=1, layoutParam=.true.) !$ call get_param(param_file, mdl, "OCEAN_OMP_HYPER_THREAD", ocean_omp_hyper_thread, & !$ "If True, use hyper-threading.", default=.false., layoutParam=.true.) !$ call set_MOM_thread_affinity(ocean_nthreads, ocean_omp_hyper_thread) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index bde8e3e219..aaa53dee59 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1468,7 +1468,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, & "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80, scale=US%m_s_to_L_T**2*US%Z_to_m) + units="m s-2", default=9.80, scale=US%m_s_to_L_T**2*US%Z_to_m) call get_param(param_file, mdl, "C_P", CS%Cp, & "The heat capacity of sea water, approximated as a constant. "//& "The default value is from the TEOS-10 definition of conservative temperature.", & diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 552216f41d..6691095b08 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -407,7 +407,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ fail_if_missing=.true.) call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, & "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80, scale=US%m_s_to_L_T**2*US%Z_to_m) + units="m s-2", default=9.80, scale=US%m_s_to_L_T**2*US%Z_to_m) call get_param(param_file, mdl, "GLEN_EXPONENT", CS%n_glen, & "nonlinearity exponent in Glen's Law", & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index c0d7e6c50c..b4092d3d43 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -2040,7 +2040,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "is permitted for the thickness diffusivity. 1.0 is the "//& "marginally unstable value in a pure layered model, but "//& "much smaller numbers (e.g. 0.1) seem to work better for "//& - "ALE-based models.", units = "nondimensional", default=0.8) + "ALE-based models.", units="nondimensional", default=0.8) call get_param(param_file, mdl, "KH_ETA_CONST", CS%Kh_eta_bg, & "The background horizontal diffusivity of the interface heights (without "//& diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 520c70172f..63b0ced556 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -393,7 +393,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS) if (CS%tidal_sal_from_file .or. CS%use_prev_tides) then call get_param(param_file, mdl, "TIDAL_INPUT_FILE", tidal_input_files, & "A list of input files for tidal information.", & - default = "", fail_if_missing=.true.) + default="", fail_if_missing=.true.) endif call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, & diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index b69cd2daae..4f13cf5793 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -277,7 +277,7 @@ logical function CVMix_shear_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "SMOOTH_RI", CS%smooth_ri, & "If true, vertically smooth the Richardson "// & "number by applying a 1-2-1 filter once.", & - default = .false.) + default=.false.) call cvmix_init_shear(mix_scheme=CS%Mix_Scheme, & KPP_nu_zero=US%Z2_T_to_m2_s*CS%Nu_Zero, & KPP_Ri_zero=CS%Ri_zero, & diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index cd29e9a536..a34c2a2e58 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -167,7 +167,7 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, call get_param(param_file, mdl, "NDIFF_REF_PRES", CS%ref_pres, & "The reference pressure (Pa) used for the derivatives of "//& "the equation of state. If negative (default), local pressure is used.", & - units="Pa", default = -1., scale=US%kg_m3_to_R*US%m_s_to_L_T**2) + units="Pa", default=-1., scale=US%kg_m3_to_R*US%m_s_to_L_T**2) call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", CS%interior_only, & "If true, only applies neutral diffusion in the ocean interior."//& "That is, the algorithm will exclude the surface and bottom"//& @@ -245,10 +245,10 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, call get_param(param_file, mdl, "NDIFF_DEBUG", CS%debug, & "Turns on verbose output for discontinuous neutral "//& "diffusion routines.", & - default = .false.) + default=.false.) call get_param(param_file, mdl, "HARD_FAIL_HEFF", CS%hard_fail_heff, & "Bring down the model if a problem with heff is detected",& - default = .true.) + default=.true.) endif if (CS%interior_only) then diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 818fa63659..f3d04980f6 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -197,7 +197,7 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "G_EARTH", CS%G_Earth, & "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80, scale=US%m_to_L**2*US%Z_to_m*US%T_to_s**2) + units="m s-2", default=9.80, scale=US%m_to_L**2*US%Z_to_m*US%T_to_s**2) call get_param(param_file, mdl, "RHO_0", CS%Rho0, & "The mean ocean density used with BOUSSINESQ true to "//& "calculate accelerations and the mass for conservation "//& diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index f8e9b342ac..07fc539979 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -64,7 +64,7 @@ subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, jus call get_param(param_file, mdl, "DISK_X_OFFSET", xOffset, & "The x-offset of the initially elevated disk in the "//& "circle_obcs test case.", units=G%x_ax_unit_short, & - default = 0.0, do_not_log=just_read) + default=0.0, do_not_log=just_read) call get_param(param_file, mdl, "DISK_IC_AMPLITUDE", IC_amp, & "Initial amplitude of interface height displacements "//& "in the circle_obcs test case.", & diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 685ffc4bee..0b8f59a6e8 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -201,7 +201,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "G_EARTH", CS%G_Earth, & "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80, scale=US%m_to_L**2*US%Z_to_m*US%T_to_s**2) + units="m s-2", default=9.80, scale=US%m_to_L**2*US%Z_to_m*US%T_to_s**2) call get_param(param_file, mdl, "RHO_0", CS%Rho0, & "The mean ocean density used with BOUSSINESQ true to "//& "calculate accelerations and the mass for conservation "//& @@ -210,7 +210,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "DUMBBELL_SLP_AMP", CS%slp_amplitude, & "Amplitude of SLP forcing in reservoirs.", & - units="Pa", default = 10000.0, scale=US%kg_m3_to_R*US%m_s_to_L_T**2) + units="Pa", default=10000.0, scale=US%kg_m3_to_R*US%m_s_to_L_T**2) call get_param(param_file, mdl, "DUMBBELL_SLP_PERIOD", CS%slp_period, & "Periodicity of SLP forcing in reservoirs.", & units="days", default=1.0) From 3a9f6c731cdd1081adf6987cf0cee64ed3ad9d32 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 4 Jan 2023 17:27:40 -0500 Subject: [PATCH 5/9] +Make z_tol arg non-optional to cut_off_column_top Made the previously optional z_tol argument to find_depth_of_pressure_in_cell and cut_off_column_top non-optional. It was already being provided except in a call for unit testing, so adding it as a parameter there led to a simple change and the elimination of a hard-coded dimensional parameter. Also replaced the hard-coded fill values over land in MOM_temp_salt_initialize_from_Z for temperatures and salinities before regridding with the runtime variables temp_land_fill and salt_land_fill that are already being used in the same routine. This does not change any answers, probably because these values are not actually used. All answers are bitwise identical, but some subroutine arguments have been made non-optional. --- src/core/MOM_density_integrals.F90 | 7 +++---- src/initialization/MOM_state_initialization.F90 | 12 +++++++----- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 1e51612e6d..6cffea5c75 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -1550,10 +1550,10 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t !! are anomalous to [R ~> kg m-3] real, intent(in) :: G_e !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2] type(EOS_type), intent(in) :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, intent(out) :: P_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa] real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m] - real, optional, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m] + real, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m] ! Local variables real :: dp ! Pressure thickness of the layer [R L2 T-2 ~> Pa] @@ -1583,8 +1583,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t Pa_left = P_t - P_tgt ! Pa_left < 0 F_r = 1. Pa_right = P_b - P_tgt ! Pa_right > 0 - Pa_tol = GxRho * 1.0e-5*US%m_to_Z - if (present(z_tol)) Pa_tol = GxRho * z_tol + Pa_tol = GxRho * z_tol F_guess = F_l - Pa_left / (Pa_right - Pa_left) * (F_r - F_l) Pa = Pa_right - Pa_left ! To get into iterative loop diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 14459f7d0a..0504994a30 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1393,7 +1393,7 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, real, dimension(nk), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] type(remapping_CS), pointer :: remap_CS !< Remapping structure for remapping T and S, !! if associated - real, optional, intent(in) :: z_tol !< The tolerance with which to find the depth + real, intent(in) :: z_tol !< The tolerance with which to find the depth !! matching the specified pressure [Z ~> m]. integer, optional, intent(in) :: remap_answer_date !< The vintage of the order of arithmetic and !! expressions to use for remapping. Values below 20190101 @@ -2809,8 +2809,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just tmpT1dIn(i,j,k) = tmpT1dIn(i,j,k-1) tmpS1dIn(i,j,k) = tmpS1dIn(i,j,k-1) else ! This next block should only ever be reached over land - tmpT1dIn(i,j,k) = -99.9*US%degC_to_C ! Change to temp_land_fill - tmpS1dIn(i,j,k) = -99.9*US%ppt_to_S ! Change to salt_land_fill + tmpT1dIn(i,j,k) = temp_land_fill + tmpS1dIn(i,j,k) = salt_land_fill endif h1(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k @@ -3129,6 +3129,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) real :: P_tot, P_t, P_b ! Pressures [R L2 T-2 ~> Pa] real :: z_out ! Output height [Z ~> m] real :: I_z_scale ! The inverse of the height scale for prescribed gradients [Z-1 ~> m-1] + real :: z_tol ! The tolerance with which to find the depth matching a specified pressure [Z ~> m]. integer :: k type(remapping_CS), pointer :: remap_CS => NULL() @@ -3143,6 +3144,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) P_tot = 0. T_ref = 20.0*US%degC_to_C S_ref = 35.0*US%ppt_to_S + z_tol = 1.0e-5*US%m_to_Z do k = 1, nk z(k) = 0.5 * ( e(K) + e(K+1) ) T_t(k) = T_ref + (0. * I_z_scale) * e(k) @@ -3159,7 +3161,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) P_t = 0. do k = 1, nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), P_t, 0.5*P_tot, & - GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out) + GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol) write(0,*) k, US%RL2_T2_to_Pa*P_t, US%RL2_T2_to_Pa*P_b, 0.5*US%RL2_T2_to_Pa*P_tot, & US%Z_to_m*e(K), US%Z_to_m*e(K+1), US%Z_to_m*z_out P_t = P_b @@ -3171,7 +3173,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) write(0,*) '' write(0,*) GV%H_to_m*h(:) call cut_off_column_top(nk, tv, GV, US, GV%g_Earth, -e(nk+1), GV%Angstrom_Z, & - T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS) + T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol) write(0,*) GV%H_to_m*h(:) end subroutine MOM_state_init_tests From 68aefe134b853b751593267646a8f9dcfb6d6015 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 5 Jan 2023 08:57:02 -0500 Subject: [PATCH 6/9] Revise tc4/MOM_input with updated parameter names Updated tc4/MOM_input to reflect the newer parameter settings, including replacing the obsolete NEW_SPONGES parameter with INTERPOLATE_SPONGE_TIME_SPACE. Without this change, NEW_SPONGES can not be formally and properly obsoleted without breaking the TC testing. All answers and testing are identical with this change. --- .testing/tc4/MOM_input | 142 +++++++++++++++++++---------------------- 1 file changed, 65 insertions(+), 77 deletions(-) diff --git a/.testing/tc4/MOM_input b/.testing/tc4/MOM_input index e33bf40bf6..591ed4c788 100644 --- a/.testing/tc4/MOM_input +++ b/.testing/tc4/MOM_input @@ -1,25 +1,25 @@ ! This file was written by the model and records the non-default parameters used at run-time. ! === module MOM === - -! === module MOM_unit_scaling === -! Parameters for doing unit scaling of variables. USE_REGRIDDING = True ! [Boolean] default = False ! If True, use the ALE algorithm (regridding/remapping). If False, use the ! layered isopycnal algorithm. -DT = 1200.0 ! [s] +DT = 1200.0 ! [s] ! The (baroclinic) dynamics time step. The time-step that is actually used will ! be an integer fraction of the forcing time-step (DT_FORCING in ocean-only mode ! or the coupling timestep in coupled mode.) -DT_THERM = 3600.0 ! [s] default = 300.0 +DT_THERM = 3600.0 ! [s] default = 1200.0 ! The thermodynamic and tracer advection time step. Ideally DT_THERM should be ! an integer multiple of DT and less than the forcing or coupling time-step, ! unless THERMO_SPANS_COUPLING is true, in which case DT_THERM can be an integer - ! multiple of the coupling timestep. By default DT_THERM is set to DT. + ! multiple of the coupling timestep. By default DT_THERM is set to DT. C_P = 3925.0 ! [J kg-1 K-1] default = 3991.86795711963 ! The heat capacity of sea water, approximated as a constant. This is only used ! if ENABLE_THERMODYNAMICS is true. The default value is from the TEOS-10 ! definition of conservative temperature. +USE_PSURF_IN_EOS = False ! [Boolean] default = True + ! If true, always include the surface pressure contributions in equation of + ! state calculations. SAVE_INITIAL_CONDS = False ! [Boolean] default = False ! If true, write the initial conditions to a file given by IC_OUTPUT_FILE. @@ -33,9 +33,6 @@ NJGLOBAL = 10 ! ! The total number of thickness grid points in the y-direction in the physical ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. -! === module MOM_hor_index === -! Sets the horizontal array index types. - ! === module MOM_verticalGrid === ! Parameters providing information about the vertical grid. NK = 2 ! [nondim] @@ -65,8 +62,9 @@ TOPO_CONFIG = "file" ! ! wall at the southern face. ! halfpipe - a zonally uniform channel with a half-sine ! profile in the meridional direction. + ! bbuilder - build topography from list of functions. ! benchmark - use the benchmark test case topography. - ! Neverland - use the Neverland test case topography. + ! Neverworld - use the Neverworld test case topography. ! DOME - use a slope and channel configuration for the ! DOME sill-overflow test case. ! ISOMIP - use a slope and channel configuration for the @@ -83,9 +81,6 @@ TOPO_CONFIG = "file" ! !MAXIMUM_DEPTH = 100.0 ! [m] ! The (diagnosed) maximum depth of the ocean. -! === module MOM_open_boundary === -! Controls where open boundaries are located, what kind of boundary condition to impose, and what data to apply, -! if any. ROTATION = "betaplane" ! default = "2omegasinlat" ! This specifies how the Coriolis parameter is specified: ! 2omegasinlat - Use twice the planetary rotation rate @@ -94,6 +89,10 @@ ROTATION = "betaplane" ! default = "2omegasinlat" ! USER - call a user modified routine. F_0 = 1.0E-04 ! [s-1] default = 0.0 ! The reference value of the Coriolis parameter with the betaplane option. +GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = False + ! If true, use an older algorithm to calculate the sine and cosines needed + ! rotate between grid-oriented directions and true north and east. Differences + ! arise at the tripolar fold. ! === module MOM_tracer_registry === @@ -106,12 +105,10 @@ DRHO_DS = 0.0 ! [kg m-3 PSU-1] default = 0.8 ! When EQN_OF_STATE=LINEAR, this is the partial derivative of density with ! salinity. -! === module MOM_restart === - ! === module MOM_tracer_flow_control === ! === module MOM_coord_initialization === -COORD_CONFIG = "linear" ! +COORD_CONFIG = "linear" ! default = "none" ! This specifies how layers are to be defined: ! ALE or none - used to avoid defining layers in ALE mode ! file - read coordinate information from the file @@ -129,6 +126,10 @@ COORD_CONFIG = "linear" ! ! ts_profile - use temperature and salinity profiles ! (read from COORD_FILE) to set layer densities. ! USER - call a user modified routine. +REMAP_UV_USING_OLD_ALG = True ! [Boolean] default = False + ! If true, uses the old remapping-via-a-delta-z method for remapping u and v. If + ! false, uses the new method that remaps between grids described by an old and + ! new thickness. REGRIDDING_COORDINATE_MODE = "Z*" ! default = "LAYER" ! Coordinate mode for vertical regridding. Choose among the following ! possibilities: LAYER - Isopycnal or stacked shallow water layers @@ -137,6 +138,7 @@ REGRIDDING_COORDINATE_MODE = "Z*" ! default = "LAYER" ! SIGMA - terrain following coordinates ! RHO - continuous isopycnal ! HYCOM1 - HyCOM-like hybrid coordinate + ! HYBGEN - Hybrid coordinate from the Hycom hybgen code ! SLIGHT - stretched coordinates above continuous isopycnal ! ADAPTIVE - optimize for smooth neutral density surfaces !ALE_RESOLUTION = 2*50.0 ! [m] @@ -150,14 +152,14 @@ REMAPPING_SCHEME = "PPM_IH4" ! default = "PLM" ! variables. It can be one of the following schemes: PCM (1st-order ! accurate) ! PLM (2nd-order accurate) + ! PLM_HYBGEN (2nd-order accurate) ! PPM_H4 (3rd-order accurate) ! PPM_IH4 (3rd-order accurate) + ! PPM_HYBGEN (3rd-order accurate) + ! WENO_HYBGEN (3rd-order accurate) ! PQM_IH4IH3 (4th-order accurate) ! PQM_IH6IH5 (5th-order accurate) -! === module MOM_grid === -! Parameters providing information about the lateral grid. - ! === module MOM_state_initialization === INIT_LAYERS_FROM_Z_FILE = True ! [Boolean] default = False ! If true, initialize the layer thicknesses, temperatures, and salinities from a @@ -181,9 +183,9 @@ SPONGE_PTEMP_VAR = "ptemp" ! default = "PTEMP" ! The name of the potential temperature variable in SPONGE_STATE_FILE. SPONGE_SALT_VAR = "salt" ! default = "SALT" ! The name of the salinity variable in SPONGE_STATE_FILE. -NEW_SPONGES = True ! [of sponge restoring data.] default = False - ! Set True if using the newer sponging code which performs on-the-fly regridding - ! in lat-lon-time. +INTERPOLATE_SPONGE_TIME_SPACE = True ! [Boolean] default = False + ! If True, perform on-the-fly regridding in lat-lon-time of sponge restoring + ! data. ! === module MOM_sponge === SPONGE_DATA_ONGRID = True ! [Boolean] default = False @@ -192,8 +194,9 @@ SPONGE_DATA_ONGRID = True ! [Boolean] default = False ! The total number of columns where sponges are applied at h points. ! === module MOM_diag_mediator === - -! === module MOM_MEKE === +DIAG_AS_CHKSUM = True ! [Boolean] default = False + ! Instead of writing diagnostics to the diag manager, write a text file + ! containing the checksum (bitcount) of the array. ! === module MOM_lateral_mixing_coeffs === @@ -202,10 +205,10 @@ LINEAR_DRAG = True ! [Boolean] default = False ! If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag law is ! cdrag*DRAG_BG_VEL*u. HBBL = 10.0 ! [m] - ! The thickness of a bottom boundary layer with a viscosity of KVBBL if - ! BOTTOMDRAGLAW is not defined, or the thickness over which near-bottom - ! velocities are averaged for the drag law if BOTTOMDRAGLAW is defined but - ! LINEAR_DRAG is not. + ! The thickness of a bottom boundary layer with a viscosity increased by + ! KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which + ! near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is + ! defined but LINEAR_DRAG is not. CDRAG = 0.002 ! [nondim] default = 0.003 ! CDRAG is the drag coefficient relating the magnitude of the velocity field to ! the bottom stress. CDRAG is only used if BOTTOMDRAGLAW is defined. @@ -214,7 +217,7 @@ DRAG_BG_VEL = 0.05 ! [m s-1] default = 0.0 ! unresolved velocity that is combined with the resolved velocity to estimate ! the velocity magnitude. DRAG_BG_VEL is only used when BOTTOMDRAGLAW is ! defined. -BBL_USE_EOS = True ! [Boolean] default = False +BBL_USE_EOS = True ! [Boolean] default = True ! If true, use the equation of state in determining the properties of the bottom ! boundary layer. Otherwise use the layer target potential densities. BBL_THICK_MIN = 0.1 ! [m] default = 0.0 @@ -228,6 +231,13 @@ KV = 1.0E-04 ! [m2 s-1] ! === module MOM_thickness_diffuse === KHTH = 500.0 ! [m2 s-1] default = 0.0 ! The background horizontal thickness diffusivity. +USE_GM_WORK_BUG = True ! [Boolean] default = False + ! If true, compute the top-layer work tendency on the u-grid with the incorrect + ! sign, for legacy reproducibility. + +! === module MOM_porous_barriers === + +! === module MOM_dynamics_split_RK2 === BE = 0.7 ! [nondim] default = 0.6 ! If SPLIT is true, BE determines the relative weighting of a 2nd-order ! Runga-Kutta baroclinic time stepping scheme (0.5) and a backward Euler scheme @@ -258,7 +268,7 @@ BOUND_CORIOLIS = True ! [Boolean] default = False ! === module MOM_PressureForce === -! === module MOM_PressureForce_AFV === +! === module MOM_PressureForce_FV === RECONSTRUCT_FOR_PRESSURE = False ! [Boolean] default = True ! If True, use vertical reconstruction of T & S within the integrals of the FV ! pressure gradient calculation. If False, use the constant-by-layer algorithm. @@ -269,17 +279,25 @@ SMAGORINSKY_AH = True ! [Boolean] default = False ! If true, use a biharmonic Smagorinsky nonlinear eddy viscosity. SMAG_BI_CONST = 0.03 ! [nondim] default = 0.0 ! The nondimensional biharmonic Smagorinsky constant, typically 0.015 - 0.06. +USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = True + ! If true, use the land mask for the computation of thicknesses at velocity + ! locations. This eliminates the dependence on arbitrary values over land or + ! outside of the domain. ! === module MOM_vert_friction === DIRECT_STRESS = True ! [Boolean] default = False ! If true, the wind stress is distributed over the topmost HMIX_STRESS of fluid - ! (like in HYCOM), and KVML may be set to a very small value. + ! (like in HYCOM), and an added mixed layer viscosity or a physically based + ! boundary layer turbulence parameterization is not needed for stability. HMIX_FIXED = 20.0 ! [m] ! The prescribed depth over which the near-surface viscosity and diffusivity are ! elevated when the bulk mixed layer is not used. -KVML = 0.01 ! [m2 s-1] default = 1.0E-04 - ! The kinematic viscosity in the mixed layer. A typical value is ~1e-2 m2 s-1. - ! KVML is not used if BULKMIXEDLAYER is true. The default is set by KV. +KV_ML_INVZ2 = 0.01 ! [m2 s-1] default = 0.0 + ! An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, with + ! the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the distance + ! from the surface, to allow for finite wind stresses to be transmitted through + ! infinitesimally thin surface layers. This is an older option for numerical + ! convenience without a strong physical basis, and its use is now discouraged. MAXVEL = 10.0 ! [m s-1] default = 3.0E+08 ! The maximum velocity allowed before the velocity components are truncated. @@ -304,23 +322,11 @@ DTBT = 10.0 ! [s or nondim] default = -0.98 ! DTBT to 0 is the same as setting it to -0.98. The value of DTBT that will ! actually be used is an integer fraction of DT, rounding down. -! === module MOM_mixed_layer_restrat === +! === module MOM_diagnostics === ! === module MOM_diabatic_driver === ! The following parameters are used for diabatic processes. -! === module MOM_CVMix_KPP === -! This is the MOM wrapper to CVMix:KPP -! See http://cvmix.github.io/ - -! === module MOM_tidal_mixing === -! Vertical Tidal Mixing Parameterization - -! === module MOM_CVMix_conv === -! Parameterization of enhanced mixing due to convection via CVMix - -! === module MOM_entrain_diffusive === - ! === module MOM_set_diffusivity === BBL_EFFIC = 0.0 ! [nondim] default = 0.2 ! The efficiency with which the energy extracted by bottom drag drives BBL @@ -332,29 +338,18 @@ KD = 0.0 ! [m2 s-1] ! The background diapycnal diffusivity of density in the interior. Zero or the ! molecular value, ~1e-7 m2 s-1, may be used. -! === module MOM_kappa_shear === -! Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008 - -! === module MOM_CVMix_shear === -! Parameterization of shear-driven turbulence via CVMix (various options) - -! === module MOM_CVMix_ddiff === -! Parameterization of mixing due to double diffusion processes via CVMix - ! === module MOM_diabatic_aux === ! The following parameters are used for auxiliary diabatic processes. -! === module MOM_regularize_layers === - ! === module MOM_opacity === +PEN_SW_ABSORB_MINTHICK = 0.001 ! [m] default = 1.0 + ! A thickness that is used to absorb the remaining penetrating shortwave heat + ! flux when it drops below PEN_SW_FLUX_ABSORB. ! === module MOM_tracer_advect === ! === module MOM_tracer_hor_diff === -! === module MOM_neutral_diffusion === -! This module implements neutral diffusion of tracers - ! === module MOM_sum_output === MAXTRUNC = 5000 ! [truncations save_interval-1] default = 0 ! The run will be stopped, and the day set to a very large value if the velocity @@ -362,6 +357,9 @@ MAXTRUNC = 5000 ! [truncations save_interval-1] default = 0 ! to stop if there is any truncation of velocities. DATE_STAMPED_STDOUT = False ! [Boolean] default = True ! If true, use dates (not times) in messages to stdout +ENERGYSAVEDAYS = 0.125 ! [days] default = 1.0 + ! The interval in units of TIMEUNIT between saves of the energies of the run and + ! other globally summed diagnostics. ! === module MOM_surface_forcing === VARIABLE_WINDS = False ! [Boolean] default = True @@ -375,19 +373,17 @@ BUOY_CONFIG = "zero" ! WIND_CONFIG = "zero" ! ! The character string that indicates how wind forcing is specified. Valid ! options include (file), (2gyre), (1gyre), (gyres), (zero), and (USER). - -! === module MOM_restart === +GUST_CONST = 0.02 ! [Pa] default = 0.0 + ! The background gustiness in the winds. +FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = True + ! If true correct a bug in the time-averaging of the gustless wind friction + ! velocity ! === module MOM_main (MOM_driver) === -DAYMAX = 0.25 ! [days] +DAYMAX = 0.25 ! [days] ! The final time of the whole simulation, in units of TIMEUNIT seconds. This ! also sets the potential end time of the present run segment if the end time is ! not set via ocean_solo_nml in input.nml. - -ENERGYSAVEDAYS = 0.125 ! [days] default = 1.44E+04 - ! The interval in units of TIMEUNIT between saves of the - ! energies of the run and other globally summed diagnostics. - RESTART_CONTROL = 3 ! default = 1 ! An integer whose bits encode which restart files are written. Add 2 (bit 1) ! for a time-stamped file, and odd (bit 0) for a non-time-stamped file. A @@ -405,21 +401,13 @@ MAXCPU = 2.88E+04 ! [wall-clock seconds] default = -1.0 ! === module MOM_file_parser === -DIAG_AS_CHKSUM = True DEBUG = True -USE_PSURF_IN_EOS = False ! [Boolean] default = False -GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True INTERPOLATE_RES_FN = True ! [Boolean] default = True GILL_EQUATORIAL_LD = False ! [Boolean] default = False -USE_GM_WORK_BUG = True ! [Boolean] default = True FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False -REMAP_UV_USING_OLD_ALG = True ! [Boolean] default = True USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = False KAPPA_SHEAR_ITER_BUG = True ! [Boolean] default = True KAPPA_SHEAR_ALL_LAYER_TKE_BUG = True ! [Boolean] default = True USE_MLD_ITERATION = False ! [Boolean] default = False -PEN_SW_ABSORB_MINTHICK = 0.001 ! [m] default = 0.001 -GUST_CONST = 0.02 ! [Pa] default = 0.02 -FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = False From a3ef1ac8ec12d2154b5eb5f0fb3f704e401f0f31 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 5 Jan 2023 09:03:55 -0500 Subject: [PATCH 7/9] +Obsoleted the runtime parameter NEW_SPONGES Formally obsoleted the runtime parameter NEW_SPONGES. The agreed upon replacement is INTERPOLATE_SPONGE_TIME_SPACE, which has been available for almost a year. There is a warning message rather than a fatal error if NEW_SPONGES is used and both are set consistently, and a hint if they are not. Also added or amended comments describing a number of the internal variables or their units in MOM_state_initialization. All answers and output are bitwise identical. --- src/diagnostics/MOM_obsolete_params.F90 | 5 +- .../MOM_state_initialization.F90 | 170 ++++++++---------- 2 files changed, 78 insertions(+), 97 deletions(-) diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index 19f3d87429..cea0c82bcf 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -95,8 +95,9 @@ subroutine find_obsolete_params(param_file) call obsolete_real(param_file, "MIN_Z_DIAG_INTERVAL") call obsolete_char(param_file, "Z_OUTPUT_GRID_FILE") - ! This parameter is on the to-do list to be obsoleted. - ! call obsolete_logical(param_file, "NEW_SPONGES", hint="Use INTERPOLATE_SPONGE_TIME_SPACE instead.") + call read_param(param_file, "INTERPOLATE_SPONGE_TIME_SPACE", test_logic) + call obsolete_logical(param_file, "NEW_SPONGES", warning_val=test_logic, & + hint="Use INTERPOLATE_SPONGE_TIME_SPACE instead.") ! Write the file version number to the model log. call log_version(param_file, mdl, version) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0504994a30..49002d4846 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -154,9 +154,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & character(len=200) :: inputdir ! The directory where NetCDF input files are. character(len=200) :: config real :: H_rescale ! A rescaling factor for thicknesses from the representation in - ! a restart file to the internal representation in this run. + ! a restart file to the internal representation in this run [various units ~> 1] real :: vel_rescale ! A rescaling factor for velocities from the representation in - ! a restart file to the internal representation in this run. + ! a restart file to the internal representation in this run [various units ~> 1] real :: dt ! The baroclinic dynamics timestep for this run [T ~> s]. logical :: from_Z_file, useALE @@ -741,8 +741,8 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f units="m", default=0.1, scale=US%m_to_Z, do_not_log=just_read) endif call get_param(param_file, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, & - "A tolerance for detecting inconsist topography and input layer "//& - "ticknesses when ADJUST_THICKNESS is false.", & + "A tolerance for detecting inconsistent topography and input layer "//& + "thicknesses when ADJUST_THICKNESS is false.", & units="m", default=1.0, scale=US%m_to_Z, & do_not_log=(just_read.or.correct_thickness)) call get_param(param_file, mdl, "INTERFACE_IC_VAR", eta_var, & @@ -1099,8 +1099,8 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read, z_top_shelf) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & eta ! The free surface height that the model should use [Z ~> m]. real :: dilate ! A ratio by which layers are dilated [nondim]. - real :: scale_factor ! A scaling factor for the eta_sfc values that are read - ! in, which can be used to change units, for example. + real :: scale_factor ! A scaling factor for the eta_sfc values that are read in, + ! which can be used to change units, for example, often [Z m-1 ~> 1]. character(len=40) :: mdl = "depress_surface" ! This subroutine's name. character(len=200) :: inputdir, eta_srf_file ! Strings for file/path character(len=200) :: filename, eta_srf_var ! Strings for file/path @@ -1194,7 +1194,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: T_t, T_b ! Top and bottom edge values for reconstructions ! of temperature within each layer [T ~> degC] character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path - real :: scale_factor ! A file-dependent scaling factor for the input pressure. + real :: scale_factor ! A file-dependent scaling factor for the input pressure [various]. real :: min_thickness ! The minimum layer thickness, recast into Z units [Z ~> m]. real :: z_tolerance ! The tolerance with which to find the depth matching a specified pressure [Z ~> m]. integer :: i, j, k @@ -1305,7 +1305,7 @@ subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & eta ! The free surface height that the model should use [Z ~> m]. ! temporary arrays - real, dimension(SZK_(GV)) :: rho_col ! potential density in the column for use in ice + real, dimension(SZK_(GV)) :: rho_col ! potential density in the column for use in ice [R ~> kg m-3] real, dimension(SZK_(GV)) :: rho_h ! potential density multiplied by thickness [R Z ~> kg m-2] real, dimension(SZK_(GV)) :: h_tmp ! temporary storage for thicknesses [H ~> m] real, dimension(SZK_(GV)) :: p_ref ! pressure for density [R Z ~> kg m-2] @@ -1401,10 +1401,12 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, !! values use more robust forms of the same remapping expressions. ! Local variables - real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions [Z ~> m] - real, dimension(nk) :: h0, S0, T0, h1, S1, T1 + real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m] + real, dimension(nk) :: h0, h1 ! Initial and remapped layer thicknesses [H ~> m or kg m-2] + real, dimension(nk) :: S0, S1 ! Initial and remapped layer salinities [S ~> ppt] + real, dimension(nk) :: T0, T1 ! Initial and remapped layer temperatures [C ~> degC] real :: P_t, P_b ! Top and bottom pressures [R L2 T-2 ~> Pa] - real :: z_out, e_top + real :: z_out, e_top ! Interface height positions [Z ~> m] logical :: answers_2018 integer :: k @@ -1568,7 +1570,7 @@ subroutine initialize_velocity_uniform(u, v, G, GV, US, param_file, just_read) !! parameters without changing u or v. ! Local variables integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - real :: initial_u_const, initial_v_const + real :: initial_u_const, initial_v_const ! Constant initial velocities [L T-1 ~> m s-1] character(len=200) :: mdl = "initialize_velocity_uniform" ! This subroutine's name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1609,7 +1611,7 @@ subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read) ! Local variables character(len=200) :: mdl = "initialize_velocity_circular" real :: circular_max_u ! The amplitude of the zonal flow [L T-1 ~> m s-1] - real :: dpi ! A local variable storing pi = 3.14159265358979... + real :: dpi ! A local variable storing pi = 3.14159265358979... [nondim] real :: psi1, psi2 ! Values of the streamfunction at two points [L2 T-1 ~> m2 s-1] 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 ; nz = GV%ke @@ -1722,7 +1724,8 @@ subroutine initialize_temp_salt_from_profile(T, S, G, GV, US, param_file, just_r logical, intent(in) :: just_read !< If true, this call will only read !! parameters without changing T or S. ! Local variables - real, dimension(SZK_(GV)) :: T0, S0 + real, dimension(SZK_(GV)) :: T0 ! The profile of temperatures [C ~> degC] + real, dimension(SZK_(GV)) :: S0 ! The profile of salinities [S ~> ppt] integer :: i, j, k character(len=200) :: filename, ts_file, inputdir ! Strings for file/path character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files @@ -1866,12 +1869,11 @@ subroutine initialize_temp_salt_linear(T, S, G, GV, US, param_file, just_read) !! this call will only read parameters !! without changing T or S. - integer :: k - real :: S_top, T_top ! Reference salinity [S ~> ppt] and temperature [C ~> degC] within surface layer - real :: S_range, T_range ! Range of salinities [S ~> ppt] and temperatures [C ~> degC] over the vertical - !real :: delta_S, delta_T - !real :: delta + ! Local variables + real :: S_top, S_range ! Reference salinity in the surface layer and its vertical range [S ~> ppt] + real :: T_top, T_range ! Reference temperature in the surface layer and its vertical range [C ~> degC] character(len=40) :: mdl = "initialize_temp_salt_linear" ! This subroutine's name. + integer :: k if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90") call get_param(param_file, mdl, "T_TOP", T_top, & @@ -1889,25 +1891,18 @@ subroutine initialize_temp_salt_linear(T, S, G, GV, US, param_file, just_read) if (just_read) return ! All run-time parameters have been read, so return. - ! Prescribe salinity - !delta_S = S_range / ( GV%ke - 1.0 ) - !S(:,:,1) = S_top - !do k=2,GV%ke - ! S(:,:,k) = S(:,:,k-1) + delta_S - !enddo + ! Prescribe salinity and temperature, with the extrapolated top interface value prescribed. do k=1,GV%ke S(:,:,k) = S_top - S_range*((real(k)-0.5)/real(GV%ke)) T(:,:,k) = T_top - T_range*((real(k)-0.5)/real(GV%ke)) enddo - ! Prescribe temperature - !delta_T = T_range / ( GV%ke - 1.0 ) - !T(:,:,1) = T_top - !do k=2,GV%ke - ! T(:,:,k) = T(:,:,k-1) + delta_T - !enddo - !delta = 1 - !T(:,:,GV%ke/2 - (delta-1):GV%ke/2 + delta) = 1.0 + ! Prescribe salinity and temperature, but with the top layer value matching the surface value. + ! S(:,:,1) = S_top ; T(:,:,1) = T_top + ! do k=2,GV%ke + ! S(:,:,k) = S_top - S_range * (real(k-1) / real(GV%ke-1)) + ! T(:,:,k) = T_top - T_range * (real(k-1) / real(GV%ke-1)) + ! enddo call callTree_leave(trim(mdl)//'()') end subroutine initialize_temp_salt_linear @@ -1945,11 +1940,17 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t real, allocatable, dimension(:,:,:) :: h ! The target interface thicknesses [H ~> m or kg m-2]. real, dimension (SZI_(G),SZJ_(G),SZK_(GV)) :: & - tmp, tmp2 ! A temporary array for tracers. + tmp, & ! A temporary array for temperatures [C ~> degC] or other tracers. + tmp2 ! A temporary array for salinities [S ~> ppt] real, dimension (SZI_(G),SZJ_(G)) :: & - tmp_2d ! A temporary array for tracers. - real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading sponge fields - real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading sponge fields + tmp_2d ! A temporary array for mixed layer densities [R ~> kg m-3] + real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading sponge target fields + ! on the vertical grid of the input file, used for both + ! temperatures [C ~> degC] and salinities [S ~> ppt] + real, allocatable, dimension(:,:,:) :: tmp_u ! Temporary array for reading sponge target zonal + ! velocities on the vertical grid of the input file [L T-1 ~> m s-1] + real, allocatable, dimension(:,:,:) :: tmp_v ! Temporary array for reading sponge target meridional + ! velocities on the vertical grid of the input file [L T-1 ~> m s-1] real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: Idamp_u(SZIB_(G),SZJ_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1] @@ -1968,7 +1969,6 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t character(len=200) :: filename, inputdir ! Strings for file/path and path. logical :: use_ALE ! True if ALE is being used, False if in layered mode - logical :: new_sponge_param ! The value of a deprecated parameter. logical :: time_space_interp_sponge ! If true use sponge data that need to be interpolated in both ! the horizontal dimension and in time prior to vertical remapping. @@ -2024,34 +2024,9 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t endif call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.) - !### NEW_SPONGES should be obsoleted properly, rather than merely deprecated, at which - ! point only the else branch of the new_sponge_param block would be retained. - call get_param(param_file, mdl, "NEW_SPONGES", new_sponge_param, & - "Set True if using the newer sponging code which "//& - "performs on-the-fly regridding in lat-lon-time"//& - "of sponge restoring data.", default=.false., do_not_log=.true.) - if (new_sponge_param) then - call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, & - "If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", & - default=.true., do_not_log=.true.) - if (.not.time_space_interp_sponge) then - call MOM_error(FATAL, " initialize_sponges: NEW_SPONGES has been deprecated, "//& - "but is set to true inconsistently with INTERPOLATE_SPONGE_TIME_SPACE. "//& - "Remove the NEW_SPONGES input line.") - else - call MOM_error(WARNING, " initialize_sponges: NEW_SPONGES has been deprecated. "//& - "Please use INTERPOLATE_SPONGE_TIME_SPACE instead. Setting "//& - "INTERPOLATE_SPONGE_TIME_SPACE = True.") - endif - call log_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, & - "If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", & - default=.true.) - else - call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, & + call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, & "If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", & default=.false.) - endif - ! Read in sponge damping rate for tracers filename = trim(inputdir)//trim(damping_file) @@ -2143,8 +2118,8 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t if ( use_temperature) then call MOM_read_data(filename, potemp_var, tmp(:,:,:), G%Domain, scale=US%degC_to_C) call set_up_sponge_field(tmp, tv%T, G, GV, nz, Layer_CSp) - call MOM_read_data(filename, salin_var, tmp(:,:,:), G%Domain, scale=US%ppt_to_S) - call set_up_sponge_field(tmp, tv%S, G, GV, nz, Layer_CSp) + call MOM_read_data(filename, salin_var, tmp2(:,:,:), G%Domain, scale=US%ppt_to_S) + call set_up_sponge_field(tmp2, tv%S, G, GV, nz, Layer_CSp) endif ! else @@ -2247,30 +2222,34 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p oda_incupd_CSp, restart_CS, Time) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< If true, T & S are state variables. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic - !! variables. + !! variables. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: u !< The zonal velocity that is being + intent(in) :: u !< The zonal velocity that is being !! initialized [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: v !< The meridional velocity that is being - !! initialized [L T-1 ~> m s-1] + intent(in) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(oda_incupd_CS), pointer :: oda_incupd_CSp !< A pointer that is set to point to the control - !! structure for this module. + !! structure for this module. type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control structure - type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in - !! overrides any value set for - !Time. + type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in + !! overrides any value set for Time. ! Local variables - real, allocatable, dimension(:,:,:) :: hoda ! The layer thk inc. and oda layer thk [H ~> m or kg m-2]. - real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda fields - real, allocatable, dimension(:,:,:) :: tmp_u, tmp_v ! Temporary arrays for reading oda fields + real, allocatable, dimension(:,:,:) :: hoda ! The layer thickness increment and oda layer thickness [H ~> m or kg m-2] + real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda tracer increments + ! on the vertical grid of the input file, used for both + ! temperatures [C ~> degC] and salinities [S ~> ppt] + real, allocatable, dimension(:,:,:) :: tmp_u ! Temporary array for reading oda zonal velocity + ! increments on the vertical grid of the input file [L T-1 ~> m s-1] + real, allocatable, dimension(:,:,:) :: tmp_v ! Temporary array for reading oda meridional velocity + ! increments on the vertical grid of the input file [L T-1 ~> m s-1] integer :: is, ie, js, je, nz integer :: isd, ied, jsd, jed @@ -2375,7 +2354,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data), source=0.0) call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T) call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, oda_incupd_CSp) - deallocate(tmp_u,tmp_v) + deallocate(tmp_u, tmp_v) endif ! calculate increments if input are full fields @@ -2415,8 +2394,8 @@ subroutine compute_global_grid_integrals(G, US) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming - real :: area_scale + real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming ! Masked and unscaled areas for sums [m2] + real :: area_scale ! A conversion factor to prepare for reproducing sums [m2 L-2 ~> 1] integer :: i,j area_scale = US%L_to_m**2 @@ -2536,7 +2515,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just real, dimension(:,:,:), allocatable, target :: tmpS1dIn ! Input salinities on a model-sized grid [S ~> ppt] real, dimension(:,:,:), allocatable :: tmp_mask_in ! The valid data mask on a model-sized grid [nondim] real, dimension(:,:,:), allocatable :: h1 ! Thicknesses [H ~> m or kg m-2]. - real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to regridding + real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to + ! regridding [H ~> m or kg m-2] real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m]. type(regridding_CS) :: regridCS ! Regridding parameters and work arrays type(remapping_CS) :: remapCS ! Remapping parameters and work arrays @@ -2685,8 +2665,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just units="m", default=0.1, scale=US%m_to_Z, & do_not_log=(just_read.or..not.correct_thickness)) call get_param(PF, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, & - "A tolerance for detecting inconsist topography and input layer "//& - "ticknesses when ADJUST_THICKNESS is false.", & + "A tolerance for detecting inconsistent topography and input layer "//& + "thicknesses when ADJUST_THICKNESS is false.", & units="m", default=1.0, scale=US%m_to_Z, & do_not_log=(just_read.or.correct_thickness)) @@ -2727,12 +2707,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just units="1e-3", default=35.0, scale=US%ppt_to_S, do_not_log=just_read) call get_param(PF, mdl, "HORIZ_INTERP_TOL_TEMP", tol_temp, & "The tolerance in temperature changes between iterations when interpolating "//& - "ifrom an nput dataset using horiz_interp_and_extrap_tracer. This routine "//& + "from an input dataset using horiz_interp_and_extrap_tracer. This routine "//& "converges slowly, so an overly small tolerance can get expensive.", & units="degC", default=1.0e-3, scale=US%degC_to_C, do_not_log=just_read) call get_param(PF, mdl, "HORIZ_INTERP_TOL_SALIN", tol_sal, & "The tolerance in salinity changes between iterations when interpolating "//& - "ifrom an nput dataset using horiz_interp_and_extrap_tracer. This routine "//& + "from an input dataset using horiz_interp_and_extrap_tracer. This routine "//& "converges slowly, so an overly small tolerance can get expensive.", & units="1e-3", default=1.0e-3, scale=US%ppt_to_S, do_not_log=just_read) @@ -2796,11 +2776,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just allocate( tmpT1dIn(isd:ied,jsd:jed,nkd), source=0.0 ) allocate( tmpS1dIn(isd:ied,jsd:jed,nkd), source=0.0 ) do j = js, je ; do i = is, ie - if (G%mask2dT(i,j)>0.) then + if (G%mask2dT(i,j) > 0.) then zTopOfCell = 0. ; zBottomOfCell = 0. tmp_mask_in(i,j,1:kd) = mask_z(i,j,:) do k = 1, nkd - if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then + if ((tmp_mask_in(i,j,k) > 0.) .and. (k <= kd)) then zBottomOfCell = max( z_edges_in(k+1), Z_bottom(i,j)) tmpT1dIn(i,j,k) = temp_z(i,j,k) tmpS1dIn(i,j,k) = salt_z(i,j,k) @@ -2831,7 +2811,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just hTarget = getCoordinateResolution( regridCS ) do j = js, je ; do i = is, ie h(i,j,:) = 0. - if (G%mask2dT(i,j)>0.) then + if (G%mask2dT(i,j) > 0.) then ! Build the target grid combining hTarget and topography zTopOfCell = 0. ; zBottomOfCell = 0. do k = 1, nz @@ -3013,7 +2993,7 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, n real, dimension(SZK_(GV)+1) :: zi_ ! A column interface heights (negative downward) [Z ~> m]. real :: slope ! The rate of change of height with density [Z R-1 ~> m4 kg-1] real :: drhodz ! A local vertical density gradient [R Z-1 ~> kg m-4] - real, parameter :: zoff=0.999 + real, parameter :: zoff = 0.999 ! A small fractional adjustment to the density differences [nondim] logical :: unstable ! True if the column is statically unstable anywhere. integer :: nlevs_data ! The number of data values in a column. logical :: work_down ! This indicates whether this pass goes up or down the water column. @@ -3028,18 +3008,18 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, n nlevs_data = nlevs(i,j) do k=1,nlevs_data ; rho_(k) = rho(i,j,k) ; enddo - unstable=.true. + unstable = .true. work_down = .true. do while (unstable) ! Modify the input profile until it no longer has densities that decrease with depth. - unstable=.false. + unstable = .false. if (work_down) then - do k=2,nlevs_data-1 ; if (rho_(k) - rho_(k-1) < 0.0 ) then + do k=2,nlevs_data-1 ; if (rho_(k) - rho_(k-1) < 0.0) then if (k == 2) then rho_(k-1) = rho_(k) - eps_rho else drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1)) - if (drhodz < 0.0) unstable=.true. + if (drhodz < 0.0) unstable = .true. rho_(k) = rho_(k-1) + drhodz*zoff*(zin(k)-zin(k-1)) endif endif ; enddo @@ -3054,7 +3034,7 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, n endif else drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1)) - if (drhodz < 0.0) unstable=.true. + if (drhodz < 0.0) unstable = .true. rho_(k) = rho_(k+1) - drhodz*(zin(k+1)-zin(k)) endif endif ; enddo From 3f82a92bbe682a6d7dfafbf2d2de822940843389 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 5 Jan 2023 11:29:33 -0500 Subject: [PATCH 8/9] Updated comments in MOM_initialize_topography Updated the comments in MOM_initialize_topography to reflect the fact that US is no longer an optional argument. Only comments are changed, and all answers are bitwise identical. --- src/initialization/MOM_fixed_initialization.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 16702b6901..0cc3794543 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -174,14 +174,13 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) end subroutine MOM_initialize_fixed -!> MOM_initialize_topography makes the appropriate call to set up the bathymetry. At this -!! point the topography is in units of [Z ~> m] or [m], depending on the presence of US. +!> MOM_initialize_topography makes the appropriate call to set up the bathymetry in units of [Z ~> m]. subroutine MOM_initialize_topography(D, max_depth, G, PF, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth [Z ~> m] or [m] + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: PF !< Parameter file structure - real, intent(out) :: max_depth !< Maximum depth of model [Z ~> m] or [m] + real, intent(out) :: max_depth !< Maximum depth of model [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! This subroutine makes the appropriate call to set up the bottom depth. From 4be437a0bbc44111fed7ae7fcde8635bddc1e9d3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 5 Jan 2023 11:31:16 -0500 Subject: [PATCH 9/9] Added units to comments in MOM_grid_initialization Added a description of the units to the comments for each of the real variables in MOM_grid_initialization and MOM_shared_initialization. Only comments are changed, and all answers are bitwise identical. --- src/initialization/MOM_grid_initialize.F90 | 139 ++++++++++-------- .../MOM_shared_initialization.F90 | 2 +- 2 files changed, 75 insertions(+), 66 deletions(-) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 964007c663..e622b11805 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -177,7 +177,7 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpU ! East face supergrid spacing [L ~> m] real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV ! North face supergrid spacing [L ~> m] real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ ! Corner latitudes or longitudes [degN] or [degE] - real, dimension(:,:), allocatable :: tmpGlbl ! A global array of axis labels + real, dimension(:,:), allocatable :: tmpGlbl ! A global array of axis labels [degrees_N] or [km] or [m] character(len=200) :: filename, grid_file, inputdir character(len=64) :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic" type(MOM_domain_type), pointer :: SGdom => NULL() ! Supergrid domain @@ -361,8 +361,8 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) ! Local variables integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, I1off, J1off integer :: niglobal, njglobal - real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) - real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) + real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) ! Axis labels [degrees_N] or [km] or [m] + real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E] or [km] or [m] real :: dx_everywhere, dy_everywhere ! Grid spacings [L ~> m]. real :: I_dx, I_dy ! Inverse grid spacings [L-1 ~> m-1]. real :: PI @@ -498,13 +498,17 @@ subroutine set_grid_metrics_spherical(G, param_file, US) type(param_file_type), intent(in) :: param_file !< Parameter file structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: PI, PI_180! PI = 3.1415926... as 4*atan(1) + real :: PI ! PI = 3.1415926... as 4*atan(1) [nondim] + real :: PI_180 ! The conversion factor from degrees to radians [radians degree-1] integer :: i, j, isd, ied, jsd, jed integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB integer :: i_offset, j_offset - real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) - real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) - real :: dLon, dLat, latitude, dL_di + real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) ! Axis labels [degrees_N] + real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E] + real :: dLon ! The change in longitude between successive grid points [degrees_E] + real :: dLat ! The change in latitude between successive grid points [degrees_N] + real :: dL_di ! dLon rescaled from degrees to radians [radians] + real :: latitude ! The latitude of a grid point [degrees_N] character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_spherical" is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -517,7 +521,7 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! Calculate the values of the metric terms that might be used ! and save them in arrays. - PI = 4.0*atan(1.0); PI_180 = atan(1.0)/45. + PI = 4.0*atan(1.0) ; PI_180 = atan(1.0)/45. call get_param(param_file, mdl, "SOUTHLAT", G%south_lat, & "The southern latitude of the domain.", units="degrees_N", & @@ -639,19 +643,23 @@ subroutine set_grid_metrics_mercator(G, param_file, US) integer :: I_off, J_off type(GPS) :: GP character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_mercator" - real :: PI, PI_2! PI = 3.1415926... as 4*atan(1), PI_2 = (PI) /2.0 - real :: y_q, y_h, jd, x_q, x_h, id + real :: PI, PI_2 ! PI = 3.1415926... as 4*atan(1), PI_2 = (PI) /2.0 [nondim] + real :: y_q, y_h ! Latitudes of a point [radians] + real :: id ! The i-grid space positions whose longitude is being sought [gridpoints] + real :: jd ! The j-grid space positions whose latitude is being sought [gridpoints] + real :: x_q, x_h ! Longitudes of a point [radians] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: & - xh, yh ! Latitude and longitude of h points in radians. + xh, yh ! Latitude and longitude of h points in radians [radians] real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: & - xu, yu ! Latitude and longitude of u points in radians. + xu, yu ! Latitude and longitude of u points in radians [radians] real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: & - xv, yv ! Latitude and longitude of v points in radians. + xv, yv ! Latitude and longitude of v points in radians [radians] real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: & - xq, yq ! Latitude and longitude of q points in radians. + xq, yq ! Latitude and longitude of q points in radians [radians] real :: fnRef ! fnRef is the value of Int_dj_dy or ! Int_dj_dy at a latitude or longitude that is - real :: jRef, iRef ! being set to be at grid index jRef or iRef. + ! being set to be at grid index jRef or iRef [gridpoints] + real :: jRef, iRef ! The grid index at which fnRef is evaluated [gridpoints] integer :: itt1, itt2 logical, parameter :: simple_area = .true. integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB @@ -860,8 +868,8 @@ end subroutine set_grid_metrics_mercator !> This function returns the grid spacing in the logical x direction in [L ~> m]. function ds_di(x, y, GP) - real, intent(in) :: x !< The longitude in question - real, intent(in) :: y !< The latitude in question + real, intent(in) :: x !< The longitude in question [radians] + real, intent(in) :: y !< The latitude in question [radians] type(GPS), intent(in) :: GP !< A structure of grid parameters real :: ds_di ! The returned grid spacing [L ~> m] @@ -874,8 +882,8 @@ end function ds_di !> This function returns the grid spacing in the logical y direction in [L ~> m]. function ds_dj(x, y, GP) - real, intent(in) :: x !< The longitude in question - real, intent(in) :: y !< The latitude in question + real, intent(in) :: x !< The longitude in question [radians] + real, intent(in) :: y !< The latitude in question [radians] type(GPS), intent(in) :: GP !< A structure of grid parameters real :: ds_dj ! The returned grid spacing [L ~> m] @@ -887,16 +895,17 @@ function ds_dj(x, y, GP) end function ds_dj !> This function returns the contribution from the line integral along one of the four sides of a -!! cell face to the area of a cell, assuming that the sides follow a linear path in latitude and -!! longitude (i.e., on a Mercator grid). +!! cell face to the area of a cell, in [radians2], assuming that the sides follow a linear path in +!! latitude and longitude (i.e., on a Mercator grid). function dL(x1, x2, y1, y2) - real, intent(in) :: x1 !< Segment starting longitude, in degrees E. - real, intent(in) :: x2 !< Segment ending longitude, in degrees E. - real, intent(in) :: y1 !< Segment ending latitude, in degrees N. - real, intent(in) :: y2 !< Segment ending latitude, in degrees N. + real, intent(in) :: x1 !< Segment starting longitude [radians] + real, intent(in) :: x2 !< Segment ending longitude [radians] + real, intent(in) :: y1 !< Segment starting latitude [radians] + real, intent(in) :: y2 !< Segment ending latitude [radians] ! Local variables - real :: dL - real :: r, dy + real :: dL ! A contribution to the spanned area the surface of the sphere [radian2] + real :: r ! A contribution from the range of latitudes, including trigonometric factors [radians] + real :: dy ! The spanned range of latitudes [radians] dy = y2 - y1 @@ -914,22 +923,24 @@ end function dL !! Newton's method that were used to polish the root. function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax) real :: find_root !< The value of y where fn(y) = fnval that will be returned - real, external :: fn !< The external function whose root is being sought - real, external :: dy_df !< The inverse of the derivative of that function - type(GPS), intent(in) :: GP !< A structure of grid parameters - real, intent(in) :: fnval !< The value of fn being sought - real, intent(in) :: y1 !< A first guess for y - real, intent(in) :: ymin !< The minimum permitted value of y - real, intent(in) :: ymax !< The maximum permitted value of y + real, external :: fn !< The external function whose root is being sought [gridpoints] + real, external :: dy_df !< The inverse of the derivative of that function [radian gridpoint-1] + type(GPS), intent(in) :: GP !< A structure of grid parameters + real, intent(in) :: fnval !< The value of fn being sought [gridpoints] + real, intent(in) :: y1 !< A first guess for y [radians] + real, intent(in) :: ymin !< The minimum permitted value of y [radians] + real, intent(in) :: ymax !< The maximum permitted value of y [radians] integer, intent(out) :: ittmax !< The number of iterations used to polish the root ! Local variables - real :: y, y_next - real :: ybot, ytop, fnbot, fntop + real :: y, y_next ! Successive guesses at the root position [radians] + real :: ybot, ytop ! Brackets bounding the root [radians] + real :: fnbot, fntop ! Values of fn at the bounding values of y [gridpoints] + real :: dy_dfn ! The inverse of the local derivative of fn with y [radian gridpoint-1] + real :: dy ! The jump to the next guess of y [radians] + real :: fny ! The difference between fn(y) and the target value [gridpoints] integer :: itt character(len=256) :: warnmesg - real :: dy_dfn, dy, fny - ! Bracket the root. Do not use the bounding values because the value at the ! function at the bounds could be infinite, as is the case for the Mercator ! grid recursion relation. (I.e., this is a search on an open interval.) @@ -1022,40 +1033,40 @@ function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax) find_root = y end function find_root -!> This function calculates and returns the value of dx/di, where x is the -!! longitude in Radians, and i is the integral north-south grid index. +!> This function calculates and returns the value of dx/di in [radian gridpoint-1], +!! where x is the longitude in Radians, and i is the integral east-west grid index. function dx_di(x, GP) - real, intent(in) :: x !< The longitude in question + real, intent(in) :: x !< The longitude in question [radians] type(GPS), intent(in) :: GP !< A structure of grid parameters - real :: dx_di + real :: dx_di ! The derivative of zonal position with the grid index [radian gridpoint-1] dx_di = (GP%len_lon * 4.0*atan(1.0)) / (180.0 * GP%niglobal) end function dx_di !> This function calculates and returns the integral of the inverse -!! of dx/di to the point x, in radians. +!! of dx/di to the point x, in radians [gridpoints] function Int_di_dx(x, GP) - real, intent(in) :: x !< The longitude in question + real, intent(in) :: x !< The longitude in question [radians] type(GPS), intent(in) :: GP !< A structure of grid parameters - real :: Int_di_dx + real :: Int_di_dx ! A position in the global i-index space [gridpoints] Int_di_dx = x * ((180.0 * GP%niglobal) / (GP%len_lon * 4.0*atan(1.0))) end function Int_di_dx -!> This subroutine calculates and returns the value of dy/dj, where y is the -!! latitude in Radians, and j is the integral north-south grid index. +!> This subroutine calculates and returns the value of dy/dj in [radian gridpoint-1], +!! where y is the latitude in Radians, and j is the integral north-south grid index. function dy_dj(y, GP) - real, intent(in) :: y !< The latitude in question + real, intent(in) :: y !< The latitude in question [radians] type(GPS), intent(in) :: GP !< A structure of grid parameters - real :: dy_dj + real :: dy_dj ! The derivative of meridional position with the grid index [radian gridpoint-1] ! Local variables - real :: PI ! 3.1415926... calculated as 4*atan(1) + real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] real :: C0 ! The constant that converts the nominal y-spacing in - ! gridpoints to the nominal spacing in Radians. + ! gridpoints to the nominal spacing in Radians [radian gridpoint-1] real :: y_eq_enhance ! The latitude in radians within which the resolution - ! is enhanced. + ! is enhanced [radians] PI = 4.0*atan(1.0) if (GP%isotropic) then C0 = (GP%len_lon * PI) / (180.0 * GP%niglobal) @@ -1074,21 +1085,19 @@ function dy_dj(y, GP) end function dy_dj !> This subroutine calculates and returns the integral of the inverse -!! of dy/dj to the point y, in radians. +!! of dy/dj to the point y in radians [gridpoints] function Int_dj_dy(y, GP) - real, intent(in) :: y !< The latitude in question + real, intent(in) :: y !< The latitude in question [radians] type(GPS), intent(in) :: GP !< A structure of grid parameters - real :: Int_dj_dy + real :: Int_dj_dy ! The grid position of latitude y [gridpoints] ! Local variables - real :: I_C0 = 0.0 ! The inverse of the constant that converts the + real :: I_C0 ! The inverse of the constant that converts the ! nominal spacing in gridpoints to the nominal - ! spacing in Radians. - real :: PI ! 3.1415926... calculated as 4*atan(1) - real :: y_eq_enhance ! The latitude in radians from - ! from the equator within which the - ! meridional grid spacing is enhanced by - ! a factor of GP%lat_enhance_factor. - real :: r + ! spacing in Radians [gridpoint radian-1] + real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] + real :: y_eq_enhance ! The latitude in radians from from the equator within which the meridional + ! grid spacing is enhanced by a factor of GP%lat_enhance_factor [radians] + real :: r ! The y grid position in the global index space [gridpoints] PI = 4.0*atan(1.0) if (GP%isotropic) then @@ -1119,9 +1128,9 @@ end function Int_dj_dy !> Extrapolates missing metric data into all the halo regions. subroutine extrapolate_metric(var, jh, missing) - real, dimension(:,:), intent(inout) :: var !< The array in which to fill in halos + real, dimension(:,:), intent(inout) :: var !< The array in which to fill in halos [A] integer, intent(in) :: jh !< The size of the halos to be filled - real, optional, intent(in) :: missing !< The missing data fill value, 0 by default. + real, optional, intent(in) :: missing !< The missing data fill value, 0 by default [A] ! Local variables real :: badval integer :: i,j diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 73be3f5843..acffc9c927 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1309,7 +1309,7 @@ subroutine compute_global_grid_integrals(G, US) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming + real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming ! Masked and unscaled cell areas [m2] real :: area_scale ! A scaling factor for area into MKS units [m2 L-2 ~> 1] integer :: i,j