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.
zdfmxl.F90 in branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 5249

Last change on this file since 5249 was 5249, checked in by davestorkey, 7 years ago

UKMO mld_zint branch: implement vertically-interpolated MLD.

File size: 16.7 KB
Line 
1MODULE zdfmxl
2   !!======================================================================
3   !!                       ***  MODULE  zdfmxl  ***
4   !! Ocean physics: mixed layer depth
5   !!======================================================================
6   !! History :  1.0  ! 2003-08  (G. Madec)  original code
7   !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop
8   !!            3.7  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl
9   !!             -   ! 2014-02  (F. Roquet)  mixed layer depth calculated using N2 instead of rhop
10   !!----------------------------------------------------------------------
11   !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and tracers variables
14   USE dom_oce         ! ocean space and time domain variables
15   USE zdf_oce         ! ocean vertical physics
16   USE in_out_manager  ! I/O manager
17   USE prtctl          ! Print control
18   USE phycst          ! physical constants
19   USE iom             ! I/O library
20   USE eosbn2          ! for zdf_mxl_zint
21   USE lib_mpp         ! MPP library
22   USE wrk_nemo        ! work arrays
23   USE timing          ! Timing
24   USE trc_oce, ONLY : lk_offline ! offline flag
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   zdf_mxl       ! called by step.F90
30
31   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]
36   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?
37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F
38
39   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
40   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
41
42   ! Namelist variables for  namzdf_mldzint
43   INTEGER          :: nn_mld_type         ! mixed layer type           
44   REAL(wp)         :: rn_zref            ! depth of initial T_ref
45   REAL(wp)         :: rn_dT_crit          ! Critical temp diff
46   REAL(wp)         :: rn_iso_frac         ! Fraction of rn_dT_crit used
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   INTEGER FUNCTION zdf_mxl_alloc()
58      !!----------------------------------------------------------------------
59      !!               ***  FUNCTION zdf_mxl_alloc  ***
60      !!----------------------------------------------------------------------
61      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated
62      IF( .NOT. ALLOCATED( nmln ) ) THEN
63         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       &
64        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
65         !
66         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
67         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.')
68         !
69      ENDIF
70   END FUNCTION zdf_mxl_alloc
71
72
73   SUBROUTINE zdf_mxl( kt )
74      !!----------------------------------------------------------------------
75      !!                  ***  ROUTINE zdfmxl  ***
76      !!                   
77      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
78      !!              with density criteria.
79      !!
80      !! ** Method  :   The mixed layer depth is the shallowest W depth with
81      !!      the density of the corresponding T point (just bellow) bellow a
82      !!      given value defined locally as rho(10m) + rho_c
83      !!               The turbocline depth is the depth at which the vertical
84      !!      eddy diffusivity coefficient (resulting from the vertical physics
85      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
86      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
87      !!
88      !! ** Action  :   nmln, hmld, hmlp, hmlpt
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
91      !
92      INTEGER  ::   ji, jj, jk   ! dummy loop indices
93      INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer
94      REAL(wp) ::   zN2_c        ! local scalar
95      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
96      !!----------------------------------------------------------------------
97      !
98      IF( nn_timing == 1 )  CALL timing_start('zdf_mxl')
99      !
100      CALL wrk_alloc( jpi,jpj, imld )
101
102      IF( kt == nit000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
105         IF(lwp) WRITE(numout,*) '~~~~~~~ '
106         !                             ! allocate zdfmxl arrays
107         IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
108      ENDIF
109
110      ! w-level of the mixing and mixed layers
111      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
112      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
113      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
114      DO jk = nlb10, jpkm1
115         DO jj = 1, jpj                ! Mixed layer level: w-level
116            DO ji = 1, jpi
117               ikt = mbkt(ji,jj)
118               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk)
119               IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
120            END DO
121         END DO
122      END DO
123      !
124      ! w-level of the turbocline
125      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
126      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
127         DO jj = 1, jpj
128            DO ji = 1, jpi
129               imkt = mikt(ji,jj)
130               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline
131            END DO
132         END DO
133      END DO
134      ! depth of the mixing and mixed layers
135      DO jj = 1, jpj
136         DO ji = 1, jpi
137            iiki = imld(ji,jj)
138            iikn = nmln(ji,jj)
139            imkt = mikt(ji,jj)
140            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth
141            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth
142            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
143         END DO
144      END DO
145      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
146         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
147         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
148      ENDIF
149     
150      ! Vertically-interpolated mixed-layer depth diagnostic
151      IF( iom_use( "mldzint" ) ) THEN
152         CALL zdf_mxl_zint( kt )
153         CALL iom_put( "mldzint" , hmld_zint )
154      ENDIF
155
156      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
157      !
158      CALL wrk_dealloc( jpi,jpj, imld )
159      !
160      IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl')
161      !
162   END SUBROUTINE zdf_mxl
163
164   SUBROUTINE zdf_mxl_zint( kt ) 
165      !!----------------------------------------------------------------------------------
166      !!                    ***  ROUTINE zdf_mxl_zint  ***
167      !                                                                       
168      !   Calculate vertically-interpolated mixed layer depth diagnostic.
169      !           
170      !   This routine can calculate the mixed layer depth diagnostic suggested by
171      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
172      !   vertically-interpolated mixed-layer depth diagnostics with other parameter
173      !   settings set in the namzdf_mldzint namelist. 
174      !
175      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the 
176      !   density has increased by an amount equivalent to a temperature difference of 
177      !   0.8C at the surface.
178      !
179      !   For other values of mld_type the mixed layer is calculated as the depth at 
180      !   which the temperature differs by 0.8C from the surface temperature. 
181      !                                                                       
182      !   David Acreman, Daley Calvert                                     
183      !
184      !!-----------------------------------------------------------------------------------
185
186      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
187      !
188      ! Local variables
189      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels
190      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level
191      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level
192      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density
193      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)
194      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature
195      REAL    :: zT_b                                   ! base temperature
196      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT
197      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference
198      REAL    :: zdz                                    ! depth difference
199      REAL    :: zdT                                    ! temperature difference
200      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon
201      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities
202      INTEGER :: ji, jj, jk                             ! loop counter
203      INTEGER :: ios
204
205      NAMELIST/namzdf_mldzint/ nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac
206
207      !!-------------------------------------------------------------------------------------
208     
209      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
210      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
211      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 
212 
213      IF( kt == nit000 ) THEN
214         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist
215         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
216901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp )
217
218         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist
219         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
220902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp )
221         IF(lwm) WRITE ( numond, namzdf_mldzint )
222
223         WRITE(numout,*) '===== Vertically-interpolated mixed layer ====='
224         WRITE(numout,*) 'nn_mld_type = ',nn_mld_type
225         WRITE(numout,*) 'rn_zref = ',rn_zref
226         WRITE(numout,*) 'rn_dT_crit = ',rn_dT_crit
227         WRITE(numout,*) 'rn_iso_frac = ',rn_iso_frac
228         WRITE(numout,*) '==============================================='
229      ENDIF
230 
231      ! Set the mixed layer depth criterion at each grid point
232      IF (nn_mld_type == 1) THEN                                         
233         ppzdep(:,:)=0.0 
234         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 
235! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
236! [assumes number of tracers less than number of vertical levels]
237         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 
238         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 
239         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 
240         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
241         ! RHO from eos (2d version) doesn't calculate north or east halo:
242         CALL lbc_lnk( zdelta_T, 'T', 1. ) 
243         zT(:,:,:) = rhop(:,:,:) 
244      ELSE
245         zdelta_T(:,:) = rn_dT_crit                     
246         zT(:,:,:) = tsn(:,:,:,jp_tem)                           
247      END IF 
248
249      ! Calculate the gradient of zT and absolute difference for use later
250      DO jk = 1 ,jpk-2 
251         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 
252         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
253      END DO 
254
255      ! Find density/temperature at the reference level (Kara et al use 10m).         
256      ! ik_ref is the index of the box centre immediately above or at the reference level
257      ! Find rn_zref in the array of model level depths and find the ref   
258      ! density/temperature by linear interpolation.                                   
259      DO jk = jpkm1, 2, -1 
260         WHERE ( fsdept(:,:,jk) > rn_zref ) 
261           ik_ref(:,:) = jk - 1 
262           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) ) 
263         END WHERE
264      END DO 
265
266      ! If the first grid box centre is below the reference level then use the
267      ! top model level to get zT_ref
268      WHERE ( fsdept(:,:,1) > rn_zref ) 
269         zT_ref = zT(:,:,1) 
270         ik_ref = 1 
271      END WHERE 
272
273      ! The number of active tracer levels is 1 less than the number of active w levels
274      ikmt(:,:) = mbathy(:,:) - 1 
275
276      ! Search for a uniform density/temperature region where adjacent levels         
277      ! differ by less than rn_iso_frac * deltaT.                                     
278      ! ik_iso is the index of the last level in the uniform layer 
279      ! ll_found indicates whether the mixed layer depth can be found by interpolation
280      ik_iso(:,:)   = ik_ref(:,:) 
281      ll_found(:,:) = .false. 
282      DO jj = 1, nlcj 
283         DO ji = 1, nlci 
284!CDIR NOVECTOR
285            DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 
286               IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
287                  ik_iso(ji,jj)   = jk 
288                  ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
289                  EXIT
290               END IF
291            END DO
292         END DO
293      END DO 
294
295      ! Use linear interpolation to find depth of mixed layer base where possible
296      hmld_zint(:,:) = rn_zref 
297      DO jj = 1, jpj 
298         DO ji = 1, jpi 
299            IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
300               zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
301               hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 
302            END IF
303         END DO
304      END DO 
305
306      ! If ll_found = .false. then calculate MLD using difference of zdelta_T   
307      ! from the reference density/temperature
308 
309! Prevent this section from working on land points
310      WHERE ( tmask(:,:,1) /= 1.0 ) 
311         ll_found = .true. 
312      END WHERE
313 
314      DO jk=1, jpk 
315         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 
316      END DO 
317 
318! Set default value where interpolation cannot be used (ll_found=false) 
319      DO jj = 1, jpj 
320         DO ji = 1, jpi 
321            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj)) 
322         END DO
323      END DO
324
325      DO jj = 1, jpj 
326         DO ji = 1, jpi 
327!CDIR NOVECTOR
328            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 
329               IF ( ll_found(ji,jj) ) EXIT
330               IF ( ll_belowml(ji,jj,jk) ) THEN               
331                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
332                  zdT  = zT_b - zT(ji,jj,jk-1)                                     
333                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
334                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz 
335                  EXIT                                                   
336               END IF
337            END DO
338         END DO
339      END DO
340
341      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 
342     
343      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
344      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
345      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 
346      !
347   END SUBROUTINE zdf_mxl_zint
348
349   !!======================================================================
350END MODULE zdfmxl
Note: See TracBrowser for help on using the repository browser.