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

Context Navigation

source: branches/2011/dev_UKM0_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3062

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