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/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90 @ 3820

Last change on this file since 3820 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

Rev	Line
[3]	1	MODULE zdfbfr
	2	!!======================================================================
	3	!! * MODULE zdfbfr *
	4	!! Ocean physics: Bottom friction
	5	!!======================================================================
[1662]	6	!! History : OPA ! 1997-06 (G. Madec, A.-M. Treguier) Original code
[1601]	7	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
[1708]	8	!! 3.2 ! 2009-09 (A.C.Coward) Correction to include barotropic contribution
[2528]	9	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[3633]	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
[1601]	12	!!----------------------------------------------------------------------
[3]	13
	14	!!----------------------------------------------------------------------
[1601]	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.
[2528]	17	!! zdf_bfr_2d : read in namelist and control the bottom friction parameters.
[3]	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and tracers variables
[3764]	20	USE dom_oce ! ocean space and time domain variables
[3]	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)
[1671]	24	USE lib_mpp ! distributed memory computing
[258]	25	USE prtctl ! Print control
[3294]	26	USE timing ! Timing
[3]	27
[3633]	28	USE phycst, ONLY: vkarmn
	29
[3]	30	IMPLICIT NONE
	31	PRIVATE
	32
[2528]	33	PUBLIC zdf_bfr ! called by step.F90
	34	PUBLIC zdf_bfr_init ! called by opa.F90
[2715]	35
[1708]	36	! !!* Namelist nambfr: bottom friction namelist *
[3764]	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)
[3]	47
	48	!! * Substitutions
[2528]	49	# include "vectopt_loop_substitute.h90"
[3]	50	# include "domzgr_substitute.h90"
	51	!!----------------------------------------------------------------------
[2715]	52	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]	53	!! $Id$
	54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	55	!!----------------------------------------------------------------------
	56	CONTAINS
	57
[2715]	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
[3]	69	SUBROUTINE zdf_bfr( kt )
	70	!!----------------------------------------------------------------------
	71	!! * ROUTINE zdf_bfr *
[3764]	72	!!
[1708]	73	!! ** Purpose : compute the bottom friction coefficient.
[3]	74	!!
[3764]	75	!! ** Method : Calculate and store part of the momentum trend due
	76	!! to bottom friction following the chosen friction type
[1662]	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.
[1708]	80	!! The coefficients are updated at each time step only
	81	!! in the quadratic case.
[1662]	82	!!
	83	!! ** Action : bfrua , bfrva bottom friction coefficients
[3]	84	!!----------------------------------------------------------------------
	85	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1601]	86	!!
[3633]	87	INTEGER :: ji, jj ! dummy loop indices
	88	INTEGER :: ikbu, ikbv ! local integers
[3820]	89	REAL(wp) :: zvu, zuv, zecu, zecv ! temporary scalars
	90	REAL(wp) :: ztmp, ztmp1 ! temporary scalars
[3]	91	!!----------------------------------------------------------------------
[3294]	92	!
	93	IF( nn_timing == 1 ) CALL timing_start('zdf_bfr')
	94	!
[1662]	95	IF( nn_bfr == 2 ) THEN ! quadratic botton friction
[1708]	96	! Calculate and store the quadratic bottom friction coefficient bfrua and bfrva
[1662]	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	!
[3633]	101
	102	IF(ln_loglayer) THEN ! "log layer" bottom friction coefficient
[3820]	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
[3633]	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))
[3820]	120	ztmp = max(ztmp, rn_bfrz0 + 1.e-10)
	121	bfrcoef2d(ji,jj) = bfrcoef2d(ji,jj) * ztmp1 * &
	122	& ( log( ztmp / rn_bfrz0 ) / vkarmn ) ** (-2)
[3633]	123	END DO
	124	END DO
	125	ENDIF
	126
[789]	127	# if defined key_vectopt_loop
[1601]	128	DO jj = 1, 1
[3]	129	!CDIR NOVERRCHK
[1601]	130	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3]	131	# else
	132	!CDIR NOVERRCHK
	133	DO jj = 2, jpjm1
	134	!CDIR NOVERRCHK
	135	DO ji = 2, jpim1
	136	# endif
[3764]	137	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
[2528]	138	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
[1662]	139	!
[2528]	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) )
[1662]	144	!
[2528]	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 )
[1662]	147	!
[3764]	148	bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj ) ) * zecu
[2528]	149	bfrva(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji ,jj+1) ) * zecv
[3]	150	END DO
	151	END DO
[3633]	152
[1601]	153	!
[1708]	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 )
[1662]	158	ENDIF
[3294]	159
[1662]	160	!
[3294]	161	IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr')
	162	!
