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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 13 years ago

set proper svn properties to all files...

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