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

Context Navigation

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 2888

Last change on this file since 2888 was 2888, checked in by davestorkey, 13 years ago
Move changes into updated BDY module and restore old OBC code. (Full merge to take place next year).
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: