Context Navigation

ablmod.F90 @ 15554

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

Download in other formats: