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/trunk/src/OCE/DIA – NEMO

source: NEMO/trunk/src/OCE/DIA/diahth.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

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