source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 10759

Last change on this file since 10759 was 10759, checked in by andmirek, 20 months ago

GMED 450 write output.namelist.dyn only for nprint > 2

File size: 22.6 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   PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init
31
32   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
36   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]
37   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint
38   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?
39   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F
40
41   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
42   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
43
44   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs
45      INTEGER   :: mld_type   ! mixed layer type     
46      REAL(wp)  :: zref       ! depth of initial T_ref
47      REAL(wp)  :: dT_crit    ! Critical temp diff
48      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit used
49   END TYPE MXL_ZINT
50
51   !! * Substitutions
52#  include "domzgr_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
55   !! $Id$
56   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   INTEGER FUNCTION zdf_mxl_alloc()
61      !!----------------------------------------------------------------------
62      !!               ***  FUNCTION zdf_mxl_alloc  ***
63      !!----------------------------------------------------------------------
64      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated
65      IF( .NOT. ALLOCATED( nmln ) ) THEN
66         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       &
67        &          htc_mld(jpi,jpj),                                                                    &
68        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
69         !
70         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
71         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.')
72         !
73      ENDIF
74   END FUNCTION zdf_mxl_alloc
75
76
77   SUBROUTINE zdf_mxl( kt )
78      !!----------------------------------------------------------------------
79      !!                  ***  ROUTINE zdfmxl  ***
80      !!                   
81      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
82      !!              with density criteria.
83      !!
84      !! ** Method  :   The mixed layer depth is the shallowest W depth with
85      !!      the density of the corresponding T point (just bellow) bellow a
86      !!      given value defined locally as rho(10m) + rho_c
87      !!               The turbocline depth is the depth at which the vertical
88      !!      eddy diffusivity coefficient (resulting from the vertical physics
89      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
90      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
91      !!
92      !! ** Action  :   nmln, hmld, hmlp, hmlpt
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
95      !
96      INTEGER  ::   ji, jj, jk      ! dummy loop indices
97      INTEGER  ::   iikn, iiki, ikt ! local integer
98      REAL(wp) ::   zN2_c           ! local scalar
99      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
100      !!----------------------------------------------------------------------
101      !
102      IF( nn_timing == 1 )  CALL timing_start('zdf_mxl')
103      !
104      CALL wrk_alloc( jpi,jpj, imld )
105
106      IF( kt <= nit000 ) THEN
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
109         IF(lwp) WRITE(numout,*) '~~~~~~~ '
110         !                             ! allocate zdfmxl arrays
111         IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
112      ENDIF
113
114      ! w-level of the mixing and mixed layers
115      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
116      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
117      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
118      DO jk = nlb10, jpkm1
119         DO jj = 1, jpj                ! Mixed layer level: w-level
120            DO ji = 1, jpi
121               ikt = mbkt(ji,jj)
122               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk)
123               IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
124            END DO
125         END DO
126      END DO
127      !
128      ! w-level of the turbocline
129      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
130      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
131         DO jj = 1, jpj
132            DO ji = 1, jpi
133               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline
134            END DO
135         END DO
136      END DO
137      ! depth of the mixing and mixed layers
138      DO jj = 1, jpj
139         DO ji = 1, jpi
140            iiki = imld(ji,jj)
141            iikn = nmln(ji,jj)
142            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth
143            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth
144            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
145         END DO
146      END DO
147      ! no need to output in offline mode
148      IF( .NOT.lk_offline ) THEN   
149         IF ( iom_use("mldr10_1") ) THEN
150            IF( ln_isfcav ) THEN
151               CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
152            ELSE
153               CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
154            END IF
155         END IF
156         IF ( iom_use("mldkz5") ) THEN
157            IF( ln_isfcav ) THEN
158               CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
159            ELSE
160               CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
161            END IF
162         END IF
163      ENDIF
164     
165      ! Vertically-interpolated mixed-layer depth diagnostic
166      CALL zdf_mxl_zint( kt )
167
168      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
169      !
170      CALL wrk_dealloc( jpi,jpj, imld )
171      !
172      IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl')
173      !
174   END SUBROUTINE zdf_mxl
175
176   SUBROUTINE zdf_mxl_zint_mld( sf ) 
177      !!----------------------------------------------------------------------------------
178      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***
179      !                                                                       
180      !   Calculate vertically-interpolated mixed layer depth diagnostic.
181      !           
182      !   This routine can calculate the mixed layer depth diagnostic suggested by
183      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
184      !   vertically-interpolated mixed-layer depth diagnostics with other parameter
185      !   settings set in the namzdf_mldzint namelist. 
186      !
187      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the 
188      !   density has increased by an amount equivalent to a temperature difference of 
189      !   0.8C at the surface.
190      !
191      !   For other values of mld_type the mixed layer is calculated as the depth at 
192      !   which the temperature differs by 0.8C from the surface temperature. 
193      !                                                                       
194      !   David Acreman, Daley Calvert                                     
195      !
196      !!-----------------------------------------------------------------------------------
197
198      TYPE(MXL_ZINT), INTENT(in)  :: sf
199
200      ! Diagnostic criteria
201      INTEGER   :: nn_mld_type   ! mixed layer type     
202      REAL(wp)  :: rn_zref       ! depth of initial T_ref
203      REAL(wp)  :: rn_dT_crit    ! Critical temp diff
204      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used
205
206      ! Local variables
207      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value
208      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels
209      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level
210      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level
211      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density
212      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)
213      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature
214      REAL    :: zT_b                                   ! base temperature
215      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT
216      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference
217      REAL    :: zdz                                    ! depth difference
218      REAL    :: zdT                                    ! temperature difference
219      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon
220      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities
221      INTEGER :: ji, jj, jk                             ! loop counter
222
223      !!-------------------------------------------------------------------------------------
224     
225      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
226      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
227      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 
228
229      ! Unpack structure
230      nn_mld_type = sf%mld_type
231      rn_zref     = sf%zref
232      rn_dT_crit  = sf%dT_crit
233      rn_iso_frac = sf%iso_frac
234
235      ! Set the mixed layer depth criterion at each grid point
236      IF( nn_mld_type == 0 ) THEN
237         zdelta_T(:,:) = rn_dT_crit
238         zT(:,:,:) = rhop(:,:,:)
239      ELSE IF( nn_mld_type == 1 ) THEN
240         ppzdep(:,:)=0.0 
241         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 
242! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
243! [assumes number of tracers less than number of vertical levels]
244         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 
245         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 
246         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 
247         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
248         ! RHO from eos (2d version) doesn't calculate north or east halo:
249         CALL lbc_lnk( zdelta_T, 'T', 1. ) 
250         zT(:,:,:) = rhop(:,:,:) 
251      ELSE
252         zdelta_T(:,:) = rn_dT_crit                     
253         zT(:,:,:) = tsn(:,:,:,jp_tem)                           
254      END IF 
255
256      ! Calculate the gradient of zT and absolute difference for use later
257      DO jk = 1 ,jpk-2 
258         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 
259         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
260      END DO 
261
262      ! Find density/temperature at the reference level (Kara et al use 10m).         
263      ! ik_ref is the index of the box centre immediately above or at the reference level
264      ! Find rn_zref in the array of model level depths and find the ref   
265      ! density/temperature by linear interpolation.                                   
266      DO jk = jpkm1, 2, -1 
267         WHERE ( fsdept(:,:,jk) > rn_zref ) 
268           ik_ref(:,:) = jk - 1 
269           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) ) 
270         END WHERE
271      END DO 
272
273      ! If the first grid box centre is below the reference level then use the
274      ! top model level to get zT_ref
275      WHERE ( fsdept(:,:,1) > rn_zref ) 
276         zT_ref = zT(:,:,1) 
277         ik_ref = 1 
278      END WHERE 
279
280      ! The number of active tracer levels is 1 less than the number of active w levels
281      ikmt(:,:) = mbathy(:,:) - 1 
282
283      ! Initialize / reset
284      ll_found(:,:) = .false.
285
286      IF ( rn_iso_frac - zepsilon > 0. ) THEN
287         ! Search for a uniform density/temperature region where adjacent levels         
288         ! differ by less than rn_iso_frac * deltaT.                                     
289         ! ik_iso is the index of the last level in the uniform layer 
290         ! ll_found indicates whether the mixed layer depth can be found by interpolation
291         ik_iso(:,:)   = ik_ref(:,:) 
292         DO jj = 1, nlcj 
293            DO ji = 1, nlci 
294!CDIR NOVECTOR
295               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 
296                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
297                     ik_iso(ji,jj)   = jk 
298                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
299                     EXIT
300                  END IF
301               END DO
302            END DO
303         END DO 
304
305         ! Use linear interpolation to find depth of mixed layer base where possible
306         hmld_zint(:,:) = rn_zref 
307         DO jj = 1, jpj 
308            DO ji = 1, jpi 
309               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
310                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
311                  hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 
312               END IF
313            END DO
314         END DO
315      END IF
316
317      ! If ll_found = .false. then calculate MLD using difference of zdelta_T   
318      ! from the reference density/temperature
319 
320! Prevent this section from working on land points
321      WHERE ( tmask(:,:,1) /= 1.0 ) 
322         ll_found = .true. 
323      END WHERE
324 
325      DO jk=1, jpk 
326         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 
327      END DO 
328 
329! Set default value where interpolation cannot be used (ll_found=false) 
330      DO jj = 1, jpj 
331         DO ji = 1, jpi 
332            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj)) 
333         END DO
334      END DO
335
336      DO jj = 1, jpj 
337         DO ji = 1, jpi 
338!CDIR NOVECTOR
339            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 
340               IF ( ll_found(ji,jj) ) EXIT
341               IF ( ll_belowml(ji,jj,jk) ) THEN               
342                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
343                  zdT  = zT_b - zT(ji,jj,jk-1)                                     
344                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
345                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz 
346                  EXIT                                                   
347               END IF
348            END DO
349         END DO
350      END DO
351
352      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 
353     
354      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
355      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
356      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 
357      !
358   END SUBROUTINE zdf_mxl_zint_mld
359
360   SUBROUTINE zdf_mxl_zint_htc( kt )
361      !!----------------------------------------------------------------------
362      !!                  ***  ROUTINE zdf_mxl_zint_htc  ***
363      !!
364      !! ** Purpose :   
365      !!
366      !! ** Method  :   
367      !!----------------------------------------------------------------------
368
369      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
370
371      INTEGER :: ji, jj, jk
372      INTEGER :: ikmax
373      REAL(wp) :: zc, zcoef
374      !
375      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel
376      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick
377
378      !!----------------------------------------------------------------------
379
380      IF( .NOT. ALLOCATED(ilevel) ) THEN
381         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
382         &         zthick(jpi,jpj), STAT=ji )
383         IF( lk_mpp  )   CALL mpp_sum(ji)
384         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
385      ENDIF
386
387      ! Find last whole model T level above the MLD
388      ilevel(:,:)   = 0
389      zthick_0(:,:) = 0._wp
390
391      DO jk = 1, jpkm1 
392         DO jj = 1, jpj
393            DO ji = 1, jpi                   
394               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk)
395               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk
396            END DO
397         END DO
398         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
399         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1)
400      END DO
401
402      ! Surface boundary condition
403      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                   
404      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
405      ENDIF
406
407      ! Deepest whole T level above the MLD
408      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
409
410      ! Integration down to last whole model T level
411      DO jk = 1, ikmax
412         DO jj = 1, jpj
413            DO ji = 1, jpi
414               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel
415               zthick(ji,jj) = zthick(ji,jj) + zc
416               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
417            END DO
418         END DO
419      END DO
420
421      ! Subsequent partial T level
422      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD
423
424      DO jj = 1, jpj
425         DO ji = 1, jpi
426            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  & 
427      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
428         END DO
429      END DO
430
431      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
432
433      ! Convert to heat content
434      zcoef = rau0 * rcp
435      htc_mld(:,:) = zcoef * htc_mld(:,:)
436
437   END SUBROUTINE zdf_mxl_zint_htc
438
439   SUBROUTINE zdf_mxl_zint( kt )
440      !!----------------------------------------------------------------------
441      !!                  ***  ROUTINE zdf_mxl_zint  ***
442      !!
443      !! ** Purpose :   
444      !!
445      !! ** Method  :   
446      !!----------------------------------------------------------------------
447
448      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
449
450      INTEGER :: ios
451      INTEGER :: jn
452
453      INTEGER :: nn_mld_diag = 0    ! number of diagnostics
454
455      CHARACTER(len=1) :: cmld
456
457      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
458      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags
459
460      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
461
462      !!----------------------------------------------------------------------
463     
464      IF( kt == nit000 ) THEN
465         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist
466         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
467901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp )
468
469         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist
470         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
471902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp )
472         IF(lwm .AND. nprint > 2) WRITE ( numond, namzdf_mldzint )
473
474         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
475
476         mld_diags(1) = sn_mld1
477         mld_diags(2) = sn_mld2
478         mld_diags(3) = sn_mld3
479         mld_diags(4) = sn_mld4
480         mld_diags(5) = sn_mld5
481
482         IF( lwp .AND. (nn_mld_diag > 0) ) THEN
483            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
484            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
485            DO jn = 1, nn_mld_diag
486               WRITE(numout,*) 'MLD criterion',jn,':'
487               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type
488               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref
489               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit
490               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac
491            END DO
492            WRITE(numout,*) '===================================================================='
493         ENDIF
494      ENDIF
495
496      IF( nn_mld_diag > 0 ) THEN
497         DO jn = 1, nn_mld_diag
498            WRITE(cmld,'(I1)') jn
499            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN
500               CALL zdf_mxl_zint_mld( mld_diags(jn) )
501
502               IF( iom_use( "mldzint_"//cmld ) ) THEN
503                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
504               ENDIF
505
506               IF( iom_use( "mldhtc_"//cmld ) )  THEN
507                  CALL zdf_mxl_zint_htc( kt )
508                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   )
509               ENDIF
510            ENDIF
511         END DO
512      ENDIF
513
514   END SUBROUTINE zdf_mxl_zint
515
516   !!======================================================================
517END MODULE zdfmxl
Note: See TracBrowser for help on using the repository browser.