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

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

GMED 450 add flush after prints

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