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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

bdyini.F90 in branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

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

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

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

Download in other formats: