Context Navigation

zdfmxl.F90 @ 6491

Last change on this file since 6491 was 6491, checked in by davestorkey, 8 years ago

Commit remaining changes to UKMO/r5518_GO6_package branch from following branches:
UKMO/dev_r5021_nn_etau_revision@6238
UKMO/dev_r5107_mld_zint@5534
UKMO/dev_r5107_eorca025_closea@6390
UKMO/restart_datestamp@5539
UKMO/icebergs_latent_heat@5821
UKMO/icebergs_restart_single_file_corrected@6480
UKMO/product_diagnostics@5971
UKMO/antarctic_partial_slip@5961
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5961 cf. r5958 of /branches/UKMO/antarctic_partial_slip/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6349 cf. r5962 of /branches/UKMO/product_diagnostics/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6480 cf. r6479 of /branches/UKMO/icebergs_restart_single_file_corrected/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5986 cf. r5852 of /branches/UKMO/icebergs_restart_single_file/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5821 cf. r5808 of /branches/UKMO/icebergs_latent_heat/NEMOGCM@6490

File size: 16.8 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, imkt ! 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	imkt = mikt(ji,jj)
131	IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = MAX( imkt, jk ) ! Turbocline
132	END DO
133	END DO
134	END DO
135	! depth of the mixing and mixed layers
136	DO jj = 1, jpj
137	DO ji = 1, jpi
138	iiki = imld(ji,jj)
139	iikn = nmln(ji,jj)
140	imkt = mikt(ji,jj)
141	hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! Turbocline depth
142	hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! Mixed layer depth
143	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
144	END DO
145	END DO
146	IF( .NOT.lk_offline ) THEN ! no need to output in offline mode
147	IF( kt >= nit000 ) THEN ! workaround for calls before SOMETHING reads the XIOS namelist
148	CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth
149	CALL iom_put( "mldkz5" , hmld ) ! turbocline depth
150	ENDIF
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) ) / fse3w(:,:,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 ( fsdept(:,:,jk) > rn_zref )
264	ik_ref(:,:) = jk - 1
265	zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,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 ( fsdept(:,:,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) = fsdept(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) = fsdept(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) = fsdept(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	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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 @ 6491

Download in other formats: