New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diahth.F90 in NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA/diahth.F90 @ 12517

Last change on this file since 12517 was 12517, checked in by cguiavarch, 4 years ago

add modernised diahth w/o key
File size: 16.9 KB
Line 
1MODULE diahth
2   !!======================================================================
3   !!                       ***  MODULE  diahth  ***
4   !! Ocean diagnostics: thermocline and 20 degree depth
5   !!======================================================================
6   !! History :  OPA  !  1994-09  (J.-P. Boulanger)  Original code
7   !!                 !  1996-11  (E. Guilyardi)  OPA8
8   !!                 !  1997-08  (G. Madec)  optimization
9   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content
10   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
11   !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag
12   !!            4.0.1!  2020-01  (G. Nurser) remove need for key lk_diahth
13   !!----------------------------------------------------------------------
14   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   !
20   USE in_out_manager  ! I/O manager
21   USE lib_mpp         ! MPP library
22   USE iom             ! I/O library
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   dia_hth       ! routine called by step.F90
28   PUBLIC   dia_hth_init  ! routine called by nemogcm.F90
29
30   LOGICAL, PUBLIC ::   ll_diahth   !: Compute further diagnostics of ML and thermocline depth
31
32   !!----------------------------------------------------------------------
33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39
40   SUBROUTINE dia_hth( kt )
41     !!---------------------------------------------------------------------
42     !!                  ***  ROUTINE dia_hth  ***
43     !!
44     !! ** Purpose : Computes
45     !!      the mixing layer depth (turbocline): avt = 5.e-4
46     !!      the depth of strongest vertical temperature gradient
47     !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01)
48     !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2
49     !!      the top of the thermochine: tn = tn(10m) - ztem2
50     !!      the pycnocline depth with density criteria equivalent to a temperature variation
51     !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
52     !!      the barrier layer thickness
53     !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
54     !!      the depth of the 20 degree isotherm (linear interpolation)
55     !!      the depth of the 28 degree isotherm (linear interpolation)
56     !!      the heat content of first 300 m
57     !!
58     !! ** Method :
59     !!-------------------------------------------------------------------
60     INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
61     !!
62     INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
63     INTEGER                          ::   iid, ilevel           ! temporary integers
64     INTEGER, DIMENSION(jpi,jpj) ::   ik20, ik28  ! levels
65     REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth
66     REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
67     REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
68     REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
69     REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars
70     REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
71     REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
72     REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2
73     REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2
74     REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3
75     REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
76     REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion
77     REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion
78     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
79     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
80     REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz
81     REAL(wp), DIMENSION(jpi,jpj) ::   zthick     ! vertical integration thickness
82     REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
83   ! note: following variables should move to local variables once iom_put is always used
84     REAL(wp), DIMENSION(jpi,jpj) ::   zhth    !: depth of the max vertical temperature gradient [m]
85     REAL(wp), DIMENSION(jpi,jpj) ::   zhd20   !: depth of 20 C isotherm                         [m]
86     REAL(wp), DIMENSION(jpi,jpj) ::   zhd28   !: depth of 28 C isotherm                         [m]
87     REAL(wp), DIMENSION(jpi,jpj) ::   zhtc3   !: heat content of first 300 m                    [W]
88
89     IF (iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1")) THEN
90        ! ------------------------------------------------------------- !
91        ! thermocline depth: strongest vertical gradient of temperature !
92        ! turbocline depth (mixing layer depth): avt = zavt5            !
93        ! MLD: rho = rho(1) + zrho3                                     !
94        ! MLD: rho = rho(1) + zrho1                                     !
95        ! ------------------------------------------------------------- !
96        zmaxdzT(:,:) = 0._wp
97        IF( nla10 > 1 ) THEN
98           DO jj = 1, jpj
99              DO ji = 1, jpi
100                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)
101                 zrho0_3(ji,jj) = zztmp
102                 zrho0_1(ji,jj) = zztmp
103                 zhth(ji,jj) = zztmp
104              END DO
105           END DO
106        ELSE IF (iom_use("mlddzt")) THEN
107           DO jj = 1, jpj
108              DO ji = 1, jpi
109                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)
110                 zhth(ji,jj) = zztmp
111              END DO
112           END DO
113        ELSE
114           zhth(:,:) = 0._wp
115
116        ENDIF
117
118        DO jk = jpkm1, 2, -1   ! loop from bottom to 2
119           DO jj = 1, jpj
120              DO ji = 1, jpi
121                 !
122                 zzdep = gdepw_n(ji,jj,jk)
123                 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)
124                 zzdep = zzdep * tmask(ji,jj,1)
125
126                 IF( zztmp > zmaxdzT(ji,jj) ) THEN
127                    zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz
128                 ENDIF
129
130                 IF( nla10 > 1 ) THEN
131                    zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
132                    IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
133                    IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
134                 ENDIF
135
136              END DO
137           END DO
138        END DO
139
140        IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline
141        IF( nla10 > 1 ) THEN
142           IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03
143           IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01
144        ENDIF
145     ENDIF
146
147     IF (iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  &
148          & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti")) THEN
149        DO jj = 1, jpj
150           DO ji = 1, jpi
151              zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)
152              zabs2   (ji,jj) = zztmp
153              ztm2    (ji,jj) = zztmp
154              zrho10_3(ji,jj) = zztmp
155              zpycn   (ji,jj) = zztmp
156           END DO
157        END DO
158        ztinv  (:,:) = 0._wp
159        zdepinv(:,:) = 0._wp
160
161        IF (iom_use("pycndep")) THEN
162           ! Preliminary computation
163           ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
164           DO jj = 1, jpj
165              DO ji = 1, jpi
166                 IF( tmask(ji,jj,nla10) == 1. ) THEN
167                    zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
168                         &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
169                         &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
170                    zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
171                         &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
172                    zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
173                    zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
174                    zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
175                    zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
176                 ELSE
177                    zdelr(ji,jj) = 0._wp
178                 ENDIF
179              END DO
180           END DO
181        ELSE
182           zdelr(:,:) = 0._wp
183        ENDIF
184
185        ! ------------------------------------------------------------- !
186        ! MLD: abs( tn - tn(10m) ) = ztem2                              !
187        ! Top of thermocline: tn = tn(10m) - ztem2                      !
188        ! MLD: rho = rho10m + zrho3                                     !
189        ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
190        ! temperature inversion: max( 0, max of tn - tn(10m) )          !
191        ! depth of temperature inversion                                !
192        ! ------------------------------------------------------------- !
193        DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
194           DO jj = 1, jpj
195              DO ji = 1, jpi
196                 !
197                 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
198                 !
199                 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
200                 IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
201                 IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
202                 zztmp = -zztmp                                          ! delta T(10m)
203                 IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
204                    ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
205                 ENDIF
206
207                 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
208                 IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
209                 IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
210                 !
211              END DO
212           END DO
213        END DO
214
215        IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1)        )   ! MLD abs(delta t) - 0.2
216        IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1)         )   ! T(10) - 0.2
217        IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03
218        IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1)        )   ! MLD delta rho equi. delta T(10m) = 0.2
219        IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref)
220        IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref)
221     ENDIF
222
223     IF(iom_use("20d") .OR. iom_use("28d")) THEN
224        ! ----------------------------------- !
225        ! search deepest level above 20C/28C  !
226        ! ----------------------------------- !
227        ik20(:,:) = 1
228        ik28(:,:) = 1
229        DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
230           DO jj = 1, jpj
231              DO ji = 1, jpi
232                 zztmp = tsn(ji,jj,jk,jp_tem)
233                 IF( zztmp >= 20. )   ik20(ji,jj) = jk
234                 IF( zztmp >= 28. )   ik28(ji,jj) = jk
235              END DO
236           END DO
237        END DO
238
239        ! --------------------------- !
240        !  Depth of 20C/28C isotherm  !
241        ! --------------------------- !
242        DO jj = 1, jpj
243           DO ji = 1, jpi
244              !
245              zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
246              !
247              iid = ik20(ji,jj)
248              IF( iid /= 1 ) THEN
249                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
250                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
251                      &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
252                      &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
253                 zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
254              ELSE
255                 zhd20(ji,jj) = 0._wp
256              ENDIF
257              !
258              iid = ik28(ji,jj)
259              IF( iid /= 1 ) THEN
260                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
261                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
262                      &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
263                      &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
264                 zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
265              ELSE
266                 zhd28(ji,jj) = 0._wp
267              ENDIF
268
269           END DO
270        END DO
271        CALL iom_put( "20d", zhd20 )   ! depth of the 20 isotherm
272        CALL iom_put( "28d", zhd28 )   ! depth of the 28 isotherm
273     ENDIF
274
275     ! ----------------------------- !
276     !  Heat content of first 300 m  !
277     ! ----------------------------- !
278     IF (iom_use("hc300")) THEN
279        ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...)
280        ilevel   = 0
281        zthick_0 = 0._wp
282        DO jk = 1, jpkm1
283           zthick_0 = zthick_0 + e3t_1d(jk)
284           IF( zthick_0 < 300. )   ilevel = jk
285        END DO
286        ! surface boundary condition
287        IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)
288        ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp
289        ENDIF
290        ! integration down to ilevel
291        DO jk = 1, ilevel
292           zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk)
293           zhtc3  (:,:) = zhtc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
294        END DO
295        ! deepest layer
296        zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
297        DO jj = 1, jpj
298           DO ji = 1, jpi
299              zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  &
300                   &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
301           END DO
302        END DO
303        ! from temperature to heat contain
304        zcoef = rau0 * rcp
305        zhtc3(:,:) = zcoef * zhtc3(:,:)
306        CALL iom_put( "hc300", zhtc3*tmask(:,:,1) )      ! first 300m heat content
307     ENDIF
308     !
309   END SUBROUTINE dia_hth
310
311
312   SUBROUTINE dia_hth_init
313      !!---------------------------------------------------------------------------
314      !!                  ***  ROUTINE dia_hth_init  ***
315      !!
316      !! ** Purpose: Initialization for ML and thermocline depths
317      !!
318      !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth
319      !!---------------------------------------------------------------------------
320       !
321      IF(lwp) THEN
322         WRITE(numout,*)
323         WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics'
324         WRITE(numout,*) '~~~~~~~~~~~~ '
325      ENDIF
326      ll_diahth = iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1") .OR.  &
327           & iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  &
328           & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti").OR. &
329           & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300")
330      IF(lwp) THEN
331         WRITE(numout,*) '   output upper ocean diagnostics (T) or not (F)       ll_diahth = ', ll_diahth
332      ENDIF
333      !
334   END SUBROUTINE dia_hth_init
335END MODULE diahth
Note: See TracBrowser for help on using the repository browser.