Spatially varying bottom drag coefficient

[c2xu]

A spatially varying bottom drag coefficient can be specified by providing a map of the spatially varying scaling factor.

[herrwang0]

I think some descriptions on the background would help explain why we need this addition. For example, we have the option to read a spatially varying background velocity for ustar with a tideamp file. The maths are obviously different, but I assume that must be a scenario that this new parameterization works better and/or is physically more justifiable than using a tideamp?

[awallcraft]

Scaling ustar does not work for DRAG_AS_BODY_FORCE, and it isn't clear to me what drag terms are scaled. CDRAG is the coefficient of quadratic bottom drag, but there is also LINEAR_DRAG and DRAG_BG_VEL and BBL_USE_TIDAL_BG. It would be better to apply the scale factor to CDRAG directly, and then there is no ambiguity (if implemented correctly including for DRAG_AS_BODY_FORCE). Then CDRAG(i,j) is just CDRAG * SCALE(i,j). Alan.

…

On Sun, Nov 9, 2025 at 6:30 AM Chengzhu Xu ***@***.***> wrote: A spatially varying bottom drag coefficient can be specified by providing a map of the spatially varying scaling factor. Note that the scaling factor is applied to the friction velocity, ustar, i.e., sqrt(CDRAG), instead of CDRAG. ------------------------------ You can view, comment on, or merge this pull request online at: #983 Commit Summary - 920f087 <920f087> Spatially varying bottom drag coefficient File Changes (1 file <https://github.com/NOAA-GFDL/MOM6/pull/983/files>) - *M* src/parameterizations/vertical/MOM_set_viscosity.F90 <https://github.com/NOAA-GFDL/MOM6/pull/983/files#diff-10cada0d2f48acc55380a1286d7bd1d075ad0f998ab9829f0ce278c5dceffca4> (40) Patch Links: - https://github.com/NOAA-GFDL/MOM6/pull/983.patch - https://github.com/NOAA-GFDL/MOM6/pull/983.diff — Reply to this email directly, view it on GitHub <#983>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AHBPTUPFIS5WE7DVU4WTX4T334QVDAVCNFSM6AAAAACLSML456VHI2DSMVQWIX3LMV43ASLTON2WKOZTGYYDIOBUGUZTAMQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[c2xu]

Thanks @herrwang0 and @awallcraft. I agree that applying the scaling factor to CDRAG is the most straightforward way, and this is what's been done in the revised commit. In theory, CDRAG should be a function of the bottom roughness, and the ability to adjust CDRAG regionally could be helpful for coastal and estuarine models where the bottom roughness is spatially varying.

[awallcraft]

This looks good now, except that using a constant cdrag_2d array does not give the same answers as using CS%cdrag with the same constant value at open boundaries, as it should:

  if (.not.CS%bottomdragmap) then
    cdrag_sqrt = sqrt(CS%cdrag)
  endif

        if (CS%bottomdragmap) then
          if (m==1) then         
            cdrag_sqrt = sqrt(0.5 * (G%mask2dT(i,j) * CS%cdrag_2d(i,j) + &
                                     G%mask2dT(i+1,j) * CS%cdrag_2d(i+1,j)))
          else
            cdrag_sqrt = sqrt(0.5 * (G%mask2dT(i,j) * CS%cdrag_2d(i,j) + &
                                     G%mask2dT(i,j+1) * CS%cdrag_2d(i,j+1)))
          endif
        endif

This is the same approach used for tidal background, but it may be wrong there as well:

    ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant  
    if (CS%BBL_use_tidal_bg) then
      do i=is,ie ; if (do_i(i)) then ; if (m==1) then
        u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
                         G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) )
      else
        u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
                         G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) )
      endif ; endif ; enddo 
    else
      do i=is,ie ; if (do_i(i)) then
        u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel
      endif ; enddo    
    endif

