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.

ablmod.F90 in NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90 @ 12592

Last change on this file since 12592 was 12592, checked in by gsamson, 4 years ago
correct mxld_abl & mxlm_abl init to avoid division by zero (ticket #2419)
File size: 60.4 KB

Line
1	MODULE ablmod
2	!!======================================================================
3	!! * MODULE ablmod *
4	!! Surface module : ABL computation to provide atmospheric data
5	!! for surface fluxes computation
6	!!======================================================================
7	!! History : 3.6 ! 2019-03 (F. Lemarié & G. Samson) Original code
8	!!----------------------------------------------------------------------
9
10	!!----------------------------------------------------------------------
11	!! abl_stp : ABL single column model
12	!! abl_zdf_tke : atmospheric vertical closure
13	!!----------------------------------------------------------------------
14	USE abl ! ABL dynamics and tracers
15	USE par_abl ! ABL constants
16
17	USE phycst ! physical constants
18	USE dom_oce, ONLY : tmask
19	USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa
20	USE sbcblk ! use rn_?fac
21	USE sbcblk_phy ! use some physical constants for flux computation
22	!
23	USE prtctl ! Print control (prt_ctl routine)
24	USE iom ! IOM library
25	USE in_out_manager ! I/O manager
26	USE lib_mpp ! MPP library
27	USE timing ! Timing
28
29	IMPLICIT NONE
30
31	PUBLIC abl_stp ! called by sbcabl.F90
32	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2, zrough
33	!! * Substitutions
34	# include "do_loop_substitute.h90"
35
36	CONTAINS
37
38
39	!===================================================================================================
40	SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, & ! in
41	& pu_dta, pv_dta, pt_dta, pq_dta, &
42	& pslp_dta, pgu_dta, pgv_dta, &
43	& pcd_du, psen, pevp, & ! in/out
44	& pwndm, ptaui, ptauj, ptaum &
45	#if defined key_si3
46	& , ptm_su, pssu_ice, pssv_ice &
47	& , pssq_ice, pcd_du_ice, psen_ice &
48	& , pevp_ice, pwndm_ice, pfrac_oce &
49	& , ptaui_ice, ptauj_ice &
50	#endif
51	& )
52	!---------------------------------------------------------------------------------------------------
53
54	!!---------------------------------------------------------------------
55	!! * ROUTINE abl_stp *
56	!!
57	!! ** Purpose : Time-integration of the ABL model
58	!!
59	!! ** Method : Compute atmospheric variables : vertical turbulence
60	!! + Coriolis term + newtonian relaxation
61	!!
62	!! ** Action : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
63	!! - Advance tracers to time n+1 (Euler backward scheme)
64	!! - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
65	!! - Advance u,v to time n+1 (Euler backward scheme)
66	!! - Apply newtonian relaxation on the dynamics and the tracers
67	!! - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
68	!!
69	!!----------------------------------------------------------------------
70	INTEGER , INTENT(in ) :: kt ! time step index
71	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psst ! sea-surface temperature [Celsius]
72	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu ! sea-surface u (U-point)
73	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point)
74	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq ! sea-surface humidity
75	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pu_dta ! large-scale windi
76	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pv_dta ! large-scale windj
77	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgu_dta ! large-scale hpgi
78	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgv_dta ! large-scale hpgj
79	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pt_dta ! large-scale pot. temp.
80	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pq_dta ! large-scale humidity
81	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp_dta ! sea-level pressure
82	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du ! Cd x Du (T-point)
83	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du
84	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du
85	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! \|\|uwnd\|\|
86	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux
87	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy
88	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! \|\|tau\|\|
89	!
90	#if defined key_si3
91	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptm_su ! ice-surface temperature [K]
92	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu_ice ! ice-surface u (U-point)
93	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point)
94	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity
95	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du_ice ! Cd x Du over ice (T-point)
96	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psen_ice ! Ch x Du over ice (T-point)
97	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pevp_ice ! Ce x Du over ice (T-point)
98	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! \|\|uwnd - uice\|\|
99	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce ! ocean fraction
100	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui_ice ! ice-surface taux stress (U-point)
101	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point)
102	#endif
103	!
104	REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j
105	REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF
106	!
107	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
108	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
109	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
110	!
111	INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices
112	REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
113	REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
114	LOGICAL :: SemiImp_Cor = .FALSE.
115	!
116	!!---------------------------------------------------------------------
117	!
118	IF(lwp .AND. kt == nit000) THEN ! control print
119	WRITE(numout,*)
120	WRITE(numout,*) 'abl_stp : ABL time stepping'
121	WRITE(numout,*) '~~~~~~'
122	ENDIF
123	!
124	IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
125	IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) )
126	!! Compute ustar squared as Cd \|\| Uatm-Uoce \|\|^2
127	!! needed for surface boundary condition of TKE
128	!! pwndm contains \| U10m - U_oce \| (see blk_oce_1 in sbcblk)
129	DO_2D_11_11
130	zzoce = pCd_du (ji,jj) * pwndm (ji,jj)
131	#if defined key_si3
132	zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)
133	ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice
134	#else
135	ustar2(ji,jj) = zzoce
136	#endif
137	zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT(Cdn_oce(ji,jj)) ) !<-- recover the value of z0 from Cdn_oce
138	END_2D
139	!
140	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
141	! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
142	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
143
144	CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
145
146	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
147	! ! 2 *** Advance tracers to time n+1
148	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
149
150	!-------------
151	DO jj = 1, jpj ! outer loop !--> tq_abl computed on (1:jpi) x (1:jpj)
152	!-------------
153	! Compute matrix elements for interior points
154	DO jk = 3, jpkam1
155	DO ji = 1, jpi ! vector opt.
156	z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
157	z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
158	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
159	END DO
160	END DO
161	! Boundary conditions
162	DO ji = 1, jpi ! vector opt.
163	! Neumann at the bottom
164	z_elem_a( ji, 2 ) = 0._wp
165	z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )
166	! Homogeneous Neumann at the top
167	z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
168	z_elem_c( ji, jpka ) = 0._wp
169	z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
170	END DO
171
172	DO jtra = 1,jptq ! loop on active tracers
173
174	DO jk = 3, jpkam1
175	DO ji = 2, jpim1
176	!DO ji = 1,jpi !!GS: to be checked
177	tq_abl ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl ( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side
178	END DO
179	END DO
180
181	IF(jtra == jp_ta) THEN
182	DO ji = 1,jpi ! surface boundary condition for temperature
183	zztmp1 = psen(ji, jj)
184	zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
185	#if defined key_si3
186	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
187	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
188	#endif
189	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
190	tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2
191	tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
192	END DO
193	ELSE ! jp_qa
194	DO ji = 1,jpi ! surface boundary condition for humidity
195	zztmp1 = pevp(ji, jj)
196	zztmp2 = pevp(ji, jj) * pssq(ji, jj)
197	#if defined key_si3
198	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj)
199	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)
200	#endif
201	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
202	tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2
203	tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
204	END DO
205	END IF
206	!!
207	!! Matrix inversion
208	!! ----------------------------------------------------------
209	DO ji = 1,jpi
210	zcff = 1._wp / z_elem_b( ji, 2 )
211	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
212	tq_abl(ji,jj,2,nt_a,jtra) = zcff * tq_abl(ji,jj,2,nt_a,jtra)
213	END DO
214
215	DO jk = 3, jpka
216	DO ji = 1,jpi
217	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
218	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
219	tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) &
220	& - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
221	END DO
222	END DO
223	!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
224	DO jk = jpkam1,2,-1
225	DO ji = 1,jpi
226	tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) + &
227	& zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
228	END DO
229	END DO
230
231	END DO !<-- loop on tracers
232	!!
233	!-------------
234	END DO ! end outer loop
235	!-------------
236
237	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
238	! ! 3 *** Compute Coriolis term with geostrophic guide
239	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240	IF( SemiImp_Cor ) THEN
241
242	!-------------
243	DO jk = 2, jpka ! outer loop
244	!-------------
245	!
246	! Advance u_abl & v_abl to time n+1
247	DO_2D_11_11
248	zcff = ( fft_abl(ji,jj) * rDt_abl )( fft_abl(ji,jj) rDt_abl ) ! (f dt)**2
249
250	u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
251	& (1._wp-gamma_Cor(1._wp-gamma_Cor)zcff)*u_abl( ji, jj, jk, nt_n ) &
252	& + rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) &
253	& / (1._wp + gamma_Corgamma_Corzcff)
254
255	v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
256	& (1._wp-gamma_Cor(1._wp-gamma_Cor)zcff)*v_abl( ji, jj, jk, nt_n ) &
257	& - rDt_abl * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) &
258	& / (1._wp + gamma_Corgamma_Corzcff)
259	END_2D
260	!
261	!-------------
262	END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
263	!-------------
264	!
265	!IF( ln_geos_winds ) THEN
266	! DO jj = 1, jpj ! outer loop
267	! DO jk = 1, jpka
268	! DO ji = 1, jpi
269	! u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) &
270	! & - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk)
271	! v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) &
272	! & + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk)
273	! END DO
274	! END DO
275	! END DO
276	!END IF
277
278	ELSE
279
280	!-------------
281	DO jk = 2, jpka ! outer loop
282	!-------------
283	!
284	IF( MOD( kt, 2 ) == 0 ) then
285	! Advance u_abl & v_abl to time n+1
286	DO jj = 1, jpj
287	DO ji = 1, jpi
288	zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) )
289	u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
290	zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) )
291	v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) ( v_abl( ji, jj, jk, nt_n ) - rDt_abl zcff )
292	u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a )
293	END DO
294	END DO
295	!
296	ELSE
297	!
298	DO jj = 1, jpj
299	DO ji = 1, jpi
300	zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) )
301	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
302	zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) )
303	u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) ( u_abl( ji, jj, jk, nt_n ) + rDt_abl zcff )
304	v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a )
305	END DO
306	END DO
307	!
308	ENDIF
309	!
310	!-------------
311	END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
312	!-------------
313
314	ENDIF
315
316	!-------------
317	!
318	IF( ln_hpgls_frc ) THEN
319	DO jj = 1, jpj ! outer loop
320	DO jk = 1, jpka
321	DO ji = 1, jpi
322	u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
323	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
324	ENDDO
325	ENDDO
326	ENDDO
327	END IF
328
329	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
330	! ! 4 *** Advance u,v to time n+1
331	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
332	!
333	! Vertical diffusion for u_abl
334	!-------------
335	DO jj = 1, jpj ! outer loop
336	!-------------
337
338	DO jk = 3, jpkam1
339	DO ji = 1, jpi
340	z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
341	z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
342	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
343	END DO
344	END DO
345
346	DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
347	!++ Surface boundary condition
348	z_elem_a( ji, 2 ) = 0._wp
349	z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
350	!
351	zztmp1 = pcd_du(ji, jj)
352	zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac
353	#if defined key_si3
354	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
355	zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) * rn_vfac
356	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
357	#endif
358	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
359	u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
360
361	IF( ln_topbc_neumann ) THEN
362	!++ Top Neumann B.C.
363	z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
364	z_elem_c( ji, jpka ) = 0._wp
365	z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
366	!u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl ( ji, jj, jpka, nt_a )
367	ELSE
368	!++ Top Dirichlet B.C.
369	z_elem_a( ji, jpka ) = 0._wp
370	z_elem_c( ji, jpka ) = 0._wp
371	z_elem_b( ji, jpka ) = e3t_abl( jpka )
372	u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
373	ENDIF
374
375	END DO
376	!!
377	!! Matrix inversion
378	!! ----------------------------------------------------------
379	DO ji = 2, jpi
380	zcff = 1._wp / z_elem_b( ji, 2 )
381	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
382	u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a)
383	END DO
384
385	DO jk = 3, jpka
386	DO ji = 2, jpi
387	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
388	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
389	u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk ,nt_a) &
390	& - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
391	END DO
392	END DO
393
394	DO jk = jpkam1,2,-1
395	DO ji = 2, jpi
396	u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a)
397	END DO
398	END DO
399
400	!-------------
401	END DO ! end outer loop
402	!-------------
403
404	!
405	! Vertical diffusion for v_abl
406	!-------------
407	DO jj = 2, jpj ! outer loop
408	!-------------
409	!
410	DO jk = 3, jpkam1
411	DO ji = 1, jpi
412	z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
413	z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
414	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
415	END DO
416	END DO
417
418	DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
419	!++ Surface boundary condition
420	z_elem_a( ji, 2 ) = 0._wp
421	z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
422	!
423	zztmp1 = pcd_du(ji, jj)
424	zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac
425	#if defined key_si3
426	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
427	zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) ) * rn_vfac
428	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
429	#endif
430	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
431	v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
432
433	IF( ln_topbc_neumann ) THEN
434	!++ Top Neumann B.C.
435	z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
436	z_elem_c( ji, jpka ) = 0._wp
437	z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
438	!v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl ( ji, jj, jpka, nt_a )
439	ELSE
440	!++ Top Dirichlet B.C.
441	z_elem_a( ji, jpka ) = 0._wp
442	z_elem_c( ji, jpka ) = 0._wp
443	z_elem_b( ji, jpka ) = e3t_abl( jpka )
444	v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
445	ENDIF
446
447	END DO
448	!!
449	!! Matrix inversion
450	!! ----------------------------------------------------------
451	DO ji = 1, jpi
452	zcff = 1._wp / z_elem_b( ji, 2 )
453	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
454	v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
455	END DO
456
457	DO jk = 3, jpka
458	DO ji = 1, jpi
459	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
460	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
461	v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk ,nt_a) &
462	& - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
463	END DO
464	END DO
465
466	DO jk = jpkam1,2,-1
467	DO ji = 1, jpi
468	v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a)
469	END DO
470	END DO
471	!
472	!-------------
473	END DO ! end outer loop
474	!-------------
475
476	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
477	! ! 5 *** Apply nudging on the dynamics and the tracers
478	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
479
480	IF( nn_dyn_restore > 0 ) THEN
481	!-------------
482	DO jk = 2, jpka ! outer loop
483	!-------------
484	DO_2D_01_01
485	zcff1 = pblh( ji, jj )
486	zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
487	zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
488	zmsk = msk_abl(ji,jj)
489	zcff2 = jp_alp3_dyn * zsig*3 + jp_alp2_dyn zsig**2 &
490	& + jp_alp1_dyn * zsig + jp_alp0_dyn
491	zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
492	! rn_Dt = rDt_abl / nn_fsbc
493	zcff = zcff * rest_eq(ji,jj)
494	u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) &
495	& + zcff * pu_dta( ji, jj, jk )
496	v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) &
497	& + zcff * pv_dta( ji, jj, jk )
498	END_2D
499	!-------------
500	END DO ! end outer loop
501	!-------------
502	END IF
503
504	!-------------
505	DO jk = 2, jpka ! outer loop
506	!-------------
507	DO_2D_11_11
508	zcff1 = pblh( ji, jj )
509	zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
510	zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
511	zmsk = msk_abl(ji,jj)
512	zcff2 = jp_alp3_tra * zsig*3 + jp_alp2_tra zsig**2 &
513	& + jp_alp1_tra * zsig + jp_alp0_tra
514	zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
515	! rn_Dt = rDt_abl / nn_fsbc
516	tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) &
517	& + zcff * pt_dta( ji, jj, jk )
518
519	tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) &
520	& + zcff * pq_dta( ji, jj, jk )
521
522	END_2D
523	!-------------
524	END DO ! end outer loop
525	!-------------
526	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
527	! ! 6 *** MPI exchanges
528	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
529	!
530	CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a) , 'T', -1., kfillmode = jpfillnothing )
531	CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing ) ! ++++ this should not be needed...
532	!
533	#if defined key_iomput
534	! first ABL level
535	IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) )
536	IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) )
537	IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl", tq_abl(:,:,2,nt_a,jp_ta) )
538	IF ( iom_use("qz1_abl") ) CALL iom_put ( "qz1_abl", tq_abl(:,:,2,nt_a,jp_qa) )
539	IF ( iom_use("uz1_dta") ) CALL iom_put ( "uz1_dta", pu_dta(:,:,2 ) )
540	IF ( iom_use("vz1_dta") ) CALL iom_put ( "vz1_dta", pv_dta(:,:,2 ) )
541	IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta", pt_dta(:,:,2 ) )
542	IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta", pq_dta(:,:,2 ) )
543	! all ABL levels
544	IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) )
545	IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) )
546	IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) )
547	IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) )
548	IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a ) )
549	IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka ) )
550	IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka ) )
551	IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxlm_abl(:,:,2:jpka ) )
552	IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh" , pblh(:,: ) )
553	! debug (to be removed)
554	IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta", pu_dta(:,:,2:jpka) )
555	IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta", pv_dta(:,:,2:jpka) )
556	IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta", pt_dta(:,:,2:jpka) )
557	IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta", pq_dta(:,:,2:jpka) )
558	IF( ln_geos_winds ) THEN
559	IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2 ) )
560	IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2 ) )
561	END IF
562	IF( ln_hpgls_frc ) THEN
563	IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
564	IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
565	END IF
566	#endif
567	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
568	! ! 7 *** Finalize flux computation
569	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
570	!
571	DO_2D_11_11
572	ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta )
573	zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa )
574	zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
575	psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention
576	pevp( ji, jj ) = rn_efacMAX( 0._wp, zcff pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
577	rhoa( ji, jj ) = zcff
578	END_2D
579
580	DO_2D_01_01
581	zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj ) )
582	zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji ,jj-1) )
583	END_2D
584	!
585	CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
586	!
587	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
588	DO_2D_11_11
589	zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) &
590	& + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj)
591	zztmp = rhoa(ji,jj) * pcd_du(ji,jj)
592	pwndm (ji,jj) = zcff
593	ptaum (ji,jj) = zztmp * zcff
594	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
595	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
596	taux_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_i(ji,jj)
597	tauy_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_j(ji,jj)
598	END_2D
599	! ... utau, vtau at U- and V_points, resp.
600	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
601	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
602	DO_2D_00_00
603	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
604	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
605	ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
606	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
607	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
608	ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
609	END_2D
610	!
611	CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
612
613	CALL iom_put( "taum_oce", ptaum )
614
615	IF(sn_cfctl%l_prtctl) THEN
616	CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' )
617	CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ' )
618	CALL prt_ctl( tab2d_2=ptauj , clinfo2= 'vtau : ' )
619	ENDIF
620
621	#if defined key_si3
622	! ------------------------------------------------------------ !
623	! Wind stress relative to the moving ice ( U10m - U_ice ) !
624	! ------------------------------------------------------------ !
625	DO_2D_00_00
626
627	zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) )
628	zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
629
630	ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &
631	& + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &
632	& * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
633	ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &
634	& + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &
635	& * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
636	END_2D
637	CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
638	!
639	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
640	& , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
641	#endif
642	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
643	! ! 8 *** Swap time indices for the next timestep
644	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
645	nt_n = 1 + MOD( kt , 2)
646	nt_a = 1 + MOD( kt+1, 2)
647	!
648	!---------------------------------------------------------------------------------------------------
649	END SUBROUTINE abl_stp
650	!===================================================================================================
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668	!===================================================================================================
669	SUBROUTINE abl_zdf_tke( )
670	!---------------------------------------------------------------------------------------------------
671
672	!!----------------------------------------------------------------------
673	!! * ROUTINE abl_zdf_tke *
674	!!
675	!! ** Purpose : Time-step Turbulente Kinetic Energy (TKE) equation
676	!!
677	!! ** Method : - source term due to shear
678	!! - source term due to stratification
679	!! - resolution of the TKE equation by inverting
680	!! a tridiagonal linear system
681	!!
682	!! ** Action : - en : now turbulent kinetic energy)
683	!! - avmu, avmv : production of TKE by shear at u and v-points
684	!! (= Kz dz[Ub] * dz[Un] )
685	!! ---------------------------------------------------------------------
686	INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
687	INTEGER, DIMENSION(1:jpi ) :: ikbl
688	REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
689	REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
690	REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2
691	REAL(wp) :: zwndi, zwndj
692	REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2
693	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2
694	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF
695	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
696	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
697	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
698	LOGICAL :: ln_Patankar = .FALSE.
699	LOGICAL :: ln_dumpvar = .FALSE.
700	LOGICAL , DIMENSION(1:jpi ) :: ln_foundl
701	!
702	tind = nt_n
703	ziRic = 1._wp / rn_Ric
704	! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
705	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
706	! ! Advance TKE equation to time n+1
707	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
708	!-------------
709	DO jj = 1, jpj ! outer loop
710	!-------------
711	!
712	! Compute vertical shear
713	DO jk = 2, jpkam1
714	DO ji = 1,jpi
715	zcff = 1.0_wp / e3w_abl( jk )**2
716	zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
717	zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
718	zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 )
719	END DO
720	END DO
721	!
722	! Compute brunt-vaisala frequency
723	DO jk = 2, jpkam1
724	DO ji = 1,jpi
725	zcff = grav * itvref / e3w_abl( jk )
726	zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta)
727	zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) &
728	& - tq_abl( ji, jj, jk , tind, jp_ta) * tq_abl( ji, jj, jk , tind, jp_qa)
729	zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi)
730	END DO
731	END DO
732	!
733	! Terms for the tridiagonal problem
734	DO jk = 2, jpkam1
735	DO ji = 1,jpi
736	zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point
737	zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
738	zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
739
740	z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal
741	z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal
742	IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE
743	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
744	& + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal
745	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side
746	ELSE
747	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
748	& + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal
749	& - e3w_abl(jk) * rDt_abl * zbuoy
750	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side
751	END IF
752	END DO
753	END DO
754
755	DO ji = 1,jpi ! vector opt.
756	zesrf = MAX( rn_esfc * ustar2(ji,jj), tke_min )
757	zetop = tke_min
758
759	z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp
760	z_elem_a ( ji, jpka ) = - 0.5 * rDt_abl * rn_Sch * (Avm_abl(ji,jj, jpka-1 )+Avm_abl(ji,jj, jpka )) / e3t_abl( jpka )
761	z_elem_c ( ji, jpka ) = 0._wp
762	z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka )
763
764	tke_abl ( ji, jj, 1 , nt_a ) = zesrf
765	tke_abl ( ji, jj, jpka , nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n )
766
767	zbn2(ji,jj, 1) = zbn2( ji,jj, 2)
768	zsh2(ji, 1) = zsh2( ji, 2)
769	zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)
770	zsh2(ji, jpka) = zsh2( ji , jpkam1)
771	END DO
772	!!
773	!! Matrix inversion
774	!! ----------------------------------------------------------
775	DO ji = 1,jpi
776	zcff = 1._wp / z_elem_b( ji, 1 )
777	zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 )
778	tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a )
779	END DO
780
781	DO jk = 2, jpka
782	DO ji = 1,jpi
783	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
784	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
785	tke_abl(ji,jj,jk,nt_a) = zcff * ( tke_abl(ji,jj,jk ,nt_a) &
786	& - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
787	END DO
788	END DO
789
790	DO jk = jpkam1,1,-1
791	DO ji = 1,jpi
792	tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a)
793	END DO
794	END DO
795
796	!!FL should not be needed because of Patankar procedure
797	tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
798
799	!!
800	!! Diagnose PBL height
801	!! ----------------------------------------------------------
802
803
804	!
805	! arrays zRH, zFC and zCF are available at this point
806	! and zFC(:, 1 ) = 0.
807	! diagnose PBL height based on zsh2 and zbn2
808	zFC ( : ,1) = 0._wp
809	ikbl( 1:jpi ) = 0
810
811	DO jk = 2,jpka
812	DO ji = 1, jpi
813	zcff = ghw_abl( jk-1 )
814	zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
815	zcff = ghw_abl( jk )
816	zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
817	zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*( &
818	zcff2 * ( zsh2( ji, jk ) - ziRic * zbn2( ji, jj, jk ) &
819	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
820	+ zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
821	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
822	& )
823	IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
824	END DO
825	END DO
826	!
827	! finalize the computation of the PBL height
828	DO ji = 1, jpi
829	jk = ikbl(ji)
830	IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
831	pblh( ji, jj ) = ( ghw_abl( jk-1 ) * zFC( ji, jk ) &
832	& - ghw_abl( jk ) * zFC( ji, jk-1 ) &
833	& ) / ( zFC( ji, jk ) - zFC( ji, jk-1 ) )
834	ELSE IF( jk==2 ) THEN
835	pblh( ji, jj ) = ghw_abl(2 )
836	ELSE
837	pblh( ji, jj ) = ghw_abl(jpka)
838	END IF
839	END DO
840	!-------------
841	END DO
842	!-------------
843	!
844	! Optional : could add pblh smoothing if pblh is noisy horizontally ...
845	IF(ln_smth_pblh) THEN
846	CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing)
847	CALL smooth_pblh( pblh, msk_abl )
848	CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing)
849	ENDIF
850	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
851	! ! Diagnostic mixing length computation
852	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
853	!
854	SELECT CASE ( nn_amxl )
855	!
856	CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom
857	# define zlup zRH
858	# define zldw zFC
859	DO jj = 1, jpj ! outer loop
860	!
861	DO ji = 1, jpi
862	mxld_abl ( ji, jj, 1 ) = mxl_min
863	mxld_abl ( ji, jj, jpka ) = mxl_min
864	mxlm_abl ( ji, jj, 1 ) = mxl_min
865	mxlm_abl ( ji, jj, jpka ) = mxl_min
866	zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc
867	zlup ( ji, jpka ) = mxl_min
868	END DO
869	!
870	DO jk = 2, jpkam1
871	DO ji = 1, jpi
872	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
873	mxlm_abl( ji, jj, jk ) = MAX( mxl_min, &
874	& SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
875	END DO
876	END DO
877	!
878	! Limit mxl
879	DO jk = jpkam1,1,-1
880	DO ji = 1, jpi
881	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
882	END DO
883	END DO
884	!
885	DO jk = 2, jpka
886	DO ji = 1, jpi
887	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
888	END DO
889	END DO
890	!
891	! DO jk = 1, jpka
892	! DO ji = 1, jpi
893	! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
894	! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) )
895	! END DO
896	! END DO
897	!
898	DO jk = 1, jpka
899	DO ji = 1, jpi
900	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
901	! zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
902	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
903	!mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
904	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
905	END DO
906	END DO
907	!
908	END DO
909	# undef zlup
910	# undef zldw
911	!
912	!
913	CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom
914	# define zlup zRH
915	# define zldw zFC
916	DO jj = 1, jpj ! outer loop
917	!
918	DO jk = 2, jpkam1
919	DO ji = 1,jpi
920	zcff = 1.0_wp / e3w_abl( jk )**2
921	zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
922	zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
923	zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
924	ENDDO
925	ENDDO
926	!
927	DO ji = 1, jpi
928	zcff = zrough(ji,jj) * rn_Lsfc
929	mxld_abl ( ji, jj, 1 ) = zcff
930	mxld_abl ( ji, jj, jpka ) = mxl_min
931	mxlm_abl ( ji, jj, 1 ) = zcff
932	mxlm_abl ( ji, jj, jpka ) = mxl_min
933	zldw ( ji, 1 ) = zcff
934	zlup ( ji, jpka ) = mxl_min
935	END DO
936	!
937	DO jk = 2, jpkam1
938	DO ji = 1, jpi
939	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
940	zcff = 2.SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rodzsh2(ji,jk) &
941	& + SQRT( rn_Rodrn_Rodzsh2(ji,jk)zsh2(ji,jk)+2.zbuoy ) )
942	mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff )
943	END DO
944	END DO
945	!
946	! Limit mxl
947	DO jk = jpkam1,1,-1
948	DO ji = 1, jpi
949	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
950	END DO
951	END DO
952	!
953	DO jk = 2, jpka
954	DO ji = 1, jpi
955	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
956	END DO
957	END DO
958	!
959	DO jk = 1, jpka
960	DO ji = 1, jpi
961	!mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
962	!zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
963	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
964	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
965	!mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
966	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
967	END DO
968	END DO
969	!
970	END DO
971	# undef zlup
972	# undef zldw
973	! !
974	CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale
975	!
976	# define zlup zRH
977	# define zldw zFC
978	! zCF is used for matrix inversion
979	!
980	DO jj = 1, jpj ! outer loop
981
982	DO ji = 1, jpi
983	zcff = zrough(ji,jj) * rn_Lsfc
984	zlup( ji, 1 ) = zcff
985	zldw( ji, 1 ) = zcff
986	zlup( ji, jpka ) = mxl_min
987	zldw( ji, jpka ) = mxl_min
988	END DO
989
990	DO jk = 2,jpka-1
991	DO ji = 1, jpi
992	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
993	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
994	END DO
995	END DO
996	!!
997	!! BL89 search for lup
998	!! ----------------------------------------------------------
999	DO jk=2,jpka-1
1000	!
1001	DO ji = 1, jpi
1002	zCF(ji,1:jpka) = 0._wp
1003	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1004	ln_foundl(ji ) = .false.
1005	END DO
1006	!
1007	DO jkup=jk+1,jpka-1
1008	DO ji = 1, jpi
1009	zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall )
1010	zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
1011	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
1012	& ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) &
1013	& + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )
1014	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1015	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
1016	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
1017	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
1018	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
1019	zlup(ji,jk) = zcff
1020	zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk)
1021	ln_foundl(ji) = .true.
1022	END IF
1023	END DO
1024	END DO
1025	!
1026	END DO
1027	!!
1028	!! BL89 search for ldwn
1029	!! ----------------------------------------------------------
1030	DO jk=2,jpka-1
1031	!
1032	DO ji = 1, jpi
1033	zCF(ji,1:jpka) = 0._wp
1034	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1035	ln_foundl(ji ) = .false.
1036	END DO
1037	!
1038	DO jkdwn=jk-1,1,-1
1039	DO ji = 1, jpi
1040	zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
1041	zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall )
1042	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
1043	& * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
1044	+ zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) )
1045	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1046	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
1047	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
1048	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
1049	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
1050	zldw(ji,jk) = zcff
1051	zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn )
1052	ln_foundl(ji) = .true.
1053	END IF
1054	END DO
1055	END DO
1056	!
1057	END DO
1058
1059	DO jk = 1, jpka
1060	DO ji = 1, jpi
1061	!zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
1062	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1063	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1064	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
1065	END DO
1066	END DO
1067
1068	END DO
1069	# undef zlup
1070	# undef zldw
1071	!
1072	CASE ( 3 ) ! Bougeault & Lacarrere 89 length-scale
1073	!
1074	# define zlup zRH
1075	# define zldw zFC
1076	! zCF is used for matrix inversion
1077	!
1078	DO jj = 1, jpj ! outer loop
1079	!
1080	DO jk = 2, jpkam1
1081	DO ji = 1,jpi
1082	zcff = 1.0_wp / e3w_abl( jk )**2
1083	zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
1084	zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
1085	zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
1086	ENDDO
1087	ENDDO
1088	zsh2(:, 1) = zsh2( :, 2)
1089	zsh2(:, jpka) = zsh2( :, jpkam1)
1090
1091	DO ji = 1, jpi
1092	zcff = zrough(ji,jj) * rn_Lsfc
1093	zlup( ji, 1 ) = zcff
1094	zldw( ji, 1 ) = zcff
1095	zlup( ji, jpka ) = mxl_min
1096	zldw( ji, jpka ) = mxl_min
1097	END DO
1098
1099	DO jk = 2,jpka-1
1100	DO ji = 1, jpi
1101	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
1102	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
1103	END DO
1104	END DO
1105	!!
1106	!! BL89 search for lup
1107	!! ----------------------------------------------------------
1108	DO jk=2,jpka-1
1109	!
1110	DO ji = 1, jpi
1111	zCF(ji,1:jpka) = 0._wp
1112	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1113	ln_foundl(ji ) = .false.
1114	END DO
1115	!
1116	DO jkup=jk+1,jpka-1
1117	DO ji = 1, jpi
1118	zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall )
1119	zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
1120	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
1121	& ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) &
1122	& + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) &
1123	& + 0.5_wp * e3t_abl(jkup) * rn_Rod * &
1124	& ( SQRT(tke_abl( ji, jj, jkup , nt_a ))*zsh2(ji,jkup ) &
1125	& + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) )
1126
1127	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1128	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
1129	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
1130	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
1131	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
1132	zlup(ji,jk) = zcff
1133	zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk)
1134	ln_foundl(ji) = .true.
1135	END IF
1136	END DO
1137	END DO
1138	!
1139	END DO
1140	!!
1141	!! BL89 search for ldwn
1142	!! ----------------------------------------------------------
1143	DO jk=2,jpka-1
1144	!
1145	DO ji = 1, jpi
1146	zCF(ji,1:jpka) = 0._wp
1147	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1148	ln_foundl(ji ) = .false.
1149	END DO
1150	!
1151	DO jkdwn=jk-1,1,-1
1152	DO ji = 1, jpi
1153	zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
1154	zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall )
1155	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
1156	& * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
1157	& +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) &
1158	& + 0.5_wp * e3t_abl(jkup) * rn_Rod * &
1159	& ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) &
1160	& + SQRT(tke_abl( ji, jj, jkdwn , nt_a ))*zsh2(ji,jkdwn ) )
1161
1162	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1163	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
1164	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
1165	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
1166	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
1167	zldw(ji,jk) = zcff
1168	zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn )
1169	ln_foundl(ji) = .true.
1170	END IF
1171	END DO
1172	END DO
1173	!
1174	END DO
1175
1176	DO jk = 1, jpka
1177	DO ji = 1, jpi
1178	!zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
1179	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1180	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1181	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
1182	END DO
1183	END DO
1184
1185	END DO
1186	# undef zlup
1187	# undef zldw
1188	!
1189	!
1190	!
1191	END SELECT
1192	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1193	! ! Finalize the computation of turbulent visc./diff.
1194	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1195
1196	!-------------
1197	DO jj = 1, jpj ! outer loop
1198	!-------------
1199	DO jk = 1, jpka
1200	DO ji = 1, jpi ! vector opt.
1201	zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) &
1202	& * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
1203	zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z)
1204	zcff = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
1205	!!FL: MAX function probably useless because of the definition of mxl_min
1206	Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak )
1207	Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak )
1208	END DO
1209	END DO
1210	!-------------
1211	END DO
1212	!-------------
1213
1214	!---------------------------------------------------------------------------------------------------
1215	END SUBROUTINE abl_zdf_tke
1216	!===================================================================================================
1217
1218
1219	!===================================================================================================
1220	SUBROUTINE smooth_pblh( pvar2d, msk )
1221	!---------------------------------------------------------------------------------------------------
1222
1223	!!----------------------------------------------------------------------
1224	!! * ROUTINE smooth_pblh *
1225	!!
1226	!! ** Purpose : 2D Hanning filter on atmospheric PBL height
1227	!!
1228	!! ---------------------------------------------------------------------
1229	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk
1230	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
1231	INTEGER :: ji,jj
1232	REAL(wp) :: smth_a, smth_b
1233	REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY
1234	REAL(wp) :: zumsk,zvmsk
1235	!!
1236	!!=========================================================
1237	!!
1238	!! Hanning filter
1239	smth_a = 1._wp / 8._wp
1240	smth_b = 1._wp / 4._wp
1241	!
1242	DO_2D_11_10
1243	zumsk = msk(ji,jj) * msk(ji+1,jj)
1244	zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk
1245	END_2D
1246
1247	DO_2D_10_11
1248	zvmsk = msk(ji,jj) * msk(ji,jj+1)
1249	zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk
1250	END_2D
1251
1252	DO_2D_10_00
1253	zFY ( ji, jj ) = zdY ( ji, jj ) &
1254	& + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) &
1255	& - (zdX ( ji, jj ) - zdX( ji-1, jj )) )
1256	END_2D
1257
1258	DO_2D_00_10
1259	zFX( ji, jj ) = zdX( ji, jj ) &
1260	& + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) &
1261	& - (zdY( ji , jj ) - zdY( ji , jj-1)) )
1262	END_2D
1263
1264	DO_2D_00_00
1265	pvar2d( ji ,jj ) = pvar2d( ji ,jj ) &
1266	& + msk(ji,jj) * smth_b * ( &
1267	& zFX( ji, jj ) - zFX( ji-1, jj ) &
1268	& +zFY( ji, jj ) - zFY( ji, jj-1 ) )
1269	END_2D
1270
1271	!---------------------------------------------------------------------------------------------------
1272	END SUBROUTINE smooth_pblh
1273	!===================================================================================================
1274
1275	!!======================================================================
1276	END MODULE ablmod

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

Download in other formats: