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 @ 11348

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

dev_r11265_ABL :

merge HPC-13_IRRMANN_BDY_optimization branch @ rev11332 with dev_r11265_ABL branch @ rev11334
allow ln_dm2dc option with ABL
cosmetic change in sbcabl.F90

identical results with rev11334 for bulk and abl orca2

File size: 50.5 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 * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
156	z_elem_c( ji, jk ) = - rdt * 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 * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )
165	! Homogeneous Neumann at the top
166	z_elem_a( ji, jpka ) = - rdt * 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 * zztmp1
188	tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rdt * 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 * zztmp1
200	tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rdt * 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 )( fft_abl(ji,jj) rdt ) ! (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 * 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 * 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 * 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 * 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 * e3t_abl(jk) * pgu_dta(ji,jj,jk)
283	v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * 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 * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
301	z_elem_c( ji, jk ) = - rdt * 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 * 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 * zztmp1
319	u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rdt * zztmp2
320
321	!++ Top Neumann B.C.
322	!z_elem_a( ji, jpka ) = - 0.5_wp * rdt * ( 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 * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
368	z_elem_c( ji, jk ) = -rdt * 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 * 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 * zztmp1
386	v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rdt * zztmp2
387	!++ Top Neumann B.C.
388	!z_elem_a( ji, jpka ) = -rdt * 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 * rdt * zcff2 ! zcff = 1 for masked points
443
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 * rdt * zcff2 ! zcff = 1 for masked points
469	tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) &
470	& + zcff * pt_dta( ji, jj, jk )
471
472	tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) &
473	& + zcff * pq_dta( ji, jj, jk )
474
475	! z_cft( ji, jj, jk ) = zcff
476	END DO
477	END DO
478	!-------------
479	END DO ! end outer loop
480	!-------------
481	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
482	! ! 6 *** MPI exchanges
483	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
484	!
485	CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), 'T', -1., kfillmode = jpfillnothing )
486	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...
487	!
488	CALL iom_put ( 'u_abl', u_abl(:,:,2:jpka,nt_a ))
489	CALL iom_put ( 'v_abl', v_abl(:,:,2:jpka,nt_a ))
490	CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta))
491	CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa))
492	! CALL iom_put ( 'coeft', z_cft(:,:,2:jpka ))
493	CALL iom_put ( 'pblh', pblh(:,: ))
494	!
495	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
496	! ! 7 *** Finalize flux computation
497	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
498
499	DO jj = 1, jpj
500	DO ji = 1, jpi
501	ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta )
502	zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa )
503	zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ...
504	& ( R_dryztemp ( 1._wp + rctv0*zhumi ) )
505	psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
506	pevp ( ji, jj ) = rn_efacMAX( 0._wp, zcff pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
507	zrhoa( ji, jj ) = zcff
508	END DO
509	END DO
510
511	DO jj = 2, jpj
512	DO ji = 2, jpi ! vect. opt.
513	zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji ,jj) + pssu(ji-1,jj) )
514	zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj ) + pssv(ji,jj-1) )
515	END DO
516	END DO
517	!
518	CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
519	!
520	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
521	DO jj = 1, jpj
522	DO ji = 1, jpi
523	zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) &
524	& + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj)
525	zztmp = zrhoa(ji,jj) * pcd_du(ji,jj)
526
527	pwndm (ji,jj) = zcff
528	ptaum (ji,jj) = zztmp * zcff
529	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
530	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
531	END DO
532	END DO
533	! ... utau, vtau at U- and V_points, resp.
534	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
535	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
536	DO jj = 2, jpjm1
537	DO ji = 2, jpim1
538	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
539	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
540	ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
541	zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
542	zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
543	ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
544	END DO
545	END DO
546	!
547	CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
548
549	CALL iom_put( "taum_oce", ptaum )
550
551	IF(ln_ctl) THEN
552	CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' )
553	CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ' )
554	CALL prt_ctl( tab2d_2=ptauj , clinfo2= 'vtau : ' )
555	ENDIF
556
557	#if defined key_si3
558	! ------------------------------------------------------------ !
559	! Wind stress relative to the moving ice ( U10m - U_ice ) !
560	! ------------------------------------------------------------ !
561	DO jj = 2, jpjm1
562	DO ji = 2, jpim1
563
564	zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
565	zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
566
567	ptaui_ice(ji,jj) = 0.5_wp * ( zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &
568	& + zrhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &
569	& * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
570	ptauj_ice(ji,jj) = 0.5_wp * ( zrhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &
571	& + zrhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &
572	& * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
573	END DO
574	END DO
575	CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
576	!
577	IF(ln_ctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
578	& , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
579	#endif
580	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
581	! ! 8 *** Swap time indices for the next timestep
582	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
583	nt_n = 1 + MOD( kt , 2)
584	nt_a = 1 + MOD( kt+1, 2)
585	!
586	!---------------------------------------------------------------------------------------------------
587	END SUBROUTINE abl_stp
588	!===================================================================================================
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606	!===================================================================================================
607	SUBROUTINE abl_zdf_tke( )
608	!---------------------------------------------------------------------------------------------------
609
610	!!----------------------------------------------------------------------
611	!! * ROUTINE abl_zdf_tke *
612	!!
613	!! ** Purpose : Time-step Turbulente Kinetic Energy (TKE) equation
614	!!
615	!! ** Method : - source term due to shear
616	!! - source term due to stratification
617	!! - resolution of the TKE equation by inverting
618	!! a tridiagonal linear system
619	!!
620	!! ** Action : - en : now turbulent kinetic energy)
621	!! - avmu, avmv : production of TKE by shear at u and v-points
622	!! (= Kz dz[Ub] * dz[Un] )
623	!! ---------------------------------------------------------------------
624	INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
625	INTEGER, DIMENSION(1:jpi ) :: ikbl
626	REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
627	REAL(wp) :: zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
628	REAL(wp) :: zdU2,zdV2
629	REAL(wp) :: zwndi,zwndj
630	REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2
631	REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2
632	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF
633	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
634	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
635	REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
636	LOGICAL :: ln_Patankar = .FALSE.
637	LOGICAL :: ln_dumpvar = .FALSE.
638	LOGICAL , DIMENSION(1:jpi ) :: ln_foundl
639	!
640	tind = nt_n
641	ziRic = 1._wp / rn_Ric
642	! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
643	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
644	! ! Advance TKE equation to time n+1
645	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
646	!-------------
647	DO jj = 1, jpj ! outer loop
648	!-------------
649	!
650	! Compute vertical shear
651	DO jk = 2, jpkam1
652	DO ji = 1,jpi
653	zcff = 1.0_wp / e3w_abl( jk )**2
654	zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
655	zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
656	zsh2(ji,jk) = zdU+zdV
657	END DO
658	END DO
659	!
660	! Compute brunt-vaisala frequency
661	DO jk = 2, jpkam1
662	DO ji = 1,jpi
663	zcff = grav * itvref / e3w_abl( jk )
664	zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta)
665	zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) &
666	& - tq_abl( ji, jj, jk , tind, jp_ta) * tq_abl( ji, jj, jk , tind, jp_qa)
667	zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi)
668	END DO
669	END DO
670	!
671	! Terms for the tridiagonal problem
672	DO jk = 2, jpkam1
673	DO ji = 1,jpi
674	zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point
675	zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
676	zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
677
678	z_elem_a( ji, jk ) = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal
679	z_elem_c( ji, jk ) = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal
680	IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE
681	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
682	& + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal
683	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) ) ! right-hand-side
684	ELSE
685	z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &
686	& + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal
687	& - e3w_abl(jk) * rdt * zbuoy
688	tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * zshear ) ! right-hand-side
689	END IF
690	END DO
691	END DO
692
693	DO ji = 1,jpi ! vector opt.
694	zesrf = MAX( 4.63_wp * ustar2(ji,jj), tke_min )
695	zetop = tke_min
696	z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp
697	z_elem_a ( ji, jpka ) = 0._wp; z_elem_c ( ji, jpka ) = 0._wp; z_elem_b ( ji, jpka ) = 1._wp
698	tke_abl( ji, jj, 1, nt_a ) = zesrf
699	tke_abl( ji, jj, jpka, nt_a ) = zetop
700	zbn2(ji,jj, 1) = zbn2( ji,jj, 2)
701	zsh2(ji, 1) = zsh2( ji, 2)
702	zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)
703	zsh2(ji, jpka) = zsh2( ji , jpkam1)
704	END DO
705	!!
706	!! Matrix inversion
707	!! ----------------------------------------------------------
708	DO ji = 1,jpi
709	zcff = 1._wp / z_elem_b( ji, 1 )
710	zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 )
711	tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a )
712	END DO
713
714	DO jk = 2, jpka
715	DO ji = 1,jpi
716	zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
717	zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
718	tke_abl(ji,jj,jk,nt_a) = zcff * ( tke_abl(ji,jj,jk ,nt_a) &
719	& - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
720	END DO
721	END DO
722
723	DO jk = jpkam1,1,-1
724	DO ji = 1,jpi
725	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)
726	END DO
727	END DO
728
729	!!FL should not be needed because of Patankar procedure
730	tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
731
732	!!
733	!! Diagnose PBL height
734	!! ----------------------------------------------------------
735
736
737	!
738	! arrays zRH, zFC and zCF are available at this point
739	! and zFC(:, 1 ) = 0.
740	! diagnose PBL height based on zsh2 and zbn2
741	zFC ( : ,1) = 0._wp
742	ikbl( 1:jpi ) = 0
743
744	DO jk = 2,jpka
745	DO ji = 1, jpi
746	zcff = ghw_abl( jk-1 )
747	zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
748	zcff = ghw_abl( jk )
749	zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
750	zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*( &
751	zcff2 * ( zsh2( ji, jk ) - ziRic * zbn2( ji, jj, jk ) &
752	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
753	+ zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
754	- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
755	& )
756	IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
757	END DO
758	END DO
759	!
760	! finalize the computation of the PBL height
761	DO ji = 1, jpi
762	jk = ikbl(ji)
763	IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
764	pblh( ji, jj ) = ( ghw_abl( jk-1 ) * zFC( ji, jk ) &
765	& - ghw_abl( jk ) * zFC( ji, jk-1 ) &
766	& ) / ( zFC( ji, jk ) - zFC( ji, jk-1 ) )
767	ELSE IF( jk==2 ) THEN
768	pblh( ji, jj ) = ghw_abl(2 )
769	ELSE
770	pblh( ji, jj ) = ghw_abl(jpka)
771	END IF
772	END DO
773	! Optional : could add pblh smoothing if pblh is noisy horizontally ...
774	!
775	!-------------
776	END DO
777	!-------------
778	IF(ln_smth_pblh) THEN
779	CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
780	CALL smooth_pblh( pblh, msk_abl )
781	CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
782	ENDIF
783	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
784	! ! Diagnostic mixing length computation
785	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
786	!
787	SELECT CASE ( nn_amxl )
788	!
789	CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom
790	# define zlup zRH
791	# define zldw zFC
792	DO jj = 1, jpj ! outer loop
793	!
794	DO ji = 1, jpi
795	mxl_abl ( ji, jj, 1 ) = 0._wp
796	mxl_abl ( ji, jj, jpka ) = mxl_min
797	zldw( ji, 1 ) = 0._wp
798	zlup( ji, jpka ) = 0._wp
799	END DO
800	!
801	DO jk = 2, jpkam1
802	DO ji = 1, jpi
803	zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
804	mxl_abl( ji, jj, jk ) = MAX( mxl_min, &
805	& SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
806	END DO
807	END DO
808	!
809	! Limit mxl
810	DO jk = jpkam1,1,-1
811	DO ji = 1, jpi
812	zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
813	END DO
814	END DO
815	!
816	DO jk = 2, jpka
817	DO ji = 1, jpi
818	zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )
819	END DO
820	END DO
821	!
822	DO jk = 1, jpka
823	DO ji = 1, jpi
824	mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
825	END DO
826	END DO
827	!
828	END DO
829	# undef zlup
830	# undef zldw
831	!
832	!
833	CASE ( 1 ) ! length-scale computed as the distance to the PBL height
834	DO jj = 1,jpj ! outer loop
835	!
836	DO ji = 1, jpi ! vector opt.
837	zcff = 1._wp / pblh( ji, jj ) ! inverse of hbl
838	DO jk = 1, jpka
839	zsig = MIN( zcff * ghw_abl( jk ), 1. )
840	zcff1 = pblh( ji, jj )
841	mxl_abl( ji, jj, jk ) = mxl_min &
842	& + zsig * ( amx1zcff1 + bmx1mxl_min ) &
843	& + zsig * zsig * ( amx2zcff1 + bmx2mxl_min ) &
844	& + zsig*3 ( amx3zcff1 + bmx3mxl_min ) &
845	& + zsig*4 ( amx4zcff1 + bmx4mxl_min ) &
846	& + zsig*5 ( amx5zcff1 + bmx5mxl_min )
847	END DO
848	END DO
849	!
850	END DO
851	!
852	CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale
853	!
854	# define zlup zRH
855	# define zldw zFC
856	! zCF is used for matrix inversion
857	!
858	DO jj = 1, jpj ! outer loop
859
860	DO ji = 1, jpi
861	zlup( ji, 1 ) = mxl_min
862	zldw( ji, 1 ) = mxl_min
863	zlup( ji, jpka ) = mxl_min
864	zldw( ji, jpka ) = mxl_min
865	END DO
866
867	DO jk = 2,jpka-1
868	DO ji = 1, jpi
869	zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
870	zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
871	END DO
872	END DO
873	!!
874	!! BL89 search for lup
875	!! ----------------------------------------------------------
876	DO jk=2,jpka-1
877	!
878	DO ji = 1, jpi
879	zCF(ji,1:jpka) = 0._wp
880	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
881	ln_foundl(ji ) = .false.
882	END DO
883	!
884	DO jkup=jk+1,jpka-1
885	DO ji = 1, jpi
886	zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
887	& ( zbn2(ji,jj,jkup )*(ghw_abl(jkup )-ghw_abl(jk)) &
888	& + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
889	IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
890	zcff2 = ghw_abl(jkup ) - ghw_abl(jk)
891	zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
892	zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / &
893	& ( zCF(ji,jkup) - zCF(ji,jkup-1) )
894	zlup(ji,jk) = zcff
895	ln_foundl(ji) = .true.
896	END IF
897	END DO
898	END DO
899	!
900	END DO
901	!!
902	!! BL89 search for ldwn
903	!! ----------------------------------------------------------
904	DO jk=2,jpka-1
905	!
906	DO ji = 1, jpi
907	zCF(ji,1:jpka) = 0._wp
908	zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
909	ln_foundl(ji ) = .false.
910	END DO
911	!
912	DO jkdwn=jk-1,1,-1
913	DO ji = 1, jpi
914	zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
915	& * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
916	+ zbn2(ji,jj,jkdwn )*(ghw_abl(jk)-ghw_abl(jkdwn )) )
917	IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
918	zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
919	zcff1 = ghw_abl(jk) - ghw_abl(jkdwn )
920	zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / &
921	& ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) )
922	zldw(ji,jk) = zcff
923	ln_foundl(ji) = .true.
924	END IF
925	END DO
926	END DO
927	!
928	END DO
929
930	DO jk = 1, jpka
931	DO ji = 1, jpi
932	mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )
933	END DO
934	END DO
935
936	END DO
937	# undef zlup
938	# undef zldw
939	!
940	END SELECT
941	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
942	! ! Finalize the computation of turbulent visc./diff.
943	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
944
945	!-------------
946	DO jj = 1, jpj ! outer loop
947	!-------------
948	DO jk = 1, jpka
949	DO ji = 1, jpi ! vector opt.
950	zcff = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk ) &
951	& * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )
952	zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z)
953	zcff = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
954	!!FL: MAX function probably useless because of the definition of mxl_min
955	Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak )
956	Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak )
957	END DO
958	END DO
959	!-------------
960	END DO
961	!-------------
962
963	!---------------------------------------------------------------------------------------------------
964	END SUBROUTINE abl_zdf_tke
965	!===================================================================================================
966
967
968	!===================================================================================================
969	SUBROUTINE smooth_pblh( pvar2d, msk )
970	!---------------------------------------------------------------------------------------------------
971
972	!!----------------------------------------------------------------------
973	!! * ROUTINE smooth_pblh *
974	!!
975	!! ** Purpose : 2D Hanning filter on atmospheric PBL height
976	!!
977	!! ---------------------------------------------------------------------
978	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk
979	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
980	INTEGER :: ji,jj
981	REAL(wp) :: smth_a, smth_b
982	REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY
983	REAL(wp) :: zumsk,zvmsk
984	!!
985	!!=========================================================
986	!!
987	!! Hanning filter
988	smth_a = 1._wp / 8._wp
989	smth_b = 1._wp / 4._wp
990	!
991	DO jj=1,jpj
992	DO ji=1,jpi-1
993	zumsk = msk(ji,jj) * msk(ji+1,jj)
994	zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk
995	END DO
996	END DO
997
998	DO jj=1,jpj-1
999	DO ji=1,jpi
1000	zvmsk = msk(ji,jj) * msk(ji,jj+1)
1001	zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk
1002	END DO
1003	END DO
1004
1005	DO jj=1,jpj-1
1006	DO ji=2,jpi-1
1007	zFY ( ji, jj ) = zdY ( ji, jj ) &
1008	& + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) &
1009	& - (zdX ( ji, jj ) - zdX( ji-1, jj )) )
1010	END DO
1011	END DO
1012
1013	DO jj=2,jpj-1
1014	DO ji=1,jpi-1
1015	zFX( ji, jj ) = zdX( ji, jj ) &
1016	& + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) &
1017	& - (zdY( ji , jj ) - zdY( ji , jj-1)) )
1018	END DO
1019	END DO
1020
1021	DO jj = 2, jpj-1
1022	DO ji = 2,jpi-1
1023	pvar2d( ji ,jj ) = pvar2d( ji ,jj ) &
1024	& + msk(ji,jj) * smth_b * ( &
1025	& zFX( ji, jj ) - zFX( ji-1, jj ) &
1026	& +zFY( ji, jj ) - zFY( ji, jj-1 ) )
1027	END DO
1028	END DO
1029	!!
1030	!---------------------------------------------------------------------------------------------------
1031	END SUBROUTINE smooth_pblh
1032	!===================================================================================================
1033
1034	!!======================================================================
1035	END MODULE ablmod

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

Download in other formats: