Context Navigation

zdfmxl.F90 @ 11134

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

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

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

Context Navigation

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 11134

Download in other formats: