New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfmxl.F90 in NEMO/branches/UKMO/r12083_mix-lyr_diag/src/OCE/ZDF – NEMO

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

Last change on this file since 12479 was 12479, checked in by jcastill, 2 years ago

Changes as in the original branch but updated to vn4.1

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