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 NEMO/branches/UKMO/r12083_mix-lyr_diag/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/r12083_mix-lyr_diag/src/OCE/ZDF/zdfmxl.F90 @ 12585

Last change on this file since 12585 was 12585, checked in by jcastill, 4 years ago

Update to CO8 version - the previous one crashed after a divide by zero error

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