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/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 7367

Last change on this file since 7367 was 7367, checked in by deazer, 8 years ago

Contains merged code for CO5 reference.
File size: 18.9 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   !!----------------------------------------------------------------------
9   !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
10   !!----------------------------------------------------------------------
11   USE oce             ! ocean dynamics and tracers variables
12   USE dom_oce         ! ocean space and time domain variables
13   USE zdf_oce         ! ocean vertical physics
14   USE in_out_manager  ! I/O manager
15   USE prtctl          ! Print control
16   USE iom             ! I/O library
17   USE eosbn2          ! Equation of state
18   USE phycst, ONLY : rau0 ! reference density
19   USE lbclnk 
20   USE lib_mpp         ! MPP library
21   USE wrk_nemo        ! work arrays
22   USE timing          ! Timing
23   USE trc_oce, ONLY : lk_offline ! offline flag
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   zdf_mxl       ! called by step.F90
29   
30   ! Namelist variables for  namzdf_karaml
31 
32   LOGICAL   :: ln_kara            ! Logical switch for calculating kara mixed
33                                     ! layer
34   LOGICAL   :: ln_kara_write25h   ! Logical switch for 25 hour outputs
35   INTEGER   :: jpmld_type         ! mixed layer type             
36   REAL(wp)  :: ppz_ref            ! depth of initial T_ref
37   REAL(wp)  :: ppdT_crit          ! Critical temp diff 
38   REAL(wp)  :: ppiso_frac         ! Fraction of ppdT_crit used
39   
40   !Used for 25h mean
41   LOGICAL, PRIVATE :: kara_25h_init = .TRUE.   !Logical used to initalise 25h
42                                                !outputs. Necissary, because we need to
43                                                !initalise the kara_25h on the zeroth
44                                                !timestep (i.e in the nemogcm_init call)
45   REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h
46   
47   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_kara  !: mixed layer depth of Kara et al   [m]
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m]
53   
54   !! * Substitutions
55#  include "domzgr_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   INTEGER FUNCTION zdf_mxl_alloc()
64      !!----------------------------------------------------------------------
65      !!               ***  FUNCTION zdf_mxl_alloc  ***
66      !!----------------------------------------------------------------------
67      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated
68      IF( .NOT. ALLOCATED( nmln ) ) THEN
69         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), & 
70                                            hmld_tref(jpi,jpj), STAT= zdf_mxl_alloc ) 
71         !
72         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
73         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.')
74         !
75      ENDIF
76   END FUNCTION zdf_mxl_alloc
77
78
79   SUBROUTINE zdf_mxl( kt )
80      !!----------------------------------------------------------------------
81      !!                  ***  ROUTINE zdfmxl  ***
82      !!                   
83      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
84      !!      with a simple density criteria. Also calculates the mixed layer
85      !!      depth of Kara et al by calling zdf_mxl_kara.
86      !!
87      !! ** Method  :   The mixed layer depth is the shallowest W depth with
88      !!      the density of the corresponding T point (just bellow) bellow a
89      !!      given value defined locally as rho(10m) + zrho_c
90      !!               The turbocline depth is the depth at which the vertical
91      !!      eddy diffusivity coefficient (resulting from the vertical physics
92      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
93      !!      value defined locally (avt_c here taken equal to 5 cm/s2)
94      !!
95      !! ** Action  :   nmln, hmld, hmlp, hmlpt
96      !!----------------------------------------------------------------------
97      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
98      !!
99      INTEGER  ::   ji, jj, jk          ! dummy loop indices
100      INTEGER  ::   iikn, iiki          ! temporary integer within a do loop
101      INTEGER, POINTER, DIMENSION(:,:) ::   imld                ! temporary workspace
102      REAL(wp) ::   zrho_c = 0.01_wp    ! density criterion for mixed layer depth
103      REAL(wp) ::   zavt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
104      REAL(wp) ::   t_ref               ! Reference temperature 
105      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth 
106      !!----------------------------------------------------------------------
107      !
108      IF( nn_timing == 1 )  CALL timing_start('zdf_mxl')
109      !
110      CALL wrk_alloc( jpi,jpj, imld )
111
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
115         IF(lwp) WRITE(numout,*) '~~~~~~~ '
116         !                             ! allocate zdfmxl arrays
117         IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
118      ENDIF
119
120      ! w-level of the mixing and mixed layers
121      nmln(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
122      imld(:,:) = mbkt(:,:) + 1
123      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
124         DO jj = 1, jpj
125            DO ji = 1, jpi
126               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + zrho_c )   nmln(ji,jj) = jk      ! Mixed layer
127               IF( avt (ji,jj,jk) < zavt_c                     )   imld(ji,jj) = jk      ! Turbocline
128            END DO
129         END DO
130      END DO
131      ! depth of the mixing and mixed layers
132     
133      CALL zdf_mxl_kara( kt ) 
134     
135      DO jj = 1, jpj
136         DO ji = 1, jpi
137            iiki = imld(ji,jj)
138            iikn = nmln(ji,jj)
139            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * tmask(ji,jj,1)    ! Turbocline depth
140            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * tmask(ji,jj,1)    ! Mixed layer depth
141            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1)                     ! depth of the last T-point inside the mixed layer
142         END DO
143      END DO
144#if defined key_iomput 
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#endif
150
151      !For the AMM model assimiation uses a temperature based mixed layer depth 
152      !This is defined here 
153      DO jj = 1, jpj 
154         DO ji = 1, jpi 
155           hmld_tref(ji,jj)=fsdept(ji,jj,1  )   
156           IF(tmask(ji,jj,1) > 0.)THEN 
157             t_ref=tsn(ji,jj,1,jp_tem) 
158             DO jk=2,jpk 
159               IF(tmask(ji,jj,jk)==0.)THEN 
160                  hmld_tref(ji,jj)=fsdept(ji,jj,jk ) 
161                  EXIT 
162               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN 
163                  hmld_tref(ji,jj)=fsdept(ji,jj,jk ) 
164               ELSE 
165                  EXIT 
166               ENDIF 
167             ENDDO 
168           ENDIF 
169         ENDDO 
170      ENDDO
171           
172      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
173      !
174      CALL wrk_dealloc( jpi,jpj, imld )
175      !
176      IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl')
177      !
178   END SUBROUTINE zdf_mxl
179   
180   
181   SUBROUTINE zdf_mxl_kara( kt ) 
182      !!----------------------------------------------------------------------------------
183      !!                    ***  ROUTINE zdf_mxl_kara  ***
184      !                                                                       
185      !   Calculate mixed layer depth according to the definition of         
186      !   Kara et al, 2000, JGR, 105, 16803. 
187      !
188      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the 
189      !   density has increased by an amount equivalent to a temperature difference of 
190      !   0.8C at the surface.
191      !
192      !   For other values of mld_type the mixed layer is calculated as the depth at 
193      !   which the temperature differs by 0.8C from the surface temperature. 
194      !                                                                       
195      !   Original version: David Acreman                                     
196      !
197      !!-----------------------------------------------------------------------------------
198     
199      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
200 
201      NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, &
202      &                       ppiso_frac, ln_kara_write25h 
203 
204      ! Local variables                                                       
205      REAL, DIMENSION(jpi,jpk) :: ppzdep      ! depth for use in calculating d(rho)
206      REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d  !Local version of tsn
207      LOGICAL :: ll_found(jpi,jpj)              ! Is T_b to be found by interpolation ?
208      LOGICAL :: ll_belowml(jpi,jpj,jpk)        ! Flag points below mixed layer when ll_found=F
209      INTEGER :: ji, jj, jk                     ! loop counter
210      INTEGER :: ik_ref(jpi,jpj)                ! index of reference level
211      INTEGER :: ik_iso(jpi,jpj)                ! index of last uniform temp level
212      REAL    :: zT(jpi,jpj,jpk)                ! Temperature or denisty
213      REAL    :: zT_ref(jpi,jpj)                ! reference temperature
214      REAL    :: zT_b                           ! base temperature
215      REAL    :: zdTdz(jpi,jpj,jpk-2)           ! gradient of zT
216      REAL    :: zmoddT(jpi,jpj,jpk-2)          ! Absolute temperature difference
217      REAL    :: zdz                            ! depth difference
218      REAL    :: zdT                            ! temperature difference
219      REAL    :: zdelta_T(jpi,jpj)              ! difference critereon
220      REAL    :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities
221      INTEGER, SAVE :: idt, i_steps
222      INTEGER, SAVE :: i_cnt_25h     ! Counter for 25 hour means
223     
224 
225      !!-------------------------------------------------------------------------------------
226 
227      IF( kt == nit000 ) THEN 
228         ! Default values
229         ln_kara      = .FALSE.
230         ln_kara_write25h = .FALSE. 
231         jpmld_type   = 1     
232         ppz_ref      = 10.0 
233         ppdT_crit    = 0.2 
234         ppiso_frac   = 0.1   
235         ! Read namelist
236         REWIND ( numnam )   
237         READ   ( numnam, namzdf_karaml ) 
238         WRITE(numout,*) '===== Kara mixed layer =====' 
239         WRITE(numout,*) 'ln_kara = ',    ln_kara
240         WRITE(numout,*) 'jpmld_type = ', jpmld_type 
241         WRITE(numout,*) 'ppz_ref = ',    ppz_ref 
242         WRITE(numout,*) 'ppdT_crit = ',  ppdT_crit 
243         WRITE(numout,*) 'ppiso_frac = ', ppiso_frac
244         WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h
245         WRITE(numout,*) '============================' 
246     
247         IF ( .NOT.ln_kara ) THEN
248            WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated" 
249         ELSE
250            IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) )
251            IF ( ln_kara_write25h .AND. kara_25h_init ) THEN
252               i_cnt_25h = 0
253               IF (.NOT. ALLOCATED(hmld_kara_25h) ) &
254               &   ALLOCATE( hmld_kara_25h(jpi,jpj) )
255               hmld_kara_25h = 0._wp
256               IF( nacc == 1 ) THEN
257                  idt = INT(rdtmin)
258               ELSE
259                  idt = INT(rdt)
260               ENDIF
261               IF( MOD( 3600,idt) == 0 ) THEN
262                  i_steps = 3600 / idt 
263               ELSE
264                  CALL ctl_stop('STOP', &
265                  & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// &
266                  & ' = 0 otherwise no hourly values are possible') 
267               ENDIF 
268            ENDIF
269         ENDIF
270      ENDIF
271     
272      IF ( ln_kara ) THEN
273         
274         !set critical ln_kara
275         ztsn_2d = tsn(:,:,1,:)
276         ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit
277     
278         ! Set the mixed layer depth criterion at each grid point
279         ppzdep = 0._wp
280         IF( jpmld_type == 1 ) THEN                                         
281            CALL eos ( tsn(:,:,1,:), &
282            &          ppzdep(:,:), zRHO1(:,:) ) 
283            CALL eos ( ztsn_2d(:,:,:), &
284            &           ppzdep(:,:), zRHO2(:,:) ) 
285            zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
286            ! RHO from eos (2d version) doesn't calculate north or east halo:
287            CALL lbc_lnk( zdelta_T, 'T', 1. ) 
288            zT(:,:,:) = rhop(:,:,:) 
289         ELSE
290            zdelta_T(:,:) = ppdT_crit                     
291            zT(:,:,:) = tsn(:,:,:,jp_tem)                           
292         ENDIF 
293     
294         ! Calculate the gradient of zT and absolute difference for use later
295         DO jk = 1 ,jpk-2 
296            zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 
297            zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
298         END DO 
299     
300         ! Find density/temperature at the reference level (Kara et al use 10m).         
301         ! ik_ref is the index of the box centre immediately above or at the reference level
302         ! Find ppz_ref in the array of model level depths and find the ref   
303         ! density/temperature by linear interpolation.                                   
304         ik_ref = -1
305         DO jk = jpkm1, 2, -1 
306            WHERE( fsdept(:,:,jk) > ppz_ref ) 
307               ik_ref(:,:) = jk - 1 
308               zT_ref(:,:) = zT(:,:,jk-1) + &
309               &             zdTdz(:,:,jk-1) * ( ppz_ref - fsdept(:,:,jk-1) ) 
310            ENDWHERE 
311         END DO
312         IF ( ANY( ik_ref  < 0 ) .OR. ANY( ik_ref  > jpkm1 ) ) THEN
313            CALL ctl_stop( "STOP", &
314            & "zdf_mxl_kara: unable to find reference level for kara ML" ) 
315         ELSE
316            ! If the first grid box centre is below the reference level then use the
317            ! top model level to get zT_ref
318            WHERE( fsdept(:,:,1) > ppz_ref ) 
319               zT_ref = zT(:,:,1) 
320               ik_ref = 1 
321            ENDWHERE 
322     
323            ! Search for a uniform density/temperature region where adjacent levels         
324            ! differ by less than ppiso_frac * deltaT.                                     
325            ! ik_iso is the index of the last level in the uniform layer 
326            ! ll_found indicates whether the mixed layer depth can be found by interpolation
327            ik_iso(:,:)   = ik_ref(:,:) 
328            ll_found(:,:) = .false. 
329            DO jj = 1, nlcj 
330               DO ji = 1, nlci 
331                 !CDIR NOVECTOR
332                  DO jk = ik_ref(ji,jj), mbathy(ji,jj)-1 
333                     IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN
334                        ik_iso(ji,jj)   = jk 
335                        ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
336                        EXIT
337                     ENDIF
338                  END DO
339               END DO
340            END DO 
341     
342            ! Use linear interpolation to find depth of mixed layer base where possible
343            hmld_kara(:,:) = ppz_ref 
344            DO jj = 1, jpj 
345               DO ji = 1, jpi 
346                  IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN
347                     zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
348                     hmld_kara(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 
349                  ENDIF
350               END DO
351            END DO 
352     
353            ! If ll_found = .false. then calculate MLD using difference of zdelta_T   
354            ! from the reference density/temperature
355     
356            ! Prevent this section from working on land points
357            WHERE( tmask(:,:,1) /= 1.0 ) 
358               ll_found = .true. 
359            ENDWHERE 
360     
361            DO jk = 1, jpk 
362               ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= &
363               & zdelta_T(:,:)
364            END DO 
365     
366            ! Set default value where interpolation cannot be used (ll_found=false) 
367            DO jj = 1, jpj 
368               DO ji = 1, jpi 
369                  IF( .NOT. ll_found(ji,jj) )  &
370                  &   hmld_kara(ji,jj) = fsdept(ji,jj,mbathy(ji,jj)) 
371               END DO
372            END DO
373     
374            DO jj = 1, jpj 
375               DO ji = 1, jpi 
376                  !CDIR NOVECTOR
377                  DO jk = ik_ref(ji,jj)+1, mbathy(ji,jj) 
378                     IF( ll_found(ji,jj) ) EXIT
379                     IF( ll_belowml(ji,jj,jk) ) THEN               
380                        zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * &
381                        &      SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
382                        zdT  = zT_b - zT(ji,jj,jk-1)                                     
383                        zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
384                        hmld_kara(ji,jj) = fsdept(ji,jj,jk-1) + zdz 
385                        EXIT                                                   
386                     ENDIF
387                  END DO
388               END DO
389            END DO
390     
391            hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1) 
392 
393            IF(  ln_kara_write25h  ) THEN
394               !Double IF required as i_steps not defined if ln_kara_write25h =
395               ! FALSE
396               IF ( ( MOD( kt, i_steps ) == 0 ) .OR.  kara_25h_init ) THEN
397                  hmld_kara_25h = hmld_kara_25h + hmld_kara
398                  i_cnt_25h = i_cnt_25h + 1
399                  IF ( kara_25h_init ) kara_25h_init = .FALSE.
400               ENDIF
401            ENDIF
402 
403#if defined key_iomput 
404            IF( kt /= nit000 ) THEN
405               CALL iom_put( "mldkara"  , hmld_kara )   
406               IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND.  ln_kara_write25h ) &
407                  CALL iom_put( "kara25h"  , ( hmld_kara_25h / 25._wp ) )
408            ENDIF
409#endif
410 
411         ENDIF
412      ENDIF
413       
414   END SUBROUTINE zdf_mxl_kara 
415   
416   !!======================================================================
417END MODULE zdfmxl
Note: See TracBrowser for help on using the repository browser.