This edge case is handled in MOM_kappa_shear.F90:

  if (CS%vertex_shear_OBC_bug) then
    !$OMP parallel do default(shared)
    do k=1,nz
      do j=JsB,JeB+1 ; do I=IsB,IeB
        h_at_u(I,j,k) = G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k)) * 0.5
      enddo ; enddo
      do J=JsB,JeB ; do i=IsB,IeB+1
        h_at_v(i,J,k) = G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k)) * 0.5
      enddo ; enddo
    enddo
  else
    ! Because G%mask2dCu(I,j) is zero if either G%mask2dT(i,j) or G%mask2dT(i+1,j) except at OBC
    ! faces, the following form give equivalent answers to those above unless OBCs are in use,
    ! although the former is clearly less complicated and costly.
    !$OMP parallel do default(shared)
      do j=JsB,JeB+1 ; do I=IsB,IeB
        h_at_u(I,j,k) = G%mask2dCu(I,j) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j) * h(i+1,j,k)) / &
                                          (G%mask2dT(i,j) + G%mask2dT(i+1,j) + 1.0e-36)
      enddo ; enddo
      do J=JsB,JeB ; do i=IsB,IeB+1
        h_at_v(i,J,k) = G%mask2dCv(i,J) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) / &
                                          (G%mask2dT(i,j) + G%mask2dT(i,j+1) + 1.0e-36)
      enddo ; enddo
    enddo
  endif

Instead of doing this, I suggest requiring cdrag_2d be defined over land.

        if (CS%bottomdragmap) then
          if (m==1) then         
            cdrag_sqrt = sqrt(0.5 * (CS%cdrag_2d(i,j) + CS%cdrag_2d(i+1,j)))
          else
            cdrag_sqrt = sqrt(0.5 * (CS%cdrag_2d(i,j) + CS%cdrag_2d(i,j+1)))
          endif
        endif

    call get_param(param_file, mdl, "CDRAG_VAR", cdrag_var, &
                 "The name of the variable in CDRAG_FILE with the spatially "//&
                 "varying bottom drag scaling factor at h points."//&
                 "Must be defined over land to allow averaging to velocity "//&
                 "points at open boundaries.", default="", &
                 do_not_log=.not.CS%bottomdragmap)

[c2xu]

This looks good now, except that using a constant cdrag_2d array does not give the same answers as using CS%cdrag with the same constant value at open boundaries, as it should:

  if (.not.CS%bottomdragmap) then
    cdrag_sqrt = sqrt(CS%cdrag)
  endif

        if (CS%bottomdragmap) then
          if (m==1) then         
            cdrag_sqrt = sqrt(0.5 * (G%mask2dT(i,j) * CS%cdrag_2d(i,j) + &
                                     G%mask2dT(i+1,j) * CS%cdrag_2d(i+1,j)))
          else
            cdrag_sqrt = sqrt(0.5 * (G%mask2dT(i,j) * CS%cdrag_2d(i,j) + &
                                     G%mask2dT(i,j+1) * CS%cdrag_2d(i,j+1)))
          endif
        endif

This is the same approach used for tidal background, but it may be wrong there as well:

    ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant  
    if (CS%BBL_use_tidal_bg) then
      do i=is,ie ; if (do_i(i)) then ; if (m==1) then
        u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
                         G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) )
      else
        u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
                         G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) )
      endif ; endif ; enddo 
    else
      do i=is,ie ; if (do_i(i)) then
        u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel
      endif ; enddo    
    endif

This edge case is handled in MOM_kappa_shear.F90:

  if (CS%vertex_shear_OBC_bug) then
    !$OMP parallel do default(shared)
    do k=1,nz
      do j=JsB,JeB+1 ; do I=IsB,IeB
        h_at_u(I,j,k) = G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k)) * 0.5
      enddo ; enddo
      do J=JsB,JeB ; do i=IsB,IeB+1
        h_at_v(i,J,k) = G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k)) * 0.5
      enddo ; enddo
    enddo
  else
    ! Because G%mask2dCu(I,j) is zero if either G%mask2dT(i,j) or G%mask2dT(i+1,j) except at OBC
    ! faces, the following form give equivalent answers to those above unless OBCs are in use,
    ! although the former is clearly less complicated and costly.
    !$OMP parallel do default(shared)
      do j=JsB,JeB+1 ; do I=IsB,IeB
        h_at_u(I,j,k) = G%mask2dCu(I,j) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j) * h(i+1,j,k)) / &
                                          (G%mask2dT(i,j) + G%mask2dT(i+1,j) + 1.0e-36)
      enddo ; enddo
      do J=JsB,JeB ; do i=IsB,IeB+1
        h_at_v(i,J,k) = G%mask2dCv(i,J) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) / &
                                          (G%mask2dT(i,j) + G%mask2dT(i,j+1) + 1.0e-36)
      enddo ; enddo
    enddo
  endif

Instead of doing this, I suggest requiring cdrag_2d be defined over land.

        if (CS%bottomdragmap) then
          if (m==1) then         
            cdrag_sqrt = sqrt(0.5 * (CS%cdrag_2d(i,j) + CS%cdrag_2d(i+1,j)))
          else
            cdrag_sqrt = sqrt(0.5 * (CS%cdrag_2d(i,j) + CS%cdrag_2d(i,j+1)))
          endif
        endif

   call get_param(param_file, mdl, "CDRAG_VAR", cdrag_var, &
                "The name of the variable in CDRAG_FILE with the spatially "//&
                "varying bottom drag scaling factor at h points."//&
                "Must be defined over land to allow averaging to velocity "//&
                "points at open boundaries.", default="", &
                do_not_log=.not.CS%bottomdragmap)

Thanks, I did not know these details. I've made changes as suggested in a separate commit, and also fixed the same issue in tidal background velocity

[herrwang0]

Thanks for the clarification and I agree this addition is a useful and valuable feature.

I agree with Alan that OBC mask needs to be handled properly. I like the simplicity and generality of the current solution, but I was wondering if it would make sense for a more defensive approach, which does not rely on the user's due diligence to include land points in input files.

I think something like the following should do the trick, and you can move it to set_visc_init to precalculate cdrag_2d_[uv].

if (G%mask2du(I,j) > 0) then
  CS%c_drag_u(I,j) = 1.0 / (G%mask2dT(i,j) + G%mask2dT(i+1,j)) &
      * (G%mask2dT(i,j) * cdrag_2d(i,j) + G%mask2dT(i+1,j) * cdrag_2d(i+1,j))
endif

Note that there would still be an issue if OBC_PROJECTION_BUG = True (tracer point outside of OBC would be unmasked), but I am not sure if it worths the trouble to do cover that scenairo.

[herrwang0]

Looks good to me!

[Hallberg-NOAA]

Overall this PR is looking good, but do we need to add error handling for the case where CDRAG_MAP is set to true but no values are provided for CDRAG_FILE or CDRAG_VAR? As written, I think that the model will fail in that case but there will not necessarily be a very useful error message that points back to the unspecified input parameters.

[Hallberg-NOAA]

Please test what happens in the case when CDRAG_MAP is set to true but no values are provided for CDRAG_FILE or CDRAG_VAR, and if the resulting error messages are not very informative, please add appropriate error handling for these cases.

[c2xu]

Please test what happens in the case when CDRAG_MAP is set to true but no values are provided for CDRAG_FILE or CDRAG_VAR, and if the resulting error messages are not very informative, please add appropriate error handling for these cases.

Thanks, an error message has been added.

[Hallberg-NOAA]

This PR is now looking good to me.

[Hallberg-NOAA]

This PR has passed pipeline testing at https://gitlab.gfdl.noaa.gov/ogrp/mom6ci/MOM6/-/pipelines/29412 with the expected warnings about a new runtime parameter.