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

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

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

Context Navigation

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

Last change on this file since 13176 was 13176, checked in by smasson, 4 years ago
Extra_Halo: rewrite prtctl, supress nn_print, see #2366
File size: 52.2 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_?fac
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
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,pssq_ice,pcd_du_ice &
47	& , psen_ice, pevp_ice, pwndm_ice, pfrac_oce &
48	& , ptaui_ice, ptauj_ice &
49	#endif
50	& )
51	!---------------------------------------------------------------------------------------------------
52
53	!!---------------------------------------------------------------------
54	!! * ROUTINE abl_stp *
55	!!
56	!! ** Purpose : Time-integration of the ABL model
57	!!
58	!! ** Method : Compute atmospheric variables : vertical turbulence
59	!! + Coriolis term + newtonian relaxation
60	!!
61	!! ** Action : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
62	!! - Advance tracers to time n+1 (Euler backward scheme)
63	!! - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
64	!! - Advance u,v to time n+1 (Euler backward scheme)
65	!! - Apply newtonian relaxation on the dynamics and the tracers
66	!! - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
67	!!
68	!!----------------------------------------------------------------------
69	INTEGER , INTENT(in ) :: kt ! time step index
70	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psst ! sea-surface temperature [Celsius]
71	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu ! sea-surface u (U-point)
72	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point)
73	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq ! sea-surface humidity
74	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pu_dta ! large-scale windi
75	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pv_dta ! large-scale windj
76	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgu_dta ! large-scale hpgi
77	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgv_dta ! large-scale hpgj
78	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pt_dta ! large-scale pot. temp.
79	REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pq_dta ! large-scale humidity
80	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp_dta ! sea-level pressure
81	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du ! Cd x Du (T-point)
82	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du
83	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du
84	REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! \|\|uwnd\|\|
85	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux
86	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy
87	REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! \|\|tau\|\|
88	!
89	#if defined key_si3
90	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptm_su ! ice-surface temperature [K]
91	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu_ice ! ice-surface u (U-point)
92	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point)
93	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity
94	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du_ice ! Cd x Du over ice (T-point)
95	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psen_ice ! Ch x Du over ice (T-point)
96	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pevp_ice ! Ce x Du over ice (T-point)
97	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! \|\|uwnd - uice\|\|
98	!REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pfrac_oce !!GS: out useless ?
99	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce !
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	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: z_cft !--FL--to be removed after the test phase
107	!
108	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
109	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
110	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
111	!
112	INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices
113	REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
114	REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
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	!! Compute ustar squared as Cd \|\| Uatm-Uoce \|\|^2
126	!! needed for surface boundary condition of TKE
127	!! pwndm contains \| U10m - U_oce \| (see blk_oce_1 in sbcblk)
128	DO_2D_11_11
129	zzoce = pCd_du (ji,jj) * pwndm (ji,jj)
130	#if defined key_si3
131	zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)
132	ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice
133	#else
134	ustar2(ji,jj) = zzoce
135	#endif
136	END_2D
137	!
138	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
139	! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
140	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
141
142	CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
143
144	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
145	! ! 2 *** Advance tracers to time n+1
146	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
147
148	!-------------
149	DO jj = 1, jpj ! outer loop !--> tq_abl computed on (1:jpi) x (1:jpj)
150	!-------------
151	! Compute matrix elements for interior points
152	DO jk = 3, jpkam1
153	DO ji = 1, jpi ! vector opt.
154	z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
155	z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
156	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
157	END DO
158	END DO
159	! Boundary conditions
160	DO ji = 1, jpi ! vector opt.
161	! Neumann at the bottom
162	z_elem_a( ji, 2 ) = 0._wp
163	z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )
164	! Homogeneous Neumann at the top
165	z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
166	z_elem_c( ji, jpka ) = 0._wp
167	z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
168	END DO
169
170	DO jtra = 1,jptq ! loop on active tracers
171
172	DO jk = 3, jpkam1
173	DO ji = 1,jpi
174	tq_abl ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl ( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side
175	END DO
176	END DO
177
178	IF(jtra == jp_ta) THEN
179	DO ji = 1,jpi ! boundary conditions for temperature
180	zztmp1 = psen(ji, jj)
181	zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
182	#if defined key_si3
183	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
184	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
185	#endif
186	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
187	tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2
188	tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
189	END DO
190	ELSE
191	DO ji = 1,jpi ! boundary conditions for humidity
192	zztmp1 = pevp(ji, jj)
193	zztmp2 = pevp(ji, jj) * pssq(ji, jj)
194	#if defined key_si3
195	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj)
196	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)
197	#endif
198	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
199	tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2
200	tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
201	END DO
202	END IF
203	!!
204	!! Matrix inversion
205	!! ----------------------------------------------------------
206	DO ji = 1,jpi
207	zcff = 1._wp / z_elem_b( ji, 2 )
208	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
209	tq_abl(ji,jj,2,nt_a,jtra) = zcff * tq_abl(ji,jj,2,nt_a,jtra)
210	END DO
211
212	DO jk = 3, jpka
213	DO ji = 1,jpi
214	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
215	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
216	tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) &
217	& - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
218	END DO
219	END DO
220	!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
221	DO jk = jpkam1,2,-1
222	DO ji = 1,jpi
223	tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) + &
224	& zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
225	END DO
226	END DO
227
228	END DO !<-- loop on tracers
229	!!
230	!-------------
231	END DO ! end outer loop
232	!-------------
233
234
235	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
236	! ! 3 *** Compute Coriolis term with geostrophic guide
237	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
238	!-------------
239	DO jk = 2, jpka ! outer loop
240	!-------------
241	!
242	! Advance u_abl & v_abl to time n+1
243	DO_2D_11_11
244	zcff = ( fft_abl(ji,jj) * rDt_abl )( fft_abl(ji,jj) rDt_abl ) ! (f dt)**2
245
246	u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
247	& (1._wp-gamma_Cor(1._wp-gamma_Cor)zcff)*u_abl( ji, jj, jk, nt_n ) &
248	& + rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) &
249	& / (1._wp + gamma_Corgamma_Corzcff)
250
251	v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
252	& (1._wp-gamma_Cor(1._wp-gamma_Cor)zcff)*v_abl( ji, jj, jk, nt_n ) &
253	& - rDt_abl * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) &
254	& / (1._wp + gamma_Corgamma_Corzcff)
255	END_2D
256	!
257	!-------------
258	END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
259	!-------------
260	!
261	IF( ln_geos_winds ) THEN
262	DO jj = 1, jpj ! outer loop
263	DO jk = 1, jpka
264	DO ji = 1, jpi
265	u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) &
266	& - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk)
267	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) &
268	& + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk)
269	END DO
270	END DO
271	END DO
272	END IF
273	!-------------
274	!
275	IF( ln_hpgls_frc ) 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 ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
280	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
281	ENDDO
282	ENDDO
283	ENDDO
284	END IF
285
286	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
287	! ! 4 *** Advance u,v to time n+1
288	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
289	!
290	! Vertical diffusion for u_abl
291	!-------------
292	DO jj = 1, jpj ! outer loop
293	!-------------
294
295	DO jk = 3, jpkam1
296	DO ji = 1, jpi
297	z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
298	z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
299	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
300	END DO
301	END DO
302
303	DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
304	!++ Surface boundary condition
305	z_elem_a( ji, 2 ) = 0._wp
306	z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
307	!
308	zztmp1 = pcd_du(ji, jj)
309	zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )
310	#if defined key_si3
311	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
312	zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) )
313	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
314	#endif
315	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
316	u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
317
318	!++ Top Neumann B.C.
319	!z_elem_a( ji, jpka ) = - 0.5_wp * rDt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )
320	!z_elem_c( ji, jpka ) = 0._wp
321	!z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
322	!++ Top Dirichlet B.C.
323	z_elem_a( ji, jpka ) = 0._wp
324	z_elem_c( ji, jpka ) = 0._wp
325	z_elem_b( ji, jpka ) = e3t_abl( jpka )
326	u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
327	END DO
328	!!
329	!! Matrix inversion
330	!! ----------------------------------------------------------
331	DO ji = 2, jpi
332	zcff = 1._wp / z_elem_b( ji, 2 )
333	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
334	u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a)
335	END DO
336
337	DO jk = 3, jpka
338	DO ji = 2, jpi
339	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
340	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
341	u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk ,nt_a) &
342	& - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
343	END DO
344	END DO
345
346	DO jk = jpkam1,2,-1
347	DO ji = 2, jpi
348	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)
349	END DO
350	END DO
351
352	!-------------
353	END DO ! end outer loop
354	!-------------
355
356	!
357	! Vertical diffusion for v_abl
358	!-------------
359	DO jj = 2, jpj ! outer loop
360	!-------------
361	!
362	DO jk = 3, jpkam1
363	DO ji = 1, jpi
364	z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
365	z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
366	z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
367	END DO
368	END DO
369
370	DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
371	!++ Surface boundary condition
372	z_elem_a( ji, 2 ) = 0._wp
373	z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
374	!
375	zztmp1 = pcd_du(ji, jj)
376	zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )
377	#if defined key_si3
378	zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
379	zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) )
380	zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
381	#endif
382	z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
383	v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
384	!++ Top Neumann B.C.
385	!z_elem_a( ji, jpka ) = -rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
386	!z_elem_c( ji, jpka ) = 0._wp
387	!z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
388	!++ Top Dirichlet B.C.
389	z_elem_a( ji, jpka ) = 0._wp
390	z_elem_c( ji, jpka ) = 0._wp
391	z_elem_b( ji, jpka ) = e3t_abl( jpka )
392	v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
393	END DO
394	!!
395	!! Matrix inversion
396	!! ----------------------------------------------------------
397	DO ji = 1, jpi
398	zcff = 1._wp / z_elem_b( ji, 2 )
399	zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
400	v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
401	END DO
402
403	DO jk = 3, jpka
404	DO ji = 1, jpi
405	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
406	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
407	v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk ,nt_a) &
408	& - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
409	END DO
410	END DO
411
412	DO jk = jpkam1,2,-1
413	DO ji = 1, jpi
414	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)
415	END DO
416	END DO
417	!
418	!-------------
419	END DO ! end outer loop
420	!-------------
421
422	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
423	! ! 5 *** Apply nudging on the dynamics and the tracers
424	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
425	z_cft(:,:,:) = 0._wp
426
427	IF( nn_dyn_restore > 0 ) THEN
428	!-------------
429	DO jk = 2, jpka ! outer loop
430	!-------------
431	DO_2D_01_01
432	zcff1 = pblh( ji, jj )
433	zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
434	zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
435	zmsk = msk_abl(ji,jj)
436	zcff2 = jp_alp3_dyn * zsig*3 + jp_alp2_dyn zsig**2 &
437	& + jp_alp1_dyn * zsig + jp_alp0_dyn
438	zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
439	! rn_Dt = rDt_abl / nn_fsbc
440	zcff = zcff * rest_eq(ji,jj)
441	z_cft( ji, jj, jk ) = zcff
442	u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) &
443	& + zcff * pu_dta( ji, jj, jk )
444	v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) &
445	& + zcff * pv_dta( ji, jj, jk )
446	END_2D
447	!-------------
448	END DO ! end outer loop
449	!-------------
450	END IF
451
452	!-------------
453	DO jk = 2, jpka ! outer loop
454	!-------------
455	DO_2D_11_11
456	zcff1 = pblh( ji, jj )
457	zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
458	zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
459	zmsk = msk_abl(ji,jj)
460	zcff2 = jp_alp3_tra * zsig*3 + jp_alp2_tra zsig**2 &
461	& + jp_alp1_tra * zsig + jp_alp0_tra
462	zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
463	! rn_Dt = rDt_abl / nn_fsbc
464	!z_cft( ji, jj, jk ) = zcff
465	tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) &
466	& + zcff * pt_dta( ji, jj, jk )
467
468	tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) &
469	& + zcff * pq_dta( ji, jj, jk )
470
471	END_2D
472	!-------------
473	END DO ! end outer loop
474	!-------------
475	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
476	! ! 6 *** MPI exchanges
477	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
478	!
479	CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a ), 'T', -1. )
480	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...
481	!
482	! first ABL level
483	IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) )
484	IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) )
485	IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl", tq_abl(:,:,2,nt_a,jp_ta) )
486	IF ( iom_use("qz1_abl") ) CALL iom_put ( "qz1_abl", tq_abl(:,:,2,nt_a,jp_qa) )
487	IF ( iom_use("uz1_dta") ) CALL iom_put ( "uz1_dta", pu_dta(:,:,2 ) )
488	IF ( iom_use("vz1_dta") ) CALL iom_put ( "vz1_dta", pv_dta(:,:,2 ) )
489	IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta", pt_dta(:,:,2 ) )
490	IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta", pq_dta(:,:,2 ) )
491	! all ABL levels
492	IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) )
493	IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) )
494	IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) )
495	IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) )
496	IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a ) )
497	IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka ) )
498	IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka ) )
499	IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka ) )
500	IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh" , pblh(:,: ) )
501	! debug (to be removed)
502	IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta", pu_dta(:,:,2:jpka) )
503	IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta", pv_dta(:,:,2:jpka) )
504	IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta", pt_dta(:,:,2:jpka) )
505	IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta", pq_dta(:,:,2:jpka) )
506	IF ( iom_use("coeft") ) CALL iom_put ( "coeft", z_cft(:,:,2:jpka) )
507	IF( ln_geos_winds ) THEN
508	IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2 ) )
509	IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2 ) )
510	END IF
511	IF( ln_hpgls_frc ) THEN
512	IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
513	IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
514	END IF
515	!
516	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
517	! ! 7 *** Finalize flux computation
518	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
519
520	DO_2D_11_11
521	ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta )
522	zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa )
523	!zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ...
524	! & ( R_dryztemp ( 1._wp + rctv0*zhumi ) )
525	zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
526	psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
527	pevp ( ji, jj ) = rn_efacMAX( 0._wp, zcff pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
528	rhoa( ji, jj ) = zcff
529	END_2D
530
531	DO_2D_01_01
532	zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji ,jj) + pssu(ji-1,jj) )
533	zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj ) + pssv(ji,jj-1) )
534	END_2D
535	!
536	CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
537	!
538	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
539	DO_2D_11_11
540	zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) &
541	& + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj)
542	zztmp = rhoa(ji,jj) * pcd_du(ji,jj)
543
544	pwndm (ji,jj) = zcff
545	ptaum (ji,jj) = zztmp * zcff
546	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
547	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
548	END_2D
549	! ... utau, vtau at U- and V_points, resp.
550	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
551	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
552	DO_2D_00_00
553	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
554	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
555	ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
556	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
557	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
558	ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
559	END_2D
560	!
561	CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
562
563	CALL iom_put( "taum_oce", ptaum )
564
565	IF(sn_cfctl%l_prtctl) THEN
566	CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' )
567	CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', &
568	& tab2d_2=ptauj , clinfo2= 'vtau : ' )
569	ENDIF
570
571	#if defined key_si3
572	! ------------------------------------------------------------ !
573	! Wind stress relative to the moving ice ( U10m - U_ice ) !
574	! ------------------------------------------------------------ !
575	DO_2D_00_00
576
577	zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
578	zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
579
580	ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &
581	& + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &
582	& * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
583	ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &
584	& + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &
585	& * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
586	END_2D
587	CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
588	!
589	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
590	& , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
591	#endif
592	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
593	! ! 8 *** Swap time indices for the next timestep
594	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
595	nt_n = 1 + MOD( nt_n, 2)
596	nt_a = 1 + MOD( nt_a, 2)
597	!
598	!---------------------------------------------------------------------------------------------------
599	END SUBROUTINE abl_stp
600	!===================================================================================================
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618	!===================================================================================================
619	SUBROUTINE abl_zdf_tke( )
620	!---------------------------------------------------------------------------------------------------
621
622	!!----------------------------------------------------------------------
623	!! * ROUTINE abl_zdf_tke *
624	!!
625	!! ** Purpose : Time-step Turbulente Kinetic Energy (TKE) equation
626	!!
627	!! ** Method : - source term due to shear
628	!! - source term due to stratification
629	!! - resolution of the TKE equation by inverting
630	!! a tridiagonal linear system
631	!!
632	!! ** Action : - en : now turbulent kinetic energy)
633	!! - avmu, avmv : production of TKE by shear at u and v-points
634	!! (= Kz dz[Ub] * dz[Un] )
635	!! ---------------------------------------------------------------------
636	INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
637	INTEGER, DIMENSION(1:jpi ) :: ikbl
638	REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
639	REAL(wp) :: zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
640	REAL(wp) :: zdU2,zdV2
641	REAL(wp) :: zwndi,zwndj
642	REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2
643	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2
644	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF
645	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
646	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
647	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
648	LOGICAL :: ln_Patankar = .FALSE.
649	LOGICAL :: ln_dumpvar = .FALSE.
650	LOGICAL , DIMENSION(1:jpi ) :: ln_foundl
651	!
652	tind = nt_n
653	ziRic = 1._wp / rn_Ric
654	! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
655	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
656	! ! Advance TKE equation to time n+1
657	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
658	!-------------
659	DO jj = 1, jpj ! outer loop
660	!-------------
661	!
662	! Compute vertical shear
663	DO jk = 2, jpkam1
664	DO ji = 1,jpi
665	zcff = 1.0_wp / e3w_abl( jk )**2
666	zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
667	zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
668	zsh2(ji,jk) = zdU+zdV
669	END DO
670	END DO
671	!
672	! Compute brunt-vaisala frequency
673	DO jk = 2, jpkam1
674	DO ji = 1,jpi
675	zcff = grav * itvref / e3w_abl( jk )
676	zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta)
677	zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) &
678	& - tq_abl( ji, jj, jk , tind, jp_ta) * tq_abl( ji, jj, jk , tind, jp_qa)
679	zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi)
680	END DO
681	END DO
682	!
683	! Terms for the tridiagonal problem
684	DO jk = 2, jpkam1
685	DO ji = 1,jpi
686	zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point
687	zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
688	zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
689
690	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
691	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
692	IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE
693	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
694	& + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal
695	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side
696	ELSE
697	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
698	& + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal
699	& - e3w_abl(jk) * rDt_abl * zbuoy
700	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side
701	END IF
702	END DO
703	END DO
704
705	DO ji = 1,jpi ! vector opt.
706	zesrf = MAX( 4.63_wp * ustar2(ji,jj), tke_min )
707	zetop = tke_min
708	z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp
709	z_elem_a ( ji, jpka ) = 0._wp; z_elem_c ( ji, jpka ) = 0._wp; z_elem_b ( ji, jpka ) = 1._wp
710	tke_abl( ji, jj, 1, nt_a ) = zesrf
711	tke_abl( ji, jj, jpka, nt_a ) = zetop
712	zbn2(ji,jj, 1) = zbn2( ji,jj, 2)
713	zsh2(ji, 1) = zsh2( ji, 2)
714	zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)
715	zsh2(ji, jpka) = zsh2( ji , jpkam1)
716	END DO
717	!!
718	!! Matrix inversion
719	!! ----------------------------------------------------------
720	DO ji = 1,jpi
721	zcff = 1._wp / z_elem_b( ji, 1 )
722	zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 )
723	tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a )
724	END DO
725
726	DO jk = 2, jpka
727	DO ji = 1,jpi
728	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
729	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
730	tke_abl(ji,jj,jk,nt_a) = zcff * ( tke_abl(ji,jj,jk ,nt_a) &
731	& - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
732	END DO
733	END DO
734
735	DO jk = jpkam1,1,-1
736	DO ji = 1,jpi
737	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)
738	END DO
739	END DO
740
741	!!FL should not be needed because of Patankar procedure
742	tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
743
744	!!
745	!! Diagnose PBL height
746	!! ----------------------------------------------------------
747
748
749	!
750	! arrays zRH, zFC and zCF are available at this point
751	! and zFC(:, 1 ) = 0.
752	! diagnose PBL height based on zsh2 and zbn2
753	zFC ( : ,1) = 0._wp
754	ikbl( 1:jpi ) = 0
755
756	DO jk = 2,jpka
757	DO ji = 1, jpi
758	zcff = ghw_abl( jk-1 )
759	zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
760	zcff = ghw_abl( jk )
761	zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
762	zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*( &
763	zcff2 * ( zsh2( ji, jk ) - ziRic * zbn2( ji, jj, jk ) &
764	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
765	+ zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
766	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
767	& )
768	IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
769	END DO
770	END DO
771	!
772	! finalize the computation of the PBL height
773	DO ji = 1, jpi
774	jk = ikbl(ji)
775	IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
776	pblh( ji, jj ) = ( ghw_abl( jk-1 ) * zFC( ji, jk ) &
777	& - ghw_abl( jk ) * zFC( ji, jk-1 ) &
778	& ) / ( zFC( ji, jk ) - zFC( ji, jk-1 ) )
779	ELSE IF( jk==2 ) THEN
780	pblh( ji, jj ) = ghw_abl(2 )
781	ELSE
782	pblh( ji, jj ) = ghw_abl(jpka)
783	END IF
784	END DO
785	!-------------
786	END DO
787	!-------------
788	!
789	! Optional : could add pblh smoothing if pblh is noisy horizontally ...
790	IF(ln_smth_pblh) THEN
791	CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
792	CALL smooth_pblh( pblh, msk_abl )
793	CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
794	ENDIF
795	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
796	! ! Diagnostic mixing length computation
797	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
798	!
799	SELECT CASE ( nn_amxl )
800	!
801	CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom
802	# define zlup zRH
803	# define zldw zFC
804	DO jj = 1, jpj ! outer loop
805	!
806	DO ji = 1, jpi
807	mxl_abl ( ji, jj, 1 ) = 0._wp
808	mxl_abl ( ji, jj, jpka ) = mxl_min
809	zldw( ji, 1 ) = 0._wp
810	zlup( ji, jpka ) = 0._wp
811	END DO
812	!
813	DO jk = 2, jpkam1
814	DO ji = 1, jpi
815	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
816	mxl_abl( ji, jj, jk ) = MAX( mxl_min, &
817	& SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
818	END DO
819	END DO
820	!
821	! Limit mxl
822	DO jk = jpkam1,1,-1
823	DO ji = 1, jpi
824	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
825	END DO
826	END DO
827	!
828	DO jk = 2, jpka
829	DO ji = 1, jpi
830	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )
831	END DO
832	END DO
833	!
834	DO jk = 1, jpka
835	DO ji = 1, jpi
836	mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
837	END DO
838	END DO
839	!
840	END DO
841	# undef zlup
842	# undef zldw
843	!
844	!
845	CASE ( 1 ) ! length-scale computed as the distance to the PBL height
846	DO jj = 1,jpj ! outer loop
847	!
848	DO ji = 1, jpi ! vector opt.
849	zcff = 1._wp / pblh( ji, jj ) ! inverse of hbl
850	DO jk = 1, jpka
851	zsig = MIN( zcff * ghw_abl( jk ), 1. )
852	zcff1 = pblh( ji, jj )
853	mxl_abl( ji, jj, jk ) = mxl_min &
854	& + zsig * ( amx1zcff1 + bmx1mxl_min ) &
855	& + zsig * zsig * ( amx2zcff1 + bmx2mxl_min ) &
856	& + zsig*3 ( amx3zcff1 + bmx3mxl_min ) &
857	& + zsig*4 ( amx4zcff1 + bmx4mxl_min ) &
858	& + zsig*5 ( amx5zcff1 + bmx5mxl_min )
859	END DO
860	END DO
861	!
862	END DO
863	!
864	CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale
865	!
866	# define zlup zRH
867	# define zldw zFC
868	! zCF is used for matrix inversion
869	!
870	DO jj = 1, jpj ! outer loop
871
872	DO ji = 1, jpi
873	zlup( ji, 1 ) = mxl_min
874	zldw( ji, 1 ) = mxl_min
875	zlup( ji, jpka ) = mxl_min
876	zldw( ji, jpka ) = mxl_min
877	END DO
878
879	DO jk = 2,jpka-1
880	DO ji = 1, jpi
881	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
882	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
883	END DO
884	END DO
885	!!
886	!! BL89 search for lup
887	!! ----------------------------------------------------------
888	DO jk=2,jpka-1
889	!
890	DO ji = 1, jpi
891	zCF(ji,1:jpka) = 0._wp
892	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
893	ln_foundl(ji ) = .false.
894	END DO
895	!
896	DO jkup=jk+1,jpka-1
897	DO ji = 1, jpi
898	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
899	& ( zbn2(ji,jj,jkup )*(ghw_abl(jkup )-ghw_abl(jk)) &
900	& + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
901	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
902	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
903	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
904	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
905	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
906	zlup(ji,jk) = zcff
907	ln_foundl(ji) = .true.
908	END IF
909	END DO
910	END DO
911	!
912	END DO
913	!!
914	!! BL89 search for ldwn
915	!! ----------------------------------------------------------
916	DO jk=2,jpka-1
917	!
918	DO ji = 1, jpi
919	zCF(ji,1:jpka) = 0._wp
920	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
921	ln_foundl(ji ) = .false.
922	END DO
923	!
924	DO jkdwn=jk-1,1,-1
925	DO ji = 1, jpi
926	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
927	& * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
928	+ zbn2(ji,jj,jkdwn )*(ghw_abl(jk)-ghw_abl(jkdwn )) )
929	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
930	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
931	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
932	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
933	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
934	zldw(ji,jk) = zcff
935	ln_foundl(ji) = .true.
936	END IF
937	END DO
938	END DO
939	!
940	END DO
941
942	DO jk = 1, jpka
943	DO ji = 1, jpi
944	mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )
945	END DO
946	END DO
947
948	END DO
949	# undef zlup
950	# undef zldw
951	!
952	END SELECT
953	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
954	! ! Finalize the computation of turbulent visc./diff.
955	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
956
957	!-------------
958	DO jj = 1, jpj ! outer loop
959	!-------------
960	DO jk = 1, jpka
961	DO ji = 1, jpi ! vector opt.
962	zcff = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk ) &
963	& * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )
964	zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z)
965	zcff = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
966	!!FL: MAX function probably useless because of the definition of mxl_min
967	Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak )
968	Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak )
969	END DO
970	END DO
971	!-------------
972	END DO
973	!-------------
974
975	!---------------------------------------------------------------------------------------------------
976	END SUBROUTINE abl_zdf_tke
977	!===================================================================================================
978
979
980	!===================================================================================================
981	SUBROUTINE smooth_pblh( pvar2d, msk )
982	!---------------------------------------------------------------------------------------------------
983
984	!!----------------------------------------------------------------------
985	!! * ROUTINE smooth_pblh *
986	!!
987	!! ** Purpose : 2D Hanning filter on atmospheric PBL height
988	!!
989	!! ---------------------------------------------------------------------
990	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk
991	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
992	INTEGER :: ji,jj
993	REAL(wp) :: smth_a, smth_b
994	REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY
995	REAL(wp) :: zumsk,zvmsk
996	!!
997	!!=========================================================
998	!!
999	!! Hanning filter
1000	smth_a = 1._wp / 8._wp
1001	smth_b = 1._wp / 4._wp
1002	!
1003	DO_2D_11_10
1004	zumsk = msk(ji,jj) * msk(ji+1,jj)
1005	zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk
1006	END_2D
1007
1008	DO_2D_10_11
1009	zvmsk = msk(ji,jj) * msk(ji,jj+1)
1010	zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk
1011	END_2D
1012
1013	DO_2D_10_00
1014	zFY ( ji, jj ) = zdY ( ji, jj ) &
1015	& + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) &
1016	& - (zdX ( ji, jj ) - zdX( ji-1, jj )) )
1017	END_2D
1018
1019	DO_2D_00_10
1020	zFX( ji, jj ) = zdX( ji, jj ) &
1021	& + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) &
1022	& - (zdY( ji , jj ) - zdY( ji , jj-1)) )
1023	END_2D
1024
1025	DO_2D_00_00
1026	pvar2d( ji ,jj ) = pvar2d( ji ,jj ) &
1027	& + msk(ji,jj) * smth_b * ( &
1028	& zFX( ji, jj ) - zFX( ji-1, jj ) &
1029	& +zFY( ji, jj ) - zFY( ji, jj-1 ) )
1030	END_2D
1031	!!
1032	!---------------------------------------------------------------------------------------------------
1033	END SUBROUTINE smooth_pblh
1034	!===================================================================================================
1035
1036	!!======================================================================
1037	END MODULE ablmod

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

Download in other formats: