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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3191

Last change on this file since 3191 was 3191, checked in by davestorkey, 12 years ago
Bug fix for BDY and fldread.F90. Update history comments for BDY. Remove redundant namelist variables in BDY.
Property svn:keywords set to `Id`
File size: 34.8 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	DO ib_bdy = 1, nb_bdy
250
251	jpbdta = 1
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 = 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	jpbdta = 1
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
333	ENDIF
334
335	ENDDO ! ib_bdy
336
337	! Allocate arrays
338	!---------------
339	ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &
340	& nbrdta(jpbdta, jpbgrd, nb_bdy) )
341
342	ALLOCATE( dta_global(jpbdta, 1, jpk) )
343
344	! Calculate global boundary index arrays or read in from file
345	!------------------------------------------------------------
346	REWIND( numnam )
347	DO ib_bdy = 1, nb_bdy
348
349	IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters
350
351	! No REWIND here because may need to read more than one nambdy_index namelist.
352	READ ( numnam, nambdy_index )
353
354	! Automatic boundary definition: if nbdysegX = -1
355	! set boundary to whole side of model domain.
356	IF( nbdysege == -1 ) THEN
357	nbdysege = 1
358	jpieob(1) = jpiglo - 1
359	jpjedt(1) = 2
360	jpjeft(1) = jpjglo - 1
361	ENDIF
362	IF( nbdysegw == -1 ) THEN
363	nbdysegw = 1
364	jpiwob(1) = 2
365	jpjwdt(1) = 2
366	jpjwft(1) = jpjglo - 1
367	ENDIF
368	IF( nbdysegn == -1 ) THEN
369	nbdysegn = 1
370	jpjnob(1) = jpjglo - 1
371	jpindt(1) = 2
372	jpinft(1) = jpiglo - 1
373	ENDIF
374	IF( nbdysegs == -1 ) THEN
375	nbdysegs = 1
376	jpjsob(1) = 2
377	jpisdt(1) = 2
378	jpisft(1) = jpiglo - 1
379	ENDIF
380
381	! ------------ T points -------------
382	igrd = 1
383	icount = 0
384	DO ir = 1, nn_rimwidth(ib_bdy)
385	! east
386	DO iseg = 1, nbdysege
387	DO ij = jpjedt(iseg), jpjeft(iseg)
388	icount = icount + 1
389	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
390	nbjdta(icount, igrd, ib_bdy) = ij
391	nbrdta(icount, igrd, ib_bdy) = ir
392	ENDDO
393	ENDDO
394	! west
395	DO iseg = 1, nbdysegw
396	DO ij = jpjwdt(iseg), jpjwft(iseg)
397	icount = icount + 1
398	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
399	nbjdta(icount, igrd, ib_bdy) = ij
400	nbrdta(icount, igrd, ib_bdy) = ir
401	ENDDO
402	ENDDO
403	! north
404	DO iseg = 1, nbdysegn
405	DO ii = jpindt(iseg), jpinft(iseg)
406	icount = icount + 1
407	nbidta(icount, igrd, ib_bdy) = ii
408	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
409	nbrdta(icount, igrd, ib_bdy) = ir
410	ENDDO
411	ENDDO
412	! south
413	DO iseg = 1, nbdysegs
414	DO ii = jpisdt(iseg), jpisft(iseg)
415	icount = icount + 1
416	nbidta(icount, igrd, ib_bdy) = ii
417	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1
418	nbrdta(icount, igrd, ib_bdy) = ir
419	ENDDO
420	ENDDO
421	ENDDO
422
423	! ------------ U points -------------
424	igrd = 2
425	icount = 0
426	DO ir = 1, nn_rimwidth(ib_bdy)
427	! east
428	DO iseg = 1, nbdysege
429	DO ij = jpjedt(iseg), jpjeft(iseg)
430	icount = icount + 1
431	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
432	nbjdta(icount, igrd, ib_bdy) = ij
433	nbrdta(icount, igrd, ib_bdy) = ir
434	ENDDO
435	ENDDO
436	! west
437	DO iseg = 1, nbdysegw
438	DO ij = jpjwdt(iseg), jpjwft(iseg)
439	icount = icount + 1
440	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
441	nbjdta(icount, igrd, ib_bdy) = ij
442	nbrdta(icount, igrd, ib_bdy) = ir
443	ENDDO
444	ENDDO
445	! north
446	DO iseg = 1, nbdysegn
447	DO ii = jpindt(iseg), jpinft(iseg) - 1
448	icount = icount + 1
449	nbidta(icount, igrd, ib_bdy) = ii
450	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
451	nbrdta(icount, igrd, ib_bdy) = ir
452	ENDDO
453	ENDDO
454	! south
455	DO iseg = 1, nbdysegs
456	DO ii = jpisdt(iseg), jpisft(iseg) - 1
457	icount = icount + 1
458	nbidta(icount, igrd, ib_bdy) = ii
459	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1
460	nbrdta(icount, igrd, ib_bdy) = ir
461	ENDDO
462	ENDDO
463	ENDDO
464
465	! ------------ V points -------------
466	igrd = 3
467	icount = 0
468	DO ir = 1, nn_rimwidth(ib_bdy)
469	! east
470	DO iseg = 1, nbdysege
471	DO ij = jpjedt(iseg), jpjeft(iseg) - 1
472	icount = icount + 1
473	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
474	nbjdta(icount, igrd, ib_bdy) = ij
475	nbrdta(icount, igrd, ib_bdy) = ir
476	ENDDO
477	ENDDO
478	! west
479	DO iseg = 1, nbdysegw
480	DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
481	icount = icount + 1
482	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
483	nbjdta(icount, igrd, ib_bdy) = ij
484	nbrdta(icount, igrd, ib_bdy) = ir
485	ENDDO
486	ENDDO
487	! north
488	DO iseg = 1, nbdysegn
489	DO ii = jpindt(iseg), jpinft(iseg)
490	icount = icount + 1
491	nbidta(icount, igrd, ib_bdy) = ii
492	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
493	nbrdta(icount, igrd, ib_bdy) = ir
494	ENDDO
495	ENDDO
496	! south
497	DO iseg = 1, nbdysegs
498	DO ii = jpisdt(iseg), jpisft(iseg)
499	icount = icount + 1
500	nbidta(icount, igrd, ib_bdy) = ii
501	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1
502	nbrdta(icount, igrd, ib_bdy) = ir
503	ENDDO
504	ENDDO
505	ENDDO
506
507	ELSE ! Read global index arrays from boundary coordinates file.
508
509	DO igrd = 1, jpbgrd
510	CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
511	DO ii = 1,nblendta(igrd,ib_bdy)
512	nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
513	END DO
514	CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
515	DO ii = 1,nblendta(igrd,ib_bdy)
516	nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
517	END DO
518	CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
519	DO ii = 1,nblendta(igrd,ib_bdy)
520	nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
521	END DO
522
523	ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
524	IF(lwp) WRITE(numout,*)
525	IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
526	IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
527	IF (ibr_max < nn_rimwidth(ib_bdy)) &
528	CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
529
530	END DO
531	CALL iom_close( inum )
532
533	ENDIF
534
535	ENDDO
536
537	! Work out dimensions of boundary data on each processor
538	! ------------------------------------------------------
539
540	iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
541	ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
542	is = mjg(1) + 1 ! if monotasking and no zoom, is=2
543	in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
544
545	DO ib_bdy = 1, nb_bdy
546	DO igrd = 1, jpbgrd
547	icount = 0
548	icountr = 0
549	idx_bdy(ib_bdy)%nblen(igrd) = 0
550	idx_bdy(ib_bdy)%nblenrim(igrd) = 0
551	DO ib = 1, nblendta(igrd,ib_bdy)
552	! check that data is in correct order in file
553	ibm1 = MAX(1,ib-1)
554	IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc...
555	IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
556	CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
557	'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
558	ENDIF
559	ENDIF
560	! check if point is in local domain
561	IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. &
562	& nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in ) THEN
563	!
564	icount = icount + 1
565	!
566	IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1
567	ENDIF
568	ENDDO
569	idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
570	idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
571	ENDDO ! igrd
572
573	! Allocate index arrays for this boundary set
574	!--------------------------------------------
575	ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
576	ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
577	ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
578	ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
579	ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
580	ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
581	ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
582	ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
583
584	! Dispatch mapping indices and discrete distances on each processor
585	! -----------------------------------------------------------------
586
587	DO igrd = 1, jpbgrd
588	icount = 0
589	! Loop on rimwidth to ensure outermost points come first in the local arrays.
590	DO ir=1, nn_rimwidth(ib_bdy)
591	DO ib = 1, nblendta(igrd,ib_bdy)
592	! check if point is in local domain and equals ir
593	IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. &
594	& nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND. &
595	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
596	!
597	icount = icount + 1
598	idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
599	idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
600	idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
601	idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
602	ENDIF
603	ENDDO
604	ENDDO
605	ENDDO
606
607	! Compute rim weights for FRS scheme
608	! ----------------------------------
609	DO igrd = 1, jpbgrd
610	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
611	nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
612	idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation
613	! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic
614	! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear
615	END DO
616	END DO
617
618	ENDDO
619
620	! ------------------------------------------------------
621	! Initialise masks and find normal/tangential directions
622	! ------------------------------------------------------
623
624	! Read global 2D mask at T-points: bdytmask
625	! -----------------------------------------
626	! bdytmask = 1 on the computational domain AND on open boundaries
627	! = 0 elsewhere
628
629	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution
630	zmask( : ,:) = 0.e0
631	zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0
632	ELSE IF( ln_mask_file ) THEN
633	CALL iom_open( cn_mask_file, inum )
634	CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
635	CALL iom_close( inum )
636	ELSE
637	zmask(:,:) = 1.e0
638	ENDIF
639
640	DO ij = 1, nlcj ! Save mask over local domain
641	DO ii = 1, nlci
642	bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
643	END DO
644	END DO
645
646	! Derive mask on U and V grid from mask on T grid
647	bdyumask(:,:) = 0.e0
648	bdyvmask(:,:) = 0.e0
649	DO ij=1, jpjm1
650	DO ii=1, jpim1
651	bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
652	bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1)
653	END DO
654	END DO
655	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
656
657
658	! Mask corrections
659	! ----------------
660	DO ik = 1, jpkm1
661	DO ij = 1, jpj
662	DO ii = 1, jpi
663	tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
664	umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
665	vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
666	bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij)
667	END DO
668	END DO
669	END DO
670
671	DO ik = 1, jpkm1
672	DO ij = 2, jpjm1
673	DO ii = 2, jpim1
674	fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) &
675	& * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
676	END DO
677	END DO
678	END DO
679
680	tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
681	bdytmask(:,:) = tmask(:,:,1)
682
683	! bdy masks and bmask are now set to zero on boundary points:
684	igrd = 1 ! In the free surface case, bmask is at T-points
685	DO ib_bdy = 1, nb_bdy
686	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
687	bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
688	ENDDO
689	ENDDO
690	!
691	igrd = 1
692	DO ib_bdy = 1, nb_bdy
693	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
694	bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
695	ENDDO
696	ENDDO
697	!
698	igrd = 2
699	DO ib_bdy = 1, nb_bdy
700	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
701	bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
702	ENDDO
703	ENDDO
704	!
705	igrd = 3
706	DO ib_bdy = 1, nb_bdy
707	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
708	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
709	ENDDO
710	ENDDO
711
712	! Lateral boundary conditions
713	CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
714	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
715
716	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
717
718	idx_bdy(ib_bdy)%flagu(:) = 0.e0
719	idx_bdy(ib_bdy)%flagv(:) = 0.e0
720	icount = 0
721
722	!flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
723	!flagu = 0 : u is tangential
724	!flagu = 1 : u is normal to the boundary and is direction is inward
725
726	igrd = 2 ! u-component
727	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
728	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
729	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
730	zefl = bdytmask(nbi ,nbj)
731	zwfl = bdytmask(nbi+1,nbj)
732	IF( zefl + zwfl == 2 ) THEN
733	icount = icount + 1
734	ELSE
735	idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
736	ENDIF
737	END DO
738
739	!flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
740	!flagv = 0 : u is tangential
741	!flagv = 1 : u is normal to the boundary and is direction is inward
742
743	igrd = 3 ! v-component
744	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
745	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
746	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
747	znfl = bdytmask(nbi,nbj )
748	zsfl = bdytmask(nbi,nbj+1)
749	IF( znfl + zsfl == 2 ) THEN
750	icount = icount + 1
751	ELSE
752	idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
753	END IF
754	END DO
755
756	IF( icount /= 0 ) THEN
757	IF(lwp) WRITE(numout,*)
758	IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', &
759	' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
760	IF(lwp) WRITE(numout,*) ' ========== '
761	IF(lwp) WRITE(numout,*)
762	nstop = nstop + 1
763	ENDIF
764
765	ENDDO
766
767	! Compute total lateral surface for volume correction:
768	! ----------------------------------------------------
769	bdysurftot = 0.e0
770	IF( ln_vol ) THEN
771	igrd = 2 ! Lateral surface at U-points
772	DO ib_bdy = 1, nb_bdy
773	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
774	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
775	nbj => idx_bdy(ib_bdy)%nbi(ib,igrd)
776	flagu => idx_bdy(ib_bdy)%flagu(ib)
777	bdysurftot = bdysurftot + hu (nbi , nbj) &
778	& * e2u (nbi , nbj) * ABS( flagu ) &
779	& * tmask_i(nbi , nbj) &
780	& * tmask_i(nbi+1, nbj)
781	ENDDO
782	ENDDO
783
784	igrd=3 ! Add lateral surface at V-points
785	DO ib_bdy = 1, nb_bdy
786	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
787	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
788	nbj => idx_bdy(ib_bdy)%nbi(ib,igrd)
789	flagv => idx_bdy(ib_bdy)%flagv(ib)
790	bdysurftot = bdysurftot + hv (nbi, nbj ) &
791	& * e1v (nbi, nbj ) * ABS( flagv ) &
792	& * tmask_i(nbi, nbj ) &
793	& * tmask_i(nbi, nbj+1)
794	ENDDO
795	ENDDO
796	!
797	IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
798	END IF
799	!
800	! Tidy up
801	!--------
802	DEALLOCATE(nbidta, nbjdta, nbrdta)
803
804	IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
805
806	END SUBROUTINE bdy_init
807
808	#else
809	!!---------------------------------------------------------------------------------
810	!! Dummy module NO open boundaries
811	!!---------------------------------------------------------------------------------
812	CONTAINS
813	SUBROUTINE bdy_init ! Dummy routine
814	END SUBROUTINE bdy_init
815	#endif
816
817	!!=================================================================================
818	END MODULE bdyini

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

Download in other formats: