Reproducibility issue with ln_traqsr=.true. and ln_sco=.true. configurations
Context
Please provide informations on how to reproduce the bug:
-
Branches impacted: branch_4.2 and version 4 releases -
Issue is evident in the reproducibility SETTE test of AMM12 with the addition of ln_traqsr=.true. -
Seems to have been recognised and fixed already on main
Analysis
There is a reproducibility issue with configurations using terrain following vertical coordinates and the penetration of solar radiation. The issue arises from the trc_oce_ext_lev
function (trc_oce.F90
) which calculates the extinction depth as:
! branch_4.2 and earlier
! Level of light extinction
pjl = jpkm1
DO jk = jpkm1, 1, -1
IF(SUM(tmask(:,:,jk)) > 0 ) THEN
zem = MAXVAL( gdepw_0(:,:,jk+1) * tmask(:,:,jk) )
IF( zem >= zhext ) pjl = jk ! last T-level reached by Qsr
ELSE
pjl = jk ! or regional sea-bed depth
ENDIF
END DO
In S- or S-Z- coordinates the decision made is going to depend on the domain decomposition with each patch, potentially, having a different range of depths. This leads to small tracer differences (after 38 timesteps of the standard SETTE test + ln_traqsr=.true.) which grow thereafter. The issue has already been addressed in the main with a new qsr_ext_lev
function (traqsr.F90
) that finds the domain max and is therefore independent of decomposition:
! main version (proto v5.)
#if defined key_vco_1d || defined key_vco_1d3d
! z- or zps coordinate (use 1D ref vertcial coordinate)
klev = jpkm1 ! Level of light extinction zco / zps
DO jk = jpkm1, 1, -1
zdw = gdepw_1d(jk+1) ! max w-depth at jk+1 level
ze3t = e3t_1d(jk ) ! minimum e3t at jk level
zhext = - pL * LOG( zcoef * ze3t ) ! extinction depth
IF( zdw >= zhext ) klev = jk ! last T-level reached by Qsr
END DO
#else
! s- or s-z- coordinate (use 3D vertical coordinate)
klev = jpkm1 ! Level of light extinction
DO jk = jpkm1, 1, -1 !
IF( SUM( tmask(:,:,jk) ) > 0 ) THEN ! ocean point at that level
zdw = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) ) ! max w-depth at jk+1 level
ze3t = MINVAL( e3t_0(:,:,jk ) , mask=(wmask(:,:,jk+1)==1) ) ! minimum e3t at jk level
zhext = - pL * LOG( zcoef * ze3t ) ! extinction depth
IF( zdw >= zhext ) klev = jk ! last T-level reached by Qsr
ELSE ! only land point at level jk
klev = jk ! local domain sea-bed level
ENDIF
END DO
CALL mpp_max('tra_qsr', klev) ! needed for reproducibility !!st may be modified to avoid this comm.
! !!st use ssmask to remove the comm ?
#endif
Fix
The nearest equivalent fix for v4.x code would be to modify trc_oce_ext_lev
with an additional IF block:
! Level of light extinction
pjl = jpkm1
DO jk = jpkm1, 1, -1
IF(SUM(tmask(:,:,jk)) > 0 ) THEN
zem = MAXVAL( gdepw_0(:,:,jk+1) * tmask(:,:,jk) )
IF( zem >= zhext ) pjl = jk ! last T-level reached by Qsr
ELSE
pjl = jk ! or regional sea-bed depth
ENDIF
END DO
IF ( ln_sco ) THEN
! Need to take the maximum over the whole model domain when
! terrain following to avoid reproducibility issues. This is
! likely to give jpkm1 if there are any very shallow regions.
zhext = pjl
CALL mpp_max( 'trc_oce_ext_lev', zhext )
pjl = INT( zhext )
ENDIF
This fixes the reproducibility issue.