Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

bdydyn2d.F90

Last change on this file was 15368, checked in by smasson, 3 years ago
trunk: final version (hopefully) for ticket #2731
Property svn:keywords set to `Id`
File size: 18.4 KB

Rev	Line
[3117]	1	MODULE bdydyn2d
	2	!!======================================================================
	3	!! * MODULE bdydyn *
[3191]	4	!! Unstructured Open Boundary Cond. : Apply boundary conditions to barotropic solution
[3117]	5	!!======================================================================
[3191]	6	!! History : 3.4 ! 2011 (D. Storkey) new module as part of BDY rewrite
[3680]	7	!! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Optimization of BDY communications
[4292]	8	!! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes
[3117]	9	!!----------------------------------------------------------------------
[4292]	10	!! bdy_dyn2d : Apply open boundary conditions to barotropic variables.
	11	!! bdy_dyn2d_frs : Apply Flow Relaxation Scheme
	12	!! bdy_dyn2d_fla : Apply Flather condition
	13	!! bdy_dyn2d_orlanski : Orlanski Radiation
	14	!! bdy_ssh : Duplicate sea level across open boundaries
[3117]	15	!!----------------------------------------------------------------------
	16	USE dom_oce ! ocean space and time domain
	17	USE bdy_oce ! ocean open boundary conditions
[4292]	18	USE bdylib ! BDY library routines
[3117]	19	USE phycst ! physical constants
[15363]	20	USE lib_mpp
[3117]	21	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[9023]	22	USE wet_dry ! Use wet dry to get reference ssh level
[3117]	23	USE in_out_manager !
	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
[4292]	28	PUBLIC bdy_dyn2d ! routine called in dynspg_ts and bdy_dyn
	29	PUBLIC bdy_ssh ! routine called in dynspg_ts or sshwzv
[3117]	30
	31	!!----------------------------------------------------------------------
[9598]	32	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]	33	!! $Id$
[10068]	34	!! Software governed by the CeCILL license (see ./LICENSE)
[3117]	35	!!----------------------------------------------------------------------
	36	CONTAINS
	37
[4354]	38	SUBROUTINE bdy_dyn2d( kt, pua2d, pva2d, pub2d, pvb2d, phur, phvr, pssh )
[3117]	39	!!----------------------------------------------------------------------
	40	!! * SUBROUTINE bdy_dyn2d *
	41	!!
	42	!! ** Purpose : - Apply open boundary conditions for barotropic variables
	43	!!
	44	!!----------------------------------------------------------------------
	45	INTEGER, INTENT(in) :: kt ! Main time step counter
[4354]	46	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
	47	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pub2d, pvb2d
	48	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: phur, phvr
	49	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh
[3117]	50	!!
[15354]	51	INTEGER :: ib_bdy, ir ! BDY set index, rim index
[15368]	52	INTEGER, DIMENSION(3) :: idir3
	53	INTEGER, DIMENSION(6) :: idir6
[15354]	54	LOGICAL :: llrim0 ! indicate if rim 0 is treated
	55	LOGICAL, DIMENSION(8) :: llsend2, llrecv2, llsend3, llrecv3 ! indicate how communications are to be carried out
	56	!!----------------------------------------------------------------------
[11536]	57
	58	llsend2(:) = .false. ; llrecv2(:) = .false.
	59	llsend3(:) = .false. ; llrecv3(:) = .false.
	60	DO ir = 1, 0, -1 ! treat rim 1 before rim 0
	61	IF( ir == 0 ) THEN ; llrim0 = .TRUE.
	62	ELSE ; llrim0 = .FALSE.
	63	END IF
	64	DO ib_bdy=1, nb_bdy
	65	SELECT CASE( cn_dyn2d(ib_bdy) )
	66	CASE('none')
	67	CYCLE
	68	CASE('frs') ! treat the whole boundary at once
	69	IF( llrim0 ) CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d )
	70	CASE('flather')
	71	CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr, llrim0 )
	72	CASE('orlanski')
	73	CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
	74	& pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo=.false. )
	75	CASE('orlanski_npo')
	76	CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
	77	& pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo=.true. )
	78	CASE DEFAULT
	79	CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
	80	END SELECT
	81	ENDDO
	82	!
	83	IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
	84	IF( nn_hls == 1 ) THEN
	85	llsend2(:) = .false. ; llrecv2(:) = .false.
	86	llsend3(:) = .false. ; llrecv3(:) = .false.
	87	END IF
	88	DO ib_bdy=1, nb_bdy
	89	SELECT CASE( cn_dyn2d(ib_bdy) )
	90	CASE('flather')
[15368]	91	idir6 = (/ jpwe, jpea, jpsw, jpse, jpnw, jpne /)
	92	llsend2(idir6) = llsend2(idir6) .OR. lsend_bdyint(ib_bdy,2,idir6,ir) ! west/east, U points
	93	idir3 = (/ jpwe, jpsw, jpnw /)
	94	llsend2(idir3) = llsend2(idir3) .OR. lsend_bdyext(ib_bdy,2,idir3,ir) ! nei might search point towards its east bdy
	95	llrecv2(idir6) = llrecv2(idir6) .OR. lrecv_bdyint(ib_bdy,2,idir6,ir) ! west/east, U points
	96	idir3 = (/ jpea, jpse, jpne /)
	97	llrecv2(idir3) = llrecv2(idir3) .OR. lrecv_bdyext(ib_bdy,2,idir3,ir) ! might search point towards bdy on the east
	98	idir6 = (/ jpso, jpno, jpsw, jpse, jpnw, jpne /)
	99	llsend3(idir6) = llsend3(idir6) .OR. lsend_bdyint(ib_bdy,3,idir6,ir) ! north/south, V points
	100	idir3 = (/ jpso, jpsw, jpse /)
	101	llsend3(idir3) = llsend3(idir3) .OR. lsend_bdyext(ib_bdy,3,idir3,ir) ! nei might search point towards its north bdy
	102	llrecv3(idir6) = llrecv3(idir6) .OR. lrecv_bdyint(ib_bdy,3,idir6,ir) ! north/south, V points
	103	idir3 = (/ jpno, jpnw, jpne /)
	104	llrecv3(idir3) = llrecv3(idir3) .OR. lrecv_bdyext(ib_bdy,3,idir3,ir) ! might search point towards bdy on the north
[11536]	105	CASE('orlanski', 'orlanski_npo')
[15354]	106	llsend2(:) = llsend2(:) .OR. lsend_bdyolr(ib_bdy,2,:,ir) ! possibly every direction, U points
	107	llrecv2(:) = llrecv2(:) .OR. lrecv_bdyolr(ib_bdy,2,:,ir) ! possibly every direction, U points
	108	llsend3(:) = llsend3(:) .OR. lsend_bdyolr(ib_bdy,3,:,ir) ! possibly every direction, V points
	109	llrecv3(:) = llrecv3(:) .OR. lrecv_bdyolr(ib_bdy,3,:,ir) ! possibly every direction, V points
[11536]	110	END SELECT
	111	END DO
	112	IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN ! if need to send/recv in at least one direction
[13226]	113	CALL lbc_lnk( 'bdydyn2d', pua2d, 'U', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
[11536]	114	END IF
	115	IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN ! if need to send/recv in at least one direction
[13226]	116	CALL lbc_lnk( 'bdydyn2d', pva2d, 'V', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
[11536]	117	END IF
	118	!
	119	END DO ! ir
	120	!
[3117]	121	END SUBROUTINE bdy_dyn2d
	122
[4354]	123	SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy, pua2d, pva2d )
[3117]	124	!!----------------------------------------------------------------------
	125	!! * SUBROUTINE bdy_dyn2d_frs *
	126	!!
	127	!! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
	128	!! at open boundaries.
	129	!!
	130	!! References :- Engedahl H., 1995: Use of the flow relaxation scheme in
	131	!! a three-dimensional baroclinic ocean model with realistic
	132	!! topography. Tellus, 365-382.
	133	!!----------------------------------------------------------------------
	134	TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices
	135	TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data
[3680]	136	INTEGER, INTENT(in) :: ib_bdy ! BDY set index
[4354]	137	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
[3117]	138	!!
[11536]	139	INTEGER :: jb ! dummy loop indices
[3117]	140	INTEGER :: ii, ij, igrd ! local integers
	141	REAL(wp) :: zwgt ! boundary weight
	142	!!----------------------------------------------------------------------
	143	!
	144	igrd = 2 ! Relaxation of zonal velocity
	145	DO jb = 1, idx%nblen(igrd)
	146	ii = idx%nbi(jb,igrd)
	147	ij = idx%nbj(jb,igrd)
	148	zwgt = idx%nbw(jb,igrd)
[4292]	149	pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
[3117]	150	END DO
	151	!
	152	igrd = 3 ! Relaxation of meridional velocity
	153	DO jb = 1, idx%nblen(igrd)
	154	ii = idx%nbi(jb,igrd)
	155	ij = idx%nbj(jb,igrd)
	156	zwgt = idx%nbw(jb,igrd)
[4292]	157	pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
[3117]	158	END DO
	159	!
	160	END SUBROUTINE bdy_dyn2d_frs
	161
	162
[11536]	163	SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr, llrim0 )
[3117]	164	!!----------------------------------------------------------------------
	165	!! * SUBROUTINE bdy_dyn2d_fla *
	166	!!
	167	!! - Apply Flather boundary conditions on normal barotropic velocities
	168	!!
	169	!! ** WARNINGS about FLATHER implementation:
	170	!!1. According to Palma and Matano, 1998 "after ssh" is used.
	171	!! In ROMS and POM implementations, it is "now ssh". In the current
	172	!! implementation (tested only in the EEL-R5 conf.), both cases were unstable.
	173	!! So I use "before ssh" in the following.
	174	!!
	175	!!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of
	176	!! fact, the model ssh just inside the dynamical boundary is used (the outside
	177	!! ssh in the code is not updated).
	178	!!
	179	!! References: Flather, R. A., 1976: A tidal model of the northwest European
	180	!! continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.
	181	!!----------------------------------------------------------------------
	182	TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices
	183	TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data
