ecearth3/sources/nemo-3.6/CONFIG/ORCA1L75_LIM3/MY_SRC/diahth.F90

MODULE diahth
   !!======================================================================
   !!                       ***  MODULE  diahth  ***
   !! Ocean diagnostics: thermocline and 20 degree depth
   !!======================================================================
   !! History :  OPA  !  1994-09  (J.-P. Boulanger)  Original code
   !!                 !  1996-11  (E. Guilyardi)  OPA8 
   !!                 !  1997-08  (G. Madec)  optimization
   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content 
   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
   !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag
   !!   NEMO     3.6  !  2015_07  (J.F. Gueremy) hd26
   !!   NEMO     3.6  !  2016_03  (M. Chevallier ) hc700 + hc2000
   !!   NEMO     3.6  !  2017_09  (C. Ethe ) taking vvl into account properly + simplyfing + cleaning 
   !!----------------------------------------------------------------------
#if   defined key_diahth   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_diahth' :                              thermocline depth diag.
   !!----------------------------------------------------------------------
   !!   dia_hth      : Compute varius diagnostics associated with the thermal vertical structure
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE phycst          ! physical constants
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! MPP library
   USE iom             ! I/O library
   USE timing          ! preformance summary

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dia_hth       ! routine called by step.F90
   PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90

   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
   ! note: following variables should move to local variables once iom_put is always used 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd26   !: depth of 26 C isotherm
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [J/m2]
   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc7   !: heat content of first 700 m                    [J/m2]
   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc2   !: heat content of first 2000 m                   [J/m2]

   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
   !! $Id$ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   FUNCTION dia_hth_alloc()
      !!---------------------------------------------------------------------
      INTEGER :: dia_hth_alloc
      !!---------------------------------------------------------------------
      !
      ALLOCATE(hth(jpi,jpj), hd20(jpi,jpj), hd26(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), &
                &  htc7(jpi,jpj),  htc2(jpi,jpj), STAT=dia_hth_alloc)
      !
      IF( lk_mpp           )   CALL mpp_sum ( dia_hth_alloc )
      IF(dia_hth_alloc /= 0)   CALL ctl_warn('dia_hth_alloc: failed to allocate arrays.')
      !
   END FUNCTION dia_hth_alloc


   SUBROUTINE dia_hth( kt )
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE dia_hth  ***
      !!
      !! ** Purpose : Computes
      !!      the mixing layer depth (turbocline): avt = 5.e-4
      !!      the depth of strongest vertical temperature gradient
      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01)
      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2       
      !!      the top of the thermochine: tn = tn(10m) - ztem2 
      !!      the pycnocline depth with density criteria equivalent to a temperature variation 
      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
      !!      the barrier layer thickness
      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
      !!      the depth of the 20 degree isotherm (linear interpolation)
      !!      the depth of the 26 degree isotherm (linear interpolation)
      !!      the depth of the 28 degree isotherm (linear interpolation)
      !!      the heat content of first 300 m (+ 700 m and 2000 m for CMIP6 :MCH)
      !!
      !! ** Method : 
      !!-------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
      !!
      INTEGER                      ::   ji, jj, jk            ! dummy loop arguments
      REAL(wp)                     ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
      REAL(wp)                     ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
      REAL(wp)                     ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
      REAL(wp)                     ::   zztmp, zzdep          ! temporary scalars inside do loop
      REAL(wp)                     ::   zu, zv, zw, zut, zvt  ! temporary workspace
      REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
      REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
      REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3      
      REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
      REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion
      REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion
      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
      REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz
      REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztemp      ! temperature
      !!----------------------------------------------------------------------
      IF( nn_timing == 1 )   CALL timing_start('dia_hth')

      IF( kt == nit000 ) THEN
         !                                      ! allocate dia_hth array
         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
         IF(lwp) WRITE(numout,*) '~~~~~~~ '
         IF(lwp) WRITE(numout,*)
         !
      ENDIF

      ! initialization
      ztinv  (:,:) = 0._wp  
      zdepinv(:,:) = 0._wp  
      zmaxdzT(:,:) = 0._wp  
      DO jj = 1, jpj
         DO ji = 1, jpi
            zztmp = bathy(ji,jj)
            hth     (ji,jj) = zztmp
            zabs2   (ji,jj) = zztmp
            ztm2    (ji,jj) = zztmp
            zrho10_3(ji,jj) = zztmp
            zpycn   (ji,jj) = zztmp
        END DO
      END DO
      IF( nla10 > 1 ) THEN 
         DO jj = 1, jpj
            DO ji = 1, jpi
               zztmp = bathy(ji,jj)
               zrho0_3(ji,jj) = zztmp
               zrho0_1(ji,jj) = zztmp
            END DO
         END DO
      ENDIF
      
      ! Preliminary computation
      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
      DO jj = 1, jpj
         DO ji = 1, jpi
            IF( tmask(ji,jj,nla10) == 1. ) THEN
               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
            ELSE
               zdelr(ji,jj) = 0._wp
            ENDIF
         END DO
      END DO

      ! ------------------------------------------------------------- !
      ! thermocline depth: strongest vertical gradient of temperature !
      ! turbocline depth (mixing layer depth): avt = zavt5            !
      ! MLD: rho = rho(1) + zrho3                                     !
      ! MLD: rho = rho(1) + zrho1                                     !
      ! ------------------------------------------------------------- !
      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
         DO jj = 1, jpj
            DO ji = 1, jpi
               !
               zzdep = fsdepw(ji,jj,jk)
               zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
               zzdep = zzdep * tmask(ji,jj,1)

               IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
               ENDIF
               
               IF( nla10 > 1 ) THEN 
                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
               ENDIF

            END DO
         END DO
      END DO
      
      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
      IF( nla10 > 1 ) THEN 
         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
      ENDIF

      ! ------------------------------------------------------------- !
      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
      ! Top of thermocline: tn = tn(10m) - ztem2                      !
      ! MLD: rho = rho10m + zrho3                                     !
      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
      ! depth of temperature inversion                                !
      ! ------------------------------------------------------------- !
      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
         DO jj = 1, jpj
            DO ji = 1, jpi
               !
               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
               !
               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
               zztmp = -zztmp                                          ! delta T(10m)
               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
               ENDIF

               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
               !
            END DO
         END DO
      END DO

      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2
      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref) 
      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref) 


      ! ------------------------------- !
      !  Depth of 20C/26C/28C isotherm  !
      ! ------------------------------- !
      IF( iom_use ("20d") ) THEN  ! depth of the 20 isotherm
         ztem2 = 20.
         CALL dia_hth_dep( ztem2, hd20 )  
         CALL iom_put( "20d", hd20 )    
      ENDIF
      !
      IF( iom_use ("26d") ) THEN  ! depth of the 26 isotherm
         ztem2 = 26.
         CALL dia_hth_dep( ztem2, hd26 )  
         CALL iom_put( "26d", hd26 )    
      ENDIF
      !
      IF( iom_use ("28d") ) THEN  ! depth of the 28 isotherm
         ztem2 = 28.
         CALL dia_hth_dep( ztem2, hd28 )  
         CALL iom_put( "28d", hd28 )    
      ENDIF
        

      ! ----------------------------- !
      !  Heat content of first 300 m  !
      ! ----------------------------- !
      IF( iom_use ("hc300") .OR. iom_use ("hc700") .OR. iom_use ("hc2000") ) THEN  
         !ztemp(:,:,:) = tsn(:,:,:,jp_tem)  ! to get heat content in J/m2
         ztemp(:,:,:) = tsn(:,:,:,jp_tem) + rt0 ! to get heat content in K*m (CMIP6)
      ENDIF

      IF( iom_use ("hc300") ) THEN  
         zzdep = 300.
         CALL  dia_hth_htc( zzdep, ztemp, htc3 )
         !htc3(:,:) = rau0_rcp * htc3(:,:)  ! to get heat content in J/m2
         CALL iom_put( "hc300", htc3 )  ! vertically integrated heat content (in K*m)
      ENDIF
      !
      ! ----------------------------- !
      !  Heat content of first 700 m  !
      ! ----------------------------- !
      IF( iom_use ("hc700") ) THEN  
         zzdep = 700.
         CALL  dia_hth_htc( zzdep, ztemp, htc7 )
         !htc7(:,:) = rau0_rcp * htc7(:,:)  ! to get heat content in J/m2
         CALL iom_put( "hc700", htc7 )  ! vertically integrated heat content (K*m) 
      ENDIF
      !
      ! ----------------------------- !
      !  Heat content of first 2000 m  !
      ! ----------------------------- !
      IF( iom_use ("hc2000") ) THEN  
         zzdep = 2000.
         CALL  dia_hth_htc( zzdep, ztemp, htc2 )
         !htc2(:,:) = rau0_rcp * htc2(:,:)  ! to get heat content in J/m2
         CALL iom_put( "hc2000", htc2 )  ! vertically integrated heat content (K*m) 
      ENDIF
      !
      IF( nn_timing == 1 )   CALL timing_stop('dia_hth')
      !
   END SUBROUTINE dia_hth
   
   SUBROUTINE dia_hth_dep( ptem, pdept )
      !
      REAL(wp), INTENT(in) :: ptem
      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept     
      !
      INTEGER  :: ji, jj, jk, iid
      REAL(wp) :: zztmp, zzdep
      INTEGER, DIMENSION(jpi,jpj) :: iktem
      
      ! --------------------------------------- !
      ! search deepest level above ptem         !
      ! --------------------------------------- !
      iktem(:,:) = 1
      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
         DO jj = 1, jpj
            DO ji = 1, jpi
               zztmp = tsn(ji,jj,jk,jp_tem)
               IF( zztmp >= ptem )   iktem(ji,jj) = jk
            END DO
         END DO
      END DO

      ! ------------------------------- !
      !  Depth of ptem isotherm         !
      ! ------------------------------- !
      DO jj = 1, jpj
         DO ji = 1, jpi
            !
            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean bottom
            !
            iid = iktem(ji,jj)
            IF( iid /= 1 ) THEN 
                zztmp =     fsdept(ji,jj,iid  )   &                     ! linear interpolation
                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
               pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
            ELSE 
               pdept(ji,jj) = 0._wp
            ENDIF
         END DO
      END DO
      !
   END SUBROUTINE dia_hth_dep


   SUBROUTINE dia_hth_htc( pdep, ptn, phtc )
      !
      REAL(wp), INTENT(in) :: pdep     ! depth over the heat content
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptn   
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc  
      !
      INTEGER  :: ji, jj, jk, ik
      REAL(wp), DIMENSION(jpi,jpj) :: zthick
      INTEGER , DIMENSION(jpi,jpj) :: ilevel


      ! surface boundary condition
      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   phtc(:,:) = 0._wp                                   
      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   phtc(:,:) = ptn(:,:,1) * sshn(:,:) * tmask(:,:,1)   
      ENDIF
      !
      ilevel(:,:) = 1
      DO jk = 2, jpkm1
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF( ( fsdept(ji,jj,jk ) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
                   ilevel(ji,jj) = jk
                   zthick(ji,jj) = zthick(ji,jj) + fse3t(ji,jj,jk)
                   phtc  (ji,jj) = phtc  (ji,jj) + fse3t(ji,jj,jk) * ptn(ji,jj,jk)
               ENDIF
            ENDDO
         ENDDO
      ENDDO
      !
      DO jj = 1, jpj
         DO ji = 1, jpi
            ik = ilevel(ji,jj)
            zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep
            phtc(ji,jj)   = phtc(ji,jj) + ptn(ji,jj,ik+1) * MIN( fse3t(ji,jj,ik+1), zthick(ji,jj) ) &
                                                                 * tmask(ji,jj,ik+1)
         END DO
      ENDDO
      !
      !
   END SUBROUTINE dia_hth_htc
   

#else
   !!----------------------------------------------------------------------
   !!   Default option :                                       Empty module
   !!----------------------------------------------------------------------
   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
CONTAINS
   SUBROUTINE dia_hth( kt )         ! Empty routine
      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
   END SUBROUTINE dia_hth
#endif

   !!======================================================================
END MODULE diahth