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.

obcini.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

Context Navigation

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90 @ 2814

Last change on this file since 2814 was 2814, checked in by davestorkey, 13 years ago
Implement tidal harmonics forcing (UKMO version) in new structure. Other bug fixes and updates.
Property svn:keywords set to `Id`
File size: 23.8 KB

Rev	Line
[2797]	1	MODULE obcini
[1601]	2	!!======================================================================
[3]	3	!! * MODULE obcini *
[2797]	4	!! Unstructured open boundaries : initialisation
[1601]	5	!!======================================================================
[2797]	6	!! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code
	7	!! - ! 2007-01 (D. Storkey) Update to use IOM module
	8	!! - ! 2007-01 (D. Storkey) Tidal forcing
	9	!! 3.0 ! 2008-04 (NEMO team) add in the reference version
	10	!! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations
	11	!! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions
	12	!! 3.4 ! 2011 (D. Storkey, J. Chanut) OBC-BDY merge
	13	!! ! --- Renamed bdyini.F90 -> obcini.F90 ---
[1601]	14	!!----------------------------------------------------------------------
[25]	15	#if defined key_obc
[1601]	16	!!----------------------------------------------------------------------
[2797]	17	!! 'key_obc' Unstructured Open Boundary Conditions
[1601]	18	!!----------------------------------------------------------------------
[2797]	19	!! obc_init : Initialization of unstructured open boundaries
[1601]	20	!!----------------------------------------------------------------------
[3]	21	USE oce ! ocean dynamics and tracers variables
[2797]	22	USE dom_oce ! ocean space and time domain
	23	USE obc_oce ! unstructured open boundary conditions
	24	USE in_out_manager ! I/O units
[3]	25	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[2797]	26	USE lib_mpp ! for mpp_sum
	27	USE iom ! I/O
[3]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[1601]	32	PUBLIC obc_init ! routine called by opa.F90
[3]	33
[1601]	34	!!----------------------------------------------------------------------
[2797]	35	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]	36	!! $Id$
[2715]	37	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1601]	38	!!----------------------------------------------------------------------
[3]	39	CONTAINS
	40
	41	SUBROUTINE obc_init
	42	!!----------------------------------------------------------------------
	43	!! * ROUTINE obc_init *
	44	!!
[2797]	45	!! ** Purpose : Initialization of the dynamics and tracer fields with
	46	!! unstructured open boundaries.
[3]	47	!!
[2797]	48	!! ** Method : Read initialization arrays (mask, indices) to identify
	49	!! an unstructured open boundary
[3]	50	!!
[2797]	51	!! ** Input : obc_init.nc, input file for unstructured open boundaries
	52	!!----------------------------------------------------------------------
	53	INTEGER :: ib_obc, ii, ij, ik, igrd, ib, ir ! dummy loop indices
	54	INTEGER :: icount, icountr, ibr_max, ilen1 ! local integers
	55	INTEGER :: iw, ie, is, in, inum, id_dummy ! - -
	56	INTEGER :: igrd_start, igrd_end, jpbdta ! - -
	57	INTEGER, POINTER :: nbi, nbj, nbr ! short cuts
	58	REAL , POINTER :: flagu, flagv ! - -
	59	REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
	60	INTEGER, DIMENSION (2) :: kdimsz
	61	INTEGER, DIMENSION(jpbgrd,jp_obc) :: nblendta ! Length of index arrays
	62	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of obc dta
	63	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
	64	REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask
	65	CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile
	66	CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid
[1601]	67	!!
[2797]	68	NAMELIST/namobc/ nb_obc, ln_coords_file, cn_coords_file, &
	69	& ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn3d, &
	70	& nn_tra, &
	71	#if defined key_lim2
	72	& nn_ice_lim2, &
	73	#endif
[2814]	74	& nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl, &
[2797]	75	& nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out, &
	76	& nn_dmp3d_in, nn_dmp3d_out
[3]	77	!!----------------------------------------------------------------------
	78
[2797]	79	IF( obc_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'obc_init : unable to allocate oce arrays' )
[3]	80
	81	IF(lwp) WRITE(numout,*)
[1601]	82	IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries'
	83	IF(lwp) WRITE(numout,*) '~~~~~~~~'
[2797]	84	!
[3]	85
[2797]	86	IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', &
	87	& ' and general open boundary condition are not compatible' )
[1151]	88
[2797]	89	cgrid= (/'T','U','V'/)
[3]	90
[2797]	91	! -----------------------------------------
	92	! Initialise and read namelist parameters
	93	! -----------------------------------------
[3]	94
[2797]	95	nb_obc = 0
	96	ln_coords_file(:) = .false.
	97	cn_coords_file(:) = ''
	98	ln_mask_file = .false.
	99	cn_mask_file(:) = ''
	100	nn_dyn2d(:) = 0
	101	nn_dyn3d(:) = 0
	102	nn_tra(:) = 0
	103	#if defined key_lim2
	104	nn_ice_lim2(:) = 0
	105	#endif
	106	ln_vol = .false.
	107	ln_clim(:) = .false.
	108	nn_dtactl(:) = -1 ! uninitialised flag
[2814]	109	nn_tides(:) = 0 ! default to no tidal forcing
[2797]	110	nn_volctl = -1 ! uninitialised flag
	111	nn_rimwidth(:) = -1 ! uninitialised flag
	112	nn_dmp2d_in(:) = -1 ! uninitialised flag
	113	nn_dmp2d_out(:) = -1 ! uninitialised flag
	114	nn_dmp3d_in(:) = -1 ! uninitialised flag
	115	nn_dmp3d_out(:) = -1 ! uninitialised flag
[3]	116
[2797]	117	REWIND( numnam )
	118	READ ( numnam, namobc )
[3]	119
[2797]	120	! -----------------------------------------
	121	! Check and write out namelist parameters
	122	! -----------------------------------------
[1601]	123
[2797]	124	! ! control prints
	125	IF(lwp) WRITE(numout,*) ' namobc'
[1151]	126
[2797]	127	IF( nb_obc .eq. 0 ) THEN
	128	IF(lwp) WRITE(numout,*) 'nb_obc = 0, NO OPEN BOUNDARIES APPLIED.'
	129	ELSE
	130	IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_obc
	131	ENDIF
[3]	132
[2797]	133	DO ib_obc = 1,nb_obc
	134	IF(lwp) WRITE(numout,*) ' '
	135	IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------'
[3]	136
[2797]	137	! ! check type of data used (nn_dtactl value)
	138	SELECT CASE( nn_dtactl(ib_obc) ) !
	139	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for obc data'
	140	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
	141	CASE DEFAULT ; CALL ctl_stop( 'nn_dtactl must be 0 or 1' )
	142	END SELECT
	143	IF(lwp) WRITE(numout,*)
[3]	144
[2797]	145	IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '
	146	SELECT CASE( nn_dyn2d(ib_obc) )
	147	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
	148	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
	149	CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition'
	150	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
	151	END SELECT
	152	IF(lwp) WRITE(numout,*)
[3]	153
[2797]	154	IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '
	155	SELECT CASE( nn_dyn3d(ib_obc) )
	156	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
	157	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
	158	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
	159	END SELECT
	160	IF(lwp) WRITE(numout,*)
[3]	161
[2797]	162	IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '
	163	SELECT CASE( nn_tra(ib_obc) )
	164	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
	165	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
	166	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' )
	167	END SELECT
	168	IF(lwp) WRITE(numout,*)
[3]	169
[2797]	170	#if defined key_lim2
	171	IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '
	172	SELECT CASE( nn_tra(ib_obc) )
	173	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
	174	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
	175	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' )
	176	END SELECT
	177	IF(lwp) WRITE(numout,*)
	178	#endif
