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/2021/ticket2632_r14588_theta_sbcblk/src/ABL – NEMO

Context Navigation

source: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ABL/ablmod.F90 @ 15548

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

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

Download in other formats: