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

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

ablmod.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL – NEMO

Context Navigation

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

Last change on this file since 11363 was 11363, checked in by gsamson, 5 years ago

dev_r11265_ABL : see #2131

fix xml files to get xios working with ABL
introduce rdt_abl = rdt * nn_fsbc to get ABL working with nn_fsbc > 1
does not change results compared to rev11360 with nn_fsbc = 1

File size: 50.9 KB

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

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

Download in other formats: