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 @ 2818

Last change on this file since 2818 was 2818, checked in by davestorkey, 13 years ago
Bug fixes for the dynspg_ts case.
Property svn:keywords set to `Id`
File size: 24.4 KB

Line
1	MODULE obcini
2	!!======================================================================
3	!! * MODULE obcini *
4	!! Unstructured open boundaries : initialisation
5	!!======================================================================
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 ---
14	!!----------------------------------------------------------------------
15	#if defined key_obc
16	!!----------------------------------------------------------------------
17	!! 'key_obc' Unstructured Open Boundary Conditions
18	!!----------------------------------------------------------------------
19	!! obc_init : Initialization of unstructured open boundaries
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers variables
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
25	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
26	USE lib_mpp ! for mpp_sum
27	USE iom ! I/O
28
29	IMPLICIT NONE
30	PRIVATE
31
32	PUBLIC obc_init ! routine called by opa.F90
33
34	!!----------------------------------------------------------------------
35	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
36	!! $Id$
37	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38	!!----------------------------------------------------------------------
39	CONTAINS
40
41	SUBROUTINE obc_init
42	!!----------------------------------------------------------------------
43	!! * ROUTINE obc_init *
44	!!
45	!! ** Purpose : Initialization of the dynamics and tracer fields with
46	!! unstructured open boundaries.
47	!!
48	!! ** Method : Read initialization arrays (mask, indices) to identify
49	!! an unstructured open boundary
50	!!
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, ibm1 ! 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
67	!!
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
74	& nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl, &
75	& nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out, &
76	& nn_dmp3d_in, nn_dmp3d_out
77	!!----------------------------------------------------------------------
78
79	IF( obc_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'obc_init : unable to allocate oce arrays' )
80
81	IF(lwp) WRITE(numout,*)
82	IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries'
83	IF(lwp) WRITE(numout,*) '~~~~~~~~'
84	!
85
86	IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', &
87	& ' and general open boundary condition are not compatible' )
88
89	cgrid= (/'T','U','V'/)
90
91	! -----------------------------------------
92	! Initialise and read namelist parameters
93	! -----------------------------------------
94
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
109	nn_tides(:) = 0 ! default to no tidal forcing
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
116
117	REWIND( numnam )
118	READ ( numnam, namobc )
119
120	! -----------------------------------------
121	! Check and write out namelist parameters
122	! -----------------------------------------
123
124	! ! control prints
125	IF(lwp) WRITE(numout,*) ' namobc'
126
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
132
133	DO ib_obc = 1,nb_obc
134	IF(lwp) WRITE(numout,*) ' '
135	IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------'
136
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,*)
144
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,*)
153
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,*)
161
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,*)
169
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
179
180	IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth
181	IF(lwp) WRITE(numout,*)
182
183	SELECT CASE( nn_tides(ib_obc) )
184	CASE(0)
185	IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing'
186	IF(lwp) WRITE(numout,*)
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
196
197	ENDDO
198
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
212
213	! -------------------------------------------------
214	! Initialise indices arrays for open boundaries
215	! -------------------------------------------------
216
217	! Work out global dimensions of boundary data
218	! ---------------------------------------------
219	DO ib_obc = 1, nb_obc
220
221	jpbdta = 1
222	IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters
223
224
225	!! 1. Read parameters from namelist
226	!! 2. Work out global size of boundary data arrays nblendta and jpbdta
227
228
229	ELSE ! Read size of arrays in boundary coordinates file.
230
231
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
239
240	ENDIF
241
242	ENDDO
243
244	! Allocate arrays
245	!---------------
246	ALLOCATE( nbidta(jpbdta, jpbgrd, nb_obc), nbjdta(jpbdta, jpbgrd, nb_obc), &
247	& nbrdta(jpbdta, jpbgrd, nb_obc) )
248
249	ALLOCATE( dta_global(jpbdta, 1, jpk) )
250
251	! Calculate global boundary index arrays or read in from file
252	!------------------------------------------------------------
253	DO ib_obc = 1, nb_obc
254
255	IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters
256
257	!! Calculate global index arrays from namelist parameters
258
259	ELSE ! Read global index arrays from boundary coordinates file.
260
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
274
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) )
281
282	END DO
283	CALL iom_close( inum )
284
285	ENDIF
286
287	ENDDO
288
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
296
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 that data is in correct order in file
305	ibm1 = MAX(1,ib-1)
306	IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc...
307	IF( nbrdta(ib,igrd,ib_obc) < nbrdta(ibm1,igrd,ib_obc) ) THEN
308	CALL ctl_stop('obc_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
309	'A utility for re-ordering boundary coordinates and data files exists in CDFTOOLS')
310	ENDIF
311	ENDIF
312	! check if point is in local domain
313	IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. &
314	& nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in ) THEN
315	!
316	icount = icount + 1
317	!
318	IF( nbrdta(ib,igrd,ib_obc) == 1 ) icountr = icountr+1
319	ENDIF
320	ENDDO
321	idx_obc(ib_obc)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
322	idx_obc(ib_obc)%nblen (igrd) = icount !: length of boundary data on each proc
323	ENDDO ! igrd
324
325	! Allocate index arrays for this boundary set
326	!--------------------------------------------
327	ilen1 = MAXVAL(idx_obc(ib_obc)%nblen(:))
328	ALLOCATE( idx_obc(ib_obc)%nbi(ilen1,jpbgrd) )
329	ALLOCATE( idx_obc(ib_obc)%nbj(ilen1,jpbgrd) )
330	ALLOCATE( idx_obc(ib_obc)%nbr(ilen1,jpbgrd) )
331	ALLOCATE( idx_obc(ib_obc)%nbmap(ilen1,jpbgrd) )
332	ALLOCATE( idx_obc(ib_obc)%nbw(ilen1,jpbgrd) )
333	ALLOCATE( idx_obc(ib_obc)%flagu(ilen1) )
334	ALLOCATE( idx_obc(ib_obc)%flagv(ilen1) )
335
336	! Dispatch mapping indices and discrete distances on each processor
337	! -----------------------------------------------------------------
338
339	DO igrd = 1, jpbgrd
340	icount = 0
341	! Loop on rimwidth to ensure outermost points come first in the local arrays.
342	DO ir=1, nn_rimwidth(ib_obc)
343	DO ib = 1, nblendta(igrd,ib_obc)
344	! check if point is in local domain and equals ir
345	IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. &
346	& nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in .AND. &
347	& nbrdta(ib,igrd,ib_obc) == ir ) THEN
348	!
349	icount = icount + 1
350	idx_obc(ib_obc)%nbi(icount,igrd) = nbidta(ib,igrd,ib_obc)- mig(1)+1
351	idx_obc(ib_obc)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_obc)- mjg(1)+1
352	idx_obc(ib_obc)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_obc)
353	idx_obc(ib_obc)%nbmap(icount,igrd) = ib
354	ENDIF
355	ENDDO
356	ENDDO
357	ENDDO
358
359	! Compute rim weights for FRS scheme
360	! ----------------------------------
361	DO igrd = 1, jpbgrd
362	DO ib = 1, idx_obc(ib_obc)%nblen(igrd)
363	nbr => idx_obc(ib_obc)%nbr(ib,igrd)
364	idx_obc(ib_obc)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation
365	! idx_obc(ib_obc)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic
366	! idx_obc(ib_obc)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear
367	END DO
368	END DO
369
370	ENDDO
371
372	! ------------------------------------------------------
373	! Initialise masks and find normal/tangential directions
374	! ------------------------------------------------------
375
376	! Read global 2D mask at T-points: obctmask
377	! -----------------------------------------
378	! obctmask = 1 on the computational domain AND on open boundaries
379	! = 0 elsewhere
380
381	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution
382	zmask( : ,:) = 0.e0
383	zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0
384	ELSE IF( ln_mask_file ) THEN
385	CALL iom_open( cn_mask_file, inum )
386	CALL iom_get ( inum, jpdom_data, 'obc_msk', zmask(:,:) )
387	CALL iom_close( inum )
388	ELSE
389	zmask(:,:) = 1.e0
390	ENDIF
391
392	DO ij = 1, nlcj ! Save mask over local domain
393	DO ii = 1, nlci
394	obctmask(ii,ij) = zmask( mig(ii), mjg(ij) )
395	END DO
396	END DO
397
398	! Derive mask on U and V grid from mask on T grid
399	obcumask(:,:) = 0.e0
400	obcvmask(:,:) = 0.e0
401	DO ij=1, jpjm1
402	DO ii=1, jpim1
403	obcumask(ii,ij)=obctmask(ii,ij)*obctmask(ii+1, ij )
404	obcvmask(ii,ij)=obctmask(ii,ij)*obctmask(ii ,ij+1)
405	END DO
406	END DO
407	CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
408
409
410	! Mask corrections
411	! ----------------
412	DO ik = 1, jpkm1
413	DO ij = 1, jpj
414	DO ii = 1, jpi
415	tmask(ii,ij,ik) = tmask(ii,ij,ik) * obctmask(ii,ij)
416	umask(ii,ij,ik) = umask(ii,ij,ik) * obcumask(ii,ij)
417	vmask(ii,ij,ik) = vmask(ii,ij,ik) * obcvmask(ii,ij)
418	bmask(ii,ij) = bmask(ii,ij) * obctmask(ii,ij)
419	END DO
420	END DO
421	END DO
422
423	DO ik = 1, jpkm1
424	DO ij = 2, jpjm1
425	DO ii = 2, jpim1
426	fmask(ii,ij,ik) = fmask(ii,ij,ik) * obctmask(ii,ij ) * obctmask(ii+1,ij ) &
427	& * obctmask(ii,ij+1) * obctmask(ii+1,ij+1)
428	END DO
429	END DO
430	END DO
431
432	tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
433	obctmask(:,:) = tmask(:,:,1)
434
435	! obc masks and bmask are now set to zero on boundary points:
436	igrd = 1 ! In the free surface case, bmask is at T-points
437	DO ib_obc = 1, nb_obc
438	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
439	bmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
440	ENDDO
441	ENDDO
442	!
443	igrd = 1
444	DO ib_obc = 1, nb_obc
445	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
446	obctmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
447	ENDDO
448	ENDDO
449	!
450	igrd = 2
451	DO ib_obc = 1, nb_obc
452	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
453	obcumask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
454	ENDDO
455	ENDDO
456	!
457	igrd = 3
458	DO ib_obc = 1, nb_obc
459	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
460	obcvmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0
461	ENDDO
462	ENDDO
463
464	! Lateral boundary conditions
465	CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( obctmask(:,:), 'T', 1. )
466	CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. )
467
468	DO ib_obc = 1, nb_obc ! Indices and directions of rim velocity components
469
470	idx_obc(ib_obc)%flagu(:) = 0.e0
471	idx_obc(ib_obc)%flagv(:) = 0.e0
472	icount = 0
473
474	!flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
475	!flagu = 0 : u is tangential
476	!flagu = 1 : u is normal to the boundary and is direction is inward
477
478	igrd = 2 ! u-component
479	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
480	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
481	nbj => idx_obc(ib_obc)%nbj(ib,igrd)
482	zefl = obctmask(nbi ,nbj)
483	zwfl = obctmask(nbi+1,nbj)
484	IF( zefl + zwfl == 2 ) THEN
485	icount = icount + 1
486	ELSE
487	idx_obc(ib_obc)%flagu(ib)=-zefl+zwfl
488	ENDIF
489	END DO
490
491	!flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
492	!flagv = 0 : u is tangential
493	!flagv = 1 : u is normal to the boundary and is direction is inward
494
495	igrd = 3 ! v-component
496	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
497	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
498	nbj => idx_obc(ib_obc)%nbj(ib,igrd)
499	znfl = obctmask(nbi,nbj )
500	zsfl = obctmask(nbi,nbj+1)
501	IF( znfl + zsfl == 2 ) THEN
502	icount = icount + 1
503	ELSE
504	idx_obc(ib_obc)%flagv(ib) = -znfl + zsfl
505	END IF
506	END DO
507
508	IF( icount /= 0 ) THEN
509	IF(lwp) WRITE(numout,*)
510	IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', &
511	' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_obc
512	IF(lwp) WRITE(numout,*) ' ========== '
513	IF(lwp) WRITE(numout,*)
514	nstop = nstop + 1
515	ENDIF
516
517	ENDDO
518
519	! Compute total lateral surface for volume correction:
520	! ----------------------------------------------------
521	obcsurftot = 0.e0
522	IF( ln_vol ) THEN
523	igrd = 2 ! Lateral surface at U-points
524	DO ib_obc = 1, nb_obc
525	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
526	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
527	nbj => idx_obc(ib_obc)%nbi(ib,igrd)
528	flagu => idx_obc(ib_obc)%flagu(ib)
529	obcsurftot = obcsurftot + hu (nbi , nbj) &
530	& * e2u (nbi , nbj) * ABS( flagu ) &
531	& * tmask_i(nbi , nbj) &
532	& * tmask_i(nbi+1, nbj)
533	ENDDO
534	ENDDO
535
536	igrd=3 ! Add lateral surface at V-points
537	DO ib_obc = 1, nb_obc
538	DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd)
539	nbi => idx_obc(ib_obc)%nbi(ib,igrd)
540	nbj => idx_obc(ib_obc)%nbi(ib,igrd)
541	flagv => idx_obc(ib_obc)%flagv(ib)
542	obcsurftot = obcsurftot + hv (nbi, nbj ) &
543	& * e1v (nbi, nbj ) * ABS( flagv ) &
544	& * tmask_i(nbi, nbj ) &
545	& * tmask_i(nbi, nbj+1)
546	ENDDO
547	ENDDO
548	!
549	IF( lk_mpp ) CALL mpp_sum( obcsurftot ) ! sum over the global domain
550	END IF
551	!
552	! Tidy up
553	!--------
554	DEALLOCATE(nbidta, nbjdta, nbrdta)
555
556	END SUBROUTINE obc_init
557
558	#else
559	!!---------------------------------------------------------------------------------
560	!! Dummy module NO open boundaries
561	!!---------------------------------------------------------------------------------
562	CONTAINS
563	SUBROUTINE obc_init ! Dummy routine
564	END SUBROUTINE obc_init
565	#endif
566
567	!!=================================================================================
568	END MODULE obcini

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

Download in other formats: