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

Context Navigation

source: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 6825

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

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

Download in other formats: