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 trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90 @ 2561

Last change on this file since 2561 was 2561, checked in by cetlod, 13 years ago

Suppression of special characters in the name of output variables & check consistency with iodef.xml
  • Property svn:keywords set to Id
File size: 15.4 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   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
11   !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag
12   !!----------------------------------------------------------------------
13#if   defined key_diahth   ||   defined key_esopa
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   USE in_out_manager  ! I/O manager
23   USE iom             ! I/O library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   dia_hth    ! routine called by step.F90
29
30   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag
31   ! note: following variables should move to local variables once iom_put is always used
32   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hth                  !: depth of the max vertical temperature gradient [m]
33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hd20                 !: depth of 20 C isotherm                         [m]
34   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hd28                 !: depth of 28 C isotherm                         [m]
35   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   htc3                 !: heat content of first 300 m                    [W]
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE dia_hth( kt )
47      !!---------------------------------------------------------------------
48      !!                  ***  ROUTINE dia_hth  ***
49      !!
50      !! ** Purpose : Computes
51      !!      the mixing layer depth (turbocline): avt = 5.e-4
52      !!      the depth of strongest vertical temperature gradient
53      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01)
54      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2       
55      !!      the top of the thermochine: tn = tn(10m) - ztem2
56      !!      the pycnocline depth with density criteria equivalent to a temperature variation
57      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
58      !!      the barrier layer thickness
59      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
60      !!      the depth of the 20 degree isotherm (linear interpolation)
61      !!      the depth of the 28 degree isotherm (linear interpolation)
62      !!      the heat content of first 300 m
63      !!
64      !! ** Method :
65      !!-------------------------------------------------------------------
66      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
67      !!
68      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
69      INTEGER                          ::   iid, ilevel           ! temporary integers
70      INTEGER, DIMENSION(jpi,jpj)      ::   ik20, ik28            ! levels
71      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth
72      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
73      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
74      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
75      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars
76      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
77      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
78      REAL(wp), DIMENSION(jpi,jpj)     ::   zabs2                 ! MLD: abs( tn - tn(10m) ) = ztem2
79      REAL(wp), DIMENSION(jpi,jpj)     ::   ztm2                  ! Top of thermocline: tn = tn(10m) - ztem2     
80      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho10_3              ! MLD: rho = rho10m + zrho3     
81      REAL(wp), DIMENSION(jpi,jpj)     ::   zpycn                 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
82      REAL(wp), DIMENSION(jpi,jpj)     ::   ztinv                 ! max of temperature inversion
83      REAL(wp), DIMENSION(jpi,jpj)     ::   zdepinv               ! depth of temperature inversion
84      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho0_3               ! MLD rho = rho(surf) = 0.03
85      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho0_1               ! MLD rho = rho(surf) = 0.01
86      REAL(wp), DIMENSION(jpi,jpj)     ::   zmaxdzT               ! max of dT/dz
87      REAL(wp), DIMENSION(jpi,jpj)     ::   zthick                ! vertical integration thickness
88      REAL(wp), DIMENSION(jpi,jpj)     ::   zdelr                 ! delta rho equivalent to deltaT = 0.2
89      !!----------------------------------------------------------------------
90
91      IF( kt == nit000 ) THEN
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
94         IF(lwp) WRITE(numout,*) '~~~~~~~ '
95         IF(lwp) WRITE(numout,*)
96      ENDIF
97
98      ! initialization
99      ztinv  (:,:) = 0._wp 
100      zdepinv(:,:) = 0._wp 
101      zmaxdzT(:,:) = 0._wp 
102      DO jj = 1, jpj
103         DO ji = 1, jpi
104            zztmp = bathy(ji,jj)
105            hth     (ji,jj) = zztmp
106            zabs2   (ji,jj) = zztmp
107            ztm2    (ji,jj) = zztmp
108            zrho10_3(ji,jj) = zztmp
109            zpycn   (ji,jj) = zztmp
110        END DO
111      END DO
112      IF( nla10 > 1 ) THEN
113         DO jj = 1, jpj
114            DO ji = 1, jpi
115               zztmp = bathy(ji,jj)
116               zrho0_3(ji,jj) = zztmp
117               zrho0_1(ji,jj) = zztmp
118            END DO
119         END DO
120      ENDIF
121     
122      ! Preliminary computation
123      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
124      DO jj = 1, jpj
125         DO ji = 1, jpi
126            IF( tmask(ji,jj,nla10) == 1. ) THEN
127               zu  =  1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10)   &
128                  &                                                           - 0.0100*tn(ji,jj,nla10)*sn(ji,jj,nla10)
129               zv  =  5891.00 + 38.000*tn(ji,jj,nla10) + 3.00*sn(ji,jj,nla10) - 0.3750*tn(ji,jj,nla10)*tn(ji,jj,nla10)
130               zut =    11.25 -  0.149*tn(ji,jj,nla10) - 0.01*sn(ji,jj,nla10)
131               zvt =    38.00 -  0.750*tn(ji,jj,nla10)
132               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
133               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
134            ELSE
135               zdelr(ji,jj) = 0._wp
136            ENDIF
137         END DO
138      END DO
139
140      ! ------------------------------------------------------------- !
141      ! thermocline depth: strongest vertical gradient of temperature !
142      ! turbocline depth (mixing layer depth): avt = zavt5            !
143      ! MLD: rho = rho(1) + zrho3                                     !
144      ! MLD: rho = rho(1) + zrho1                                     !
145      ! ------------------------------------------------------------- !
146      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               !
150               zzdep = fsdepw(ji,jj,jk)
151               zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
152               zzdep = zzdep * tmask(ji,jj,1)
153
154               IF( zztmp > zmaxdzT(ji,jj) ) THEN                       
155                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
156               ENDIF
157               
158               IF( nla10 > 1 ) THEN
159                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
160                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
161                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
162               ENDIF
163
164            END DO
165         END DO
166      END DO
167     
168      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
169      IF( nla10 > 1 ) THEN
170         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
171         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
172      ENDIF
173
174      ! ------------------------------------------------------------- !
175      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
176      ! Top of thermocline: tn = tn(10m) - ztem2                      !
177      ! MLD: rho = rho10m + zrho3                                     !
178      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
179      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
180      ! depth of temperature inversion                                !
181      ! ------------------------------------------------------------- !
182      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
183         DO jj = 1, jpj
184            DO ji = 1, jpi
185               !
186               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
187               !
188               zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk)                  ! - delta T(10m)
189               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
190               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
191               zztmp = -zztmp                                          ! delta T(10m)
192               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
193                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
194               ENDIF
195
196               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
197               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
198               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
199               !
200            END DO
201         END DO
202      END DO
203
204      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2
205      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
206      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
207      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
208      CALL iom_put( "BLT"     , ztm2 - zpycn )   ! Barrier Layer Thickness
209      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)
210      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)
211
212
213      ! ----------------------------------- !
214      ! search deepest level above 20C/28C  !
215      ! ----------------------------------- !
216      ik20(:,:) = 1
217      ik28(:,:) = 1
218      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
219         DO jj = 1, jpj
220            DO ji = 1, jpi
221               zztmp = tn(ji,jj,jk)
222               IF( zztmp >= 20. )   ik20(ji,jj) = jk
223               IF( zztmp >= 28. )   ik28(ji,jj) = jk
224            END DO
225         END DO
226      END DO
227
228      ! --------------------------- !
229      !  Depth of 20C/28C isotherm  !
230      ! --------------------------- !
231      DO jj = 1, jpj
232         DO ji = 1, jpi
233            !
234            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
235            !
236            iid = ik20(ji,jj)
237            IF( iid /= 1 ) THEN
238               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
239                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
240                  &  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
241                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
242               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
243            ELSE
244               hd20(ji,jj) = 0._wp
245            ENDIF
246            !
247            iid = ik28(ji,jj)
248            IF( iid /= 1 ) THEN
249               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
250                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
251                  &  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
252                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
253               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
254            ELSE
255               hd28(ji,jj) = 0._wp
256            ENDIF
257
258         END DO
259      END DO
260      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
261      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
262
263      ! ----------------------------- !
264      !  Heat content of first 300 m  !
265      ! ----------------------------- !
266
267      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...)
268      ilevel   = 0
269      zthick_0 = 0._wp
270      DO jk = 1, jpkm1                     
271         zthick_0 = zthick_0 + e3t_0(jk)
272         IF( zthick_0 < 300. )   ilevel = jk
273      END DO
274      ! surface boundary condition
275      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
276      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)   
277      ENDIF
278      ! integration down to ilevel
279      DO jk = 1, ilevel
280         zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
281         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tn(:,:,jk) * tmask(:,:,jk)
282      END DO
283      ! deepest layer
284      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
285      DO jj = 1, jpj
286         DO ji = 1, jpi
287            htc3(ji,jj) = htc3(ji,jj) + tn(ji,jj,ilevel+1) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
288         END DO
289      END DO
290      ! from temperature to heat contain
291      zcoef = rau0 * rcp
292      htc3(:,:) = zcoef * htc3(:,:)
293      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
294      !
295   END SUBROUTINE dia_hth
296
297#else
298   !!----------------------------------------------------------------------
299   !!   Default option :                                       Empty module
300   !!----------------------------------------------------------------------
301   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
302CONTAINS
303   SUBROUTINE dia_hth( kt )         ! Empty routine
304      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
305   END SUBROUTINE dia_hth
306#endif
307
308   !!======================================================================
309END MODULE diahth
Note: See TracBrowser for help on using the repository browser.