[3]	179
[2797]	180	IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth
	181	IF(lwp) WRITE(numout,*)
[3]	182
[2814]	183	SELECT CASE( nn_tides(ib_obc) )
	184	CASE(0)
	185	IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing'
[2797]	186	IF(lwp) WRITE(numout,*)
[2814]	187	CASE(1)
	188	IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution'
	189	IF(lwp) WRITE(numout,*)
	190	CASE(2)
	191	IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions'
	192	IF(lwp) WRITE(numout,*)
	193	CASE DEFAULT
	194	CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' )
	195	END SELECT
[3]	196
[2797]	197	ENDDO
[3]	198
[2797]	199	IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)
	200	IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
	201	IF(lwp) WRITE(numout,*)
	202	SELECT CASE ( nn_volctl )
	203	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant'
	204	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'
	205	CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )
	206	END SELECT
	207	IF(lwp) WRITE(numout,*)
	208	ELSE
	209	IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
	210	IF(lwp) WRITE(numout,*)
	211	ENDIF
[3]	212
[2797]	213	! -------------------------------------------------
	214	! Initialise indices arrays for open boundaries
	215	! -------------------------------------------------
[3]	216
[2797]	217	! Work out global dimensions of boundary data
	218	! ---------------------------------------------
	219	DO ib_obc = 1, nb_obc
[3]	220
[2797]	221	jpbdta = 1
	222	IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters
	223
[3]	224
[2797]	225	!! 1. Read parameters from namelist
	226	!! 2. Work out global size of boundary data arrays nblendta and jpbdta
[3]	227
	228
[2797]	229	ELSE ! Read size of arrays in boundary coordinates file.
[3]	230
	231
[2797]	232	CALL iom_open( cn_coords_file(ib_obc), inum )
	233	jpbdta = 1
	234	DO igrd = 1, jpbgrd
	235	id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )
	236	nblendta(igrd,ib_obc) = kdimsz(1)
	237	jpbdta = MAX(jpbdta, kdimsz(1))
	238	ENDDO
[3]	239
[2797]	240	ENDIF
[3]	241
[2797]	242	ENDDO
[3]	243
[2797]	244	! Allocate arrays
	245	!---------------
	246	ALLOCATE( nbidta(jpbdta, jpbgrd, nb_obc), nbjdta(jpbdta, jpbgrd, nb_obc), &
	247	& nbrdta(jpbdta, jpbgrd, nb_obc) )
[3]	248
[2797]	249	ALLOCATE( dta_global(jpbdta, 1, jpk) )
[3]	250
[2797]	251	! Calculate global boundary index arrays or read in from file
	252	!------------------------------------------------------------
	253	DO ib_obc = 1, nb_obc
[3]	254
[2797]	255	IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters
[2149]	256
[2797]	257	!! Calculate global index arrays from namelist parameters
[2149]	258
[2797]	259	ELSE ! Read global index arrays from boundary coordinates file.
[2149]	260
[2797]	261	DO igrd = 1, jpbgrd
	262	CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) )
	263	DO ii = 1,nblendta(igrd,ib_obc)
	264	nbidta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) )
	265	END DO
	266	CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) )
	267	DO ii = 1,nblendta(igrd,ib_obc)
	268	nbjdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) )
	269	END DO
	270	CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) )
	271	DO ii = 1,nblendta(igrd,ib_obc)
	272	nbrdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) )
	273	END DO
[2149]	274
[2797]	275	ibr_max = MAXVAL( nbrdta(:,igrd,ib_obc) )
	276	IF(lwp) WRITE(numout,*)
	277	IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
	278	IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_obc)
	279	IF (ibr_max < nn_rimwidth(ib_obc)) &
	280	CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_obc) )
[2528]	281
[1528]	282	END DO
[2797]	283	CALL iom_close( inum )
[2149]	284
[2797]	285	ENDIF
[3]	286
[2797]	287	ENDDO
[2149]	288
[2797]	289	! Work out dimensions of boundary data on each processor
	290	! ------------------------------------------------------
	291
	292	iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
	293	ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
	294	is = mjg(1) + 1 ! if monotasking and no zoom, is=2
	295	in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
[3]	296
[2797]	297	DO ib_obc = 1, nb_obc
	298	DO igrd = 1, jpbgrd
	299	icount = 0
	300	icountr = 0
	301	idx_obc(ib_obc)%nblen(igrd) = 0
	302	idx_obc(ib_obc)%nblenrim(igrd) = 0
	303	DO ib = 1, nblendta(igrd,ib_obc)
	304	! check if point is in local domain
	305	IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. &
	306	& nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in ) THEN
	307	!
	308	icount = icount + 1
	309	!
	310	IF( nbrdta(ib,igrd,ib_obc) == 1 ) icountr = icountr+1
	311	ENDIF
	312	ENDDO
	313	idx_obc(ib_obc)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
	314	idx_obc(ib_obc)%nblen (igrd) = icount !: length of boundary data on each proc
	315	ENDDO ! igrd
[3]	316
[2797]	317	! Allocate index arrays for this boundary set
	318	!--------------------------------------------
	319	ilen1 = MAXVAL(idx_obc(ib_obc)%nblen(:))
	320	ALLOCATE( idx_obc(ib_obc)%nbi(ilen1,jpbgrd) )
	321	ALLOCATE( idx_obc(ib_obc)%nbj(ilen1,jpbgrd) )
	322	ALLOCATE( idx_obc(ib_obc)%nbr(ilen1,jpbgrd) )
	323	ALLOCATE( idx_obc(ib_obc)%nbmap(ilen1,jpbgrd) )
	324	ALLOCATE( idx_obc(ib_obc)%nbw(ilen1,jpbgrd) )
	325	ALLOCATE( idx_obc(ib_obc)%flagu(ilen1) )
	326	ALLOCATE( idx_obc(ib_obc)%flagv(ilen1) )
[25]	327
[2797]	328	! Dispatch mapping indices and discrete distances on each processor
	329	! -----------------------------------------------------------------
[3]	330
[2797]	331	DO igrd = 1, jpbgrd
	332	icount = 0
	333	! Loop on rimwidth to ensure outermost points come first in the local arrays.
	334	DO ir=1, nn_rimwidth(ib_obc)
	335	DO ib = 1, nblendta(igrd,ib_obc)
	336	! check if point is in local domain and equals ir
	337	IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. &
	338	& nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in .AND. &
	339	& nbrdta(ib,igrd,ib_obc) == ir ) THEN
	340	!
	341	icount = icount + 1
	342	idx_obc(ib_obc)%nbi(icount,igrd) = nbidta(ib,igrd,ib_obc)- mig(1)+1
	343	idx_obc(ib_obc)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_obc)- mjg(1)+1
	344	idx_obc(ib_obc)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_obc)
	345	idx_obc(ib_obc)%nbmap(icount,igrd) = ib
	346	ENDIF
	347	ENDDO
	348	ENDDO
	349	ENDDO
[25]	350
[2814]	351	! Compute rim weights for FRS scheme
	352	! ----------------------------------
[2797]	353	DO igrd = 1, jpbgrd
	354	DO ib = 1, idx_obc(ib_obc)%nblen(igrd)
	355	nbr => idx_obc(ib_obc)%nbr(ib,igrd)
	356	idx_obc(ib_obc)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation
	357	! idx_obc(ib_obc)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic
	358	! idx_obc(ib_obc)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear
	359	END DO
	360	END DO
[3]	361
[2797]	362	ENDDO
[25]	363
[2797]	364	! ------------------------------------------------------
	365	! Initialise masks and find normal/tangential directions
	366	! ------------------------------------------------------
[3]	367
[2797]	368	! Read global 2D mask at T-points: obctmask
	369	! -----------------------------------------
	370	! obctmask = 1 on the computational domain AND on open boundaries
	371	! = 0 elsewhere
	372
	373	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution
	374	zmask( : ,:) = 0.e0
	375	zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0
	376	ELSE IF( ln_mask_file ) THEN
	377	CALL iom_open( cn_mask_file, inum )
	378	CALL iom_get ( inum, jpdom_data, 'obc_msk', zmask(:,:) )
	379	CALL iom_close( inum )
	380	ELSE
	381	zmask(:,:) = 1.e0
	382	ENDIF
	383
	384	DO ij = 1, nlcj ! Save mask over local domain
	385	DO ii = 1, nlci
	386	obctmask(ii,ij) = zmask( mig(ii), mjg(ij) )
[3]	387	END DO
[2797]	388	END DO
[25]	389
[2797]	390	! Derive mask on U and V grid from mask on T grid
	391	obcumask(:,:) = 0.e0
	392	obcvmask(:,:) = 0.e0
	393	DO ij=1, jpjm1
	394	DO ii=1, jpim1
	395	obcumask(ii,ij)=obctmask(ii,ij)*obctmask(ii+1, ij )
	396	obcvmask(ii,ij)=obctmask(ii,ij)*obctmask(ii ,ij+1)
	397	END DO
	398	END DO
	399	CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
[3]	400
[474]	401
[2797]	402	! Mask corrections
	403	! ----------------
	404	DO ik = 1, jpkm1
	405	DO ij = 1, jpj
	406	DO ii = 1, jpi
	407	tmask(ii,ij,ik) = tmask(ii,ij,ik) * obctmask(ii,ij)
	408	umask(ii,ij,ik) = umask(ii,ij,ik) * obcumask(ii,ij)
	409	vmask(ii,ij,ik) = vmask(ii,ij,ik) * obcvmask(ii,ij)
	410	bmask(ii,ij) = bmask(ii,ij) * obctmask(ii,ij)
	411	END DO
	412	END DO
	413	END DO
[3]	414
[2797]	415	DO ik = 1, jpkm1
	416	DO ij = 2, jpjm1
	417	DO ii = 2, jpim1
	418	fmask(ii,ij,ik) = fmask(ii,ij,ik) * obctmask(ii,ij ) * obctmask(ii+1,ij ) &
	419	& * obctmask(ii,ij+1) * obctmask(ii+1,ij+1)
	420	END DO
	421	END DO
	422	END DO
[3]	423
[2797]	424	tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
	425	obctmask(:,:) = tmask(:,:,1)
[3]	426
[2797]	427	! obc masks and bmask are now set to zero on boundary points:
	428	igrd = 1 ! In the free surface case, bmask is at T-points
	429	DO ib_obc = 1, nb_obc
	430	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	431	bmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
	432	ENDDO
	433	ENDDO
	434	!
	435	igrd = 1
	436	DO ib_obc = 1, nb_obc
	437	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	438	obctmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
	439	ENDDO
	440	ENDDO
	441	!
	442	igrd = 2
	443	DO ib_obc = 1, nb_obc
	444	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	445	obcumask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
	446	ENDDO
	447	ENDDO
	448	!
	449	igrd = 3
	450	DO ib_obc = 1, nb_obc
	451	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	452	obcvmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
	453	ENDDO
	454	ENDDO
[3]	455
[2797]	456	! Lateral boundary conditions
	457	CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( obctmask(:,:), 'T', 1. )
	458	CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. )
[3]	459
[2797]	460	DO ib_obc = 1, nb_obc ! Indices and directions of rim velocity components
[3]	461
[2797]	462	idx_obc(ib_obc)%flagu(:) = 0.e0
	463	idx_obc(ib_obc)%flagv(:) = 0.e0
	464	icount = 0
[3]	465
[2797]	466	!flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
	467	!flagu = 0 : u is tangential
	468	!flagu = 1 : u is normal to the boundary and is direction is inward
	469
	470	igrd = 2 ! u-component
	471	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	472	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
	473	nbj => idx_obc(ib_obc)%nbj(ib,igrd)
	474	zefl = obctmask(nbi ,nbj)
	475	zwfl = obctmask(nbi+1,nbj)
	476	IF( zefl + zwfl == 2 ) THEN
	477	icount = icount + 1
	478	ELSE
	479	idx_obc(ib_obc)%flagu(ib)=-zefl+zwfl
[25]	480	ENDIF
[2797]	481	END DO
[3]	482
[2797]	483	!flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
	484	!flagv = 0 : u is tangential
	485	!flagv = 1 : u is normal to the boundary and is direction is inward
[25]	486
[2797]	487	igrd = 3 ! v-component
	488	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	489	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
	490	nbj => idx_obc(ib_obc)%nbj(ib,igrd)
	491	znfl = obctmask(nbi,nbj )
	492	zsfl = obctmask(nbi,nbj+1)
	493	IF( znfl + zsfl == 2 ) THEN
	494	icount = icount + 1
	495	ELSE
	496	idx_obc(ib_obc)%flagv(ib) = -znfl + zsfl
	497	END IF
	498	END DO
	499
	500	IF( icount /= 0 ) THEN
	501	IF(lwp) WRITE(numout,*)
	502	IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', &
	503	' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_obc
	504	IF(lwp) WRITE(numout,*) ' ========== '
	505	IF(lwp) WRITE(numout,*)
	506	nstop = nstop + 1
	507	ENDIF
	508
	509	ENDDO
[3]	510
[2797]	511	! Compute total lateral surface for volume correction:
	512	! ----------------------------------------------------
	513	obcsurftot = 0.e0
	514	IF( ln_vol ) THEN
	515	igrd = 2 ! Lateral surface at U-points
	516	DO ib_obc = 1, nb_obc
	517	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	518	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
	519	nbj => idx_obc(ib_obc)%nbi(ib,igrd)
	520	flagu => idx_obc(ib_obc)%flagu(ib)
	521	obcsurftot = obcsurftot + hu (nbi , nbj) &
	522	& * e2u (nbi , nbj) * ABS( flagu ) &
	523	& * tmask_i(nbi , nbj) &
	524	& * tmask_i(nbi+1, nbj)
	525	ENDDO
	526	ENDDO
[3]	527
[2797]	528	igrd=3 ! Add lateral surface at V-points
	529	DO ib_obc = 1, nb_obc
	530	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
	531	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
	532	nbj => idx_obc(ib_obc)%nbi(ib,igrd)
	533	flagv => idx_obc(ib_obc)%flagv(ib)
	534	obcsurftot = obcsurftot + hv (nbi, nbj ) &
	535	& * e1v (nbi, nbj ) * ABS( flagv ) &
	536	& * tmask_i(nbi, nbj ) &
	537	& * tmask_i(nbi, nbj+1)
	538	ENDDO
	539	ENDDO
	540	!
	541	IF( lk_mpp ) CALL mpp_sum( obcsurftot ) ! sum over the global domain
	542	END IF
	543	!
	544	! Tidy up
	545	!--------
	546	DEALLOCATE(nbidta, nbjdta, nbrdta)
[3]	547
	548	END SUBROUTINE obc_init
[25]	549
[3]	550	#else
[25]	551	!!---------------------------------------------------------------------------------
[2814]	552	!! Dummy module NO open boundaries
[25]	553	!!---------------------------------------------------------------------------------
[3]	554	CONTAINS
[25]	555	SUBROUTINE obc_init ! Dummy routine
[3]	556	END SUBROUTINE obc_init
	557	#endif
	558
[25]	559	!!=================================================================================
[3]	560	END MODULE obcini

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

Download in other formats: