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.

bdyini.F90 in branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3908

Last change on this file since 3908 was 3908, checked in by davestorkey, 11 years ago
Bug fixes for BDY - see ticket #1105.
Property svn:keywords set to `Id`
File size: 34.9 KB

Line
1	MODULE bdyini
2	!!======================================================================
3	!! * MODULE bdyini *
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) rewrite in preparation for OBC-BDY merge
13	!!----------------------------------------------------------------------
14	#if defined key_bdy
15	!!----------------------------------------------------------------------
16	!! 'key_bdy' Unstructured Open Boundary Conditions
17	!!----------------------------------------------------------------------
18	!! bdy_init : Initialization of unstructured open boundaries
19	!!----------------------------------------------------------------------
20	USE timing ! Timing
21	USE oce ! ocean dynamics and tracers variables
22	USE dom_oce ! ocean space and time domain
23	USE bdy_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 bdy_init ! routine called in nemo_init
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 bdy_init
42	!!----------------------------------------------------------------------
43	!! * ROUTINE bdy_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 : bdy_init.nc, input file for unstructured open boundaries
52	!!----------------------------------------------------------------------
53	! namelist variables
54	!-------------------
55	INTEGER, PARAMETER :: jp_nseg = 100
56	INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs
57	INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft
58	INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft
59	INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft
60	INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft
61
62	! local variables
63	!-------------------
64	INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
65	INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers
66	INTEGER :: iw, ie, is, in, inum, id_dummy ! - -
67	INTEGER :: igrd_start, igrd_end, jpbdta ! - -
68	INTEGER, POINTER :: nbi, nbj, nbr ! short cuts
69	REAL , POINTER :: flagu, flagv ! - -
70	REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
71	INTEGER, DIMENSION (2) :: kdimsz
72	INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays
73	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta
74	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
75	REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask
76	CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile
77	CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid
78	!!
79	NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, &
80	& ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
81	& nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta, &
82	#if defined key_lim2
83	& nn_ice_lim2, nn_ice_lim2_dta, &
84	#endif
85	& ln_vol, nn_volctl, nn_rimwidth
86	!!
87	NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft, &
88	nbdysegw, jpiwob, jpjwdt, jpjwft, &
89	nbdysegn, jpjnob, jpindt, jpinft, &
90	nbdysegs, jpjsob, jpisdt, jpisft
91
92	!!----------------------------------------------------------------------
93
94	IF( nn_timing == 1 ) CALL timing_start('bdy_init')
95
96	IF( bdy_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
97
98	IF(lwp) WRITE(numout,*)
99	IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
100	IF(lwp) WRITE(numout,*) '~~~~~~~~'
101	!
102
103	IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', &
104	& ' and general open boundary condition are not compatible' )
105
106	cgrid= (/'t','u','v'/)
107
108	! -----------------------------------------
109	! Initialise and read namelist parameters
110	! -----------------------------------------
111
112	nb_bdy = 0
113	ln_coords_file(:) = .false.
114	cn_coords_file(:) = ''
115	ln_mask_file = .false.
116	cn_mask_file(:) = ''
117	nn_dyn2d(:) = 0
118	nn_dyn2d_dta(:) = -1 ! uninitialised flag
119	nn_dyn3d(:) = 0
120	nn_dyn3d_dta(:) = -1 ! uninitialised flag
121	nn_tra(:) = 0
122	nn_tra_dta(:) = -1 ! uninitialised flag
123	#if defined key_lim2
124	nn_ice_lim2(:) = 0
125	nn_ice_lim2_dta(:)= -1 ! uninitialised flag
126	#endif
127	ln_vol = .false.
128	nn_volctl = -1 ! uninitialised flag
129	nn_rimwidth(:) = -1 ! uninitialised flag
130
131	REWIND( numnam )
132	READ ( numnam, nambdy )
133
134	! -----------------------------------------
135	! Check and write out namelist parameters
136	! -----------------------------------------
137
138	! ! control prints
139	IF(lwp) WRITE(numout,*) ' nambdy'
140
141	IF( nb_bdy .eq. 0 ) THEN
142	IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
143	ELSE
144	IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
145	ENDIF
146
147	DO ib_bdy = 1,nb_bdy
148	IF(lwp) WRITE(numout,*) ' '
149	IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'
150
151	IF( ln_coords_file(ib_bdy) ) THEN
152	IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
153	ELSE
154	IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
155	ENDIF
156	IF(lwp) WRITE(numout,*)
157
158	IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '
159	SELECT CASE( nn_dyn2d(ib_bdy) )
160	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
161	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
162	CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition'
163	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
164	END SELECT
165	IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
166	SELECT CASE( nn_dyn2d_dta(ib_bdy) ) !
167	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
168	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
169	CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' tidal harmonic forcing taken from file'
170	CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files'
171	CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
172	END SELECT
173	ENDIF
174	IF(lwp) WRITE(numout,*)
175
176	IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '
177	SELECT CASE( nn_dyn3d(ib_bdy) )
178	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
179	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
180	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
181	END SELECT
182	IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
183	SELECT CASE( nn_dyn3d_dta(ib_bdy) ) !
184	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
185	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
186	CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
187	END SELECT
188	ENDIF
189	IF(lwp) WRITE(numout,*)
190
191	IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '
192	SELECT CASE( nn_tra(ib_bdy) )
193	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
194	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
195	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' )
196	END SELECT
197	IF( nn_tra(ib_bdy) .gt. 0 ) THEN
198	SELECT CASE( nn_tra_dta(ib_bdy) ) !
199	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
200	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
201	CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
202	END SELECT
203	ENDIF
204	IF(lwp) WRITE(numout,*)
205
206	#if defined key_lim2
207	IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '
208	SELECT CASE( nn_ice_lim2(ib_bdy) )
209	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition'
210	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
211	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' )
212	END SELECT
213	IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN
214	SELECT CASE( nn_ice_lim2_dta(ib_bdy) ) !
215	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
216	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
217	CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
218	END SELECT
219	ENDIF
220	IF(lwp) WRITE(numout,*)
221	#endif
222
223	IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)
224	IF(lwp) WRITE(numout,*)
225
226	ENDDO
227
228	IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)
229	IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
230	IF(lwp) WRITE(numout,*)
231	SELECT CASE ( nn_volctl )
232	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant'
233	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'
234	CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )
235	END SELECT
236	IF(lwp) WRITE(numout,*)
237	ELSE
238	IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
239	IF(lwp) WRITE(numout,*)
240	ENDIF
241
242	! -------------------------------------------------
243	! Initialise indices arrays for open boundaries
244	! -------------------------------------------------
245
246	! Work out global dimensions of boundary data
247	! ---------------------------------------------
248	REWIND( numnam )
249	jpbdta = 1
250	DO ib_bdy = 1, nb_bdy
251
252	IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
253
254	! No REWIND here because may need to read more than one nambdy_index namelist.
255	READ ( numnam, nambdy_index )
256
257	! Automatic boundary definition: if nbdysegX = -1
258	! set boundary to whole side of model domain.
259	IF( nbdysege == -1 ) THEN
260	nbdysege = 1
261	jpieob(1) = jpiglo - 1
262	jpjedt(1) = 2
263	jpjeft(1) = jpjglo - 1
264	ENDIF
265	IF( nbdysegw == -1 ) THEN
266	nbdysegw = 1
267	jpiwob(1) = 2
268	jpjwdt(1) = 2
269	jpjwft(1) = jpjglo - 1
270	ENDIF
271	IF( nbdysegn == -1 ) THEN
272	nbdysegn = 1
273	jpjnob(1) = jpjglo - 1
274	jpindt(1) = 2
275	jpinft(1) = jpiglo - 1
276	ENDIF
277	IF( nbdysegs == -1 ) THEN
278	nbdysegs = 1
279	jpjsob(1) = 2
280	jpisdt(1) = 2
281	jpisft(1) = jpiglo - 1
282	ENDIF
283
284	nblendta(:,ib_bdy) = 0
285	DO iseg = 1, nbdysege
286	igrd = 1
287	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1
288	igrd = 2
289	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1
290	igrd = 3
291	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)
292	ENDDO
293	DO iseg = 1, nbdysegw
294	igrd = 1
295	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1
296	igrd = 2
297	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1
298	igrd = 3
299	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)
300	ENDDO
301	DO iseg = 1, nbdysegn
302	igrd = 1
303	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
304	igrd = 2
305	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)
306	igrd = 3
307	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
308	ENDDO
309	DO iseg = 1, nbdysegs
310	igrd = 1
311	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1
312	igrd = 2
313	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)
314	igrd = 3
315	nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1
316	ENDDO
317
318	nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
319	jpbdta = MAX( jpbdta, MAXVAL(nblendta(:,ib_bdy)) )
320
321
322	ELSE ! Read size of arrays in boundary coordinates file.
323
324
325	CALL iom_open( cn_coords_file(ib_bdy), inum )
326
327	DO igrd = 1, jpbgrd
328	id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )
329	nblendta(igrd,ib_bdy) = kdimsz(1)
330	jpbdta = MAX(jpbdta, kdimsz(1))
331	ENDDO
332	CALL iom_close( inum )
333
334	ENDIF
335
336	ENDDO ! ib_bdy
337
338	! Allocate arrays
339	!---------------
340	ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &
341	& nbrdta(jpbdta, jpbgrd, nb_bdy) )
342
343	ALLOCATE( dta_global(jpbdta, 1, jpk) )
344
345	! Calculate global boundary index arrays or read in from file
346	!------------------------------------------------------------
347	REWIND( numnam )
348	DO ib_bdy = 1, nb_bdy
349
350	IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters
351
352	! No REWIND here because may need to read more than one nambdy_index namelist.
353	READ ( numnam, nambdy_index )
354
355	! Automatic boundary definition: if nbdysegX = -1
356	! set boundary to whole side of model domain.
357	IF( nbdysege == -1 ) THEN
358	nbdysege = 1
359	jpieob(1) = jpiglo - 1
360	jpjedt(1) = 2
361	jpjeft(1) = jpjglo - 1
362	ENDIF
363	IF( nbdysegw == -1 ) THEN
364	nbdysegw = 1
365	jpiwob(1) = 2
366	jpjwdt(1) = 2
367	jpjwft(1) = jpjglo - 1
368	ENDIF
369	IF( nbdysegn == -1 ) THEN
370	nbdysegn = 1
371	jpjnob(1) = jpjglo - 1
372	jpindt(1) = 2
373	jpinft(1) = jpiglo - 1
374	ENDIF
375	IF( nbdysegs == -1 ) THEN
376	nbdysegs = 1
377	jpjsob(1) = 2
378	jpisdt(1) = 2
379	jpisft(1) = jpiglo - 1
380	ENDIF
381
382	! ------------ T points -------------
383	igrd = 1
384	icount = 0
385	DO ir = 1, nn_rimwidth(ib_bdy)
386	! east
387	DO iseg = 1, nbdysege
388	DO ij = jpjedt(iseg), jpjeft(iseg)
389	icount = icount + 1
390	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
391	nbjdta(icount, igrd, ib_bdy) = ij
392	nbrdta(icount, igrd, ib_bdy) = ir
393	ENDDO
394	ENDDO
395	! west
396	DO iseg = 1, nbdysegw
397	DO ij = jpjwdt(iseg), jpjwft(iseg)
398	icount = icount + 1
399	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
400	nbjdta(icount, igrd, ib_bdy) = ij
401	nbrdta(icount, igrd, ib_bdy) = ir
402	ENDDO
403	ENDDO
404	! north
405	DO iseg = 1, nbdysegn
406	DO ii = jpindt(iseg), jpinft(iseg)
407	icount = icount + 1
408	nbidta(icount, igrd, ib_bdy) = ii
409	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
410	nbrdta(icount, igrd, ib_bdy) = ir
411	ENDDO
412	ENDDO
413	! south
414	DO iseg = 1, nbdysegs
415	DO ii = jpisdt(iseg), jpisft(iseg)
416	icount = icount + 1
417	nbidta(icount, igrd, ib_bdy) = ii
418	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
419	nbrdta(icount, igrd, ib_bdy) = ir
420	ENDDO
421	ENDDO
422	ENDDO
423
424	! ------------ U points -------------
425	igrd = 2
426	icount = 0
427	DO ir = 1, nn_rimwidth(ib_bdy)
428	! east
429	DO iseg = 1, nbdysege
430	DO ij = jpjedt(iseg), jpjeft(iseg)
431	icount = icount + 1
432	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
433	nbjdta(icount, igrd, ib_bdy) = ij
434	nbrdta(icount, igrd, ib_bdy) = ir
435	ENDDO
436	ENDDO
437	! west
438	DO iseg = 1, nbdysegw
439	DO ij = jpjwdt(iseg), jpjwft(iseg)
440	icount = icount + 1
441	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
442	nbjdta(icount, igrd, ib_bdy) = ij
443	nbrdta(icount, igrd, ib_bdy) = ir
444	ENDDO
445	ENDDO
446	! north
447	DO iseg = 1, nbdysegn
448	DO ii = jpindt(iseg), jpinft(iseg) - 1
449	icount = icount + 1
450	nbidta(icount, igrd, ib_bdy) = ii
451	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
452	nbrdta(icount, igrd, ib_bdy) = ir
453	ENDDO
454	ENDDO
455	! south
456	DO iseg = 1, nbdysegs
457	DO ii = jpisdt(iseg), jpisft(iseg) - 1
458	icount = icount + 1
459	nbidta(icount, igrd, ib_bdy) = ii
460	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
461	nbrdta(icount, igrd, ib_bdy) = ir
462	ENDDO
463	ENDDO
464	ENDDO
465
466	! ------------ V points -------------
467	igrd = 3
468	icount = 0
469	DO ir = 1, nn_rimwidth(ib_bdy)
470	! east
471	DO iseg = 1, nbdysege
472	DO ij = jpjedt(iseg), jpjeft(iseg) - 1
473	icount = icount + 1
474	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
475	nbjdta(icount, igrd, ib_bdy) = ij
476	nbrdta(icount, igrd, ib_bdy) = ir
477	ENDDO
478	ENDDO
479	! west
480	DO iseg = 1, nbdysegw
481	DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
482	icount = icount + 1
483	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
484	nbjdta(icount, igrd, ib_bdy) = ij
485	nbrdta(icount, igrd, ib_bdy) = ir
486	ENDDO
487	ENDDO
488	! north
489	DO iseg = 1, nbdysegn
490	DO ii = jpindt(iseg), jpinft(iseg)
491	icount = icount + 1
492	nbidta(icount, igrd, ib_bdy) = ii
493	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
494	nbrdta(icount, igrd, ib_bdy) = ir
495	ENDDO
496	ENDDO
497	! south
498	DO iseg = 1, nbdysegs
499	DO ii = jpisdt(iseg), jpisft(iseg)
500	icount = icount + 1
501	nbidta(icount, igrd, ib_bdy) = ii
502	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
503	nbrdta(icount, igrd, ib_bdy) = ir
504	ENDDO
505	ENDDO
506	ENDDO
507
508	ELSE ! Read global index arrays from boundary coordinates file.
509
510	CALL iom_open( cn_coords_file(ib_bdy), inum )
511	DO igrd = 1, jpbgrd
512	CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
513	DO ii = 1,nblendta(igrd,ib_bdy)
514	nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
515	END DO
516	CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
517	DO ii = 1,nblendta(igrd,ib_bdy)
518	nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
519	END DO
520	CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
521	DO ii = 1,nblendta(igrd,ib_bdy)
522	nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
523	END DO
524
525	ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
526	IF(lwp) WRITE(numout,*)
527	IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
528	IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
529	IF (ibr_max < nn_rimwidth(ib_bdy)) &
530	CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
531
532	END DO
533	CALL iom_close( inum )
534
535	ENDIF
536
537	ENDDO
538
539	! Work out dimensions of boundary data on each processor
540	! ------------------------------------------------------
541
542	iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
543	ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
544	is = mjg(1) + 1 ! if monotasking and no zoom, is=2
545	in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
546
547	DO ib_bdy = 1, nb_bdy
548	DO igrd = 1, jpbgrd
549	icount = 0
550	icountr = 0
551	idx_bdy(ib_bdy)%nblen(igrd) = 0
552	idx_bdy(ib_bdy)%nblenrim(igrd) = 0
553	DO ib = 1, nblendta(igrd,ib_bdy)
554	! check that data is in correct order in file
555	ibm1 = MAX(1,ib-1)
556	IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc...
557	IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
558	CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
559	'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
560	ENDIF
561	ENDIF
562	! check if point is in local domain
563	IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. &
564	& nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in ) THEN
565	!
566	icount = icount + 1
567	!
568	IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1
569	ENDIF
570	ENDDO
571	idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
572	idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
573	ENDDO ! igrd
574
575	! Allocate index arrays for this boundary set
576	!--------------------------------------------
577	ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
578	ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
579	ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
580	ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
581	ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
582	ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
583	ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
584	ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
585
586	! Dispatch mapping indices and discrete distances on each processor
587	! -----------------------------------------------------------------
588
589	DO igrd = 1, jpbgrd
590	icount = 0
591	! Loop on rimwidth to ensure outermost points come first in the local arrays.
592	DO ir=1, nn_rimwidth(ib_bdy)
593	DO ib = 1, nblendta(igrd,ib_bdy)
594	! check if point is in local domain and equals ir
595	IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. &
596	& nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND. &
597	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
598	!
599	icount = icount + 1
600	idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
601	idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
602	idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
603	idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
604	ENDIF
605	ENDDO
606	ENDDO
607	ENDDO
608
609	! Compute rim weights for FRS scheme
610	! ----------------------------------
611	DO igrd = 1, jpbgrd
612	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
613	nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
614	idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation
615	! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic
616	! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear
617	END DO
618	END DO
619
620	ENDDO
621
622	! ------------------------------------------------------
623	! Initialise masks and find normal/tangential directions
624	! ------------------------------------------------------
625
626	! Read global 2D mask at T-points: bdytmask
627	! -----------------------------------------
628	! bdytmask = 1 on the computational domain AND on open boundaries
629	! = 0 elsewhere
630
631	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution
632	zmask( : ,:) = 0.e0
633	zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0
634	ELSE IF( ln_mask_file ) THEN
635	CALL iom_open( cn_mask_file, inum )
636	CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
637	CALL iom_close( inum )
638	ELSE
639	zmask(:,:) = 1.e0
640	ENDIF
641
642	DO ij = 1, nlcj ! Save mask over local domain
643	DO ii = 1, nlci
644	bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
645	END DO
646	END DO
647
648	! Derive mask on U and V grid from mask on T grid
649	bdyumask(:,:) = 0.e0
650	bdyvmask(:,:) = 0.e0
651	DO ij=1, jpjm1
652	DO ii=1, jpim1
653	bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
654	bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1)
655	END DO
656	END DO
657	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
658
659
660	! Mask corrections
661	! ----------------
662	DO ik = 1, jpkm1
663	DO ij = 1, jpj
664	DO ii = 1, jpi
665	tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
666	umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
667	vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
668	bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij)
669	END DO
670	END DO
671	END DO
672
673	DO ik = 1, jpkm1
674	DO ij = 2, jpjm1
675	DO ii = 2, jpim1
676	fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) &
677	& * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
678	END DO
679	END DO
680	END DO
681
682	tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
683	bdytmask(:,:) = tmask(:,:,1)
684
685	! bdy masks and bmask are now set to zero on boundary points:
686	igrd = 1 ! In the free surface case, bmask is at T-points
687	DO ib_bdy = 1, nb_bdy
688	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
689	bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
690	ENDDO
691	ENDDO
692	!
693	igrd = 1
694	DO ib_bdy = 1, nb_bdy
695	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
696	bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
697	ENDDO
698	ENDDO
699	!
700	igrd = 2
701	DO ib_bdy = 1, nb_bdy
702	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
703	bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
704	ENDDO
705	ENDDO
706	!
707	igrd = 3
708	DO ib_bdy = 1, nb_bdy
709	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
710	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
711	ENDDO
712	ENDDO
713
714	! Lateral boundary conditions
715	CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
716	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
717
718	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
719
720	idx_bdy(ib_bdy)%flagu(:) = 0.e0
721	idx_bdy(ib_bdy)%flagv(:) = 0.e0
722	icount = 0
723
724	!flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
725	!flagu = 0 : u is tangential
726	!flagu = 1 : u is normal to the boundary and is direction is inward
727
728	igrd = 2 ! u-component
729	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
730	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
731	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
732	zefl = bdytmask(nbi ,nbj)
733	zwfl = bdytmask(nbi+1,nbj)
734	IF( zefl + zwfl == 2 ) THEN
735	icount = icount + 1
736	ELSE
737	idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
738	ENDIF
739	END DO
740
741	!flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
742	!flagv = 0 : u is tangential
743	!flagv = 1 : u is normal to the boundary and is direction is inward
744
745	igrd = 3 ! v-component
746	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
747	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
748	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
749	znfl = bdytmask(nbi,nbj )
750	zsfl = bdytmask(nbi,nbj+1)
751	IF( znfl + zsfl == 2 ) THEN
752	icount = icount + 1
753	ELSE
754	idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
755	END IF
756	END DO
757
758	IF( icount /= 0 ) THEN
759	IF(lwp) WRITE(numout,*)
760	IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', &
761	' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
762	IF(lwp) WRITE(numout,*) ' ========== '
763	IF(lwp) WRITE(numout,*)
764	nstop = nstop + 1
765	ENDIF
766
767	ENDDO
768
769	! Compute total lateral surface for volume correction:
770	! ----------------------------------------------------
771	bdysurftot = 0.e0
772	IF( ln_vol ) THEN
773	igrd = 2 ! Lateral surface at U-points
774	DO ib_bdy = 1, nb_bdy
775	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
776	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
777	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
778	flagu => idx_bdy(ib_bdy)%flagu(ib)
779	bdysurftot = bdysurftot + hu (nbi , nbj) &
780	& * e2u (nbi , nbj) * ABS( flagu ) &
781	& * tmask_i(nbi , nbj) &
782	& * tmask_i(nbi+1, nbj)
783	ENDDO
784	ENDDO
785
786	igrd=3 ! Add lateral surface at V-points
787	DO ib_bdy = 1, nb_bdy
788	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
789	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
790	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
791	flagv => idx_bdy(ib_bdy)%flagv(ib)
792	bdysurftot = bdysurftot + hv (nbi, nbj ) &
793	& * e1v (nbi, nbj ) * ABS( flagv ) &
794	& * tmask_i(nbi, nbj ) &
795	& * tmask_i(nbi, nbj+1)
796	ENDDO
797	ENDDO
798	!
799	IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
800	END IF
801	!
802	! Tidy up
803	!--------
804	DEALLOCATE(nbidta, nbjdta, nbrdta)
805
806	IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
807
808	END SUBROUTINE bdy_init
809
810	#else
811	!!---------------------------------------------------------------------------------
812	!! Dummy module NO open boundaries
813	!!---------------------------------------------------------------------------------
814	CONTAINS
815	SUBROUTINE bdy_init ! Dummy routine
816	END SUBROUTINE bdy_init
817	#endif
818
819	!!=================================================================================
820	END MODULE bdyini

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

Download in other formats: