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

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

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

Download in other formats: