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_withNOOS/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_withNOOS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 9474

Last change on this file since 9474 was 9474, checked in by deazer, 4 years ago

25h mean bugfix needed for midnight crossover

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