Bug in computation of potential temp from conservative temp (`eos_pt_from_ct`)

rabbit · 11 October 2025 13:34

We are running nemo with TEOS-10 and open cavities. We came across the bug that the output variable toce_pot(thetao) is masked with surface tmask

In src/OCE/TRA/eosbn2.F90 on line 1086:

       ztm = tmask(ji,jj,1)

Which is then applied to the derived potential temperature

         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm

I think that there is no need to perform multiplication by ztm at all, which erroneously masks open cavities. ztm can be removed from the routine.

I checked that the bug is present in v5 as well.

Could you add the fix on nemo forge, @smasson so that this is not forgotten? And thank you for fixing hfbasin, we recently pulled in the patch to nemo4.2

Dorothee · 16 October 2025 20:47

Hi!
Here is a fix. I have added a 3d version of eos_pt_from_ct
In src/OCE/TRA/eosbn2.F90:

   FUNCTION eos_pt_from_ct_3d( ctmp, psal ) RESULT( ptmp )
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE eos_pt_from_ct  ***
      !!
      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celsius]
      !!
      !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm
      !!       checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
      !!
      !! Reference  :   TEOS-10, UNESCO
      !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ctmp   ! Cons. Temp   [Celsius]
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity     [psu]
      ! Leave result array automatic rather than making explicitly allocated
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ptmp   ! potential temperature [Celsius]
      !
      INTEGER  ::   ji, jj, jk           ! dummy loop indices
      REAL(wp) ::   zt , zs , ztm        ! local scalars
      REAL(wp) ::   zn , zd              ! local scalars
      REAL(wp) ::   zdeltaS , z1_S0 , z1_T0
      !!----------------------------------------------------------------------
      !
      IF( ln_timing )   CALL timing_start('eos_pt_from_ct_3d')
      !
      zdeltaS = 5._wp
      z1_S0   = 0.875_wp/35.16504_wp
      z1_T0   = 1._wp/40._wp
      !
      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
         !
         zt  = ctmp   (ji,jj,jk) * z1_T0
         zs  = SQRT( ABS( psal(ji,jj,jk) + zdeltaS ) * z1_S0 )
         !jkt = mikt(ji,jj)
         ztm = tmask(ji,jj,jk)
         !
         zn = ((((-2.1385727895e-01_wp*zt   &
            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
            !
         zd = (2.0035003456_wp*zt   &
            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
            !
         ptmp(ji,jj,jk) = ( zt / z1_T0 + zn / zd ) * ztm
            !
      END_3D
      !
      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct_3d')
      !
   END FUNCTION eos_pt_from_ct_3d

In src/OCE/DIA/diaar5.F90:

      IF( ln_teos10 ) THEN        ! ! potential temperature (TEOS-10 case)
         IF( iom_use( 'toce_pot') .OR. iom_use( 'temptot_pot' ) .OR. iom_use( 'sst_pot' )  &
                                  .OR. iom_use( 'ssttot' ) .OR.  iom_use( 'tosmint_pot' ) ) THEN
            !
            ALLOCATE( ztpot(jpi,jpj,jpk) )
            ztpot(:,:,jpk) = 0._wp
            !
            ztpot(:,:,:) =eos_pt_from_ct_3d( ts(:,:,:,jp_tem,Kmm), ts(:,:,:,jp_sal,Kmm) )            !
            CALL iom_put( 'toce_pot', ztpot(:,:,:) )  ! potential temperature (TEOS-10 case)
            CALL iom_put( 'sst_pot' , ztpot(:,:,1) )  ! surface temperature

Dorothee · 17 October 2025 20:10

More elegant way:

   FUNCTION eos_pt_from_ct( ctmp, psal, pjk ) RESULT( ptmp )
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE eos_pt_from_ct  ***
      !!
      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celsius]
      !!
      !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm
      !!       checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
      !!
      !! Reference  :   TEOS-10, UNESCO
      !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp   [Celsius]
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity     [psu]
      INTEGER, OPTIONAL, INTENT(in   )            ::   pjk   ! level for tmask
      ! Leave result array automatic rather than making explicitly allocated
      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celsius]
      !
      INTEGER  ::   ji, jj, zjk          ! dummy loop indices
      REAL(wp) ::   zt , zs , ztm        ! local scalars
      REAL(wp) ::   zn , zd              ! local scalars
      REAL(wp) ::   zdeltaS , z1_S0 , z1_T0
      !!----------------------------------------------------------------------
      !
      IF( ln_timing )   CALL timing_start('eos_pt_from_ct')
      !
      zdeltaS = 5._wp
      z1_S0   = 0.875_wp/35.16504_wp
      z1_T0   = 1._wp/40._wp
      !
      IF (PRESENT(pjk)) THEN
         zjk = pjk
      ELSE
         zjk = 1 ! default value 
      END IF
      !
      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
         !
         zt  = ctmp   (ji,jj) * z1_T0
         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * z1_S0 )
         ztm = tmask(ji,jj,zjk)
         !
         zn = ((((-2.1385727895e-01_wp*zt   &
            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
            !
         zd = (2.0035003456_wp*zt   &
            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
            !
         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
            !
      END_2D
      !
      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct')
      !
   END FUNCTION eos_pt_from_ct

Topic		Replies	Views
Open boundary tracer: TEOS10 cons. temperature, abs. salinity and conversion function to prac. salinity Regional P4Z , TRA	3	625	24 June 2021
New regional configuration setup: open boundaries Regional BDY	24	2851	21 October 2021
Abort due to unrealistic or NaN value for SSH, velocity or salinity	13	1163	30 August 2022
[v4.0.x] p4zopt.f90 crash: out of lower bound value for `rkrgb` v4.0.x P4Z , PISCES	3	423	12 May 2022
QG Leith viscosity implementation Global DYN , LDF , orca , R025	7	374	9 January 2023

Bug in computation of potential temp from conservative temp (`eos_pt_from_ct`)

Related topics