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/releases/release-4.0/src/OCE/DIA – NEMO

source: NEMO/releases/release-4.0/src/OCE/DIA/diahth.F90 @ 11090

Last change on this file since 11090 was 11090, checked in by agn, 5 years ago

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