[3680]	184	INTEGER, INTENT(in) :: ib_bdy ! BDY set index
[4354]	185	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
[11536]	186	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh, phur, phvr
	187	LOGICAL , INTENT(in) :: llrim0 ! indicate if rim 0 is treated
	188	INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1)
[3117]	189	INTEGER :: jb, igrd ! dummy loop indices
[11536]	190	INTEGER :: ii, ij ! 2D addresses
	191	INTEGER :: iiTrim, ijTrim ! T pts i/j-indice on the rim
	192	INTEGER :: iiToce, ijToce, iiUoce, ijVoce ! T, U and V pts i/j-indice of the ocean next to the rim
	193	REAL(wp) :: flagu, flagv ! short cuts
	194	REAL(wp) :: zfla ! Flather correction
	195	REAL(wp) :: z1_2 !
	196	REAL(wp), DIMENSION(jpi,jpj) :: sshdta ! 2D version of dta%ssh
[3117]	197	!!----------------------------------------------------------------------
	198
[4292]	199	z1_2 = 0.5_wp
	200
[3117]	201	! ---------------------------------!
	202	! Flather boundary conditions :!
[11536]	203	! ---------------------------------!
[3117]	204
[11536]	205	! Fill temporary array with ssh data (here we use spgu with the alias sshdta):
[3117]	206	igrd = 1
[11536]	207	IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd)
	208	ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd)
	209	END IF
	210	!
	211	DO jb = ibeg, iend
[3117]	212	ii = idx%nbi(jb,igrd)
	213	ij = idx%nbj(jb,igrd)
[11536]	214	IF( ll_wd ) THEN ; sshdta(ii, ij) = dta%ssh(jb) - ssh_ref
	215	ELSE ; sshdta(ii, ij) = dta%ssh(jb)
[9023]	216	ENDIF
[3117]	217	END DO
	218	!
[11536]	219	igrd = 2 ! Flather bc on u-velocity
[3117]	220	! ! remember that flagu=-1 if normal velocity direction is outward
	221	! ! I think we should rather use after ssh ?
[11536]	222	IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd)
	223	ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd)
	224	END IF
	225	DO jb = ibeg, iend
	226	ii = idx%nbi(jb,igrd)
	227	ij = idx%nbj(jb,igrd)
	228	flagu = idx%flagu(jb,igrd)
	229	IF( flagu == 0. ) THEN
	230	pua2d(ii,ij) = dta%u2d(jb)
	231	ELSE ! T pts j-indice on the rim on the ocean next to the rim on T and U points
	232	IF( flagu == 1. ) THEN ; iiTrim = ii ; iiToce = ii+1 ; iiUoce = ii+1 ; ENDIF
	233	IF( flagu == -1. ) THEN ; iiTrim = ii+1 ; iiToce = ii ; iiUoce = ii-1 ; ENDIF
	234	!
	235	! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received
	236	IF( iiTrim > jpi .OR. iiToce > jpi .OR. iiUoce > jpi .OR. iiUoce < 1 ) CYCLE
	237	!
	238	zfla = dta%u2d(jb) - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iiToce,ij) - sshdta(iiTrim,ij) )
	239	!
	240	! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) :
	241	! mix Flather scheme with velocity of the ocean next to the rim
	242	pua2d(ii,ij) = z1_2 * ( pua2d(iiUoce,ij) + zfla )
	243	END IF
[3117]	244	END DO
	245	!
	246	igrd = 3 ! Flather bc on v-velocity
	247	! ! remember that flagv=-1 if normal velocity direction is outward
[11536]	248	IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd)
	249	ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd)
	250	END IF
	251	DO jb = ibeg, iend
	252	ii = idx%nbi(jb,igrd)
	253	ij = idx%nbj(jb,igrd)
	254	flagv = idx%flagv(jb,igrd)
	255	IF( flagv == 0. ) THEN
	256	pva2d(ii,ij) = dta%v2d(jb)
	257	ELSE ! T pts j-indice on the rim on the ocean next to the rim on T and V points
	258	IF( flagv == 1. ) THEN ; ijTrim = ij ; ijToce = ij+1 ; ijVoce = ij+1 ; ENDIF
	259	IF( flagv == -1. ) THEN ; ijTrim = ij+1 ; ijToce = ij ; ijVoce = ij-1 ; ENDIF
	260	!
	261	! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received
	262	IF( ijTrim > jpj .OR. ijToce > jpj .OR. ijVoce > jpj .OR. ijVoce < 1 ) CYCLE
	263	!
	264	zfla = dta%v2d(jb) - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii,ijToce) - sshdta(ii,ijTrim) )
	265	!
	266	! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) :
	267	! mix Flather scheme with velocity of the ocean next to the rim
	268	pva2d(ii,ij) = z1_2 * ( pva2d(ii,ijVoce) + zfla )
	269	END IF
[3117]	270	END DO
	271	!
	272	END SUBROUTINE bdy_dyn2d_fla
[4292]	273
	274
[11536]	275	SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo )
[4292]	276	!!----------------------------------------------------------------------
	277	!! * SUBROUTINE bdy_dyn2d_orlanski *
	278	!!
	279	!! - Apply Orlanski radiation condition adaptively:
	280	!! - radiation plus weak nudging at outflow points
	281	!! - no radiation and strong nudging at inflow points
	282	!!
	283	!!
	284	!! References: Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)
	285	!!----------------------------------------------------------------------
	286	TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices
	287	TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data
	288	INTEGER, INTENT(in) :: ib_bdy ! number of current open boundary set
[4354]	289	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
	290	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d
[4292]	291	LOGICAL, INTENT(in) :: ll_npo ! flag for NPO version
[11536]	292	LOGICAL, INTENT(in) :: llrim0 ! indicate if rim 0 is treated
[4292]	293	INTEGER :: ib, igrd ! dummy loop indices
	294	INTEGER :: ii, ij, iibm1, ijbm1 ! indices
	295	!!----------------------------------------------------------------------
	296	!
	297	igrd = 2 ! Orlanski bc on u-velocity;
	298	!
[11536]	299	CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, llrim0, ll_npo )
[4292]	300
	301	igrd = 3 ! Orlanski bc on v-velocity
	302	!
[11536]	303	CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, llrim0, ll_npo )
[4292]	304	!
	305	END SUBROUTINE bdy_dyn2d_orlanski
	306
[9124]	307
[4292]	308	SUBROUTINE bdy_ssh( zssh )
	309	!!----------------------------------------------------------------------
	310	!! * SUBROUTINE bdy_ssh *
	311	!!
	312	!! ** Purpose : Duplicate sea level across open boundaries
	313	!!
	314	!!----------------------------------------------------------------------
[11536]	315	REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) :: zssh ! Sea level, need 3 dimensions to be used by bdy_nmn
[4292]	316	!!
[11536]	317	INTEGER :: ib_bdy, ir ! bdy index, rim index
	318	INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1)
	319	LOGICAL :: llrim0 ! indicate if rim 0 is treated
[15354]	320	LOGICAL, DIMENSION(8) :: llsend1, llrecv1 ! indicate how communications are to be carried out
[11536]	321	!!----------------------------------------------------------------------
	322	llsend1(:) = .false. ; llrecv1(:) = .false.
	323	DO ir = 1, 0, -1 ! treat rim 1 before rim 0
	324	IF( nn_hls == 1 ) THEN ; llsend1(:) = .false. ; llrecv1(:) = .false. ; END IF
	325	IF( ir == 0 ) THEN ; llrim0 = .TRUE.
	326	ELSE ; llrim0 = .FALSE.
	327	END IF
	328	DO ib_bdy = 1, nb_bdy
	329	CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh, llrim0 ) ! zssh is masked
	330	llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:,ir) ! possibly every direction, T points
	331	llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(ib_bdy,1,:,ir) ! possibly every direction, T points
[4292]	332	END DO
[11536]	333	IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
	334	IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN ! if need to send/recv in at least one direction
[13226]	335	CALL lbc_lnk( 'bdydyn2d', zssh(:,:,1), 'T', 1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
[11536]	336	END IF
[4292]	337	END DO
[11536]	338	!
[4292]	339	END SUBROUTINE bdy_ssh
	340
[3117]	341	!!======================================================================
	342	END MODULE bdydyn2d
[4292]	343

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/OCE/BDY/bdydyn2d.F90

Download in other formats: