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

Context Navigation

source: branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 5249

Last change on this file since 5249 was 5249, checked in by davestorkey, 9 years ago
UKMO mld_zint branch: implement vertically-interpolated MLD.
File size: 16.7 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, SAVE, DIMENSION(:,:) :: hmld_zint !: vertically-interpolated mixed layer depth [m]
36	LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ?
37	LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F
38
39	REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth
40	REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth
41
42	! Namelist variables for namzdf_mldzint
43	INTEGER :: nn_mld_type ! mixed layer type
44	REAL(wp) :: rn_zref ! depth of initial T_ref
45	REAL(wp) :: rn_dT_crit ! Critical temp diff
46	REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used
47
48	!! * Substitutions
49	# include "domzgr_substitute.h90"
50	!!----------------------------------------------------------------------
51	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
52	!! $Id$
53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54	!!----------------------------------------------------------------------
55	CONTAINS
56
57	INTEGER FUNCTION zdf_mxl_alloc()
58	!!----------------------------------------------------------------------
59	!! * FUNCTION zdf_mxl_alloc *
60	!!----------------------------------------------------------------------
61	zdf_mxl_alloc = 0 ! set to zero if no array to be allocated
62	IF( .NOT. ALLOCATED( nmln ) ) THEN
63	ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), &
64	& ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
65	!
66	IF( lk_mpp ) CALL mpp_sum ( zdf_mxl_alloc )
67	IF( zdf_mxl_alloc /= 0 ) CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.')
68	!
69	ENDIF
70	END FUNCTION zdf_mxl_alloc
71
72
73	SUBROUTINE zdf_mxl( kt )
74	!!----------------------------------------------------------------------
75	!! * ROUTINE zdfmxl *
76	!!
77	!! ** Purpose : Compute the turbocline depth and the mixed layer depth
78	!! with density criteria.
79	!!
80	!! ** Method : The mixed layer depth is the shallowest W depth with
81	!! the density of the corresponding T point (just bellow) bellow a
82	!! given value defined locally as rho(10m) + rho_c
83	!! The turbocline depth is the depth at which the vertical
84	!! eddy diffusivity coefficient (resulting from the vertical physics
85	!! alone, not the isopycnal part, see trazdf.F) fall below a given
86	!! value defined locally (avt_c here taken equal to 5 cm/s2 by default)
87	!!
88	!! ** Action : nmln, hmld, hmlp, hmlpt
89	!!----------------------------------------------------------------------
90	INTEGER, INTENT(in) :: kt ! ocean time-step index
91	!
92	INTEGER :: ji, jj, jk ! dummy loop indices
93	INTEGER :: iikn, iiki, ikt, imkt ! local integer
94	REAL(wp) :: zN2_c ! local scalar
95	INTEGER, POINTER, DIMENSION(:,:) :: imld ! 2D workspace
96	!!----------------------------------------------------------------------
97	!
98	IF( nn_timing == 1 ) CALL timing_start('zdf_mxl')
99	!
100	CALL wrk_alloc( jpi,jpj, imld )
101
102	IF( kt == nit000 ) THEN
103	IF(lwp) WRITE(numout,*)
104	IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
105	IF(lwp) WRITE(numout,*) '~~~~~~~ '
106	! ! allocate zdfmxl arrays
107	IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
108	ENDIF
109
110	! w-level of the mixing and mixed layers
111	nmln(:,:) = nlb10 ! Initialization to the number of w ocean point
112	hmlp(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2
113	zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
114	DO jk = nlb10, jpkm1
115	DO jj = 1, jpj ! Mixed layer level: w-level
116	DO ji = 1, jpi
117	ikt = mbkt(ji,jj)
118	hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk)
119	IF( hmlp(ji,jj) < zN2_c ) nmln(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
120	END DO
121	END DO
122	END DO
123	!
124	! w-level of the turbocline
125	imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point
126	DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10
127	DO jj = 1, jpj
128	DO ji = 1, jpi
129	imkt = mikt(ji,jj)
130	IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = MAX( imkt, jk ) ! Turbocline
131	END DO
132	END DO
133	END DO
134	! depth of the mixing and mixed layers
135	DO jj = 1, jpj
136	DO ji = 1, jpi
137	iiki = imld(ji,jj)
138	iikn = nmln(ji,jj)
139	imkt = mikt(ji,jj)
140	hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! Turbocline depth
141	hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj) ! Mixed layer depth
142	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
143	END DO
144	END DO
145	IF( .NOT.lk_offline ) THEN ! no need to output in offline mode
146	CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth
147	CALL iom_put( "mldkz5" , hmld ) ! turbocline depth
148	ENDIF
149
150	! Vertically-interpolated mixed-layer depth diagnostic
151	IF( iom_use( "mldzint" ) ) THEN
152	CALL zdf_mxl_zint( kt )
153	CALL iom_put( "mldzint" , hmld_zint )
154	ENDIF
155
156	IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
157	!
158	CALL wrk_dealloc( jpi,jpj, imld )
159	!
160	IF( nn_timing == 1 ) CALL timing_stop('zdf_mxl')
161	!
162	END SUBROUTINE zdf_mxl
163
164	SUBROUTINE zdf_mxl_zint( kt )
165	!!----------------------------------------------------------------------------------
166	!! * ROUTINE zdf_mxl_zint *
167	!
168	! Calculate vertically-interpolated mixed layer depth diagnostic.
169	!
170	! This routine can calculate the mixed layer depth diagnostic suggested by
171	! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
172	! vertically-interpolated mixed-layer depth diagnostics with other parameter
173	! settings set in the namzdf_mldzint namelist.
174	!
175	! If mld_type=1 the mixed layer depth is calculated as the depth at which the
176	! density has increased by an amount equivalent to a temperature difference of
177	! 0.8C at the surface.
178	!
179	! For other values of mld_type the mixed layer is calculated as the depth at
180	! which the temperature differs by 0.8C from the surface temperature.
181	!
182	! David Acreman, Daley Calvert
183	!
184	!!-----------------------------------------------------------------------------------
185
186	INTEGER, INTENT(in) :: kt ! ocean time-step index
187	!
188	! Local variables
189	INTEGER, POINTER, DIMENSION(:,:) :: ikmt ! number of active tracer levels
190	INTEGER, POINTER, DIMENSION(:,:) :: ik_ref ! index of reference level
191	INTEGER, POINTER, DIMENSION(:,:) :: ik_iso ! index of last uniform temp level
192	REAL, POINTER, DIMENSION(:,:,:) :: zT ! Temperature or density
193	REAL, POINTER, DIMENSION(:,:) :: ppzdep ! depth for use in calculating d(rho)
194	REAL, POINTER, DIMENSION(:,:) :: zT_ref ! reference temperature
195	REAL :: zT_b ! base temperature
196	REAL, POINTER, DIMENSION(:,:,:) :: zdTdz ! gradient of zT
197	REAL, POINTER, DIMENSION(:,:,:) :: zmoddT ! Absolute temperature difference
198	REAL :: zdz ! depth difference
199	REAL :: zdT ! temperature difference
200	REAL, POINTER, DIMENSION(:,:) :: zdelta_T ! difference critereon
201	REAL, POINTER, DIMENSION(:,:) :: zRHO1, zRHO2 ! Densities
202	INTEGER :: ji, jj, jk ! loop counter
203	INTEGER :: ios
204
205	NAMELIST/namzdf_mldzint/ nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac
206
207	!!-------------------------------------------------------------------------------------
208	!
209	CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso)
210	CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )
211	CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT )
212
213	IF( kt == nit000 ) THEN
214	REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist
215	READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
216	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp )
217
218	REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist
219	READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
220	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp )
221	IF(lwm) WRITE ( numond, namzdf_mldzint )
222
223	WRITE(numout,*) '===== Vertically-interpolated mixed layer ====='
224	WRITE(numout,*) 'nn_mld_type = ',nn_mld_type
225	WRITE(numout,*) 'rn_zref = ',rn_zref
226	WRITE(numout,*) 'rn_dT_crit = ',rn_dT_crit
227	WRITE(numout,*) 'rn_iso_frac = ',rn_iso_frac
228	WRITE(numout,*) '==============================================='
229	ENDIF
230
231	! Set the mixed layer depth criterion at each grid point
232	IF (nn_mld_type == 1) THEN
233	ppzdep(:,:)=0.0
234	call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )
235	! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
236	! [assumes number of tracers less than number of vertical levels]
237	zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)
238	zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit
239	CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )
240	zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0
241	! RHO from eos (2d version) doesn't calculate north or east halo:
242	CALL lbc_lnk( zdelta_T, 'T', 1. )
243	zT(:,:,:) = rhop(:,:,:)
244	ELSE
245	zdelta_T(:,:) = rn_dT_crit
246	zT(:,:,:) = tsn(:,:,:,jp_tem)
247	END IF
248
249	! Calculate the gradient of zT and absolute difference for use later
250	DO jk = 1 ,jpk-2
251	zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)
252	zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
253	END DO
254
255	! Find density/temperature at the reference level (Kara et al use 10m).
256	! ik_ref is the index of the box centre immediately above or at the reference level
257	! Find rn_zref in the array of model level depths and find the ref
258	! density/temperature by linear interpolation.
259	DO jk = jpkm1, 2, -1
260	WHERE ( fsdept(:,:,jk) > rn_zref )
261	ik_ref(:,:) = jk - 1
262	zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) )
263	END WHERE
264	END DO
265
266	! If the first grid box centre is below the reference level then use the
267	! top model level to get zT_ref
268	WHERE ( fsdept(:,:,1) > rn_zref )
269	zT_ref = zT(:,:,1)
270	ik_ref = 1
271	END WHERE
272
273	! The number of active tracer levels is 1 less than the number of active w levels
274	ikmt(:,:) = mbathy(:,:) - 1
275
276	! Search for a uniform density/temperature region where adjacent levels
277	! differ by less than rn_iso_frac * deltaT.
278	! ik_iso is the index of the last level in the uniform layer
279	! ll_found indicates whether the mixed layer depth can be found by interpolation
280	ik_iso(:,:) = ik_ref(:,:)
281	ll_found(:,:) = .false.
282	DO jj = 1, nlcj
283	DO ji = 1, nlci
284	!CDIR NOVECTOR
285	DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1
286	IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
287	ik_iso(ji,jj) = jk
288	ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
289	EXIT
290	END IF
291	END DO
292	END DO
293	END DO
294
295	! Use linear interpolation to find depth of mixed layer base where possible
296	hmld_zint(:,:) = rn_zref
297	DO jj = 1, jpj
298	DO ji = 1, jpi
299	IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
300	zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
301	hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz
302	END IF
303	END DO
304	END DO
305
306	! If ll_found = .false. then calculate MLD using difference of zdelta_T
307	! from the reference density/temperature
308
309	! Prevent this section from working on land points
310	WHERE ( tmask(:,:,1) /= 1.0 )
311	ll_found = .true.
312	END WHERE
313
314	DO jk=1, jpk
315	ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)
316	END DO
317
318	! Set default value where interpolation cannot be used (ll_found=false)
319	DO jj = 1, jpj
320	DO ji = 1, jpi
321	IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj))
322	END DO
323	END DO
324
325	DO jj = 1, jpj
326	DO ji = 1, jpi
327	!CDIR NOVECTOR
328	DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)
329	IF ( ll_found(ji,jj) ) EXIT
330	IF ( ll_belowml(ji,jj,jk) ) THEN
331	zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )
332	zdT = zT_b - zT(ji,jj,jk-1)
333	zdz = zdT / zdTdz(ji,jj,jk-1)
334	hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz
335	EXIT
336	END IF
337	END DO
338	END DO
339	END DO
340
341	hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)
342	!
343	CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso)
344	CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )
345	CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT )
346	!
347	END SUBROUTINE zdf_mxl_zint
348
349	!!======================================================================
350	END MODULE zdfmxl

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

Download in other formats: