Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

ablmod.F90

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/ABL/ablmod.F90

Download in other formats: