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_r13296_HPC-07_mocavero_mpi3/src/ABL – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ABL/ablmod.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago
Add neighborhood collectives calls in the NEMO src - ticket #2496
File size: 63.8 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( 1, 1, 1, 1 )
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( 1, 1, 1, 1 )
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( 0, 1, 0, 1 )
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( 1, 1, 1, 1 )
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	#if defined key_mpi3
532	CALL lbc_lnk_nc_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1._wp, v_abl(:,:,:,nt_a) , 'T', -1._wp )
533	CALL lbc_lnk_nc_multi( '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...
534	#else
535	CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1._wp, v_abl(:,:,:,nt_a) , 'T', -1._wp )
536	CALL lbc_lnk_multi( '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...
537	#endif
538	!
539	#if defined key_iomput
540	! 2D & first ABL level
541	IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh", pblh(:,: ) )
542	IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) )
543	IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) )
544	IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl", tq_abl(:,:,2,nt_a,jp_ta) )
545	IF ( iom_use("qz1_abl") ) CALL iom_put ( "qz1_abl", tq_abl(:,:,2,nt_a,jp_qa) )
546	IF ( iom_use("uz1_dta") ) CALL iom_put ( "uz1_dta", pu_dta(:,:,2 ) )
547	IF ( iom_use("vz1_dta") ) CALL iom_put ( "vz1_dta", pv_dta(:,:,2 ) )
548	IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta", pt_dta(:,:,2 ) )
549	IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta", pq_dta(:,:,2 ) )
550	! debug 2D
551	IF( ln_geos_winds ) THEN
552	IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) )
553	IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) )
554	END IF
555	IF( ln_hpgls_frc ) THEN
556	IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
557	IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
558	END IF
559	! 3D (all ABL levels)
560	IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) )
561	IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) )
562	IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) )
563	IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) )
564	IF ( iom_use("tke_abl" ) ) CALL iom_put ( "tke_abl" , tke_abl(:,:,2:jpka,nt_a ) )
565	IF ( iom_use("avm_abl" ) ) CALL iom_put ( "avm_abl" , avm_abl(:,:,2:jpka ) )
566	IF ( iom_use("avt_abl" ) ) CALL iom_put ( "avt_abl" , avt_abl(:,:,2:jpka ) )
567	IF ( iom_use("mxlm_abl") ) CALL iom_put ( "mxlm_abl", mxlm_abl(:,:,2:jpka ) )
568	IF ( iom_use("mxld_abl") ) CALL iom_put ( "mxld_abl", mxld_abl(:,:,2:jpka ) )
569	! debug 3D
570	IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta", pu_dta(:,:,2:jpka) )
571	IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta", pv_dta(:,:,2:jpka) )
572	IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta", pt_dta(:,:,2:jpka) )
573	IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta", pq_dta(:,:,2:jpka) )
574	IF( ln_geos_winds ) THEN
575	IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka) )
576	IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) )
577	END IF
578	IF( ln_hpgls_frc ) THEN
579	IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
580	IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
581	END IF
582	#endif
583	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
584	! ! 7 *** Finalize flux computation
585	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
586	!
587	DO_2D( 1, 1, 1, 1 )
588	ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta )
589	zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa )
590	zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
591	psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention
592	pevp( ji, jj ) = rn_efacMAX( 0._wp, zcff pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
593	rhoa( ji, jj ) = zcff
594	END_2D
595
596	DO_2D( 0, 1, 0, 1 )
597	zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) )
598	zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) )
599	END_2D
600	!
601	#if defined key_mpi3
602	CALL lbc_lnk_nc_multi( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp )
603	#else
604	CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp )
605	#endif
606	!
607	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
608	DO_2D( 1, 1, 1, 1 )
609	zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) &
610	& + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj)
611	zztmp = rhoa(ji,jj) * pcd_du(ji,jj)
612
613	pwndm (ji,jj) = zcff
614	ptaum (ji,jj) = zztmp * zcff
615	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
616	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
617	END_2D
618	! ... utau, vtau at U- and V_points, resp.
619	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
620	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
621	DO_2D( 0, 0, 0, 0 )
622	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
623	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
624	ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
625	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
626	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
627	ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
628	END_2D
629	!
630	#if defined key_mpi3
631	CALL lbc_lnk_nc_multi( 'ablmod', ptaui(:,:), 'U', -1.0_wp, ptauj(:,:), 'V', -1.0_wp )
632	#else
633	CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1.0_wp, ptauj(:,:), 'V', -1.0_wp )
634	#endif
635
636	CALL iom_put( "taum_oce", ptaum )
637
638	IF(sn_cfctl%l_prtctl) THEN
639	CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', mask1=umask, &
640	& tab2d_2=ptauj , clinfo2=' vtau : ', mask2=vmask )
641	CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' )
642	ENDIF
643
644	#if defined key_si3
645	! ------------------------------------------------------------ !
646	! Wind stress relative to the moving ice ( U10m - U_ice ) !
647	! ------------------------------------------------------------ !
648	DO_2D( 0, 0, 0, 0 )
649	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) ) &
650	& * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
651	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) ) &
652	& * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
653	END_2D
654	#if defined key_mpi3
655	CALL lbc_lnk_nc_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
656	#else
657	CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
658	#endif
659	!
660	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
661	& , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
662	! ------------------------------------------------------------ !
663	! Wind stress relative to the moving ice ( U10m - U_ice ) !
664	! ------------------------------------------------------------ !
665	DO_2D( 0, 0, 0, 0 )
666
667	zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) )
668	zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
669
670	ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &
671	& + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &
672	& * ( zztmp1 - pssu_ice(ji,jj) )
673	ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &
674	& + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &
675	& * ( zztmp2 - pssv_ice(ji,jj) )
676	END_2D
677	#if defined key_mpi3
678	CALL lbc_lnk_nc_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp )
679	#else
680	CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp )
681	#endif
682	!
683	IF(sn_cfctl%l_prtctl) THEN
684	CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask, &
685	& tab2d_2=ptauj_ice , clinfo2=' vtau_ice : ', mask2=vmask )
686	END IF
687	#endif
688	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
689	! ! 8 *** Swap time indices for the next timestep
690	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
691	nt_n = 1 + MOD( nt_n, 2)
692	nt_a = 1 + MOD( nt_a, 2)
693	!
694	!---------------------------------------------------------------------------------------------------
695	END SUBROUTINE abl_stp
696	!===================================================================================================
697
698
699
700
701	!===================================================================================================
702	SUBROUTINE abl_zdf_tke( )
703	!---------------------------------------------------------------------------------------------------
704
705	!!----------------------------------------------------------------------
706	!! * ROUTINE abl_zdf_tke *
707	!!
708	!! ** Purpose : Time-step Turbulente Kinetic Energy (TKE) equation
709	!!
710	!! ** Method : - source term due to shear
711	!! - source term due to stratification
712	!! - resolution of the TKE equation by inverting
713	!! a tridiagonal linear system
714	!!
715	!! ** Action : - en : now turbulent kinetic energy)
716	!! - avmu, avmv : production of TKE by shear at u and v-points
717	!! (= Kz dz[Ub] * dz[Un] )
718	!! ---------------------------------------------------------------------
719	INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
720	INTEGER, DIMENSION(1:jpi ) :: ikbl
721	REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
722	REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
723	REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89
724	REAL(wp) :: zwndi, zwndj
725	REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2
726	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2
727	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF
728	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
729	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
730	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
731	LOGICAL :: ln_Patankar = .FALSE.
732	LOGICAL :: ln_dumpvar = .FALSE.
733	LOGICAL , DIMENSION(1:jpi ) :: ln_foundl
734	!
735	tind = nt_n
736	ziRic = 1._wp / rn_Ric
737	! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
738	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
739	! ! Advance TKE equation to time n+1
740	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
741	!-------------
742	DO jj = 1, jpj ! outer loop
743	!-------------
744	!
745	! Compute vertical shear
746	DO jk = 2, jpkam1
747	DO ji = 1, jpi
748	zcff = 1.0_wp / e3w_abl( jk )**2
749	zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
750	zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
751	zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 )
752	END DO
753	END DO
754	!
755	! Compute brunt-vaisala frequency
756	DO jk = 2, jpkam1
757	DO ji = 1,jpi
758	zcff = grav * itvref / e3w_abl( jk )
759	zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta)
760	zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) &
761	& - tq_abl( ji, jj, jk , tind, jp_ta) * tq_abl( ji, jj, jk , tind, jp_qa)
762	zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi)
763	END DO
764	END DO
765	!
766	! Terms for the tridiagonal problem
767	DO jk = 2, jpkam1
768	DO ji = 1, jpi
769	zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point
770	zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
771	zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
772
773	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
774	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
775	IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE
776	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
777	& + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal
778	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side
779	ELSE
780	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
781	& + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal
782	& - e3w_abl(jk) * rDt_abl * zbuoy
783	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side
784	END IF
785	END DO
786	END DO
787
788	DO ji = 1,jpi ! vector opt.
789	zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min )
790	zetop = tke_min
791
792	z_elem_a ( ji, 1 ) = 0._wp
793	z_elem_c ( ji, 1 ) = 0._wp
794	z_elem_b ( ji, 1 ) = 1._wp
795	tke_abl ( ji, jj, 1, nt_a ) = zesrf
796
797	!++ Top Neumann B.C.
798	!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 )
799	!z_elem_c ( ji, jpka ) = 0._wp
800	!z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka )
801	!tke_abl ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n )
802
803	!++ Top Dirichlet B.C.
804	z_elem_a ( ji, jpka ) = 0._wp
805	z_elem_c ( ji, jpka ) = 0._wp
806	z_elem_b ( ji, jpka ) = 1._wp
807	tke_abl ( ji, jj, jpka, nt_a ) = zetop
808
809	zbn2 ( ji, jj, 1 ) = zbn2 ( ji, jj, 2 )
810	zsh2 ( ji, 1 ) = zsh2 ( ji, 2 )
811	zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 )
812	zsh2 ( ji, jpka ) = zsh2 ( ji , jpkam1 )
813	END DO
814	!!
815	!! Matrix inversion
816	!! ----------------------------------------------------------
817	DO ji = 1,jpi
818	zcff = 1._wp / z_elem_b( ji, 1 )
819	zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 )
820	tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a )
821	END DO
822
823	DO jk = 2, jpka
824	DO ji = 1,jpi
825	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
826	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
827	tke_abl(ji,jj,jk,nt_a) = zcff * ( tke_abl(ji,jj,jk ,nt_a) &
828	& - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
829	END DO
830	END DO
831
832	DO jk = jpkam1,1,-1
833	DO ji = 1,jpi
834	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)
835	END DO
836	END DO
837
838	!!FL should not be needed because of Patankar procedure
839	tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
840
841	!!
842	!! Diagnose PBL height
843	!! ----------------------------------------------------------
844
845
846	!
847	! arrays zRH, zFC and zCF are available at this point
848	! and zFC(:, 1 ) = 0.
849	! diagnose PBL height based on zsh2 and zbn2
850	zFC ( : ,1) = 0._wp
851	ikbl( 1:jpi ) = 0
852
853	DO jk = 2,jpka
854	DO ji = 1, jpi
855	zcff = ghw_abl( jk-1 )
856	zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
857	zcff = ghw_abl( jk )
858	zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
859	zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*( &
860	zcff2 * ( zsh2( ji, jk ) - ziRic * zbn2( ji, jj, jk ) &
861	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
862	+ zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
863	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
864	& )
865	IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
866	END DO
867	END DO
868	!
869	! finalize the computation of the PBL height
870	DO ji = 1, jpi
871	jk = ikbl(ji)
872	IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
873	pblh( ji, jj ) = ( ghw_abl( jk-1 ) * zFC( ji, jk ) &
874	& - ghw_abl( jk ) * zFC( ji, jk-1 ) &
875	& ) / ( zFC( ji, jk ) - zFC( ji, jk-1 ) )
876	ELSE IF( jk==2 ) THEN
877	pblh( ji, jj ) = ghw_abl(2 )
878	ELSE
879	pblh( ji, jj ) = ghw_abl(jpka)
880	END IF
881	END DO
882	!-------------
883	END DO
884	!-------------
885	!
886	! Optional : could add pblh smoothing if pblh is noisy horizontally ...
887	IF(ln_smth_pblh) THEN
888	#if defined key_mpi3
889	CALL lbc_lnk_nc_multi( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
890	#else
891	CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
892	#endif
893	CALL smooth_pblh( pblh, msk_abl )
894	#if defined key_mpi3
895	CALL lbc_lnk_nc_multi( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
896	#else
897	CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
898	#endif
899	ENDIF
900	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
901	! ! Diagnostic mixing length computation
902	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
903	!
904	SELECT CASE ( nn_amxl )
905	!
906	CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom
907	# define zlup zRH
908	# define zldw zFC
909	DO jj = 1, jpj ! outer loop
910	!
911	DO ji = 1, jpi
912	mxld_abl( ji, jj, 1 ) = mxl_min
913	mxld_abl( ji, jj, jpka ) = mxl_min
914	mxlm_abl( ji, jj, 1 ) = mxl_min
915	mxlm_abl( ji, jj, jpka ) = mxl_min
916	zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc
917	zlup ( ji, jpka ) = mxl_min
918	END DO
919	!
920	DO jk = 2, jpkam1
921	DO ji = 1, jpi
922	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
923	mxlm_abl( ji, jj, jk ) = MAX( mxl_min, &
924	& SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
925	END DO
926	END DO
927	!
928	! Limit mxl
929	DO jk = jpkam1,1,-1
930	DO ji = 1, jpi
931	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
932	END DO
933	END DO
934	!
935	DO jk = 2, jpka
936	DO ji = 1, jpi
937	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
938	END DO
939	END DO
940	!
941	! DO jk = 1, jpka
942	! DO ji = 1, jpi
943	! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
944	! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) )
945	! END DO
946	! END DO
947	!
948	DO jk = 1, jpka
949	DO ji = 1, jpi
950	! zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
951	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
952	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
953	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
954	END DO
955	END DO
956	!
957	END DO
958	# undef zlup
959	# undef zldw
960	!
961	!
962	CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom
963	# define zlup zRH
964	# define zldw zFC
965	DO jj = 1, jpj ! outer loop
966	!
967	DO jk = 2, jpkam1
968	DO ji = 1,jpi
969	zcff = 1.0_wp / e3w_abl( jk )**2
970	zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
971	zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
972	zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
973	ENDDO
974	ENDDO
975	!
976	DO ji = 1, jpi
977	zcff = zrough(ji,jj) * rn_Lsfc
978	mxld_abl ( ji, jj, 1 ) = zcff
979	mxld_abl ( ji, jj, jpka ) = mxl_min
980	mxlm_abl ( ji, jj, 1 ) = zcff
981	mxlm_abl ( ji, jj, jpka ) = mxl_min
982	zldw ( ji, 1 ) = zcff
983	zlup ( ji, jpka ) = mxl_min
984	END DO
985	!
986	DO jk = 2, jpkam1
987	DO ji = 1, jpi
988	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
989	zcff = 2.0_wpSQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rodzsh2(ji,jk) &
990	& + SQRT(rn_Rodrn_Rodzsh2(ji,jk)zsh2(ji,jk)+2.0_wpzbuoy ) )
991	mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff )
992	END DO
993	END DO
994	!
995	! Limit mxl
996	DO jk = jpkam1,1,-1
997	DO ji = 1, jpi
998	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
999	END DO
1000	END DO
1001	!
1002	DO jk = 2, jpka
1003	DO ji = 1, jpi
1004	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
1005	END DO
1006	END DO
1007	!
1008	DO jk = 1, jpka
1009	DO ji = 1, jpi
1010	!mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1011	!zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
1012	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1013	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1014	!mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
1015	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
1016	END DO
1017	END DO
1018	!
1019	END DO
1020	# undef zlup
1021	# undef zldw
1022	!
1023	CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale
1024	!
1025	# define zlup zRH
1026	# define zldw zFC
1027	! zCF is used for matrix inversion
1028	!
1029	DO jj = 1, jpj ! outer loop
1030
1031	DO ji = 1, jpi
1032	zcff = zrough(ji,jj) * rn_Lsfc
1033	zlup( ji, 1 ) = zcff
1034	zldw( ji, 1 ) = zcff
1035	zlup( ji, jpka ) = mxl_min
1036	zldw( ji, jpka ) = mxl_min
1037	END DO
1038
1039	DO jk = 2,jpka-1
1040	DO ji = 1, jpi
1041	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
1042	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
1043	END DO
1044	END DO
1045	!!
1046	!! BL89 search for lup
1047	!! ----------------------------------------------------------
1048	DO jk=2,jpka-1
1049	!
1050	DO ji = 1, jpi
1051	zCF(ji,1:jpka) = 0._wp
1052	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1053	ln_foundl(ji ) = .false.
1054	END DO
1055	!
1056	DO jkup=jk+1,jpka-1
1057	DO ji = 1, jpi
1058	zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall )
1059	zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
1060	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
1061	& ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) &
1062	& + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )
1063	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1064	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
1065	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
1066	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
1067	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
1068	zlup(ji,jk) = zcff
1069	zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk)
1070	ln_foundl(ji) = .true.
1071	END IF
1072	END DO
1073	END DO
1074	!
1075	END DO
1076	!!
1077	!! BL89 search for ldwn
1078	!! ----------------------------------------------------------
1079	DO jk=2,jpka-1
1080	!
1081	DO ji = 1, jpi
1082	zCF(ji,1:jpka) = 0._wp
1083	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1084	ln_foundl(ji ) = .false.
1085	END DO
1086	!
1087	DO jkdwn=jk-1,1,-1
1088	DO ji = 1, jpi
1089	zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
1090	zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall )
1091	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
1092	& * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
1093	+ zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) )
1094	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1095	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
1096	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
1097	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
1098	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
1099	zldw(ji,jk) = zcff
1100	zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn )
1101	ln_foundl(ji) = .true.
1102	END IF
1103	END DO
1104	END DO
1105	!
1106	END DO
1107
1108	DO jk = 1, jpka
1109	DO ji = 1, jpi
1110	!zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
1111	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1112	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1113	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
1114	END DO
1115	END DO
1116
1117	END DO
1118	# undef zlup
1119	# undef zldw
1120	!
1121	CASE ( 3 ) ! Bougeault & Lacarrere 89 length-scale
1122	!
1123	# define zlup zRH
1124	# define zldw zFC
1125	! zCF is used for matrix inversion
1126	!
1127	DO jj = 1, jpj ! outer loop
1128	!
1129	DO jk = 2, jpkam1
1130	DO ji = 1,jpi
1131	zcff = 1.0_wp / e3w_abl( jk )**2
1132	zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
1133	zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
1134	zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
1135	ENDDO
1136	ENDDO
1137	zsh2(:, 1) = zsh2( :, 2)
1138	zsh2(:, jpka) = zsh2( :, jpkam1)
1139
1140	DO ji = 1, jpi
1141	zcff = zrough(ji,jj) * rn_Lsfc
1142	zlup( ji, 1 ) = zcff
1143	zldw( ji, 1 ) = zcff
1144	zlup( ji, jpka ) = mxl_min
1145	zldw( ji, jpka ) = mxl_min
1146	END DO
1147
1148	DO jk = 2,jpka-1
1149	DO ji = 1, jpi
1150	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
1151	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
1152	END DO
1153	END DO
1154	!!
1155	!! BL89 search for lup
1156	!! ----------------------------------------------------------
1157	DO jk=2,jpka-1
1158	!
1159	DO ji = 1, jpi
1160	zCF(ji,1:jpka) = 0._wp
1161	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1162	ln_foundl(ji ) = .false.
1163	END DO
1164	!
1165	DO jkup=jk+1,jpka-1
1166	DO ji = 1, jpi
1167	zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall )
1168	zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
1169	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
1170	& ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) &
1171	& + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) &
1172	& + 0.5_wp * e3t_abl(jkup) * rn_Rod * &
1173	& ( SQRT(tke_abl( ji, jj, jkup , nt_a ))*zsh2(ji,jkup ) &
1174	& + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) )
1175
1176	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1177	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
1178	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
1179	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
1180	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
1181	zlup(ji,jk) = zcff
1182	zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk)
1183	ln_foundl(ji) = .true.
1184	END IF
1185	END DO
1186	END DO
1187	!
1188	END DO
1189	!!
1190	!! BL89 search for ldwn
1191	!! ----------------------------------------------------------
1192	DO jk=2,jpka-1
1193	!
1194	DO ji = 1, jpi
1195	zCF(ji,1:jpka) = 0._wp
1196	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
1197	ln_foundl(ji ) = .false.
1198	END DO
1199	!
1200	DO jkdwn=jk-1,1,-1
1201	DO ji = 1, jpi
1202	zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
1203	zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall )
1204	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
1205	& * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
1206	& +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) &
1207	& + 0.5_wp * e3t_abl(jkup) * rn_Rod * &
1208	& ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) &
1209	& + SQRT(tke_abl( ji, jj, jkdwn , nt_a ))*zsh2(ji,jkdwn ) )
1210
1211	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1212	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
1213	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
1214	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
1215	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
1216	zldw(ji,jk) = zcff
1217	zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn )
1218	ln_foundl(ji) = .true.
1219	END IF
1220	END DO
1221	END DO
1222	!
1223	END DO
1224
1225	DO jk = 1, jpka
1226	DO ji = 1, jpi
1227	!zcff = 2.SQRT(2.)( zldw( ji, jk )(-2._wp/3._wp) + zlup( ji, jk )(-2._wp/3._wp) )**(-3._wp/2._wp)
1228	zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1229	mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1230	mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
1231	END DO
1232	END DO
1233
1234	END DO
1235	# undef zlup
1236	# undef zldw
1237	!
1238	!
1239	END SELECT
1240	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1241	! ! Finalize the computation of turbulent visc./diff.
1242	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1243
1244	!-------------
1245	DO jj = 1, jpj ! outer loop
1246	!-------------
1247	DO jk = 1, jpka
1248	DO ji = 1, jpi ! vector opt.
1249	zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) &
1250	& * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
1251	zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z)
1252	zcff = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
1253	!!FL: MAX function probably useless because of the definition of mxl_min
1254	Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak )
1255	Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak )
1256	END DO
1257	END DO
1258	!-------------
1259	END DO
1260	!-------------
1261
1262	!---------------------------------------------------------------------------------------------------
1263	END SUBROUTINE abl_zdf_tke
1264	!===================================================================================================
1265
1266
1267	!===================================================================================================
1268	SUBROUTINE smooth_pblh( pvar2d, msk )
1269	!---------------------------------------------------------------------------------------------------
1270
1271	!!----------------------------------------------------------------------
1272	!! * ROUTINE smooth_pblh *
1273	!!
1274	!! ** Purpose : 2D Hanning filter on atmospheric PBL height
1275	!!
1276	!! ---------------------------------------------------------------------
1277	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk
1278	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
1279	INTEGER :: ji,jj
1280	REAL(wp) :: smth_a, smth_b
1281	REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY
1282	REAL(wp) :: zumsk,zvmsk
1283	!!
1284	!!=========================================================
1285	!!
1286	!! Hanning filter
1287	smth_a = 1._wp / 8._wp
1288	smth_b = 1._wp / 4._wp
1289	!
1290	DO_2D( 1, 1, 1, 0 )
1291	zumsk = msk(ji,jj) * msk(ji+1,jj)
1292	zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk
1293	END_2D
1294
1295	DO_2D( 1, 0, 1, 1 )
1296	zvmsk = msk(ji,jj) * msk(ji,jj+1)
1297	zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk
1298	END_2D
1299
1300	DO_2D( 1, 0, 0, 0 )
1301	zFY ( ji, jj ) = zdY ( ji, jj ) &
1302	& + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) &
1303	& - (zdX ( ji, jj ) - zdX( ji-1, jj )) )
1304	END_2D
1305
1306	DO_2D( 0, 0, 1, 0 )
1307	zFX( ji, jj ) = zdX( ji, jj ) &
1308	& + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) &
1309	& - (zdY( ji , jj ) - zdY( ji , jj-1)) )
1310	END_2D
1311
1312	DO_2D( 0, 0, 0, 0 )
1313	pvar2d( ji ,jj ) = pvar2d( ji ,jj ) &
1314	& + msk(ji,jj) * smth_b * ( &
1315	& zFX( ji, jj ) - zFX( ji-1, jj ) &
1316	& +zFY( ji, jj ) - zFY( ji, jj-1 ) )
1317	END_2D
1318
1319	!---------------------------------------------------------------------------------------------------
1320	END SUBROUTINE smooth_pblh
1321	!===================================================================================================
1322
1323	!!======================================================================
1324	END MODULE ablmod

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

Download in other formats: