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/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfmxl.F90 @ 15381

Last change on this file since 15381 was 15381, checked in by hadjt, 3 years ago
ZDF/zdfmxl.F90 Added kara mld, all enabled it in the regional means.
File size: 33.9 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 trc_oce , ONLY: l_offline ! ocean space and time domain variables
16	USE zdf_oce ! ocean vertical physics
17	USE eosbn2 ! for zdf_mxl_zint
18	!
19	USE in_out_manager ! I/O manager
20	USE prtctl ! Print control
21	USE phycst ! physical constants
22	USE iom ! I/O library
23	USE lib_mpp ! MPP library
24	!JT
25	USE lbclnk ! (or mpp link)
26	!JT
27
28	IMPLICIT NONE
29	PRIVATE
30
31	PUBLIC zdf_mxl ! called by zdfphy.F90
32
33	INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP)
34	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld !: mixing layer depth (turbocline) [m] (used by TOP)
35	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] (used by LDF)
36	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: depth of the last T-point inside the mixed layer [m] (used by LDF)
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	!JT
43	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld_kara !: mixed layer depth of Kara et al [m]
44
45
46
47	! Namelist variables for namzdf_karaml
48	LOGICAL :: ln_kara ! Logical switch for calculating kara mixed
49	! layer
50	LOGICAL :: ln_kara_write25h ! Logical switch for 25 hour outputs
51	INTEGER :: jpmld_type ! mixed layer type
52	REAL(wp) :: ppz_ref ! depth of initial T_ref
53	REAL(wp) :: ppdT_crit ! Critical temp diff
54	REAL(wp) :: ppiso_frac ! Fraction of ppdT_crit used
55
56	!Used for 25h mean
57	LOGICAL, PRIVATE :: kara_25h_init = .TRUE. !Logical used to initalise 25h
58	!outputs. Necissary, because we need to
59	!initalise the kara_25h on the zeroth
60	!timestep (i.e in the nemogcm_init call)
61	REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h
62
63
64
65	!JT
66
67	REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth
68	REAL(wp), PUBLIC :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth
69
70	TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs
71	INTEGER :: mld_type ! mixed layer type
72	REAL(wp) :: zref ! depth of initial T_ref
73	REAL(wp) :: dT_crit ! Critical temp diff
74	REAL(wp) :: iso_frac ! Fraction of rn_dT_crit
75	END TYPE MXL_ZINT
76
77	!!----------------------------------------------------------------------
78	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
79	!! $Id$
80	!! Software governed by the CeCILL license (see ./LICENSE)
81	!!----------------------------------------------------------------------
82	CONTAINS
83
84	INTEGER FUNCTION zdf_mxl_alloc()
85	!!----------------------------------------------------------------------
86	!! * FUNCTION zdf_mxl_alloc *
87	!!----------------------------------------------------------------------
88	zdf_mxl_alloc = 0 ! set to zero if no array to be allocated
89	IF( .NOT. ALLOCATED( nmln ) ) THEN
90	ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), &
91	& htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
92	!
93	CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc )
94	IF( zdf_mxl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_alloc: failed to allocate arrays.' )
95	!
96	ENDIF
97	END FUNCTION zdf_mxl_alloc
98
99
100	SUBROUTINE zdf_mxl( kt )
101	!!----------------------------------------------------------------------
102	!! * ROUTINE zdfmxl *
103	!!
104	!! ** Purpose : Compute the turbocline depth and the mixed layer depth
105	!! with density criteria.
106	!!
107	!! ** Method : The mixed layer depth is the shallowest W depth with
108	!! the density of the corresponding T point (just bellow) bellow a
109	!! given value defined locally as rho(10m) + rho_c
110	!! The turbocline depth is the depth at which the vertical
111	!! eddy diffusivity coefficient (resulting from the vertical physics
112	!! alone, not the isopycnal part, see trazdf.F) fall below a given
113	!! value defined locally (avt_c here taken equal to 5 cm/s2 by default)
114	!!
115	!! ** Action : nmln, hmld, hmlp, hmlpt
116	!!----------------------------------------------------------------------
117	INTEGER, INTENT(in) :: kt ! ocean time-step index
118	!
119	INTEGER :: ji, jj, jk ! dummy loop indices
120	INTEGER :: iikn, iiki, ikt ! local integer
121	REAL(wp) :: zN2_c ! local scalar
122	INTEGER, DIMENSION(jpi,jpj) :: imld ! 2D workspace
123	!!----------------------------------------------------------------------
124	!
125	IF( kt == nit000 ) THEN
126	IF(lwp) WRITE(numout,*)
127	IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
128	IF(lwp) WRITE(numout,*) '~~~~~~~ '
129	! ! allocate zdfmxl arrays
130	IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
131	ENDIF
132	!
133	! w-level of the mixing and mixed layers
134	nmln(:,:) = nlb10 ! Initialization to the number of w ocean point
135	hmlp(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2
136	zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
137	DO jk = nlb10, jpkm1
138	DO jj = 1, jpj ! Mixed layer level: w-level
139	DO ji = 1, jpi
140	ikt = mbkt(ji,jj)
141	hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
142	IF( hmlp(ji,jj) < zN2_c ) nmln(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
143	END DO
144	END DO
145	END DO
146	!
147	! w-level of the turbocline and mixing layer (iom_use)
148	imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point
149	DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10
150	DO jj = 1, jpj
151	DO ji = 1, jpi
152	IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline
153	END DO
154	END DO
155	END DO
156	! depth of the mixing and mixed layers
157	!JT
158	CALL zdf_mxl_kara( kt ) ! kara MLD
159	!JT
160
161	DO jj = 1, jpj
162	DO ji = 1, jpi
163	iiki = imld(ji,jj)
164	iikn = nmln(ji,jj)
165	hmld (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth
166	hmlp (ji,jj) = gdepw_n(ji,jj,iikn ) * ssmask(ji,jj) ! Mixed layer depth
167	hmlpt(ji,jj) = gdept_n(ji,jj,iikn-1) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer
168	END DO
169	END DO
170	!
171	IF( .NOT.l_offline ) THEN
172	IF( iom_use("mldr10_1") ) THEN
173	IF( ln_isfcav ) THEN ; CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness
174	ELSE ; CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth
175	END IF
176	END IF
177	IF( iom_use("mldkz5") ) THEN
178	IF( ln_isfcav ) THEN ; CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness
179	ELSE ; CALL iom_put( "mldkz5" , hmld ) ! turbocline depth
180	END IF
181	ENDIF
182	ENDIF
183	!
184	! Vertically-interpolated mixed-layer depth diagnostic
185	CALL zdf_mxl_zint( kt )
186	!
187	IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ' )
188	!
189	END SUBROUTINE zdf_mxl
190
191	SUBROUTINE zdf_mxl_zint_mld( sf )
192	!!----------------------------------------------------------------------------------
193	!! * ROUTINE zdf_mxl_zint_mld *
194	!
195	! Calculate vertically-interpolated mixed layer depth diagnostic.
196	!
197	! This routine can calculate the mixed layer depth diagnostic suggested by
198	! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
199	! vertically-interpolated mixed-layer depth diagnostics with other parameter
200	! settings set in the namzdf_mldzint namelist.
201	!
202	! If mld_type=1 the mixed layer depth is calculated as the depth at which the
203	! density has increased by an amount equivalent to a temperature difference of
204	! 0.8C at the surface.
205	!
206	! For other values of mld_type the mixed layer is calculated as the depth at
207	! which the temperature differs by 0.8C from the surface temperature.
208	!
209	! David Acreman, Daley Calvert
210	!
211	!!-----------------------------------------------------------------------------------
212
213	TYPE(MXL_ZINT), INTENT(in) :: sf
214
215	! Diagnostic criteria
216	INTEGER :: nn_mld_type ! mixed layer type
217	REAL(wp) :: rn_zref ! depth of initial T_ref
218	REAL(wp) :: rn_dT_crit ! Critical temp diff
219	REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used
220
221	! Local variables
222	REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value
223	INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels
224	INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level
225	INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level
226	REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density
227	REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho)
228	REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature
229	REAL :: zT_b ! base temperature
230	REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT
231	REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference
232	REAL :: zdz ! depth difference
233	REAL :: zdT ! temperature difference
234	REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon
235	REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities
236	INTEGER :: ji, jj, jk ! loop counter
237
238	!!-------------------------------------------------------------------------------------
239	!
240	! Unpack structure
241	nn_mld_type = sf%mld_type
242	rn_zref = sf%zref
243	rn_dT_crit = sf%dT_crit
244	rn_iso_frac = sf%iso_frac
245
246	! Set the mixed layer depth criterion at each grid point
247	IF( nn_mld_type == 0 ) THEN
248	zdelta_T(:,:) = rn_dT_crit
249	zT(:,:,:) = rhop(:,:,:)
250	ELSE IF( nn_mld_type == 1 ) THEN
251	ppzdep(:,:)=0.0
252	call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )
253	! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
254	! [assumes number of tracers less than number of vertical levels]
255	zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)
256	zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit
257	CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )
258	zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0
259	! RHO from eos (2d version) doesn't calculate north or east halo:
260	CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )
261	zT(:,:,:) = rhop(:,:,:)
262	ELSE
263	zdelta_T(:,:) = rn_dT_crit
264	zT(:,:,:) = tsn(:,:,:,jp_tem)
265	END IF
266
267	! Calculate the gradient of zT and absolute difference for use later
268	DO jk = 1 ,jpk-2
269	zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)
270	zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
271	END DO
272
273	! Find density/temperature at the reference level (Kara et al use 10m).
274	! ik_ref is the index of the box centre immediately above or at the reference level
275	! Find rn_zref in the array of model level depths and find the ref
276	! density/temperature by linear interpolation.
277	DO jk = jpkm1, 2, -1
278	WHERE ( gdept_n(:,:,jk) > rn_zref )
279	ik_ref(:,:) = jk - 1
280	zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )
281	END WHERE
282	END DO
283
284	! If the first grid box centre is below the reference level then use the
285	! top model level to get zT_ref
286	WHERE ( gdept_n(:,:,1) > rn_zref )
287	zT_ref = zT(:,:,1)
288	ik_ref = 1
289	END WHERE
290
291	! The number of active tracer levels is 1 less than the number of active w levels
292	ikmt(:,:) = mbkt(:,:) - 1
293
294	! Initialize / reset
295	ll_found(:,:) = .false.
296
297	IF ( rn_iso_frac - zepsilon > 0. ) THEN
298	! Search for a uniform density/temperature region where adjacent levels
299	! differ by less than rn_iso_frac * deltaT.
300	! ik_iso is the index of the last level in the uniform layer
301	! ll_found indicates whether the mixed layer depth can be found by interpolation
302	ik_iso(:,:) = ik_ref(:,:)
303	DO jj = 1, nlcj
304	DO ji = 1, nlci
305	!CDIR NOVECTOR
306	DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1
307	IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
308	ik_iso(ji,jj) = jk
309	ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
310	EXIT
311	END IF
312	END DO
313	END DO
314	END DO
315
316	! Use linear interpolation to find depth of mixed layer base where possible
317	hmld_zint(:,:) = rn_zref
318	DO jj = 1, jpj
319	DO ji = 1, jpi
320	IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
321	zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
322	hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz
323	END IF
324	END DO
325	END DO
326	END IF
327
328	! If ll_found = .false. then calculate MLD using difference of zdelta_T
329	! from the reference density/temperature
330
331	! Prevent this section from working on land points
332	WHERE ( tmask(:,:,1) /= 1.0 )
333	ll_found = .true.
334	END WHERE
335
336	DO jk=1, jpk
337	ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)
338	END DO
339
340	! Set default value where interpolation cannot be used (ll_found=false)
341	DO jj = 1, jpj
342	DO ji = 1, jpi
343	IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj))
344	END DO
345	END DO
346
347	DO jj = 1, jpj
348	DO ji = 1, jpi
349	!CDIR NOVECTOR
350	DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)
351	IF ( ll_found(ji,jj) ) EXIT
352	IF ( ll_belowml(ji,jj,jk) ) THEN
353	zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )
354	zdT = zT_b - zT(ji,jj,jk-1)
355	zdz = zdT / zdTdz(ji,jj,jk-1)
356	hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz
357	EXIT
358	END IF
359	END DO
360	END DO
361	END DO
362
363	hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)
364	!
365	END SUBROUTINE zdf_mxl_zint_mld
366
367	SUBROUTINE zdf_mxl_zint_htc( kt )
368	!!----------------------------------------------------------------------
369	!! * ROUTINE zdf_mxl_zint_htc *
370	!!
371	!! ** Purpose :
372	!!
373	!! ** Method :
374	!!----------------------------------------------------------------------
375
376	INTEGER, INTENT(in) :: kt ! ocean time-step index
377
378	INTEGER :: ji, jj, jk
379	INTEGER :: ikmax
380	REAL(wp) :: zc, zcoef
381	!
382	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel
383	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick
384
385	!!----------------------------------------------------------------------
386
387	IF( .NOT. ALLOCATED(ilevel) ) THEN
388	ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
389	& zthick(jpi,jpj), STAT=ji )
390	IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji )
391	IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
392	ENDIF
393
394	! Find last whole model T level above the MLD
395	ilevel(:,:) = 0
396	zthick_0(:,:) = 0._wp
397
398	DO jk = 1, jpkm1
399	DO jj = 1, jpj
400	DO ji = 1, jpi
401	zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk)
402	IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk
403	END DO
404	END DO
405	WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
406	WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1)
407	END DO
408
409	! Surface boundary condition
410	IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)
411	ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp
412	ENDIF
413
414	! Deepest whole T level above the MLD
415	ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
416
417	! Integration down to last whole model T level
418	DO jk = 1, ikmax
419	DO jj = 1, jpj
420	DO ji = 1, jpi
421	zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1 ) ) ! 0 below ilevel
422	zthick(ji,jj) = zthick(ji,jj) + zc
423	htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
424	END DO
425	END DO
426	END DO
427
428	! Subsequent partial T level
429	zthick(:,:) = hmld_zint(:,:) - zthick(:,:) ! remaining thickness to reach MLD
430
431	DO jj = 1, jpj
432	DO ji = 1, jpi
433	htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) &
434	& * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
435	END DO
436	END DO
437
438	WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
439
440	! Convert to heat content
441	zcoef = rau0 * rcp
442	htc_mld(:,:) = zcoef * htc_mld(:,:)
443
444	END SUBROUTINE zdf_mxl_zint_htc
445
446	SUBROUTINE zdf_mxl_zint( kt )
447	!!----------------------------------------------------------------------
448	!! * ROUTINE zdf_mxl_zint *
449	!!
450	!! ** Purpose :
451	!!
452	!! ** Method :
453	!!----------------------------------------------------------------------
454
455	INTEGER, INTENT(in) :: kt ! ocean time-step index
456
457	INTEGER :: ios
458	INTEGER :: jn
459
460	INTEGER :: nn_mld_diag = 0 ! number of diagnostics
461
462	CHARACTER(len=1) :: cmld
463
464	TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
465	TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags
466
467	NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
468
469	!!----------------------------------------------------------------------
470
471	IF( kt == nit000 ) THEN
472	REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist
473	READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
474	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' )
475
476	REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist
477	READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
478	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' )
479	IF(lwm) WRITE ( numond, namzdf_mldzint )
480
481	IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
482
483	mld_diags(1) = sn_mld1
484	mld_diags(2) = sn_mld2
485	mld_diags(3) = sn_mld3
486	mld_diags(4) = sn_mld4
487	mld_diags(5) = sn_mld5
488
489	IF( lwp .AND. (nn_mld_diag > 0) ) THEN
490	WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
491	WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
492	DO jn = 1, nn_mld_diag
493	WRITE(numout,*) 'MLD criterion',jn,':'
494	WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type
495	WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref
496	WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit
497	WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac
498	END DO
499	WRITE(numout,*) '===================================================================='
500	ENDIF
501	ENDIF
502
503	IF( nn_mld_diag > 0 ) THEN
504	DO jn = 1, nn_mld_diag
505	WRITE(cmld,'(I1)') jn
506	IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN
507	CALL zdf_mxl_zint_mld( mld_diags(jn) )
508
509	IF( iom_use( "mldzint_"//cmld ) ) THEN
510	CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
511	ENDIF
512
513	IF( iom_use( "mldhtc_"//cmld ) ) THEN
514	CALL zdf_mxl_zint_htc( kt )
515	CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) )
516	ENDIF
517	ENDIF
518	END DO
519	ENDIF
520
521	END SUBROUTINE zdf_mxl_zint
522
523
524
525	SUBROUTINE zdf_mxl_kara( kt )
526	!!----------------------------------------------------------------------------------
527	!! * ROUTINE zdf_mxl_kara *
528	!
529	! Calculate mixed layer depth according to the definition of
530	! Kara et al, 2000, JGR, 105, 16803.
531	!
532	! If mld_type=1 the mixed layer depth is calculated as the depth at which the
533	! density has increased by an amount equivalent to a temperature difference of
534	! 0.8C at the surface.
535	!
536	! For other values of mld_type the mixed layer is calculated as the depth at
537	! which the temperature differs by 0.8C from the surface temperature.
538	!
539	! Original version: David Acreman
540	!
541	!!-----------------------------------------------------------------------------------
542
543	INTEGER, INTENT( in ) :: kt ! ocean time-step index
544
545	NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, &
546	& ppiso_frac, ln_kara_write25h
547
548	! Local variables
549	REAL, DIMENSION(jpi,jpk) :: ppzdep ! depth for use in calculating d(rho)
550	REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d !Local version of tsn
551	LOGICAL :: ll_found(jpi,jpj) ! Is T_b to be found by interpolation ?
552	LOGICAL :: ll_belowml(jpi,jpj,jpk) ! Flag points below mixed layer when ll_found=F
553	INTEGER :: ji, jj, jk ! loop counter
554	INTEGER :: ik_ref(jpi,jpj) ! index of reference level
555	INTEGER :: ik_iso(jpi,jpj) ! index of last uniform temp level
556	REAL :: zT(jpi,jpj,jpk) ! Temperature or denisty
557	REAL :: zT_ref(jpi,jpj) ! reference temperature
558	REAL :: zT_b ! base temperature
559	REAL :: zdTdz(jpi,jpj,jpk-2) ! gradient of zT
560	REAL :: zmoddT(jpi,jpj,jpk-2) ! Absolute temperature difference
561	REAL :: zdz ! depth difference
562	REAL :: zdT ! temperature difference
563	REAL :: zdelta_T(jpi,jpj) ! difference critereon
564	REAL :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities
565	INTEGER, SAVE :: idt, i_steps
566	INTEGER, SAVE :: i_cnt_25h ! Counter for 25 hour means
567	INTEGER :: ios ! Local integer output status for namelist read
568
569
570
571	!!-------------------------------------------------------------------------------------
572
573	IF( kt == nit000 ) THEN
574	! Default values
575	ln_kara = .FALSE.
576	ln_kara_write25h = .FALSE.
577	jpmld_type = 1
578	ppz_ref = 10.0
579	ppdT_crit = 0.2
580	ppiso_frac = 0.1
581	! Read namelist
582	REWIND ( numnam_ref )
583	READ ( numnam_ref, namzdf_karaml, IOSTAT=ios, ERR= 901 )
584	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in reference namelist' )
585
586	REWIND( numnam_cfg ) ! Namelist nam_diadiaop in configuration namelist 3D hourly diagnostics
587	READ ( numnam_cfg, namzdf_karaml, IOSTAT = ios, ERR = 902 )
588	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist' )
589	IF(lwm) WRITE ( numond, namzdf_karaml )
590
591
592	WRITE(numout,*) '===== Kara mixed layer ====='
593	WRITE(numout,*) 'ln_kara = ', ln_kara
594	WRITE(numout,*) 'jpmld_type = ', jpmld_type
595	WRITE(numout,*) 'ppz_ref = ', ppz_ref
596	WRITE(numout,*) 'ppdT_crit = ', ppdT_crit
597	WRITE(numout,*) 'ppiso_frac = ', ppiso_frac
598	WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h
599	WRITE(numout,*) '============================'
600
601	IF ( .NOT.ln_kara ) THEN
602	WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated"
603	ELSE
604	IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) )
605	IF ( ln_kara_write25h .AND. kara_25h_init ) THEN
606	i_cnt_25h = 0
607	IF (.NOT. ALLOCATED(hmld_kara_25h) ) &
608	& ALLOCATE( hmld_kara_25h(jpi,jpj) )
609	hmld_kara_25h = 0._wp
610	!IF( nacc == 1 ) THEN
611	! idt = INT(rdtmin)
612	!ELSE
613	! idt = INT(rdt)
614	!ENDIF
615
616	idt = INT(rdt)
617	IF( MOD( 3600,idt) == 0 ) THEN
618	i_steps = 3600 / idt
619	ELSE
620	CALL ctl_stop('STOP', &
621	& 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// &
622	& ' = 0 otherwise no hourly values are possible')
623	ENDIF
624	ENDIF
625	ENDIF
626	ENDIF
627
628	IF ( ln_kara ) THEN
629
630	!set critical ln_kara
631	ztsn_2d = tsn(:,:,1,:)
632	ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit
633
634	! Set the mixed layer depth criterion at each grid point
635	ppzdep = 0._wp
636	IF( jpmld_type == 1 ) THEN
637	CALL eos ( tsn(:,:,1,:), &
638	& ppzdep(:,:), zRHO1(:,:) )
639	CALL eos ( ztsn_2d(:,:,:), &
640	& ppzdep(:,:), zRHO2(:,:) )
641	zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0
642	! RHO from eos (2d version) doesn't calculate north or east halo:
643	CALL lbc_lnk( 'zdf_mxl_kara',zdelta_T, 'T', 1. )
644	zT(:,:,:) = rhop(:,:,:)
645	ELSE
646	zdelta_T(:,:) = ppdT_crit
647	zT(:,:,:) = tsn(:,:,:,jp_tem)
648	ENDIF
649
650	! Calculate the gradient of zT and absolute difference for use later
651	DO jk = 1 ,jpk-2
652	zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)
653	zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
654	END DO
655
656	! Find density/temperature at the reference level (Kara et al use 10m).
657	! ik_ref is the index of the box centre immediately above or at the reference level
658	! Find ppz_ref in the array of model level depths and find the ref
659	! density/temperature by linear interpolation.
660	ik_ref = -1
661	DO jk = jpkm1, 2, -1
662	WHERE( gdept_n(:,:,jk) > ppz_ref )
663	ik_ref(:,:) = jk - 1
664	zT_ref(:,:) = zT(:,:,jk-1) + &
665	& zdTdz(:,:,jk-1) * ( ppz_ref - gdept_n(:,:,jk-1) )
666	ENDWHERE
667	END DO
668	IF ( ANY( ik_ref < 0 ) .OR. ANY( ik_ref > jpkm1 ) ) THEN
669	CALL ctl_stop( "STOP", &
670	& "zdf_mxl_kara: unable to find reference level for kara ML" )
671	ELSE
672	! If the first grid box centre is below the reference level then use the
673	! top model level to get zT_ref
674	WHERE( gdept_n(:,:,1) > ppz_ref )
675	zT_ref = zT(:,:,1)
676	ik_ref = 1
677	ENDWHERE
678
679	! Search for a uniform density/temperature region where adjacent levels
680	! differ by less than ppiso_frac * deltaT.
681	! ik_iso is the index of the last level in the uniform layer
682	! ll_found indicates whether the mixed layer depth can be found by interpolation
683	ik_iso(:,:) = ik_ref(:,:)
684	ll_found(:,:) = .false.
685	DO jj = 1, nlcj
686	DO ji = 1, nlci
687	!CDIR NOVECTOR
688	DO jk = ik_ref(ji,jj), mbkt(ji,jj)-1 !mbathy(ji,jj)-1
689	IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN
690	ik_iso(ji,jj) = jk
691	ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
692	EXIT
693	ENDIF
694	END DO
695	END DO
696	END DO
697
698	! Use linear interpolation to find depth of mixed layer base where possible
699	hmld_kara(:,:) = ppz_ref
700	DO jj = 1, jpj
701	DO ji = 1, jpi
702	IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN
703	zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
704	hmld_kara(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz
705	ENDIF
706	END DO
707	END DO
708
709	! If ll_found = .false. then calculate MLD using difference of zdelta_T
710	! from the reference density/temperature
711
712	! Prevent this section from working on land points
713	WHERE( tmask(:,:,1) /= 1.0 )
714	ll_found = .true.
715	ENDWHERE
716
717	DO jk = 1, jpk
718	ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= &
719	& zdelta_T(:,:)
720	END DO
721
722	! Set default value where interpolation cannot be used (ll_found=false)
723	DO jj = 1, jpj
724	DO ji = 1, jpi
725	IF( .NOT. ll_found(ji,jj) ) &
726	& hmld_kara(ji,jj) = gdept_n(ji,jj,mbkt(ji,jj))! mbathy(ji,jj))
727	END DO
728	END DO
729
730	DO jj = 1, jpj
731	DO ji = 1, jpi
732	!CDIR NOVECTOR
733	DO jk = ik_ref(ji,jj)+1, mbkt(ji,jj) !mbathy(ji,jj)
734	IF( ll_found(ji,jj) ) EXIT
735	IF( ll_belowml(ji,jj,jk) ) THEN
736	zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * &
737	& SIGN(1.0, zdTdz(ji,jj,jk-1) )
738	zdT = zT_b - zT(ji,jj,jk-1)
739	zdz = zdT / zdTdz(ji,jj,jk-1)
740	hmld_kara(ji,jj) = gdept_n(ji,jj,jk-1) + zdz
741	EXIT
742	ENDIF
743	END DO
744	END DO
745	END DO
746
747	hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1)
748
749	IF( ln_kara_write25h ) THEN
750	!Double IF required as i_steps not defined if ln_kara_write25h =
751	! FALSE
752	IF ( ( MOD( kt, i_steps ) == 0 ) .OR. kara_25h_init ) THEN
753	hmld_kara_25h = hmld_kara_25h + hmld_kara
754	i_cnt_25h = i_cnt_25h + 1
755	IF ( kara_25h_init ) kara_25h_init = .FALSE.
756	ENDIF
757	ENDIF
758
759	!#if defined key_iomput
760	IF( kt /= nit000 ) THEN
761	CALL iom_put( "mldkara" , hmld_kara )
762	IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND. ln_kara_write25h ) &
763	CALL iom_put( "kara25h" , ( hmld_kara_25h / 25._wp ) )
764	ENDIF
765	!#endif
766
767	ENDIF
768	ENDIF
769
770	END SUBROUTINE zdf_mxl_kara
771
772	!!======================================================================
773	END MODULE zdfmxl

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

Download in other formats: