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

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

Initial restructuring.
Use fldread to read open boundary data.

Property svn:keywords set to Id

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

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

Download in other formats: