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

Context Navigation

source: branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 6443

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

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

Download in other formats: