Incomplete implementation of spatially varying drag coefficient?

[andrew-c-ross]

#983 added the useful option to pass a 2D map to scale the constant CDRAG into a spatially varying drag coefficient.

In MOM_set_viscosity.F90, if the CDRAG_MAP option is on, the map is used to set the 2D cdrag_u and cdrag_v:

MOM6/src/parameterizations/vertical/MOM_set_viscosity.F90

Lines 3306 to 3324 in ba3da16

	if (CS%bottomdragmap) then
	if (len_trim(cdrag_file)==0 .or. len_trim(cdrag_var)==0) then
	call MOM_error(FATAL,"CDRAG_FILE and CDRAG_VAR are required when using CDRAG_MAP.")
	endif
	allocate(cdrag_h(isd:ied,jsd:jed), source=0.0)
	allocate(CS%cdrag_u(IsdB:IedB,jsd:jed), source=0.0)
	allocate(CS%cdrag_v(isd:ied,JsdB:JedB), source=0.0)
	filename = trim(CS%inputdir) // trim(cdrag_file)
	call log_param(param_file, mdl, "INPUTDIR/CDRAG_FILE", filename)
	call MOM_read_data(filename, cdrag_var, cdrag_h, G%domain, scale=CS%cdrag)
	call pass_var(cdrag_h, G%domain)
	do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0) then
	CS%cdrag_u(I,j) = (G%mask2dT(i,j) * cdrag_h(i,j) + G%mask2dT(i+1,j) * cdrag_h(i+1,j)) / &
	(G%mask2dT(i,j) + G%mask2dT(i+1,j))
	endif ; enddo ; enddo
	do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0) then
	CS%cdrag_v(i,J) = (G%mask2dT(i,j) * cdrag_h(i,j) + G%mask2dT(i,j+1) * cdrag_h(i,j+1)) / &
	(G%mask2dT(i,j) + G%mask2dT(i,j+1))
	endif ; enddo ; enddo

And cdrag_u and/or cdrag_v are used for the final drag coefficient:

MOM6/src/parameterizations/vertical/MOM_set_viscosity.F90

Lines 660 to 668 in ba3da16

	if (CS%bottomdragmap) then
	if (m==1) then
	cdrag_sqrt = sqrt(CS%cdrag_u(i,j))
	else
	cdrag_sqrt = sqrt(CS%cdrag_v(i,j))
	endif
	cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
	cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
	endif

in MOM_set_diffusivity.F90, only the constant CDRAG parameter is used. The CDRAG_MAP option is never read, and the set_visc_CS structure is not available to get cdrag_u/v from. For example,

MOM6/src/parameterizations/vertical/MOM_set_diffusivity.F90

Lines 2021 to 2039 in ba3da16

	cdrag_sqrt = sqrt(CS%cdrag)
	! Find the vertical distances across layers.
	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
	!$OMP parallel default(shared) private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl)
	!$OMP do
	do J=js-1,je
	! Determine ustar and the square magnitude of the velocity in the bottom boundary layer.
	! Together these give the TKE source and vertical decay scale.
	do i=is,ie
	do_i(i) = .false. ; vstar(i,J) = 0.0 ; vhtot(i) = 0.0 ; htot(i) = 0.0
	enddo
	if (allocated(visc%Kv_bbl_v)) then
	do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then
	do_i(i) = .true.
	vstar(i,J) = visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J))
	endif ; enddo
	endif

MOM_internal_tides.F90 is also missing the 2D cdrag capability.

If CDRAG_MAP is True, do we also need to use it in MOM_set_diffusivity and MOM_internal_tides?

cc @charliestock