Calculate opacity in ocean_shortwave_csiro module

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,14 +364,15 @@ end subroutine sw_source
 !
 ! </DESCRIPTION>
 
-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
   type(ocean_diag_tracer_type),   intent(inout) :: T_diag(:)
   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
 

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)"
 !
 ! </DESCRIPTION>
-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
 ! </SUBROUTINE> NAME="sw_source_csiro"