Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

bdyini.F90

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

Download in other formats: