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 @ 9181

Last change on this file since 9181 was 9181, checked in by kingr, 6 years ago
Adding changes required for surface temperature increments to be projected through mixed layer.
File size: 26.8 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_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	!!----------------------------------------------------------------------
70	CONTAINS
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)
517	901 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 )
521	902 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	ENDIF
600	ENDIF
601
602	ENDIF
603
604	END SUBROUTINE zdf_mxl_zint
605
606	!!======================================================================
607	END MODULE zdfmxl

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

Download in other formats: