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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 17.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   !!   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#if   defined key_diahth
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      ENDIF
134
135      ! initialization
136      ztinv  (:,:) = 0._wp 
137      zdepinv(:,:) = 0._wp 
138      zmaxdzT(:,:) = 0._wp 
139      DO jj = 1, jpj
140         DO ji = 1, jpi
141            zztmp = bathy(ji,jj)
142            hth     (ji,jj) = zztmp
143            zabs2   (ji,jj) = zztmp
144            ztm2    (ji,jj) = zztmp
145            zrho10_3(ji,jj) = zztmp
146            zpycn   (ji,jj) = zztmp
147        END DO
148      END DO
149      IF( nla10 > 1 ) THEN
150         DO jj = 1, jpj
151            DO ji = 1, jpi
152               zztmp = bathy(ji,jj)
153               zrho0_3(ji,jj) = zztmp
154               zrho0_1(ji,jj) = zztmp
155            END DO
156         END DO
157      ENDIF
158     
159      ! Preliminary computation
160      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
161      DO jj = 1, jpj
162         DO ji = 1, jpi
163            IF( tmask(ji,jj,nla10) == 1. ) THEN
164               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
165                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
166                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
167               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
168                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
169               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
170               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
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      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
186         DO jj = 1, jpj
187            DO ji = 1, jpi
188               !
189               zzdep = fsdepw(ji,jj,jk)
190               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)
191               zzdep = zzdep * tmask(ji,jj,1)
192
193               IF( zztmp > zmaxdzT(ji,jj) ) THEN                       
194                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
195               ENDIF
196               
197               IF( nla10 > 1 ) THEN
198                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
199                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
200                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
201               ENDIF
202
203            END DO
204         END DO
205      END DO
206     
207      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
208      IF( nla10 > 1 ) THEN
209         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
210         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
211      ENDIF
212
213      ! ------------------------------------------------------------- !
214      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
215      ! Top of thermocline: tn = tn(10m) - ztem2                      !
216      ! MLD: rho = rho10m + zrho3                                     !
217      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
218      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
219      ! depth of temperature inversion                                !
220      ! ------------------------------------------------------------- !
221      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
222         DO jj = 1, jpj
223            DO ji = 1, jpi
224               !
225               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
226               !
227               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
228               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
229               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
230               zztmp = -zztmp                                          ! delta T(10m)
231               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
232                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
233               ENDIF
234
235               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
236               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
237               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
238               !
239            END DO
240         END DO
241      END DO
242
243      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2
244      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
245      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
246      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
247      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)
248      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)
249
250
251      ! ----------------------------------- !
252      ! search deepest level above 20C/28C  !
253      ! ----------------------------------- !
254      ik20(:,:) = 1
255      ik28(:,:) = 1
256      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
257         DO jj = 1, jpj
258            DO ji = 1, jpi
259               zztmp = tsn(ji,jj,jk,jp_tem)
260               IF( zztmp >= 20. )   ik20(ji,jj) = jk
261               IF( zztmp >= 28. )   ik28(ji,jj) = jk
262            END DO
263         END DO
264      END DO
265
266      ! --------------------------- !
267      !  Depth of 20C/28C isotherm  !
268      ! --------------------------- !
269      DO jj = 1, jpj
270         DO ji = 1, jpi
271            !
272            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
273            !
274            iid = ik20(ji,jj)
275            IF( iid /= 1 ) THEN
276               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
277                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
278                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
279                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
280               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
281            ELSE
282               hd20(ji,jj) = 0._wp
283            ENDIF
284            !
285            iid = ik28(ji,jj)
286            IF( iid /= 1 ) THEN
287               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
288                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
289                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
290                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
291               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
292            ELSE
293               hd28(ji,jj) = 0._wp
294            ENDIF
295
296         END DO
297      END DO
298      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
299      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
300
301      ! ----------------------------- !
302      !  Heat content of first 300 m  !
303      ! ----------------------------- !
304
305      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...)
306      ilevel   = 0
307      zthick_0 = 0._wp
308      DO jk = 1, jpkm1                     
309         zthick_0 = zthick_0 + e3t_1d(jk)
310         IF( zthick_0 < 300. )   ilevel = jk
311      END DO
312      ! surface boundary condition
313      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
314      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
315      ENDIF
316      ! integration down to ilevel
317      DO jk = 1, ilevel
318         zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
319         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
320      END DO
321      ! deepest layer
322      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
323      DO jj = 1, jpj
324         DO ji = 1, jpi
325            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )  &
326                                                                   * tmask(ji,jj,ilevel+1)
327         END DO
328      END DO
329      ! from temperature to heat contain
330      zcoef = rau0 * rcp
331      htc3(:,:) = zcoef * htc3(:,:)
332      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
333      !
334      IF( nn_timing == 1 )   CALL timing_stop('dia_hth')
335      !
336   END SUBROUTINE dia_hth
337
338#else
339   !!----------------------------------------------------------------------
340   !!   Default option :                                       Empty module
341   !!----------------------------------------------------------------------
342   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
343CONTAINS
344   SUBROUTINE dia_hth( kt )         ! Empty routine
345      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
346   END SUBROUTINE dia_hth
347#endif
348
349   !!======================================================================
350END MODULE diahth
Note: See TracBrowser for help on using the repository browser.