Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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
!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
!! 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag
!!----------------------------------------------------------------------
!! dia_hth : Compute varius diagnostics associated with the mixed layer
!!----------------------------------------------------------------------
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
LOGICAL, SAVE :: l_hth !: thermocline-20d depths flag
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: diahth.F90 15234 2021-09-08 14:07:02Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dia_hth( kt, Kmm )
!!---------------------------------------------------------------------
!! *** 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 28 degree isotherm (linear interpolation)
!! the heat content of first 300 m
!!
!! ** Method :
!!-------------------------------------------------------------------
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level 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) :: ztemp = 0._wp
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrho10_3 ! MLD: rho = rho10m + zrho3
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztinv ! max of temperature inversion
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdepinv ! depth of temperature inversion
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmaxdzT ! max of dT/dz
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdelr ! delta rho equivalent to deltaT = 0.2
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhth ! depth of the max vertical temperature gradient [m]
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D working array
!!---------------------------------------------------------------------
IF( .NOT. l_istiled .OR. ntile == 1) THEN ! Do only for the first tile
IF( kt == nit000 ) THEN
!
l_hth = iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) .OR. &
& iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. &
& iom_use( '20d' ) .OR. iom_use( '26d' ) .OR. iom_use( '28d' ) .OR. &
& iom_use( 'hc300' ) .OR. iom_use( 'hc700' ) .OR. iom_use( 'hc2000' ) .OR. &
& iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' )
!
! ! allocate dia_hth array
IF( l_hth ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
IF(lwp) WRITE(numout,*) '~~~~~~~ '
IF(lwp) WRITE(numout,*)
ENDIF
ALLOCATE( z2d(T2D(0)), zhth(T2D(0)), zabs2(T2D(0)), ztm2(T2D(0)), zrho10_3(T2D(0)), zpycn(T2D(0)), &
& ztinv(T2D(0)), zdepinv(T2D(0)), zrho0_3(T2D(0)), zrho0_1(T2D(0)), zmaxdzT(T2D(0)), zdelr(T2D(0)))
! initialization
IF( iom_use( 'tinv' ) ) ztinv (:,:) = 0._wp
IF( iom_use( 'depti' ) ) zdepinv(:,:) = 0._wp
IF( iom_use( 'mlddzt' ) ) zmaxdzT(:,:) = 0._wp
IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) &
& .OR. iom_use( 'mldr10_3' ) .OR. iom_use( 'pycndep' ) ) THEN
zabs2 (ji,jj) = zztmp
ztm2 (ji,jj) = zztmp
zrho10_3(ji,jj) = zztmp
zpycn (ji,jj) = zztmp
END_2D
ENDIF
IF( iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
zrho0_3(ji,jj) = zztmp
zrho0_1(ji,jj) = zztmp
END_2D
ENDIF
ENDIF
IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
! ------------------------------------------------------------- !
! 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_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! loop from bottom to 2
!
zzdep = gdepw(ji,jj,jk,Kmm)
zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) &
& / 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
zhth (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_3D
CALL iom_put( 'mlddzt', zhth ) ! 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
!
ENDIF
!
IF( iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. &
& iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' ) ) THEN
!
! Preliminary computation
! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
IF( tmask(ji,jj,nla10) == 1. ) THEN
zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) &
& - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) &
& - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
zv = 5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00 * ts(ji,jj,nla10,jp_sal,Kmm) &
& - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
zut = 11.25 - 0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01 * ts(ji,jj,nla10,jp_sal,Kmm)
zvt = 38.00 - 0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
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_2D
!
! ------------------------------------------------------------- !
! 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_3DS( 0, 0, 0, 0, jpkm1, nlb10, -1 ) ! loop from bottom to nlb10
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
!
zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
!
zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ! - 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_3D
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)
!
ENDIF
! ------------------------------- !
! Depth of 20C/26C/28C isotherm !
! ------------------------------- !
IF( iom_use ('20d') ) THEN ! depth of the 20 isotherm
CALL dia_hth_dep( Kmm, 20., z2d )
CALL iom_put( '20d', z2d )
ENDIF
!
IF( iom_use ('26d') ) THEN ! depth of the 26 isotherm
CALL dia_hth_dep( Kmm, 26., z2d )
CALL iom_put( '26d', z2d )
ENDIF
!
IF( iom_use ('28d') ) THEN ! depth of the 28 isotherm
CALL dia_hth_dep( Kmm, 28., z2d )
CALL iom_put( '28d', z2d )
ENDIF
! ----------------------------- !
! Heat content of first 300 m !
! ----------------------------- !
IF( iom_use ('hc300') ) THEN
CALL dia_hth_htc( Kmm, 300., ts(:,:,:,jp_tem,Kmm), z2d )
CALL iom_put( 'hc300', rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
ENDIF
!
! ----------------------------- !
! Heat content of first 700 m !
! ----------------------------- !
IF( iom_use ('hc700') ) THEN
CALL dia_hth_htc( Kmm, 700., ts(:,:,:,jp_tem,Kmm), z2d )
CALL iom_put( 'hc700', rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
ENDIF
!
! ----------------------------- !
! Heat content of first 2000 m !
! ----------------------------- !
IF( iom_use ('hc2000') ) THEN
CALL dia_hth_htc( Kmm, 2000., ts(:,:,:,jp_tem,Kmm), z2d )
CALL iom_put( 'hc2000', rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
DEALLOCATE( z2d, zhth, zabs2, ztm2, zrho10_3, zpycn, ztinv, &
& zdepinv, zrho0_3, zrho0_1, zmaxdzT, zdelr )
ENDIF
!
IF( ln_timing ) CALL timing_stop('dia_hth')
!
END SUBROUTINE dia_hth
SUBROUTINE dia_hth_dep( Kmm, ptem, pdept )
!
INTEGER , INTENT(in) :: Kmm ! ocean time level index
REAL(wp), INTENT(in) :: ptem
REAL(wp), DIMENSION(T2D(0)), INTENT(out) :: pdept
!
INTEGER :: ji, jj, jk, iid
REAL(wp) :: zztmp, zzdep
INTEGER, DIMENSION(T2D(0)) :: iktem
! --------------------------------------- !
! search deepest level above ptem !
! --------------------------------------- !
DO_2D( 0, 0, 0, 0 )
iktem(ji,jj) = 1
END_2D
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! beware temperature is not always decreasing with depth => loop from top to bottom
zztmp = ts(ji,jj,jk,jp_tem,Kmm)
IF( zztmp >= ptem ) iktem(ji,jj) = jk
END_3D
! ------------------------------- !
! Depth of ptem isotherm !
! ------------------------------- !
!
zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the ocean bottom
!
iid = iktem(ji,jj)
IF( iid /= 1 ) THEN
zztmp = gdept(ji,jj,iid ,Kmm) & ! linear interpolation
& + ( gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm) ) &
& * ( ptem*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm) ) &
& / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (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_2D
!
END SUBROUTINE dia_hth_dep
SUBROUTINE dia_hth_htc( Kmm, pdep, pt, phtc )
!
INTEGER , INTENT(in) :: Kmm ! ocean time level index
REAL(wp), INTENT(in) :: pdep ! depth over the heat content
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt
REAL(wp), DIMENSION(T2D(0)) , INTENT(inout) :: phtc
REAL(wp), DIMENSION(T2D(0)) :: zthick
INTEGER , DIMENSION(T2D(0)) :: ilevel
IF( .NOT. ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
zthick(ji,jj) = 0._wp
phtc (ji,jj) = 0._wp
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zthick(ji,jj) = ssh(ji,jj,Kmm)
phtc (ji,jj) = pt(ji,jj,1) * ssh(ji,jj,Kmm) * tmask(ji,jj,1)
END_2D
DO_2D( 0, 0, 0, 0 )
ilevel(ji,jj) = 1
END_2D
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( ( gdepw(ji,jj,jk+1,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
ilevel(ji,jj) = jk+1
zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
phtc (ji,jj) = phtc (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk)
ENDIF
END_3D
!
ik = ilevel(ji,jj)
IF( tmask(ji,jj,ik) == 1 ) THEN
zthick(ji,jj) = MIN ( gdepw(ji,jj,ik+1,Kmm), pdep ) - zthick(ji,jj) ! remaining thickness to reach dephw pdep
phtc(ji,jj) = phtc(ji,jj) + pt(ji,jj,ik) * zthick(ji,jj)
ENDIF
END_2D
!
END SUBROUTINE dia_hth_htc
!!======================================================================
END MODULE diahth