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_ASINTER-01-05_merged/src/ABL – NEMO

Context Navigation

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ABL/ablmod.F90 @ 12015

Last change on this file since 12015 was 12015, checked in by gsamson, 4 years ago
dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK
File size: 53.1 KB

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

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

Download in other formats: