source: NEMO/branches/UKMO/2018_NEMO4_beta_mirror_10037_Benchmark/src/OCE/ZDF/zdfmxl.F90 @ 10280

Last change on this file since 10280 was 10280, checked in by andmirek, 2 years ago

GMED 425 merge dev_r9950_GO8_package

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