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 branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90 @ 4146

Last change on this file since 4146 was 3820, checked in by hliu, 11 years ago
add the 2D-enhancement for loglayer bfr coefficient. see ticket 1058
Property svn:keywords set to `Id`
File size: 15.5 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 momentum Kz at the ocean bottom due to the type of bottom friction chosen
16	!! zdf_bfr_init : read in namelist and control the bottom friction parameters.
17	!! zdf_bfr_2d : read in namelist and control the bottom friction parameters.
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers variables
20	USE dom_oce ! ocean space and time domain variables
21	USE zdf_oce ! ocean vertical physics variables
22	USE in_out_manager ! I/O manager
23	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
24	USE lib_mpp ! distributed memory computing
25	USE prtctl ! Print control
26	USE timing ! Timing
27
28	USE phycst, ONLY: vkarmn
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC zdf_bfr ! called by step.F90
34	PUBLIC zdf_bfr_init ! called by opa.F90
35
36	! !!* Namelist nambfr: bottom friction namelist *
37	INTEGER , PUBLIC :: nn_bfr = 0 ! = 0/1/2/3 type of bottom friction (PUBLIC for TAM)
38	REAL(wp), PUBLIC :: rn_bfri1 = 4.0e-4_wp ! bottom drag coefficient (linear case) (PUBLIC for TAM)
39	REAL(wp), PUBLIC :: rn_bfri2 = 1.0e-3_wp ! bottom drag coefficient (non linear case) (PUBLIC for TAM)
40	REAL(wp), PUBLIC :: rn_bfeb2 = 2.5e-3_wp ! background bottom turbulent kinetic energy [m2/s2] (PUBLIC for TAM)
41	REAL(wp), PUBLIC :: rn_bfrien = 30._wp ! local factor to enhance coefficient bfri (PUBLIC for TAM)
42	LOGICAL , PUBLIC :: ln_bfr2d = .false. ! logical switch for 2D enhancement (PUBLIC for TAM)
43	LOGICAL , PUBLIC :: ln_loglayer = .false. ! switch for log layer bfr coeff. (PUBLIC for TAM)
44	REAL(wp), PUBLIC :: rn_bfrz0 = 0.003_wp ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
45	LOGICAL , PUBLIC :: ln_bfrimp = .false. ! 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 :: ikbu, ikbv ! local integers
89	REAL(wp) :: zvu, zuv, zecu, zecv ! temporary scalars
90	REAL(wp) :: ztmp, ztmp1 ! temporary scalars
91	!!----------------------------------------------------------------------
92	!
93	IF( nn_timing == 1 ) CALL timing_start('zdf_bfr')
94	!
95	IF( nn_bfr == 2 ) THEN ! quadratic botton friction
96	! Calculate and store the quadratic bottom friction coefficient bfrua and bfrva
97	! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]}
98	! from these the trend due to bottom friction: -F_h/e3U can be calculated
99	! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]}
100	!
101
102	IF(ln_loglayer) THEN ! "log layer" bottom friction coefficient
103
104	! add 2D-enhancement bottom friction
105	ztmp1 = 1._wp
106	IF(ABS(rn_bfri2) >= 1.e-10 ) THEN
107	ztmp1 = 1._wp / rn_bfri2
108	ELSE
109	CALL ctl_stop( 'rn_bfri2 must not be less than 1.e-10')
110	END IF
111
112	# if defined key_vectopt_loop
113	DO jj = 1, 1
114	DO ji = 1, jpij ! vector opt. (forced unrolling)
115	# else
116	DO jj = 1, jpj
117	DO ji = 1, jpi
118	# endif
119	ztmp = 0.5_wp * fse3t(ji,jj,mbkt(ji,jj))
120	ztmp = max(ztmp, rn_bfrz0 + 1.e-10)
121	bfrcoef2d(ji,jj) = bfrcoef2d(ji,jj) * ztmp1 * &
122	& ( log( ztmp / rn_bfrz0 ) / vkarmn ) ** (-2)
123	END DO
124	END DO
125	ENDIF
126
127	# if defined key_vectopt_loop
128	DO jj = 1, 1
129	!CDIR NOVERRCHK
130	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
131	# else
132	!CDIR NOVERRCHK
133	DO jj = 2, jpjm1
134	!CDIR NOVERRCHK
135	DO ji = 2, jpim1
136	# endif
137	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
138	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
139	!
140	zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) &
141	& + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) )
142	zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) &
143	& + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) )
144	!
145	zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 )
146	zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 )
147	!
148	bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj ) ) * zecu
149	bfrva(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji ,jj+1) ) * zecv
150	END DO
151	END DO
152
153	!
154	CALL lbc_lnk( bfrua, 'U', 1. ) ; CALL lbc_lnk( bfrva, 'V', 1. ) ! Lateral boundary condition
155	!
156	IF(ln_ctl) CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr - u: ', mask1=umask, &
157	& tab2d_2=bfrva, clinfo2= ' v: ', mask2=vmask,ovlap=1 )
158	ENDIF
159
160	!
161	IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr')
162	!
163	END SUBROUTINE zdf_bfr
164
165
166	SUBROUTINE zdf_bfr_init
167	!!----------------------------------------------------------------------
168	!! * ROUTINE zdf_bfr_init *
169	!!
170	!! ** Purpose : Initialization of the bottom friction
171	!!
172	!! ** Method : Read the nammbf namelist and check their consistency
173	!! called at the first timestep (nit000)
174	!!----------------------------------------------------------------------
175	USE iom ! I/O module for ehanced bottom friction file
176	!!
177	INTEGER :: inum ! logical unit for enhanced bottom friction file
178	INTEGER :: ji, jj ! dummy loop indexes
179	INTEGER :: ikbu, ikbv ! temporary integers
180	INTEGER :: ictu, ictv ! - -
181	REAL(wp) :: zminbfr, zmaxbfr ! temporary scalars
182	REAL(wp) :: zfru, zfrv ! - -
183	!!
184	NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
185	& rn_bfrien, ln_bfrimp, ln_loglayer
186	!!----------------------------------------------------------------------
187	!
188	IF( nn_timing == 1 ) CALL timing_start('zdf_bfr_init')
189	!
190	REWIND ( numnam ) !* Read Namelist nam_bfr : bottom momentum boundary condition
191	READ ( numnam, nambfr )
192
193	! !* Parameter control and print
194	IF(lwp) WRITE(numout,*)
195	IF(lwp) WRITE(numout,*) 'zdf_bfr : momentum bottom friction'
196	IF(lwp) WRITE(numout,*) '~~~~~~~'
197	IF(lwp) WRITE(numout,*) ' Namelist nam_bfr : set bottom friction parameters'
198
199	! ! allocate zdfbfr arrays
200	IF( zdf_bfr_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' )
201
202	! ! Make sure ln_zdfexp=.false. when use implicit bfr
203	IF( ln_bfrimp .AND. ln_zdfexp ) THEN
204	IF(lwp) THEN
205	WRITE(numout,*)
206	WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.'
207	WRITE(numout,*) ' but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!'
208	WRITE(ctmp1,*) ' set either ln_zdfexp = .false or ln_bfrimp = .false.'
209	CALL ctl_stop( ctmp1 )
210	END IF
211	END IF
212
213	SELECT CASE (nn_bfr)
214	!
215	CASE( 0 )
216	IF(lwp) WRITE(numout,*) ' free-slip '
217	bfrua(:,:) = 0.e0
218	bfrva(:,:) = 0.e0
219	!
220	CASE( 1 )
221	IF(lwp) WRITE(numout,*) ' linear botton friction'
222	IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri1 = ', rn_bfri1
223	IF( ln_bfr2d ) THEN
224	IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_bfr2d = ', ln_bfr2d
225	IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien
226	ENDIF
227	!
228	bfrcoef2d(:,:) = rn_bfri1 ! initialize bfrcoef2d to the namelist variable
229	!
230	IF(ln_bfr2d) THEN
231	! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
232	CALL iom_open('bfr_coef.nc',inum)
233	CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
234	CALL iom_close(inum)
235	bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
236	ENDIF
237	bfrua(:,:) = - bfrcoef2d(:,:)
238	bfrva(:,:) = - bfrcoef2d(:,:)
239	!
240	CASE( 2 )
241	IF(lwp) WRITE(numout,*) ' quadratic botton friction'
242	IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri2 = ', rn_bfri2
243	IF(lwp) WRITE(numout,*) ' background tke rn_bfeb2 = ', rn_bfeb2
244	IF( ln_bfr2d ) THEN
245	IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_bfr2d = ', ln_bfr2d
246	IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien
247	ENDIF
248	bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable
249
250	!
251	IF(ln_bfr2d) THEN
252	! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
253	CALL iom_open('bfr_coef.nc',inum)
254	CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
255	CALL iom_close(inum)
256	bfrcoef2d(:,:)= rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
257	ENDIF
258	!
259	CASE DEFAULT
260	IF(lwp) WRITE(ctmp1,*) ' bad flag value for nn_bfr = ', nn_bfr
261	CALL ctl_stop( ctmp1 )
262	!
263	END SELECT
264
265	IF( nn_bfr /= 2 .AND. ln_loglayer ) THEN
266	IF(lwp) THEN
267	WRITE(numout,*)
268	WRITE(numout,*) 'Loglayer can only be by applied for quadratic bottom friction'
269	WRITE(numout,*) 'but you have set: nn_bfr /= 2 and ln_loglayer=.true.!!!!'
270	WRITE(ctmp1,*) 'check nn_bfr and ln_loglayer (should be 2 and true)'
271	CALL ctl_stop( ctmp1 )
272	END IF
273	END IF
274
275
276
277	IF(lwp) WRITE(numout,*) ' implicit bottom friction switch ln_bfrimp = ', ln_bfrimp
278	!
279	! Basic stability check on bottom friction coefficient
280	!
281	ictu = 0 ! counter for stability criterion breaches at U-pts
282	ictv = 0 ! counter for stability criterion breaches at V-pts
283	zminbfr = 1.e10_wp ! initialise tracker for minimum of bottom friction coefficient
284	zmaxbfr = -1.e10_wp ! initialise tracker for maximum of bottom friction coefficient
285	!
286	# if defined key_vectopt_loop
287	DO jj = 1, 1
288	!CDIR NOVERRCHK
289	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
290	# else
291	!CDIR NOVERRCHK
292	DO jj = 2, jpjm1
293	!CDIR NOVERRCHK
294	DO ji = 2, jpim1
295	# endif
296	ikbu = mbku(ji,jj) ! deepest ocean level at u- and v-points
297	ikbv = mbkv(ji,jj)
298	zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
299	zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
300	IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN
301	IF( ln_ctl ) THEN
302	WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
303	WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru
304	ENDIF
305	ictu = ictu + 1
306	ENDIF
307	IF( ABS( bfrcoef2d(ji,jj) ) > zfrv ) THEN
308	IF( ln_ctl ) THEN
309	WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv
310	WRITE(numout,*) 'BFR ', bfrcoef2d(ji,jj), zfrv
311	ENDIF
312	ictv = ictv + 1
313	ENDIF
314	zminbfr = MIN( zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) ) )
315	zmaxbfr = MAX( zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) ) )
316	END DO
317	END DO
318	IF( lk_mpp ) THEN
319	CALL mpp_sum( ictu )
320	CALL mpp_sum( ictv )
321	CALL mpp_min( zminbfr )
322	CALL mpp_max( zmaxbfr )
323	ENDIF
324	IF( .NOT.ln_bfrimp) THEN
325	IF( lwp .AND. ictu + ictv > 0 ) THEN
326	WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
327	WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
328	WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
329	WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
330	ENDIF
331	ENDIF
332	!
333	IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr_init')
334	!
335	END SUBROUTINE zdf_bfr_init
336
337	!!======================================================================
338	END MODULE zdfbfr

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

Download in other formats: