source: branches/UKMO/r8395_mix-lyr_diag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 11302

Last change on this file since 11302 was 11302, checked in by jcastill, 15 months ago

Addition of mixed layer diagnostics, following the changes in zdfmxl.F90 in the branch UKMO/AMM15_v3_6_STABLE_package between revisions 8058 and 10162

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