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 branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 9473

Last change on this file since 9473 was 9473, checked in by deazer, 5 years ago

Bug fix for 25h mean, needs reinitialzation of the mean variable at midnight

File size: 25.0 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
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   !: mixed layer depth at t-points        [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)         ::   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. Necassary, 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
57   INTEGER, SAVE :: i_cnt_25h                   ! Counter for 25 hour means
58
59   REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   hmld_zint_25h
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
65   !! $Id$
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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),       &
77        &          htc_mld(jpi,jpj),                                                                    &
78        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
79         !
80         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
81         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.')
82         !
83      ENDIF
84   END FUNCTION zdf_mxl_alloc
85
86
87   SUBROUTINE zdf_mxl( kt )
88      !!----------------------------------------------------------------------
89      !!                  ***  ROUTINE zdfmxl  ***
90      !!                   
91      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
92      !!              with density criteria.
93      !!
94      !! ** Method  :   The mixed layer depth is the shallowest W depth with
95      !!      the density of the corresponding T point (just bellow) bellow a
96      !!      given value defined locally as rho(10m) + rho_c
97      !!               The turbocline depth is the depth at which the vertical
98      !!      eddy diffusivity coefficient (resulting from the vertical physics
99      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
100      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
101      !!
102      !! ** Action  :   nmln, hmld, hmlp, hmlpt
103      !!----------------------------------------------------------------------
104      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
105      !
106      INTEGER  ::   ji, jj, jk   ! dummy loop indices
107      INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer
108      REAL(wp) ::   zN2_c        ! local scalar
109      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start('zdf_mxl')
113      !
114      CALL wrk_alloc( jpi,jpj, imld )
115
116      IF( kt == nit000 ) THEN
117         IF(lwp) WRITE(numout,*)
118         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
119         IF(lwp) WRITE(numout,*) '~~~~~~~ '
120         !                             ! allocate zdfmxl arrays
121         IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
122      ENDIF
123
124      ! w-level of the mixing and mixed layers
125      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
126      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
127      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
128      DO jk = nlb10, jpkm1
129         DO jj = 1, jpj                ! Mixed layer level: w-level
130            DO ji = 1, jpi
131               ikt = mbkt(ji,jj)
132               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk)
133               IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
134            END DO
135         END DO
136      END DO
137      !
138      ! w-level of the turbocline
139      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
140      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
141         DO jj = 1, jpj
142            DO ji = 1, jpi
143               imkt = mikt(ji,jj)
144               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline
145            END DO
146         END DO
147      END DO
148      ! depth of the mixing and mixed layers
149      DO jj = 1, jpj
150         DO ji = 1, jpi
151            iiki = imld(ji,jj)
152            iikn = nmln(ji,jj)
153            imkt = mikt(ji,jj)
154            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth
155            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth
156            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
157         END DO
158      END DO
159      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
160         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
161         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
162      ENDIF
163     
164      ! Vertically-interpolated mixed-layer depth diagnostic
165      CALL zdf_mxl_zint( kt )
166
167      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
168      !
169      CALL wrk_dealloc( jpi,jpj, imld )
170      !
171      IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl')
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), 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, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels
208      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level
209      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level
210      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density
211      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)
212      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature
213      REAL    :: zT_b                                   ! base temperature
214      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT
215      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference
216      REAL    :: zdz                                    ! depth difference
217      REAL    :: zdT                                    ! temperature difference
218      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon
219      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities
220      INTEGER :: ji, jj, jk                             ! loop counter
221
222      !!-------------------------------------------------------------------------------------
223     
224      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
225      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
226      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 
227
228      ! Unpack structure
229      nn_mld_type = sf%mld_type
230      rn_zref     = sf%zref
231      rn_dT_crit  = sf%dT_crit
232      rn_iso_frac = sf%iso_frac
233
234      ! Set the mixed layer depth criterion at each grid point
235      IF( nn_mld_type == 0 ) THEN
236         zdelta_T(:,:) = rn_dT_crit
237         zT(:,:,:) = rhop(:,:,:)
238      ELSE IF( nn_mld_type == 1 ) THEN
239         ppzdep(:,:)=0.0 
240         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 
241! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
242! [assumes number of tracers less than number of vertical levels]
243         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 
244         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 
245         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 
246         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
247         ! RHO from eos (2d version) doesn't calculate north or east halo:
248         CALL lbc_lnk( zdelta_T, 'T', 1. ) 
249         zT(:,:,:) = rhop(:,:,:) 
250      ELSE
251         zdelta_T(:,:) = rn_dT_crit                     
252         zT(:,:,:) = tsn(:,:,:,jp_tem)                           
253      END IF 
254
255      ! Calculate the gradient of zT and absolute difference for use later
256      DO jk = 1 ,jpk-2 
257         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 
258         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
259      END DO 
260
261      ! Find density/temperature at the reference level (Kara et al use 10m).         
262      ! ik_ref is the index of the box centre immediately above or at the reference level
263      ! Find rn_zref in the array of model level depths and find the ref   
264      ! density/temperature by linear interpolation.                                   
265      DO jk = jpkm1, 2, -1 
266         WHERE ( fsdept(:,:,jk) > rn_zref ) 
267           ik_ref(:,:) = jk - 1 
268           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) ) 
269         END WHERE
270      END DO 
271
272      ! If the first grid box centre is below the reference level then use the
273      ! top model level to get zT_ref
274      WHERE ( fsdept(:,:,1) > rn_zref ) 
275         zT_ref = zT(:,:,1) 
276         ik_ref = 1 
277      END WHERE 
278
279      ! The number of active tracer levels is 1 less than the number of active w levels
280      ikmt(:,:) = mbathy(:,:) - 1 
281
282      ! Initialize / reset
283      ll_found(:,:) = .false.
284
285      IF ( rn_iso_frac - zepsilon > 0. ) THEN
286         ! Search for a uniform density/temperature region where adjacent levels         
287         ! differ by less than rn_iso_frac * deltaT.                                     
288         ! ik_iso is the index of the last level in the uniform layer 
289         ! ll_found indicates whether the mixed layer depth can be found by interpolation
290         ik_iso(:,:)   = ik_ref(:,:) 
291         DO jj = 1, nlcj 
292            DO ji = 1, nlci 
293!CDIR NOVECTOR
294               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 
295                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
296                     ik_iso(ji,jj)   = jk 
297                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
298                     EXIT
299                  END IF
300               END DO
301            END DO
302         END DO 
303
304         ! Use linear interpolation to find depth of mixed layer base where possible
305         hmld_zint(:,:) = rn_zref 
306         DO jj = 1, jpj 
307            DO ji = 1, jpi 
308               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
309                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
310                  hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 
311               END IF
312            END DO
313         END DO
314      END IF
315
316      ! If ll_found = .false. then calculate MLD using difference of zdelta_T   
317      ! from the reference density/temperature
318 
319! Prevent this section from working on land points
320      WHERE ( tmask(:,:,1) /= 1.0 ) 
321         ll_found = .true. 
322      END WHERE
323 
324      DO jk=1, jpk 
325         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 
326      END DO 
327 
328! Set default value where interpolation cannot be used (ll_found=false) 
329      DO jj = 1, jpj 
330         DO ji = 1, jpi 
331            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj)) 
332         END DO
333      END DO
334
335      DO jj = 1, jpj 
336         DO ji = 1, jpi 
337!CDIR NOVECTOR
338            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 
339               IF ( ll_found(ji,jj) ) EXIT
340               IF ( ll_belowml(ji,jj,jk) ) THEN               
341                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
342                  zdT  = zT_b - zT(ji,jj,jk-1)                                     
343                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
344                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz 
345                  EXIT                                                   
346               END IF
347            END DO
348         END DO
349      END DO
350
351      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 
352     
353      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
354      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
355      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 
356      !
357   END SUBROUTINE zdf_mxl_zint_mld
358
359   SUBROUTINE zdf_mxl_zint_htc( kt )
360      !!----------------------------------------------------------------------
361      !!                  ***  ROUTINE zdf_mxl_zint_htc  ***
362      !!
363      !! ** Purpose :   
364      !!
365      !! ** Method  :   
366      !!----------------------------------------------------------------------
367
368      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
369
370      INTEGER :: ji, jj, jk
371      INTEGER :: ikmax
372      REAL(wp) :: zc, zcoef
373      !
374      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel
375      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick
376
377      !!----------------------------------------------------------------------
378
379      IF( .NOT. ALLOCATED(ilevel) ) THEN
380         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
381         &         zthick(jpi,jpj), STAT=ji )
382         IF( lk_mpp  )   CALL mpp_sum(ji)
383         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
384      ENDIF
385
386      ! Find last whole model T level above the MLD
387      ilevel(:,:)   = 0
388      zthick_0(:,:) = 0._wp
389
390      DO jk = 1, jpkm1 
391         DO jj = 1, jpj
392            DO ji = 1, jpi                   
393               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk)
394               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk
395            END DO
396         END DO
397         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
398         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1)
399      END DO
400
401      ! Surface boundary condition
402      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                   
403      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
404      ENDIF
405
406      ! Deepest whole T level above the MLD
407      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
408
409      ! Integration down to last whole model T level
410      DO jk = 1, ikmax
411         DO jj = 1, jpj
412            DO ji = 1, jpi
413               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel
414               zthick(ji,jj) = zthick(ji,jj) + zc
415               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
416            END DO
417         END DO
418      END DO
419
420      ! Subsequent partial T level
421      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD
422
423      DO jj = 1, jpj
424         DO ji = 1, jpi
425            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  & 
426      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
427         END DO
428      END DO
429
430      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
431
432      ! Convert to heat content
433      zcoef = rau0 * rcp
434      htc_mld(:,:) = zcoef * htc_mld(:,:)
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      INTEGER :: nn_mld_diag = 0    ! number of diagnostics
453
454      INTEGER :: i_steps          ! no of timesteps per hour
455      INTEGER :: ierror             ! logical error message   
456
457      REAL(wp) :: zdt             ! timestep variable 
458
459      CHARACTER(len=1) :: cmld
460
461      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
462      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags
463
464      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
465
466      !!----------------------------------------------------------------------
467     
468      IF( kt == nit000 ) THEN
469         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist
470         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
471901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp )
472
473         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist
474         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
475902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp )
476         IF(lwm) WRITE ( numond, namzdf_mldzint )
477
478         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
479
480         mld_diags(1) = sn_mld1
481         mld_diags(2) = sn_mld2
482         mld_diags(3) = sn_mld3
483         mld_diags(4) = sn_mld4
484         mld_diags(5) = sn_mld5
485
486         IF( nn_mld_diag > 0 ) THEN
487            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
488            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
489            DO jn = 1, nn_mld_diag
490               WRITE(numout,*) 'MLD criterion',jn,':'
491               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type
492               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref
493               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit
494               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac
495            END DO
496            WRITE(numout,*) '===================================================================='
497         ENDIF
498      ENDIF
499
500      IF( nn_mld_diag > 0 ) THEN
501         DO jn = 1, nn_mld_diag
502            WRITE(cmld,'(I1)') jn
503            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN
504               CALL zdf_mxl_zint_mld( mld_diags(jn) )
505
506               IF( iom_use( "mldzint_"//cmld ) ) THEN
507                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
508               ENDIF
509
510               IF( iom_use( "mldhtc_"//cmld ) )  THEN
511                  CALL zdf_mxl_zint_htc( kt )
512                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   )
513               ENDIF
514
515               IF( iom_use( "mldzint25h_"//cmld ) ) THEN
516                  IF( .NOT. mld_25h_write ) mld_25h_write = .TRUE.
517                  zdt = rdt
518                  IF( nacc == 1 ) zdt = rdtmin
519                  IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
520                     i_steps = 3600/INT(zdt)
521                  ELSE
522                     CALL ctl_stop('STOP', 'zdf_mxl_zint 25h: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
523                  ENDIF
524                  IF( ( mld_25h_init ) .OR. ( kt == nit000 ) ) THEN
525                     i_cnt_25h = 1
526                     IF( .NOT. ALLOCATED(hmld_zint_25h) ) THEN
527                        ALLOCATE( hmld_zint_25h(jpi,jpj,nn_mld_diag), STAT=ierror )
528                        IF( ierror > 0 )  CALL ctl_stop( 'zdf_mxl_zint 25h: unable to allocate hmld_zint_25h' )   
529                     ENDIF
530                     hmld_zint_25h(:,:,jn) = hmld_zint(:,:)
531                  ENDIF
532                  IF( MOD( kt, i_steps ) == 0 .AND.  kt .NE. nn_it000 ) THEN
533                     hmld_zint_25h(:,:,jn) = hmld_zint_25h(:,:,jn) + hmld_zint(:,:)
534                  ENDIF
535                  IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
536                     CALL iom_put( "mldzint25h_"//cmld , hmld_zint_25h(:,:,jn) / 25._wp   )
537                  ENDIF
538               ENDIF
539
540            ENDIF
541         END DO
542                 
543         IF(  mld_25h_write  ) THEN
544            IF( ( MOD( kt, i_steps ) == 0 ) .OR.  mld_25h_init ) THEN
545               IF (lwp) THEN
546                  WRITE(numout,*) 'zdf_mxl_zint (25h) : Summed the following number of hourly values so far',i_cnt_25h
547          ENDIF
548               i_cnt_25h = i_cnt_25h + 1
549               IF( mld_25h_init ) mld_25h_init = .FALSE.
550            ENDIF
551            IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
552               i_cnt_25h = 1 
553               DO jn = 1, nn_mld_diag
554                     hmld_zint_25h(:,:,jn) = hmld_zint(:,:)
555               ENDDO
556            ENDIF
557         ENDIF
558                 
559      ENDIF
560
561   END SUBROUTINE zdf_mxl_zint
562
563   !!======================================================================
564END MODULE zdfmxl
Note: See TracBrowser for help on using the repository browser.