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

source: trunk/NEMO/OPA_SRC/DIA/diahth.F90 @ 1484

Last change on this file since 1484 was 1484, checked in by smasson, 15 years ago

bugfix in 300m heat contain, see ticket:342
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.0 KB
Line 
1MODULE diahth
2   !!======================================================================
3   !!                       ***  MODULE  diahth  ***
4   !! Ocean diagnostics: thermocline and 20 degree depth
5   !!======================================================================
6#if   defined key_diahth   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_diahth' :                              thermocline depth diag.
9   !!----------------------------------------------------------------------
10   !!   dia_hth      : Compute diagnostics associated with the thermocline
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE phycst          ! physical constants
16   USE in_out_manager  ! I/O manager
17   USE iom
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC dia_hth    ! routine called by step.F90
24
25   !! * Shared module variables
26   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag
27   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
28      hth  ,      &  !: depth of the max vertical temperature gradient (m)
29      hd20 ,      &  !: depth of 20 C isotherm (m)
30      hd28 ,      &  !: depth of 28 C isotherm (m)
31      htc3           !: heat content of first 300 m
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35   !!----------------------------------------------------------------------
36   !!   OPA 9.0 , LOCEAN-IPSL (2005)
37   !! $Id$
38   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE dia_hth( kt )
44      !!---------------------------------------------------------------------
45      !!                  ***  ROUTINE dia_hth  ***
46      !!
47      !! ** Purpose :
48      !!      Computes the depth of strongest vertical temperature gradient
49      !!      Computes the depth of the 20 degree isotherm
50      !!      Computes the depth of the 28 degree isotherm
51      !!      Computes the heat content of first 300 m
52      !!
53      !! ** Method :
54      !!
55      !! History :
56      !!        !  94-09  (J.-P. Boulanger)  Original code
57      !!        !  96-11  (E. Guilyardi)  OPA8
58      !!        !  97-08  (G. Madec)  optimization
59      !!        !  99-07  (E. Guilyardi)  hd28 + heat content
60      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
61      !!-------------------------------------------------------------------
62      !! * Arguments
63      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
64
65      !! * Local declarations
66      INTEGER :: ji, jj, jk         ! dummy loop arguments
67      INTEGER :: iid, iif, ilevel   ! temporary integers
68      INTEGER, DIMENSION(jpi) ::   idepth
69      INTEGER, DIMENSION(jpi,jpj) ::   ikc
70
71      REAL(wp) :: zd, zmoy, zthick_0, zcoef              ! temporary scalars
72      REAL(wp), DIMENSION(jpi) ::   zmax
73      REAL(wp), DIMENSION(jpi,jpj) ::   zthick
74      REAL(wp), DIMENSION(jpi,jpk) ::   zdzt
75      !!----------------------------------------------------------------------
76
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
80         IF(lwp) WRITE(numout,*) '~~~~~~~ '
81         IF(lwp) WRITE(numout,*)
82      ENDIF
83
84
85      ! -------------------------- !
86      !  Depth of the thermocline  !
87      ! -------------------------- !
88      ! The depth of the thermocline is defined as the depth of the
89      ! strongest vertical temperature gradient
90     
91      DO jj = 1, jpj
92         
93         ! vertical gradient of temperature
94         DO jk = 2, jpkm1
95            zdzt(:,jk) = ( tn(:,jj,jk-1) - tn(:,jj,jk) ) / fse3w(:,jj,jk) * tmask(:,jj,jk)
96         END DO
97         
98         ! search the level of maximum vertical temperature gradient
99         zmax  (:) = 0.e0
100         idepth(:) = 1
101         DO jk = jpkm1, 2, -1
102            DO ji = 1, jpi
103               IF( zdzt(ji,jk) > zmax(ji) ) THEN
104                  zmax  (ji) = zdzt(ji,jk)
105                  idepth(ji) = jk
106               ENDIF
107            END DO
108         END DO
109
110         ! depth of the thermocline
111         DO ji = 1, jpi
112            hth(ji,jj) = fsdepw(ji,jj,idepth(ji))
113         END DO
114         
115      END DO
116      CALL iom_put( "thermod", hth )   ! depth of the thermocline
117
118
119      ! ----------------------- !
120      !  Depth of 20C isotherm  !
121      ! ----------------------- !
122
123      ! initialization to the number of ocean w-point mbathy
124      ! (cf dommsk, minimum value: 1)
125      ikc(:,:) = 1
126
127      ! search the depth of 20 degrees isotherm
128      ! ( starting from the top, last level above 20C, if not exist, = 1)
129      DO jk = 1, jpkm1
130         DO jj = 1, jpj
131            DO ji = 1, jpi
132               IF( tn(ji,jj,jk) >= 20. ) ikc(ji,jj) = jk
133            END DO
134         END DO
135      END DO
136     
137      ! Depth of 20C isotherm
138      DO jj = 1, jpj
139         DO ji = 1, jpi
140            iid = ikc(ji,jj)
141            iif = mbathy(ji,jj)
142            IF( iid /= 1 ) THEN 
143               ! linear interpolation
144               zd =  fsdept(ji,jj,iid)   &
145                  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) )   &
146                  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid) )   &
147                  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid)    &
148                  + (1.-tmask(ji,jj,1))                       )
149               ! bound by the ocean depth, minimum value, first T-point depth
150               hd20(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif))
151            ELSE
152               hd20(ji,jj)=0.
153            ENDIF
154         END DO
155      END DO
156      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
157
158      ! ----------------------- !
159      !  Depth of 28C isotherm  !
160      ! ----------------------- !
161     
162      ! initialization to the number of ocean w-point mbathy
163      ! (cf dommsk, minimum value: 1)
164      ikc(:,:) = 1
165     
166      ! search the depth of 28 degrees isotherm
167      ! ( starting from the top, last level above 28C, if not exist, = 1)
168      DO jk = 1, jpkm1
169         DO jj = 1, jpj
170            DO ji = 1, jpi
171               IF( tn(ji,jj,jk) >= 28. ) ikc(ji,jj) = jk
172            END DO
173         END DO
174      END DO
175     
176      ! Depth of 28C isotherm
177      DO jj = 1, jpj
178         DO ji = 1, jpi
179            iid = ikc(ji,jj)
180            iif = mbathy(ji,jj)
181            IF( iid /= 1 ) THEN 
182               ! linear interpolation
183               zd =  fsdept(ji,jj,iid)   &
184                  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) )   &
185                  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid) )   &
186                  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid)    &
187                  + ( 1. - tmask(ji,jj,1) )  )
188               ! bound by the ocean depth, minimum value, first T-point depth
189               hd28(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) )
190            ELSE
191               hd28(ji,jj) = 0.
192            ENDIF
193         END DO
194      END DO
195      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
196
197      ! ----------------------------------------- !
198      !  Heat content of first 300 m (18 levels)  !
199      ! ----------------------------------------- !
200
201      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...)
202      ilevel = 0
203      zthick_0 = 0.e0
204      DO jk = 1, jpk-1                       
205         zthick_0 = zthick_0 + e3t_0(jk)
206         IF( zthick_0 < 300. )   ilevel = jk
207      END DO
208      ! surface boundary condition
209      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0.e0        ;   htc3(:,:) = 0.e0                                     
210      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)   
211      ENDIF
212      ! integration down to ilevel
213      DO jk = 1, ilevel
214         zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
215         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tn(:,:,jk) * tmask(:,:,jk)
216      END DO
217      ! deepest layer
218      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
219      htc3(:,:) = htc3(:,:) + tn(:,:,ilevel+1) * MIN( fse3t(:,:,ilevel+1), zthick(:,:) ) * tmask(:,:,ilevel+1)
220      ! from temperature to heat contain
221      zcoef = rau0 * rcp
222      htc3(:,:) = zcoef * htc3(:,:)
223      CALL iom_put( "hc300", htc3 )   ! first 300m heaat content
224
225   END SUBROUTINE dia_hth
226
227#else
228   !!----------------------------------------------------------------------
229   !!   Default option :                                       Empty module
230   !!----------------------------------------------------------------------
231   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
232CONTAINS
233   SUBROUTINE dia_hth( kt )         ! Empty routine
234      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
235   END SUBROUTINE dia_hth
236#endif
237
238   !!======================================================================
239END MODULE diahth
Note: See TracBrowser for help on using the repository browser.