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/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90 @ 11305

Last change on this file since 11305 was 11305, checked in by smasson, 5 years ago
dev_r11265_ABL : first implementation of ABL, see #2131
File size: 43.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
20	USE sbcblk ! use some physical constants for flux computation
21	!
22	USE prtctl ! Print control (prt_ctl routine)
23	USE iom ! IOM library
24	USE in_out_manager ! I/O manager
25	USE lib_mpp ! MPP library
26	USE timing ! Timing
27
28	IMPLICIT NONE
29
30	PUBLIC abl_stp ! called by sbcabl.F90
31	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2
32	!! * Substitutions
33	# include "vectopt_loop_substitute.h90"
34
35	CONTAINS
36
37
38	!===================================================================================================
39	SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, & ! in
40	& pu_dta, pv_dta, pt_dta, pq_dta, &
41	& pslp_dta, pgu_dta, pgv_dta, &
42	& pcd_du, psen, pevp, & ! in/out
43	& pwndm, ptaui, ptauj, ptaum )
44	!---------------------------------------------------------------------------------------------------
45
46	!!---------------------------------------------------------------------
47	!! * ROUTINE abl_stp *
48	!!
49	!! ** Purpose : Time-integration of the ABL model
50	!!
51	!! ** Method : Compute atmospheric variables : vertical turbulence
52	!! + Coriolis term + newtonian relaxation
53	!!
54	!! ** Action : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
55	!! - Advance tracers to time n+1 (Euler backward scheme)
56	!! - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
57	!! - Advance u,v to time n+1 (Euler backward scheme)
58	!! - Apply newtonian relaxation on the dynamics and the tracers
59	!! - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
60	!!
61	!!----------------------------------------------------------------------
62	INTEGER , INTENT(in ) :: kt ! time step index
63	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psst ! sea-surface temperature [Celsius]
64	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu ! sea-surface u (U-point)
65	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point)
66	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq ! sea-surface humidity
67	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pu_dta ! large-scale windi
68	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pv_dta ! large-scale windj
69	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgu_dta ! large-scale hpgi
70	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgv_dta ! large-scale hpgj
71	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pt_dta ! large-scale pot. temp.
72	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pq_dta ! large-scale humidity
73	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp_dta ! sea-level pressure
74	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du ! Cd x Du (T-point)
75	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du
76	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du
77	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! \|\|uwnd\|\|
78	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux
79	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy
80	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! \|\|tau\|\|
81	!
82	REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zrhoa, zwnd_i, zwnd_j
83	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC
84	REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF
85	REAL(wp), DIMENSION(1:jpi, jptq ) :: zBC
86	!
87	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
88	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
89	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
90	!
91	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka ) :: z_cft !--FL--to be removed after the test phase
92	!
93	INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices
94	REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1
95	REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv
96	LOGICAL :: ln_old_coriolis = .FALSE. ! possibility to switch off Coriolis term
97	!
98	!!---------------------------------------------------------------------
99	!
100	IF(lwp .AND. kt == nit000) THEN ! control print
101	WRITE(numout,*)
102	WRITE(numout,*) 'abl_stp : ABL time stepping'
103	WRITE(numout,*) '~~~~~~'
104	ENDIF
105	!
106	IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
107	!! Compute ustar squared as Cd \|\| Uatm-Uoce \|\|^2
108	!! needed for surface boundary condition of TKE
109	!! pwndm contains \| U10m - U_oce \| (see blk_oce_1 in sbcblk)
110	DO jj = 1,jpj
111	DO ji = 1,jpi
112	ustar2(ji,jj) = pCd_du(ji,jj) * pwndm(ji,jj)
113	END DO
114	END DO
115	!
116	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
117	! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
118	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
119
120	CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
121
122	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
123	! ! 2 *** Advance tracers to time n+1
124	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
125
126	!-------------
127	DO jj = 1, jpj ! outer loop !--> tq_abl computed on (1:jpi) x (1:jpj)
128	!-------------
129	! Compute matrix elements for interior points
130	DO jk = 3, jpkam1
131	DO ji = 1, jpi ! vector opt.
132	z_elem_a( ji, jk ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
133	z_elem_c( ji, jk ) = - rdt * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
134	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
135	END DO
136	END DO
137	! Boundary conditions
138	DO ji = 1, jpi ! vector opt.
139	! Neumann at the bottom
140	z_elem_a( ji, 2 ) = 0._wp
141	z_elem_c( ji, 2 ) = - rdt * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )
142	! Homogeneous Neumann at the top
143	z_elem_a( ji, jpka ) = - rdt * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
144	z_elem_c( ji, jpka ) = 0._wp
145	z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
146	END DO
147
148	DO jtra = 1,jptq ! loop on active tracers
149
150	DO jk = 3, jpkam1
151	DO ji = fs_2, fs_jpim1
152	tq_abl ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl ( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side
153	END DO
154	END DO
155
156	IF(jtra == jp_ta) THEN
157	DO ji = 1,jpi ! surface boundary condition for temperature
158	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * psen(ji, jj)
159	tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) &
160	& + rdt * psen(ji, jj) * ( psst(ji, jj) + rt0 )
161	tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
162	END DO
163	ELSE
164	DO ji = 1,jpi ! surface boundary condition for humidity
165	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * pevp(ji, jj)
166	tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) &
167	& + rdt * pevp(ji, jj) * pssq(ji, jj)
168	tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
169	END DO
170	END IF
171	!!
172	!! Matrix inversion
173	!! ----------------------------------------------------------
174	DO ji = 1,jpi
175	zcff = 1._wp / z_elem_b( ji, 2 )
176	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
177	tq_abl(ji,jj,2,nt_a,jtra) = zcff * tq_abl(ji,jj,2,nt_a,jtra)
178	END DO
179
180	DO jk = 3, jpka
181	DO ji = 1,jpi
182	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
183	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
184	tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) &
185	& - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
186	END DO
187	END DO
188	!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
189	DO jk = jpkam1,2,-1
190	DO ji = 1,jpi
191	tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) + &
192	& zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
193	END DO
194	END DO
195
196	END DO !<-- loop on tracers
197	!!
198	!-------------
199	END DO ! end outer loop
200	!-------------
201
202
203	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
204	! ! 3 *** Compute Coriolis term with geostrophic guide
205	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
206	!-------------
207	DO jk = 2, jpka ! outer loop
208	!-------------
209	!
210	! Advance u_abl & v_abl to time n+1
211	DO jj = 1, jpj
212	DO ji = 1, jpi
213	zcff = ( ff_t(ji,jj) * rdt )( ff_t(ji,jj) rdt ) ! (f dt)**2
214
215	u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
216	& (1._wp-gamma_Cor(1._wp-gamma_Cor)zcff)*u_abl( ji, jj, jk, nt_n ) &
217	& + rdt * ff_t(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) &
218	& / (1._wp + gamma_Corgamma_Corzcff)
219
220	v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
221	& (1._wp-gamma_Cor(1._wp-gamma_Cor)zcff)*v_abl( ji, jj, jk, nt_n ) &
222	& - rdt * ff_t(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) &
223	& / (1._wp + gamma_Corgamma_Corzcff)
224	END DO
225	END DO
226	!
227	!-------------
228	END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
229	!-------------
230	!
231	IF( ln_geos_winds ) THEN
232	DO jj = 1, jpj ! outer loop
233	DO jk = 1, jpka
234	DO ji = 1, jpi
235	u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) &
236	& - rdt * e3t_abl(jk) * ff_t(ji , jj) * pgv_dta(ji ,jj ,jk)
237	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) &
238	& + rdt * e3t_abl(jk) * ff_t(ji, jj ) * pgu_dta(ji ,jj ,jk)
239	END DO
240	END DO
241	END DO
242	END IF
243	!-------------
244	!
245	IF( ln_hpgls_frc ) THEN
246	DO jj = 1, jpj ! outer loop
247	DO jk = 1, jpka
248	DO ji = 1, jpi
249	u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk)
250	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk)
251	ENDDO
252	ENDDO
253	ENDDO
254	END IF
255
256	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
257	! ! 4 *** Advance u,v to time n+1
258	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
259	!
260	! Vertical diffusion for u_abl
261	!-------------
262	DO jj = 1, jpj ! outer loop
263	!-------------
264
265	DO jk = 3, jpkam1
266	DO ji = 1, jpi
267	z_elem_a( ji, jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
268	z_elem_c( ji, jk ) = - rdt * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
269	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
270	END DO
271	END DO
272
273	DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
274	!++ Surface boundary condition
275	z_elem_a( ji, 2 ) = 0._wp
276	z_elem_c( ji, 2 ) = - rdt * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
277	zcff = rdt * pcd_du(ji, jj)
278	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff
279	u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + 0.5_wp * zcff * ( pssu(ji-1, jj) + pssu(ji,jj) )
280	!++ Top Neumann B.C.
281	!z_elem_a( ji, jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )
282	!z_elem_c( ji, jpka ) = 0._wp
283	!z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
284	!++ Top Dirichlet B.C.
285	z_elem_a( ji, jpka ) = 0._wp
286	z_elem_c( ji, jpka ) = 0._wp
287	z_elem_b( ji, jpka ) = e3t_abl( jpka )
288	u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
289	END DO
290	!!
291	!! Matrix inversion
292	!! ----------------------------------------------------------
293	DO ji = 2, jpi
294	zcff = 1._wp / z_elem_b( ji, 2 )
295	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
296	u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a)
297	END DO
298
299	DO jk = 3, jpka
300	DO ji = 2, jpi
301	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
302	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
303	u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk ,nt_a) &
304	& - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
305	END DO
306	END DO
307
308	DO jk = jpkam1,2,-1
309	DO ji = 2, jpi
310	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)
311	END DO
312	END DO
313
314	!-------------
315	END DO ! end outer loop
316	!-------------
317
318	!
319	! Vertical diffusion for v_abl
320	!-------------
321	DO jj = 2, jpj ! outer loop
322	!-------------
323	!
324	DO jk = 3, jpkam1
325	DO ji = 1, jpi
326	z_elem_a( ji, jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
327	z_elem_c( ji, jk ) = -rdt * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
328	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
329	END DO
330	END DO
331
332	DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
333	!++ Surface boundary condition
334	z_elem_a( ji, 2 ) = 0._wp
335	z_elem_c( ji, 2 ) = - rdt * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
336	zcff = rdt * pcd_du(ji, jj)
337	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff
338	v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + 0.5_wp * zcff * ( pssv(ji, jj) + pssv(ji, jj-1) )
339	!++ Top Neumann B.C.
340	!z_elem_a( ji, jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
341	!z_elem_c( ji, jpka ) = 0._wp
342	!z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
343	!++ Top Dirichlet B.C.
344	z_elem_a( ji, jpka ) = 0._wp
345	z_elem_c( ji, jpka ) = 0._wp
346	z_elem_b( ji, jpka ) = e3t_abl( jpka )
347	v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
348	END DO
349	!!
350	!! Matrix inversion
351	!! ----------------------------------------------------------
352	DO ji = 1, jpi
353	zcff = 1._wp / z_elem_b( ji, 2 )
354	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
355	v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
356	END DO
357
358	DO jk = 3, jpka
359	DO ji = 1, jpi
360	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
361	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
362	v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk ,nt_a) &
363	& - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
364	END DO
365	END DO
366
367	DO jk = jpkam1,2,-1
368	DO ji = 1, jpi
369	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)
370	END DO
371	END DO
372	!
373	!-------------
374	END DO ! end outer loop
375	!-------------
376
377	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
378	! ! 5 *** Apply nudging on the dynamics and the tracers
379	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
380	z_cft(:,:,:) = 0._wp
381
382	IF( nn_dyn_restore > 0 ) THEN
383	!-------------
384	DO jk = 2, jpka ! outer loop
385	!-------------
386	DO jj = 2, jpj
387	DO ji = 2, jpi
388	zcff1 = pblh( ji, jj )
389	zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
390	zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
391	zmsk = msk_abl(ji,jj)
392	zcff2 = jp_alp3_dyn * zsig*3 + jp_alp2_dyn zsig**2 &
393	& + jp_alp1_dyn * zsig + jp_alp0_dyn
394	zcff = (1._wp-zmsk) + zmsk * rdt * zcff2 ! zcff = 1 for masked points
395
396	u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) &
397	& + zcff * pu_dta( ji, jj, jk )
398	v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) &
399	& + zcff * pv_dta( ji, jj, jk )
400	END DO
401	END DO
402	!-------------
403	END DO ! end outer loop
404	!-------------
405	END IF
406
407	!-------------
408	DO jk = 2, jpka ! outer loop
409	!-------------
410	DO jj = 1,jpj
411	DO ji = 1,jpi
412	zcff1 = pblh( ji, jj )
413	zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
414	zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
415	zmsk = msk_abl(ji,jj)
416	zcff2 = jp_alp3_tra * zsig*3 + jp_alp2_tra zsig**2 &
417	& + jp_alp1_tra * zsig + jp_alp0_tra
418	zcff = (1._wp-zmsk) + zmsk * rdt * zcff2 ! zcff = 1 for masked points
419	tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) &
420	& + zcff * pt_dta( ji, jj, jk )
421
422	tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) &
423	& + zcff * pq_dta( ji, jj, jk )
424
425	z_cft( ji, jj, jk ) = zcff
426	END DO
427	END DO
428	!-------------
429	END DO ! end outer loop
430	!-------------
431	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
432	! ! 6 *** MPI exchanges
433	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
434	!
435	CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), kfillmode = jpfillnothing )
436	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...
437	!
438	CALL iom_put ( 'u_abl', u_abl(:,:,2:jpka,nt_a ))
439	CALL iom_put ( 'v_abl', v_abl(:,:,2:jpka,nt_a ))
440	CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta))
441	CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa))
442	! CALL iom_put ( 'coeft', z_cft(:,:,2:jpka ))
443	CALL iom_put ( 'pblh', pblh(:,: ))
444	!
445	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
446	! ! 7 *** Finalize flux computation
447	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
448
449	DO jj = 1, jpj
450	DO ji = 1, jpi
451	ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta )
452	zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa )
453	zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ...
454	& ( R_dryztemp ( 1._wp + rctv0*zhumi ) )
455	psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
456	pevp ( ji, jj ) = rn_efacMAX( 0._wp, zcff pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
457	zrhoa( ji, jj ) = zcff
458	END DO
459	END DO
460
461	DO jj = 2, jpj
462	DO ji = 2, jpi ! vect. opt.
463	zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji ,jj) + pssu(ji-1,jj) )
464	zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj ) + pssv(ji,jj-1) )
465	END DO
466	END DO
467	!
468	CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
469	!
470	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
471	DO jj = 1, jpj
472	DO ji = 1, jpi
473	zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) &
474	& + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) * msk_abl(ji,jj)
475	zztmp = zrhoa(ji,jj) * pcd_du(ji,jj)
476
477	pwndm (ji,jj) = zcff
478	ptaum (ji,jj) = zztmp * zcff
479	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
480	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
481	END DO
482	END DO
483	! ... utau, vtau at U- and V_points, resp.
484	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
485	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
486	DO jj = 2, jpjm1
487	DO ji = 2, jpim1
488	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
489	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
490	ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
491	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
492	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
493	ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
494	END DO
495	END DO
496	!
497	CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
498
499	CALL iom_put( "taum_oce", ptaum )
500
501	IF(ln_ctl) THEN
502	CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' )
503	CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ' )
504	CALL prt_ctl( tab2d_2=ptauj , clinfo2= 'vtau : ' )
505	ENDIF
506	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
507	! ! 8 *** Swap time indices for the next timestep
508	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
509	nt_n = 1 + MOD( kt , 2)
510	nt_a = 1 + MOD( kt+1, 2)
511	!
512	!---------------------------------------------------------------------------------------------------
513	END SUBROUTINE abl_stp
514	!===================================================================================================
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532	!===================================================================================================
533	SUBROUTINE abl_zdf_tke( )
534	!---------------------------------------------------------------------------------------------------
535
536	!!----------------------------------------------------------------------
537	!! * ROUTINE abl_zdf_tke *
538	!!
539	!! ** Purpose : Time-step Turbulente Kinetic Energy (TKE) equation
540	!!
541	!! ** Method : - source term due to shear
542	!! - source term due to stratification
543	!! - resolution of the TKE equation by inverting
544	!! a tridiagonal linear system
545	!!
546	!! ** Action : - en : now turbulent kinetic energy)
547	!! - avmu, avmv : production of TKE by shear at u and v-points
548	!! (= Kz dz[Ub] * dz[Un] )
549	!! ---------------------------------------------------------------------
550	INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
551	INTEGER, DIMENSION(1:jpi ) :: ikbl
552	REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
553	REAL(wp) :: zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
554	REAL(wp) :: zdU2,zdV2
555	REAL(wp) :: zwndi,zwndj
556	REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2
557	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2
558	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF
559	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
560	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
561	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
562	LOGICAL :: ln_Patankar = .FALSE.
563	LOGICAL :: ln_dumpvar = .FALSE.
564	LOGICAL , DIMENSION(1:jpi ) :: ln_foundl
565	!
566	tind = nt_n
567	ziRic = 1._wp / rn_Ric
568	! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
569	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
570	! ! Advance TKE equation to time n+1
571	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
572	!-------------
573	DO jj = 1, jpj ! outer loop
574	!-------------
575	!
576	! Compute vertical shear
577	DO jk = 2, jpkam1
578	DO ji = 1,jpi
579	zcff = 1.0_wp / e3w_abl( jk )**2
580	zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
581	zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
582	zsh2(ji,jk) = zdU+zdV
583	END DO
584	END DO
585	!
586	! Compute brunt-vaisala frequency
587	DO jk = 2, jpkam1
588	DO ji = 1,jpi
589	zcff = grav * itvref / e3w_abl( jk )
590	zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta)
591	zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) &
592	& - tq_abl( ji, jj, jk , tind, jp_ta) * tq_abl( ji, jj, jk , tind, jp_qa)
593	zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi)
594	END DO
595	END DO
596	!
597	! Terms for the tridiagonal problem
598	DO jk = 2, jpkam1
599	DO ji = 1,jpi
600	zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point
601	zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
602	zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
603
604	z_elem_a( ji, jk ) = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal
605	z_elem_c( ji, jk ) = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal
606	IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE
607	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
608	& + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal
609	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) ) ! right-hand-side
610	ELSE
611	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
612	& + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal
613	& - e3w_abl(jk) * rdt * zbuoy
614	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * zshear ) ! right-hand-side
615	END IF
616	END DO
617	END DO
618
619	DO ji = 1,jpi ! vector opt.
620	zesrf = MAX( 4.63_wp * ustar2(ji,jj), tke_min )
621	zetop = tke_min
622	z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp
623	z_elem_a ( ji, jpka ) = 0._wp; z_elem_c ( ji, jpka ) = 0._wp; z_elem_b ( ji, jpka ) = 1._wp
624	tke_abl( ji, jj, 1, nt_a ) = zesrf
625	tke_abl( ji, jj, jpka, nt_a ) = zetop
626	zbn2(ji,jj, 1) = zbn2( ji,jj, 2)
627	zsh2(ji, 1) = zsh2( ji, 2)
628	zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)
629	zsh2(ji, jpka) = zsh2( ji , jpkam1)
630	END DO
631	!!
632	!! Matrix inversion
633	!! ----------------------------------------------------------
634	DO ji = 1,jpi
635	zcff = 1._wp / z_elem_b( ji, 1 )
636	zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 )
637	tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a )
638	END DO
639
640	DO jk = 2, jpka
641	DO ji = 1,jpi
642	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
643	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
644	tke_abl(ji,jj,jk,nt_a) = zcff * ( tke_abl(ji,jj,jk ,nt_a) &
645	& - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
646	END DO
647	END DO
648
649	DO jk = jpkam1,1,-1
650	DO ji = 1,jpi
651	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)
652	END DO
653	END DO
654
655	!!FL should not be needed because of Patankar procedure
656	tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
657
658	!!
659	!! Diagnose PBL height
660	!! ----------------------------------------------------------
661
662
663	!
664	! arrays zRH, zFC and zCF are available at this point
665	! and zFC(:, 1 ) = 0.
666	! diagnose PBL height based on zsh2 and zbn2
667	zFC ( : ,1) = 0._wp
668	ikbl( 1:jpi ) = 0
669
670	DO jk = 2,jpka
671	DO ji = 1, jpi
672	zcff = ghw_abl( jk-1 )
673	zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
674	zcff = ghw_abl( jk )
675	zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*( &
676	zcff2 * ( zsh2( ji, jk ) - ziRic * zbn2( ji, jj, jk ) &
677	- rn_Cek * ( ff_t( ji, jj ) * ff_t( ji, jj ) ) ) &
678	+ zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
679	- rn_Cek * ( ff_t( ji, jj ) * ff_t( ji, jj ) ) ) &
680	& )
681	IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
682	END DO
683	END DO
684	!
685	! finalize the computation of the PBL height
686	DO ji = 1, jpi
687	jk = ikbl(ji)
688	IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
689	pblh( ji, jj ) = ( ghw_abl( jk-1 ) * zFC( ji, jk ) &
690	& - ghw_abl( jk ) * zFC( ji, jk-1 ) &
691	& ) / ( zFC( ji, jk ) - zFC( ji, jk-1 ) )
692	ELSE IF( jk==2 ) THEN
693	pblh( ji, jj ) = ghw_abl(2 )
694	ELSE
695	pblh( ji, jj ) = ghw_abl(jpka)
696	END IF
697	END DO
698	! Optional : could add pblh smoothing if pblh is noisy horizontally ...
699	!
700	!-------------
701	END DO
702	!-------------
703	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
704	! ! Diagnostic mixing length computation
705	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
706	!
707	SELECT CASE ( nn_amxl )
708	!
709	CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom
710	# define zlup zRH
711	# define zldw zFC
712	DO jj = 1, jpj ! outer loop
713	!
714	DO ji = 1, jpi
715	mxl_abl ( ji, jj, 1 ) = 0._wp
716	mxl_abl ( ji, jj, jpka ) = mxl_min
717	zldw( ji, 1 ) = 0._wp
718	zlup( ji, jpka ) = 0._wp
719	END DO
720	!
721	DO jk = 2, jpkam1
722	DO ji = 1, jpi
723	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
724	mxl_abl( ji, jj, jk ) = MAX( mxl_min, &
725	& SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
726	END DO
727	END DO
728	!
729	! Limit mxl
730	DO jk = jpkam1,1,-1
731	DO ji = 1, jpi
732	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
733	END DO
734	END DO
735	!
736	DO jk = 2, jpka
737	DO ji = 1, jpi
738	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )
739	END DO
740	END DO
741	!
742	DO jk = 1, jpka
743	DO ji = 1, jpi
744	mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
745	END DO
746	END DO
747	!
748	END DO
749	# undef zlup
750	# undef zldw
751	!
752	!
753	CASE ( 1 ) ! length-scale computed as the distance to the PBL height
754	DO jj = 1,jpj ! outer loop
755	!
756	DO ji = 1, jpi ! vector opt.
757	zcff = 1._wp / pblh( ji, jj ) ! inverse of hbl
758	DO jk = 1, jpka
759	zsig = MIN( zcff * ghw_abl( jk ), 1. )
760	zcff1 = pblh( ji, jj )
761	mxl_abl( ji, jj, jk ) = mxl_min &
762	& + zsig * ( amx1zcff1 + bmx1mxl_min ) &
763	& + zsig * zsig * ( amx2zcff1 + bmx2mxl_min ) &
764	& + zsig*3 ( amx3zcff1 + bmx3mxl_min ) &
765	& + zsig*4 ( amx4zcff1 + bmx4mxl_min ) &
766	& + zsig*5 ( amx5zcff1 + bmx5mxl_min )
767	END DO
768	END DO
769	!
770	END DO
771	!
772	CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale
773	!
774	# define zlup zRH
775	# define zldw zFC
776	! zCF is used for matrix inversion
777	!
778	DO jj = 1, jpj ! outer loop
779
780	DO ji = 1, jpi
781	zlup( ji, 1 ) = mxl_min
782	zldw( ji, 1 ) = mxl_min
783	zlup( ji, jpka ) = mxl_min
784	zldw( ji, jpka ) = mxl_min
785	END DO
786
787	DO jk = 2,jpka-1
788	DO ji = 1, jpi
789	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
790	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
791	END DO
792	END DO
793	!!
794	!! BL89 search for lup
795	!! ----------------------------------------------------------
796	DO jk=2,jpka-1
797	!
798	DO ji = 1, jpi
799	zCF(ji,1:jpka) = 0._wp
800	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
801	ln_foundl(ji ) = .false.
802	END DO
803	!
804	DO jkup=jk+1,jpka-1
805	DO ji = 1, jpi
806	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
807	& ( zbn2(ji,jj,jkup )*(ghw_abl(jkup )-ghw_abl(jk)) &
808	& + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
809	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
810	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
811	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
812	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
813	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
814	zlup(ji,jk) = zcff
815	ln_foundl(ji) = .true.
816	END IF
817	END DO
818	END DO
819	!
820	END DO
821	!!
822	!! BL89 search for ldwn
823	!! ----------------------------------------------------------
824	DO jk=2,jpka-1
825	!
826	DO ji = 1, jpi
827	zCF(ji,1:jpka) = 0._wp
828	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
829	ln_foundl(ji ) = .false.
830	END DO
831	!
832	DO jkdwn=jk-1,1,-1
833	DO ji = 1, jpi
834	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
835	& * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
836	+ zbn2(ji,jj,jkdwn )*(ghw_abl(jk)-ghw_abl(jkdwn )) )
837	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
838	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
839	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
840	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
841	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
842	zldw(ji,jk) = zcff
843	ln_foundl(ji) = .true.
844	END IF
845	END DO
846	END DO
847	!
848	END DO
849
850	DO jk = 1, jpka
851	DO ji = 1, jpi
852	mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )
853	END DO
854	END DO
855
856	END DO
857	# undef zlup
858	# undef zldw
859	!
860	END SELECT
861	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
862	! ! Finalize the computation of turbulent visc./diff.
863	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
864
865	!-------------
866	DO jj = 1, jpj ! outer loop
867	!-------------
868	DO jk = 1, jpka
869	DO ji = 1, jpi ! vector opt.
870	zcff = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk ) &
871	& * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )
872	zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z)
873	zcff = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
874	!!FL: MAX function probably useless because of the definition of mxl_min
875	Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak )
876	Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak )
877	END DO
878	END DO
879	!-------------
880	END DO
881	!-------------
882
883	!---------------------------------------------------------------------------------------------------
884	END SUBROUTINE abl_zdf_tke
885	!===================================================================================================
886
887
888	!!======================================================================
889	END MODULE ablmod

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

Download in other formats: