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

Context Navigation

source: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3367

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

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

Download in other formats: