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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 17.5 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
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   dia_hth       ! routine called by step.F90
30   PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90
31
32   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
33   ! note: following variables should move to local variables once iom_put is always used
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 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
39   !! * Control permutation of array indices
40#  include "oce_ftrans.h90"
41#  include "dom_oce_ftrans.h90"
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   FUNCTION dia_hth_alloc()
53      !!---------------------------------------------------------------------
54      INTEGER :: dia_hth_alloc
55      !!---------------------------------------------------------------------
56      !
57      ALLOCATE(hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc)
58      !
59      IF( lk_mpp           )   CALL mpp_sum ( dia_hth_alloc )
60      IF(dia_hth_alloc /= 0)   CALL ctl_warn('dia_hth_alloc: failed to allocate arrays.')
61      !
62   END FUNCTION dia_hth_alloc
63
64
65   SUBROUTINE dia_hth( kt )
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      !!
87      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
88      INTEGER                          ::   iid, ilevel           ! temporary integers
89      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels
90      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth
91      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
92      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
93      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
94      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars
95      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
96      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
97      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2
98      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
99      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3     
100      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
101      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion
102      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion
103      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
104      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
105      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz
106      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness
107      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
108      !!----------------------------------------------------------------------
109
110      IF( kt == nit000 ) THEN
111         !                                      ! allocate dia_hth array
112         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
113
114         IF(.not. ALLOCATED(ik20))THEN
115            ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), &
116               &      zabs2(jpi,jpj),   &
117               &      ztm2(jpi,jpj),    &
118               &      zrho10_3(jpi,jpj),&
119               &      zpycn(jpi,jpj),   &
120               &      ztinv(jpi,jpj),   &
121               &      zdepinv(jpi,jpj), &
122               &      zrho0_3(jpi,jpj), &
123               &      zrho0_1(jpi,jpj), &
124               &      zmaxdzT(jpi,jpj), &
125               &      zthick(jpi,jpj),  &
126               &      zdelr(jpi,jpj), STAT=ji)
127            IF( lk_mpp  )   CALL mpp_sum(ji)
128            IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' )
129         END IF
130
131         IF(lwp) WRITE(numout,*)
132         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
133         IF(lwp) WRITE(numout,*) '~~~~~~~ '
134         IF(lwp) WRITE(numout,*)
135      ENDIF
136
137      ! initialization
138      ztinv  (:,:) = 0._wp 
139      zdepinv(:,:) = 0._wp 
140      zmaxdzT(:,:) = 0._wp 
141      DO jj = 1, jpj
142         DO ji = 1, jpi
143            zztmp = bathy(ji,jj)
144            hth     (ji,jj) = zztmp
145            zabs2   (ji,jj) = zztmp
146            ztm2    (ji,jj) = zztmp
147            zrho10_3(ji,jj) = zztmp
148            zpycn   (ji,jj) = zztmp
149        END DO
150      END DO
151      IF( nla10 > 1 ) THEN
152         DO jj = 1, jpj
153            DO ji = 1, jpi
154               zztmp = bathy(ji,jj)
155               zrho0_3(ji,jj) = zztmp
156               zrho0_1(ji,jj) = zztmp
157            END DO
158         END DO
159      ENDIF
160     
161      ! Preliminary computation
162      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
163      DO jj = 1, jpj
164         DO ji = 1, jpi
165            IF( tmask(ji,jj,nla10) == 1. ) THEN
166               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)   &
167                  &                                                           - 0.0100*tn(ji,jj,nla10)*sn(ji,jj,nla10)
168               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)
169               zut =    11.25 -  0.149*tn(ji,jj,nla10) - 0.01*sn(ji,jj,nla10)
170               zvt =    38.00 -  0.750*tn(ji,jj,nla10)
171               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
172               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
173            ELSE
174               zdelr(ji,jj) = 0._wp
175            ENDIF
176         END DO
177      END DO
178
179      ! ------------------------------------------------------------- !
180      ! thermocline depth: strongest vertical gradient of temperature !
181      ! turbocline depth (mixing layer depth): avt = zavt5            !
182      ! MLD: rho = rho(1) + zrho3                                     !
183      ! MLD: rho = rho(1) + zrho1                                     !
184      ! ------------------------------------------------------------- !
185#if defined key_z_first
186      DO jj = 1, jpj
187         DO ji = 1, jpi
188            DO jk = jpkm1, 2, -1   ! loop from bottom to 2
189#else
190      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
191         DO jj = 1, jpj
192            DO ji = 1, jpi
193#endif
194               !
195               zzdep = fsdepw(ji,jj,jk)
196               zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
197               zzdep = zzdep * tmask(ji,jj,1)
198
199               IF( zztmp > zmaxdzT(ji,jj) ) THEN                       
200                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
201               ENDIF
202               
203               IF( nla10 > 1 ) THEN
204                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
205                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
206                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
207               ENDIF
208
209            END DO
210         END DO
211      END DO
212     
213      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
214      IF( nla10 > 1 ) THEN
215         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
216         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
217      ENDIF
218
219      ! ------------------------------------------------------------- !
220      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
221      ! Top of thermocline: tn = tn(10m) - ztem2                      !
222      ! MLD: rho = rho10m + zrho3                                     !
223      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
224      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
225      ! depth of temperature inversion                                !
226      ! ------------------------------------------------------------- !
227#if defined key_z_first
228      DO jj = 1, jpj
229         DO ji = 1, jpi
230            DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
231#else
232      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
233         DO jj = 1, jpj
234            DO ji = 1, jpi
235#endif
236               !
237               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
238               !
239               zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk)                  ! - delta T(10m)
240               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
241               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
242               zztmp = -zztmp                                          ! delta T(10m)
243               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
244                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
245               ENDIF
246
247               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
248               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
249               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
250               !
251            END DO
252         END DO
253      END DO
254
255      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2
256      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
257      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
258      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
259      CALL iom_put( "BLT"     , ztm2 - zpycn )   ! Barrier Layer Thickness
260      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)
261      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)
262
263
264      ! ----------------------------------- !
265      ! search deepest level above 20C/28C  !
266      ! ----------------------------------- !
267      ik20(:,:) = 1
268      ik28(:,:) = 1
269#if defined key_z_first
270      DO jj = 1, jpj
271         DO ji = 1, jpi
272            DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth =>
273               !               ! loop from top to bottom
274#else
275      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
276         DO jj = 1, jpj
277            DO ji = 1, jpi
278#endif
279               zztmp = tn(ji,jj,jk)
280               IF( zztmp >= 20. )   ik20(ji,jj) = jk
281               IF( zztmp >= 28. )   ik28(ji,jj) = jk
282            END DO
283         END DO
284      END DO
285
286      ! --------------------------- !
287      !  Depth of 20C/28C isotherm  !
288      ! --------------------------- !
289      DO jj = 1, jpj
290         DO ji = 1, jpi
291            !
292            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
293            !
294            iid = ik20(ji,jj)
295            IF( iid /= 1 ) THEN
296               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
297                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
298                  &  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
299                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
300               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
301            ELSE
302               hd20(ji,jj) = 0._wp
303            ENDIF
304            !
305            iid = ik28(ji,jj)
306            IF( iid /= 1 ) THEN
307               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
308                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
309                  &  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
310                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
311               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
312            ELSE
313               hd28(ji,jj) = 0._wp
314            ENDIF
315
316         END DO
317      END DO
318      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
319      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
320
321      ! ----------------------------- !
322      !  Heat content of first 300 m  !
323      ! ----------------------------- !
324
325      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...)
326      ilevel   = 0
327      zthick_0 = 0._wp
328      DO jk = 1, jpkm1                     
329         zthick_0 = zthick_0 + e3t_0(jk)
330         IF( zthick_0 < 300. )   ilevel = jk
331      END DO
332      ! surface boundary condition
333      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
334      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)   
335      ENDIF
336      ! integration down to ilevel
337      DO jk = 1, ilevel
338         zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
339         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tn(:,:,jk) * tmask(:,:,jk)
340      END DO
341      ! deepest layer
342      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
343      DO jj = 1, jpj
344         DO ji = 1, jpi
345            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)
346         END DO
347      END DO
348      ! from temperature to heat contain
349      zcoef = rau0 * rcp
350      htc3(:,:) = zcoef * htc3(:,:)
351      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
352      !
353   END SUBROUTINE dia_hth
354
355#else
356   !!----------------------------------------------------------------------
357   !!   Default option :                                       Empty module
358   !!----------------------------------------------------------------------
359   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
360CONTAINS
361   SUBROUTINE dia_hth( kt )         ! Empty routine
362      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
363   END SUBROUTINE dia_hth
364#endif
365
366   !!======================================================================
367END MODULE diahth
Note: See TracBrowser for help on using the repository browser.