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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

ablmod.F90 in NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/ABL – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/ABL/ablmod.F90 @ 13229

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

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

Download in other formats: