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/NEMO_4.0.4_CO9_package-mix-lyr/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_package-mix-lyr/src/OCE/ZDF/zdfmxl.F90 @ 15505

Last change on this file since 15505 was 15505, checked in by jcastill, 8 months ago

Changes in zdfmxl as in 4.0.1 version of the branch: NEMO/branches/UKMO/r12083_mix-lyr_diag@12836

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