[3]	163	END SUBROUTINE zdf_bfr
	164
	165
	166	SUBROUTINE zdf_bfr_init
	167	!!----------------------------------------------------------------------
	168	!! * ROUTINE zdf_bfr_init *
[3764]	169	!!
[3]	170	!! ** Purpose : Initialization of the bottom friction
	171	!!
	172	!! ** Method : Read the nammbf namelist and check their consistency
[1662]	173	!! called at the first timestep (nit000)
[3]	174	!!----------------------------------------------------------------------
[1662]	175	USE iom ! I/O module for ehanced bottom friction file
	176	!!
[3633]	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 ! - -
[1708]	183	!!
[3633]	184	NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
	185	& rn_bfrien, ln_bfrimp, ln_loglayer
[3]	186	!!----------------------------------------------------------------------
[3294]	187	!
	188	IF( nn_timing == 1 ) CALL timing_start('zdf_bfr_init')
	189	!
[1662]	190	REWIND ( numnam ) !* Read Namelist nam_bfr : bottom momentum boundary condition
	191	READ ( numnam, nambfr )
[3]	192
[1662]	193	! !* Parameter control and print
	194	IF(lwp) WRITE(numout,*)
[3]	195	IF(lwp) WRITE(numout,*) 'zdf_bfr : momentum bottom friction'
	196	IF(lwp) WRITE(numout,*) '~~~~~~~'
[1708]	197	IF(lwp) WRITE(numout,*) ' Namelist nam_bfr : set bottom friction parameters'
[3]	198
[2715]	199	! ! allocate zdfbfr arrays
	200	IF( zdf_bfr_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' )
	201
[3294]	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
[1662]	213	SELECT CASE (nn_bfr)
[2715]	214	!
[3]	215	CASE( 0 )
[1708]	216	IF(lwp) WRITE(numout,*) ' free-slip '
[1662]	217	bfrua(:,:) = 0.e0
	218	bfrva(:,:) = 0.e0
	219	!
[3]	220	CASE( 1 )
[1708]	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
[1662]	226	ENDIF
[1601]	227	!
[1662]	228	bfrcoef2d(:,:) = rn_bfri1 ! initialize bfrcoef2d to the namelist variable
[1708]	229	!
[3764]	230	IF(ln_bfr2d) THEN
[1662]	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)
[2715]	235	bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
[1662]	236	ENDIF
	237	bfrua(:,:) = - bfrcoef2d(:,:)
	238	bfrva(:,:) = - bfrcoef2d(:,:)
	239	!
[3]	240	CASE( 2 )
[1708]	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
[1662]	247	ENDIF
	248	bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable
[3633]	249
[1708]	250	!
[3764]	251	IF(ln_bfr2d) THEN
[1662]	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)
[1708]	256	bfrcoef2d(:,:)= rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
[1662]	257	ENDIF
[1601]	258	!
[3]	259	CASE DEFAULT
[1708]	260	IF(lwp) WRITE(ctmp1,*) ' bad flag value for nn_bfr = ', nn_bfr
[474]	261	CALL ctl_stop( ctmp1 )
[1601]	262	!
[3]	263	END SELECT
[3820]	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
[3294]	277	IF(lwp) WRITE(numout,*) ' implicit bottom friction switch ln_bfrimp = ', ln_bfrimp
[1601]	278	!
[1662]	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
[1708]	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
[1662]	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
[2528]	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
[1708]	300	IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN
	301	IF( ln_ctl ) THEN
[2528]	302	WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
	303	WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru
[1662]	304	ENDIF
	305	ictu = ictu + 1
	306	ENDIF
[1708]	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
[1662]	311	ENDIF
	312	ictv = ictv + 1
	313	ENDIF
[1708]	314	zminbfr = MIN( zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) ) )
	315	zmaxbfr = MAX( zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) ) )
[1662]	316	END DO
	317	END DO
[1708]	318	IF( lk_mpp ) THEN
[1662]	319	CALL mpp_sum( ictu )
	320	CALL mpp_sum( ictv )
[1708]	321	CALL mpp_min( zminbfr )
	322	CALL mpp_max( zmaxbfr )
[1662]	323	ENDIF
[3294]	324	IF( .NOT.ln_bfrimp) THEN
[2528]	325	IF( lwp .AND. ictu + ictv > 0 ) THEN
[1662]	326	WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
	327	WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
[1708]	328	WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
	329	WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
[1662]	330	ENDIF
[3294]	331	ENDIF
[1662]	332	!
[3294]	333	IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr_init')
	334	!
[3]	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: