From 3f3a2e2c8d3fe6351dd8cc337b4a584166d589bd Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Tue, 20 Aug 2024 16:13:06 +1000 Subject: [PATCH] Calculate opacity in ocean_shortwave_csiro module --- .../ocean_param/sources/ocean_shortwave.F90 | 9 ++++--- .../sources/ocean_shortwave_csiro.F90 | 27 ++++++++++++++++++- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/src/mom5/ocean_param/sources/ocean_shortwave.F90 b/src/mom5/ocean_param/sources/ocean_shortwave.F90 index 6b5140486f..5cf141da43 100644 --- a/src/mom5/ocean_param/sources/ocean_shortwave.F90 +++ b/src/mom5/ocean_param/sources/ocean_shortwave.F90 @@ -294,11 +294,11 @@ subroutine sw_source (Time, Thickness, Dens, T_prog, T_diag, swflx, swflx_vis, T if(use_shortwave_gfdl) then call sw_source_gfdl (Time, Thickness, T_diag(:), swflx, swflx_vis, index_irr, Temp, sw_frac_zt, opacity) elseif(use_shortwave_csiro) then - call sw_source_csiro (Time, Thickness, T_diag(:), swflx, index_irr, Temp, sw_frac_zt) + call sw_source_csiro (Time, Thickness, T_diag(:), swflx, index_irr, Temp, sw_frac_zt, opacity) elseif(use_shortwave_jerlov) then call sw_source_jerlov (Thickness, T_diag(:), swflx, swflx_vis, index_irr, Temp, sw_frac_zt, opacity) elseif(use_shortwave_ext) then - call sw_source_ext(Time, Thickness, T_diag(:), swflx, Temp, sw_frac_zt) + call sw_source_ext(Time, Thickness, T_diag(:), swflx, Temp, sw_frac_zt, opacity) endif ! add heating rate to thickness*density weighted temperature @@ -364,7 +364,7 @@ end subroutine sw_source ! ! -subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt) +subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt, opacity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness @@ -372,6 +372,7 @@ subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt) real, dimension(isd:,jsd:) , intent(in) :: swflx type(ocean_prog_tracer_type), intent(inout) :: Temp real, dimension(isd:,jsd:,:), intent(inout) :: sw_frac_zt + real, dimension(isd:,jsd:,:), intent(inout) :: opacity integer :: k, j, i real :: sw_heat_rate(isd:ied,jsd:jed,1:nk) @@ -381,7 +382,7 @@ subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt) ! ! for our testing purposes we'll just invoke one of our methods ! Note here that index_irr is "global" - call sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt) + call sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt, opacity) ! the heating rate is stored in sw_heat_rate = Temp%wrk1 diff --git a/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 b/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 index 43f4072734..a960d2cef4 100644 --- a/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 +++ b/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 @@ -100,6 +100,7 @@ module ocean_shortwave_csiro_mod use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end use mpp_mod, only: CLOCK_ROUTINE use time_interp_external_mod, only: time_interp_external, init_external_field +use constants_mod, only: epsln use ocean_domains_mod, only: get_local_indices use ocean_types_mod, only: ocean_time_type, ocean_domain_type, ocean_grid_type @@ -306,7 +307,7 @@ end subroutine ocean_shortwave_csiro_init ! surface tracer flux "stf(i,j,temp)" ! ! -subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt) +subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt, opacity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness @@ -315,6 +316,7 @@ subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_ integer, intent(in) :: index_irr type(ocean_prog_tracer_type), intent(inout) :: Temp real, dimension(isd:,jsd:,:), intent(inout) :: sw_frac_zt + real, dimension(isd:,jsd:,:), intent(inout) :: opacity real, dimension(isd:ied,jsd:jed,0:nk) :: sw_frac_zw real, dimension(isd:ied,jsd:jed) :: zt_sw @@ -422,6 +424,29 @@ subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_ enddo enddo + ! Compute opacity, which is used by biogeochemistry. + ! We split off the k=1 level, since sw_frac_zw(k=0)=sw_frac_top, + ! which is typically set to sw_frac_top=0.0 for purposes of + ! accounting (as swflx is also in stf(index_temp). A value + ! of sw_frac_zw(k=0)=0.0 to compute opacity would result in a + ! negative opacity at k=1, which is not physical. Instead, for + ! purposes of opacity calculation, we need sw_frac_zw(k=0)=1.0. + k=1 + do j=jsc,jec + do i=isc,iec + opacity(i,j,k) = -log( sw_frac_zw(i,j,k) / (F_vis + epsln) + epsln) & + / (Thickness%dzt(i,j,k) + epsln) + enddo + enddo + do k=2,nk-1 + do j=jsc,jec + do i=isc,iec + opacity(i,j,k) = -log( sw_frac_zw(i,j,k) / (sw_frac_zw(i,j,k-1) + epsln) + epsln) & + / (Thickness%dzt(i,j,k) + epsln) + enddo + enddo + enddo + end subroutine sw_source_csiro ! NAME="sw_source_csiro"