source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 18 months ago

Full set of changes as in the original branch

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