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.

zdfbfr.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90 @ 4528

Last change on this file since 4528 was 4381, checked in by jchanut, 10 years ago
bug with log bottom friction and namelist update
Property svn:keywords set to `Id`
File size: 16.6 KB

Line
1	MODULE zdfbfr
2	!!======================================================================
3	!! * MODULE zdfbfr *
4	!! Ocean physics: Bottom friction
5	!!======================================================================
6	!! History : OPA ! 1997-06 (G. Madec, A.-M. Treguier) Original code
7	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
8	!! 3.2 ! 2009-09 (A.C.Coward) Correction to include barotropic contribution
9	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
10	!! 3.4 ! 2011-11 (H. Liu) implementation of semi-implicit bottom friction option
11	!! ! 2012-06 (H. Liu) implementation of Log Layer bottom friction option
12	!!----------------------------------------------------------------------
13
14	!!----------------------------------------------------------------------
15	!! zdf_bfr : update bottom friction coefficient (non-linear bottom friction only)
16	!! zdf_bfr_init : read in namelist and control the bottom friction parameters.
17	!!----------------------------------------------------------------------
18	USE oce ! ocean dynamics and tracers variables
19	USE dom_oce ! ocean space and time domain variables
20	USE zdf_oce ! ocean vertical physics variables
21	USE in_out_manager ! I/O manager
22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
23	USE lib_mpp ! distributed memory computing
24	USE prtctl ! Print control
25	USE timing ! Timing
26	USE wrk_nemo ! Memory Allocation
27	USE phycst, ONLY: vkarmn
28
29	IMPLICIT NONE
30	PRIVATE
31
32	PUBLIC zdf_bfr ! called by step.F90
33	PUBLIC zdf_bfr_init ! called by nemogcm.F90
34
35	! !!* Namelist nambfr: bottom friction namelist *
36	INTEGER , PUBLIC :: nn_bfr ! = 0/1/2/3 type of bottom friction (PUBLIC for TAM)
37	REAL(wp), PUBLIC :: rn_bfri1 ! bottom drag coefficient (linear case) (PUBLIC for TAM)
38	REAL(wp), PUBLIC :: rn_bfri2 ! bottom drag coefficient (non linear case) (PUBLIC for TAM)
39	REAL(wp), PUBLIC :: rn_bfri2_max ! Maximum bottom drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM)
40	REAL(wp), PUBLIC :: rn_bfeb2 ! background bottom turbulent kinetic energy [m2/s2] (PUBLIC for TAM)
41	REAL(wp), PUBLIC :: rn_bfrien ! local factor to enhance coefficient bfri (PUBLIC for TAM)
42	LOGICAL , PUBLIC :: ln_bfr2d ! logical switch for 2D enhancement (PUBLIC for TAM)
43	LOGICAL , PUBLIC :: ln_loglayer ! switch for log layer bfr coeff. (PUBLIC for TAM)
44	REAL(wp), PUBLIC :: rn_bfrz0 ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
45	LOGICAL , PUBLIC :: ln_bfrimp ! logical switch for implicit bottom friction
46	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d ! 2D bottom drag coefficient (PUBLIC for TAM)
47
48	!! * Substitutions
49	# include "vectopt_loop_substitute.h90"
50	# include "domzgr_substitute.h90"
51	!!----------------------------------------------------------------------
52	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
53	!! $Id$
54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
55	!!----------------------------------------------------------------------
56	CONTAINS
57
58	INTEGER FUNCTION zdf_bfr_alloc()
59	!!----------------------------------------------------------------------
60	!! * FUNCTION zdf_bfr_alloc *
61	!!----------------------------------------------------------------------
62	ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
63	!
64	IF( lk_mpp ) CALL mpp_sum ( zdf_bfr_alloc )
65	IF( zdf_bfr_alloc /= 0 ) CALL ctl_warn('zdf_bfr_alloc: failed to allocate arrays.')
66	END FUNCTION zdf_bfr_alloc
67
68
69	SUBROUTINE zdf_bfr( kt )
70	!!----------------------------------------------------------------------
71	!! * ROUTINE zdf_bfr *
72	!!
73	!! ** Purpose : compute the bottom friction coefficient.
74	!!
75	!! ** Method : Calculate and store part of the momentum trend due
76	!! to bottom friction following the chosen friction type
77	!! (free-slip, linear, or quadratic). The component
78	!! calculated here is multiplied by the bottom velocity in
79	!! dyn_bfr to provide the trend term.
80	!! The coefficients are updated at each time step only
81	!! in the quadratic case.
82	!!
83	!! ** Action : bfrua , bfrva bottom friction coefficients
84	!!----------------------------------------------------------------------
85	INTEGER, INTENT( in ) :: kt ! ocean time-step index
86	!!
87	INTEGER :: ji, jj ! dummy loop indices
88	INTEGER :: ikbt, ikbu, ikbv ! local integers
89	REAL(wp) :: zvu, zuv, zecu, zecv, ztmp ! temporary scalars
90	REAL(wp), POINTER, DIMENSION(:,:) :: zbfrt
91	!!----------------------------------------------------------------------
92	!
93	IF( nn_timing == 1 ) CALL timing_start('zdf_bfr')
94	!
95	IF( kt == nit000 .AND. lwp ) THEN
96	WRITE(numout,*)
97	WRITE(numout,*) 'zdf_bfr : Set bottom friction coefficient (non-linear case)'
98	WRITE(numout,*) '~~~~~~~~'
99	ENDIF
100	!
101	IF( nn_bfr == 2 ) THEN ! quadratic bottom friction only
102	!
103	CALL wrk_alloc( jpi, jpj, zbfrt )
104
105	IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient
106
107	# if defined key_vectopt_loop
108	DO jj = 1, 1
109	!CDIR NOVERRCHK
110	DO ji = 1, jpij ! vector opt. (forced unrolling)
111	# else
112	!CDIR NOVERRCHK
113	DO jj = 1, jpj
114	!CDIR NOVERRCHK
115	DO ji = 1, jpi
116	# endif
117	ikbt = mbkt(ji,jj)
118	! JC: possible WAD implementation should modify line below if layers vanish
119	ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
120	zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
121	zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max)
122	END DO
123	END DO
124	!
125	ELSE
126	zbfrt(:,:) = bfrcoef2d(:,:)
127	ENDIF
128
129	# if defined key_vectopt_loop
130	DO jj = 1, 1
131	!CDIR NOVERRCHK
132	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
133	# else
134	!CDIR NOVERRCHK
135	DO jj = 2, jpjm1
136	!CDIR NOVERRCHK
137	DO ji = 2, jpim1
138	# endif
139	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
140	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
141	!
142	zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) &
143	& + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) )
144	zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) &
145	& + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) )
146	!
147	zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 )
148	zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 )
149	!
150	bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) * zecu
151	bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) * zecv
152	END DO
153	END DO
154
155	!
156	CALL lbc_lnk( bfrua, 'U', 1. ) ; CALL lbc_lnk( bfrva, 'V', 1. ) ! Lateral boundary condition
157	!
158	IF(ln_ctl) CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr - u: ', mask1=umask, &
159	& tab2d_2=bfrva, clinfo2= ' v: ', mask2=vmask,ovlap=1 )
160	CALL wrk_dealloc( jpi,jpj,zbfrt )
161	ENDIF
162	!
163	IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr')
164	!
165	END SUBROUTINE zdf_bfr
166
167
168	SUBROUTINE zdf_bfr_init
169	!!----------------------------------------------------------------------
170	!! * ROUTINE zdf_bfr_init *
171	!!
172	!! ** Purpose : Initialization of the bottom friction
173	!!
174	!! ** Method : Read the nambfr namelist and check their consistency
175	!! called at the first timestep (nit000)
176	!!----------------------------------------------------------------------
177	USE iom ! I/O module for ehanced bottom friction file
178	!!
179	INTEGER :: inum ! logical unit for enhanced bottom friction file
180	INTEGER :: ji, jj ! dummy loop indexes
181	INTEGER :: ikbt, ikbu, ikbv ! temporary integers
182	INTEGER :: ictu, ictv ! - -
183	INTEGER :: ios
184	REAL(wp) :: zminbfr, zmaxbfr ! temporary scalars
185	REAL(wp) :: ztmp, zfru, zfrv ! - -
186	!!
187	NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
188	& rn_bfrien, ln_bfrimp, ln_loglayer
189	!!----------------------------------------------------------------------
190	!
191	IF( nn_timing == 1 ) CALL timing_start('zdf_bfr_init')
192	!
193	! !* Allocate zdfbfr arrays
194	IF( zdf_bfr_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' )
195	!
196	! !* Parameter control and print
197	!
198	REWIND( numnam_ref ) ! Namelist nambfr in reference namelist : Bottom momentum boundary condition
199	READ ( numnam_ref, nambfr, IOSTAT = ios, ERR = 901)
200	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in reference namelist', lwp )
201
202	REWIND( numnam_cfg ) ! Namelist nambfr in configuration namelist : Bottom momentum boundary condition
203	READ ( numnam_cfg, nambfr, IOSTAT = ios, ERR = 902 )
204	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in configuration namelist', lwp )
205	WRITE ( numond, nambfr )
206	IF(lwp) WRITE(numout,*)
207	IF(lwp) WRITE(numout,*) 'zdf_bfr_init : momentum bottom friction'
208	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
209	IF(lwp) WRITE(numout,*) ' Namelist nam_bfr : set bottom friction parameters'
210	!
211	SELECT CASE (nn_bfr)
212	!
213	CASE( 0 )
214	IF(lwp) WRITE(numout,*) ' free-slip '
215	bfrua(:,:) = 0.e0
216	bfrva(:,:) = 0.e0
217	!
218	CASE( 1 )
219	IF(lwp) WRITE(numout,*) ' linear botton friction'
220	IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri1 = ', rn_bfri1
221	IF( ln_bfr2d ) THEN
222	IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_bfr2d = ', ln_bfr2d
223	IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien
224	ENDIF
225	!
226	IF(ln_bfr2d) THEN
227	! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
228	CALL iom_open('bfr_coef.nc',inum)
229	CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
230	CALL iom_close(inum)
231	bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
232	ELSE
233	bfrcoef2d(:,:) = rn_bfri1 ! initialize bfrcoef2d to the namelist variable
234	ENDIF
235	!
236	bfrua(:,:) = - bfrcoef2d(:,:)
237	bfrva(:,:) = - bfrcoef2d(:,:)
238	!
239	CASE( 2 )
240	IF(lwp) WRITE(numout,*) ' quadratic bottom friction'
241	IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri2 = ', rn_bfri2
242	IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_bfri2_max = ', rn_bfri2_max
243	IF(lwp) WRITE(numout,*) ' background tke rn_bfeb2 = ', rn_bfeb2
244	IF(lwp) WRITE(numout,*) ' log formulation ln_bfr2d = ', ln_loglayer
245	IF(lwp) WRITE(numout,*) ' bottom roughness rn_bfrz0 [m] = ', rn_bfrz0
246	IF( rn_bfrz0<=0.e0 ) THEN
247	WRITE(ctmp1,*) ' bottom roughness must be strictly positive'
248	CALL ctl_stop( ctmp1 )
249	ENDIF
250	IF( ln_bfr2d ) THEN
251	IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_bfr2d = ', ln_bfr2d
252	IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien
253	ENDIF
254	!
255	IF(ln_bfr2d) THEN
256	! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
257	CALL iom_open('bfr_coef.nc',inum)
258	CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
259	CALL iom_close(inum)
260	!
261	bfrcoef2d(:,:) = rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
262	ELSE
263	bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable
264	ENDIF
265	!
266	IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
267	# if defined key_vectopt_loop
268	DO jj = 1, 1
269	!CDIR NOVERRCHK
270	DO ji = 1, jpij ! vector opt. (forced unrolling)
271	# else
272	!CDIR NOVERRCHK
273	DO jj = 1, jpj
274	!CDIR NOVERRCHK
275	DO ji = 1, jpi
276	# endif
277	ikbt = mbkt(ji,jj)
278	ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
279	bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
280	bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max)
281	END DO
282	END DO
283	ENDIF
284	!
285	CASE DEFAULT
286	IF(lwp) WRITE(ctmp1,*) ' bad flag value for nn_bfr = ', nn_bfr
287	CALL ctl_stop( ctmp1 )
288	!
289	END SELECT
290	!
291	IF(lwp) WRITE(numout,*) ' implicit bottom friction switch ln_bfrimp = ', ln_bfrimp
292	!
293	! ! Make sure ln_zdfexp=.false. when use implicit bfr
294	IF( ln_bfrimp .AND. ln_zdfexp ) THEN
295	IF(lwp) THEN
296	WRITE(numout,*)
297	WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.'
298	WRITE(numout,*) ' but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!'
299	WRITE(ctmp1,*) ' set either ln_zdfexp = .false or ln_bfrimp = .false.'
300	CALL ctl_stop( ctmp1 )
301	END IF
302	END IF
303	!
304	! Basic stability check on bottom friction coefficient
305	!
306	ictu = 0 ! counter for stability criterion breaches at U-pts
307	ictv = 0 ! counter for stability criterion breaches at V-pts
308	zminbfr = 1.e10_wp ! initialise tracker for minimum of bottom friction coefficient
309	zmaxbfr = -1.e10_wp ! initialise tracker for maximum of bottom friction coefficient
310	!
311	# if defined key_vectopt_loop
312	DO jj = 1, 1
313	!CDIR NOVERRCHK
314	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
315	# else
316	!CDIR NOVERRCHK
317	DO jj = 2, jpjm1
318	!CDIR NOVERRCHK
319	DO ji = 2, jpim1
320	# endif
321	ikbu = mbku(ji,jj) ! deepest ocean level at u- and v-points
322	ikbv = mbkv(ji,jj)
323	zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
324	zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
325	IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN
326	IF( ln_ctl ) THEN
327	WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
328	WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru
329	ENDIF
330	ictu = ictu + 1
331	ENDIF
332	IF( ABS( bfrcoef2d(ji,jj) ) > zfrv ) THEN
333	IF( ln_ctl ) THEN
334	WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv
335	WRITE(numout,*) 'BFR ', bfrcoef2d(ji,jj), zfrv
336	ENDIF
337	ictv = ictv + 1
338	ENDIF
339	zminbfr = MIN( zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) ) )
340	zmaxbfr = MAX( zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) ) )
341	END DO
342	END DO
343	IF( lk_mpp ) THEN
344	CALL mpp_sum( ictu )
345	CALL mpp_sum( ictv )
346	CALL mpp_min( zminbfr )
347	CALL mpp_max( zmaxbfr )
348	ENDIF
349	IF( .NOT.ln_bfrimp) THEN
350	IF( lwp .AND. ictu + ictv > 0 ) THEN
351	WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
352	WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
353	WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
354	WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
355	ENDIF
356	ENDIF
357	!
358	IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr_init')
359	!
360	END SUBROUTINE zdf_bfr_init
361
362	!!======================================================================
363	END MODULE zdfbfr

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

Download in other formats: