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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zdfmxl.F90 in branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

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

Last change on this file since 9473 was 9473, checked in by deazer, 6 years ago
Bug fix for 25h mean, needs reinitialzation of the mean variable at midnight
File size: 25.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: