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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 17.6 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 lib_mpp         ! MPP library
24   USE iom             ! I/O library
25   USE timing          ! preformance summary
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dia_hth       ! routine called by step.F90
31   PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90
32
33   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
34   ! note: following variables should move to local variables once iom_put is always used
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m]
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m]
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W]
39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   FUNCTION dia_hth_alloc()
50      !!---------------------------------------------------------------------
51      INTEGER :: dia_hth_alloc
52      !!---------------------------------------------------------------------
53      !
54      ALLOCATE(hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc)
55      !
56      IF( lk_mpp           )   CALL mpp_sum ( dia_hth_alloc )
57      IF(dia_hth_alloc /= 0)   CALL ctl_warn('dia_hth_alloc: failed to allocate arrays.')
58      !
59   END FUNCTION dia_hth_alloc
60
61
62   SUBROUTINE dia_hth( kt )
63      !!---------------------------------------------------------------------
64      !!                  ***  ROUTINE dia_hth  ***
65      !!
66      !! ** Purpose : Computes
67      !!      the mixing layer depth (turbocline): avt = 5.e-4
68      !!      the depth of strongest vertical temperature gradient
69      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01)
70      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2       
71      !!      the top of the thermochine: tn = tn(10m) - ztem2
72      !!      the pycnocline depth with density criteria equivalent to a temperature variation
73      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
74      !!      the barrier layer thickness
75      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
76      !!      the depth of the 20 degree isotherm (linear interpolation)
77      !!      the depth of the 28 degree isotherm (linear interpolation)
78      !!      the heat content of first 300 m
79      !!
80      !! ** Method :
81      !!-------------------------------------------------------------------
82      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
83      !!
84      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
85      INTEGER                          ::   iid, ilevel           ! temporary integers
86      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels
87      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth
88      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
89      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
90      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
91      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars
92      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
93      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
94      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2
95      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
96      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3     
97      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
98      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion
99      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion
100      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
101      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
102      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz
103      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness
104      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
105      !!----------------------------------------------------------------------
106      IF( nn_timing == 1 )   CALL timing_start('dia_hth')
107
108      IF( kt == nit000 ) THEN
109         !                                      ! allocate dia_hth array
110         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
111
112         IF(.not. ALLOCATED(ik20))THEN
113            ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), &
114               &      zabs2(jpi,jpj),   &
115               &      ztm2(jpi,jpj),    &
116               &      zrho10_3(jpi,jpj),&
117               &      zpycn(jpi,jpj),   &
118               &      ztinv(jpi,jpj),   &
119               &      zdepinv(jpi,jpj), &
120               &      zrho0_3(jpi,jpj), &
121               &      zrho0_1(jpi,jpj), &
122               &      zmaxdzT(jpi,jpj), &
123               &      zthick(jpi,jpj),  &
124               &      zdelr(jpi,jpj), STAT=ji)
125            IF( lk_mpp  )   CALL mpp_sum(ji)
126            IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' )
127         END IF
128
129         IF(lwp) WRITE(numout,*)
130         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
131         IF(lwp) WRITE(numout,*) '~~~~~~~ '
132         IF(lwp) WRITE(numout,*)
133         IF(lwp .AND. lflush) CALL flush(numout)
134      ENDIF
135
136      ! initialization
137      ztinv  (:,:) = 0._wp 
138      zdepinv(:,:) = 0._wp 
139      zmaxdzT(:,:) = 0._wp 
140      DO jj = 1, jpj
141         DO ji = 1, jpi
142            zztmp = bathy(ji,jj)
143            hth     (ji,jj) = zztmp
144            zabs2   (ji,jj) = zztmp
145            ztm2    (ji,jj) = zztmp
146            zrho10_3(ji,jj) = zztmp
147            zpycn   (ji,jj) = zztmp
148        END DO
149      END DO
150      IF( nla10 > 1 ) THEN
151         DO jj = 1, jpj
152            DO ji = 1, jpi
153               zztmp = bathy(ji,jj)
154               zrho0_3(ji,jj) = zztmp
155               zrho0_1(ji,jj) = zztmp
156            END DO
157         END DO
158      ENDIF
159     
160      ! Preliminary computation
161      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
162      DO jj = 1, jpj
163         DO ji = 1, jpi
164            IF( tmask(ji,jj,nla10) == 1. ) THEN
165               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
166                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
167                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
168               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
169                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
170               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
171               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
172               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
173               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
174            ELSE
175               zdelr(ji,jj) = 0._wp
176            ENDIF
177         END DO
178      END DO
179
180      ! ------------------------------------------------------------- !
181      ! thermocline depth: strongest vertical gradient of temperature !
182      ! turbocline depth (mixing layer depth): avt = zavt5            !
183      ! MLD: rho = rho(1) + zrho3                                     !
184      ! MLD: rho = rho(1) + zrho1                                     !
185      ! ------------------------------------------------------------- !
186      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
187         DO jj = 1, jpj
188            DO ji = 1, jpi
189               !
190               zzdep = fsdepw(ji,jj,jk)
191               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)
192               zzdep = zzdep * tmask(ji,jj,1)
193
194               IF( zztmp > zmaxdzT(ji,jj) ) THEN                       
195                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
196               ENDIF
197               
198               IF( nla10 > 1 ) THEN
199                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
200                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
201                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
202               ENDIF
203
204            END DO
205         END DO
206      END DO
207     
208      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
209      IF( nla10 > 1 ) THEN
210         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
211         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
212      ENDIF
213
214      ! ------------------------------------------------------------- !
215      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
216      ! Top of thermocline: tn = tn(10m) - ztem2                      !
217      ! MLD: rho = rho10m + zrho3                                     !
218      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
219      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
220      ! depth of temperature inversion                                !
221      ! ------------------------------------------------------------- !
222      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
223         DO jj = 1, jpj
224            DO ji = 1, jpi
225               !
226               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
227               !
228               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
229               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
230               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
231               zztmp = -zztmp                                          ! delta T(10m)
232               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
233                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
234               ENDIF
235
236               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
237               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
238               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
239               !
240            END DO
241         END DO
242      END DO
243
244      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2
245      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
246      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
247      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
248      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)
249      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)
250
251
252      ! ----------------------------------- !
253      ! search deepest level above 20C/28C  !
254      ! ----------------------------------- !
255      ik20(:,:) = 1
256      ik28(:,:) = 1
257      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
258         DO jj = 1, jpj
259            DO ji = 1, jpi
260               zztmp = tsn(ji,jj,jk,jp_tem)
261               IF( zztmp >= 20. )   ik20(ji,jj) = jk
262               IF( zztmp >= 28. )   ik28(ji,jj) = jk
263            END DO
264         END DO
265      END DO
266
267      ! --------------------------- !
268      !  Depth of 20C/28C isotherm  !
269      ! --------------------------- !
270      DO jj = 1, jpj
271         DO ji = 1, jpi
272            !
273            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
274            !
275            iid = ik20(ji,jj)
276            IF( iid /= 1 ) THEN
277               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
278                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
279                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
280                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
281               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
282            ELSE
283               hd20(ji,jj) = 0._wp
284            ENDIF
285            !
286            iid = ik28(ji,jj)
287            IF( iid /= 1 ) THEN
288               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
289                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
290                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
291                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
292               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
293            ELSE
294               hd28(ji,jj) = 0._wp
295            ENDIF
296
297         END DO
298      END DO
299      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
300      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
301
302      ! ----------------------------- !
303      !  Heat content of first 300 m  !
304      ! ----------------------------- !
305
306      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...)
307      ilevel   = 0
308      zthick_0 = 0._wp
309      DO jk = 1, jpkm1                     
310         zthick_0 = zthick_0 + e3t_1d(jk)
311         IF( zthick_0 < 300. )   ilevel = jk
312      END DO
313      ! surface boundary condition
314      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
315      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
316      ENDIF
317      ! integration down to ilevel
318      DO jk = 1, ilevel
319         zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
320         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
321      END DO
322      ! deepest layer
323      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
324      DO jj = 1, jpj
325         DO ji = 1, jpi
326            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )  &
327                                                                   * tmask(ji,jj,ilevel+1)
328         END DO
329      END DO
330      ! from temperature to heat contain
331      zcoef = rau0 * rcp
332      htc3(:,:) = zcoef * htc3(:,:)
333      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
334      !
335      IF( nn_timing == 1 )   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      ! ocean time-step index
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.