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

Context Navigation

source: branches/CMIP5_IPSL/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3307

Last change on this file since 3307 was 1528, checked in by rblod, 15 years ago
Suppress rigid-lid option, see ticket #486
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 18.9 KB

Line
1	MODULE bdyini
2	!!======================================================================
3	!! * MODULE bdyini *
4	!! Unstructured open boundaries : initialisation
5	!!======================================================================
6	!! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code
7	!! - ! 2007-01 (D. Storkey) Update to use IOM module
8	!! - ! 2007-01 (D. Storkey) Tidal forcing
9	!! 3.0 ! 2008-04 (NEMO team) add in the reference version
10	!!----------------------------------------------------------------------
11	#if defined key_bdy
12	!!----------------------------------------------------------------------
13	!! 'key_bdy' Unstructured Open Boundary Conditions
14	!!----------------------------------------------------------------------
15	!! bdy_init : Initialization of unstructured open boundaries
16	!!----------------------------------------------------------------------
17	USE oce ! ocean dynamics and tracers variables
18	USE dom_oce ! ocean space and time domain
19	USE bdy_oce ! unstructured open boundary conditions
20	USE bdytides ! tides at open boundaries initialization (tide_init routine)
21	USE in_out_manager ! I/O units
22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
23	USE lib_mpp ! for mpp_sum
24	USE iom ! I/O
25
26	IMPLICIT NONE
27	PRIVATE
28
29	PUBLIC bdy_init ! routine called by opa.F90
30
31	!!----------------------------------------------------------------------
32	!! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
33	!! $Id$
34	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
35	!!---------------------------------------------------------------------------------
36
37	CONTAINS
38
39	SUBROUTINE bdy_init
40	!!----------------------------------------------------------------------
41	!! * ROUTINE bdy_init *
42	!!
43	!! ** Purpose : Initialization of the dynamics and tracer fields with
44	!! unstructured open boundaries.
45	!!
46	!! ** Method : Read initialization arrays (mask, indices) to identify
47	!! an unstructured open boundary
48	!!
49	!! ** Input : bdy_init.nc, input file for unstructured open boundaries
50	!!
51	!!----------------------------------------------------------------------
52	INTEGER :: ii, ij, ik, igrd, ib, ir ! dummy loop indices
53	INTEGER :: icount, icountr
54	INTEGER :: ib_len, ibr_max
55	INTEGER :: iw, ie, is, in
56	INTEGER :: inum ! temporary logical unit
57	INTEGER :: id_dummy ! temporary integers
58	INTEGER :: igrd_start, igrd_end ! start and end of loops on igrd
59	INTEGER, DIMENSION (2) :: kdimsz
60	INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta
61	INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbrdta ! Discrete distance from rim points
62	REAL(wp) :: zefl, zwfl, znfl, zsfl ! temporary scalars
63	REAL(wp) , DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask
64	REAL(wp) , DIMENSION(jpbdta,1) :: zdta ! temporary array
65	CHARACTER(LEN=80),DIMENSION(3) :: clfile
66	!!
67	NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, &
68	& ln_bdy_tides, ln_bdy_clim, ln_bdy_vol, ln_bdy_mask, &
69	& ln_bdy_dyn_fla, ln_bdy_dyn_frs, ln_bdy_tra_frs, &
70	& nbdy_dta , nb_rimwidth , volbdy
71	!!----------------------------------------------------------------------
72
73	IF(lwp) WRITE(numout,*)
74	IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries'
75	IF(lwp) WRITE(numout,*) '~~~~~~~~'
76	!
77	IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', &
78	' and unstructured open boundary condition are not compatible' )
79
80	#if defined key_obc
81	CALL ctl_stop( 'Straight open boundaries,', &
82	' and unstructured open boundaries are not compatible' )
83	#endif
84
85	! Read namelist parameters
86	! ---------------------------
87	REWIND( numnam )
88	READ ( numnam, nambdy )
89
90	! control prints
91	IF(lwp) WRITE(numout,*) ' nambdy'
92
93	! Check nbdy_dta value
94	IF(lwp) WRITE(numout,*) 'nbdy_dta =', nbdy_dta
95	IF(lwp) WRITE(numout,*) ' '
96	SELECT CASE( nbdy_dta )
97	CASE( 0 )
98	IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
99	CASE( 1 )
100	IF(lwp) WRITE(numout,*) ' boundary data taken from file'
101	CASE DEFAULT
102	CALL ctl_stop( 'nbdy_dta must be 0 or 1' )
103	END SELECT
104
105	IF(lwp) WRITE(numout,*) ' '
106	IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth
107
108	IF(lwp) WRITE(numout,*) ' '
109	IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy
110
111	IF (ln_bdy_vol) THEN
112	SELECT CASE ( volbdy ) ! Check volbdy value
113	CASE( 1 )
114	IF(lwp) WRITE(numout,*) ' The total volume will be constant'
115	CASE( 0 )
116	IF(lwp) WRITE(numout,*) ' The total volume will vary according'
117	IF(lwp) WRITE(numout,*) ' to the surface E-P flux'
118	CASE DEFAULT
119	CALL ctl_stop( 'volbdy must be 0 or 1' )
120	END SELECT
121	ELSE
122	IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries'
123	IF(lwp) WRITE(numout,*) ' '
124	ENDIF
125
126	IF (ln_bdy_tides) THEN
127	IF(lwp) WRITE(numout,*) ' '
128	IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries'
129	IF(lwp) WRITE(numout,*) ' '
130	ENDIF
131
132	IF (ln_bdy_dyn_fla) THEN
133	IF(lwp) WRITE(numout,*) ' '
134	IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries'
135	IF(lwp) WRITE(numout,*) ' '
136	ENDIF
137
138	IF (ln_bdy_dyn_frs) THEN
139	IF(lwp) WRITE(numout,*) ' '
140	IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries'
141	IF(lwp) WRITE(numout,*) ' '
142	ENDIF
143
144	IF (ln_bdy_tra_frs) THEN
145	IF(lwp) WRITE(numout,*) ' '
146	IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries'
147	IF(lwp) WRITE(numout,*) ' '
148	ENDIF
149
150	! Read tides namelist
151	! ------------------------
152	IF( ln_bdy_tides ) CALL tide_init
153
154	! Read arrays defining unstructured open boundaries
155	! -------------------------------------------------
156
157	! Read global 2D mask at T-points: bdytmask
158	! *****************************************
159	! bdytmask = 1 on the computational domain AND on open boundaries
160	! = 0 elsewhere
161
162	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
163	zmask( : ,:) = 0.e0
164	zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0
165	ELSE IF ( ln_bdy_mask ) THEN
166	CALL iom_open( filbdy_mask, inum )
167	CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
168	CALL iom_close( inum )
169	ELSE
170	zmask(:,:) = 1.e0
171	ENDIF
172
173	! Save mask over local domain
174	DO ij = 1, nlcj
175	DO ii = 1, nlci
176	bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
177	END DO
178	END DO
179
180	! Derive mask on U and V grid from mask on T grid
181	bdyumask(:,:) = 0.e0
182	bdyvmask(:,:) = 0.e0
183	DO ij=1, jpjm1
184	DO ii=1, jpim1
185	bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
186	bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1)
187	END DO
188	END DO
189
190	! Lateral boundary conditions
191	CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
192	CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
193
194	! Read discrete distance and mapping indices
195	! ******************************************
196	nbidta(:,:) = 0.e0
197	nbjdta(:,:) = 0.e0
198	nbrdta(:,:) = 0.e0
199
200	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
201	icount = 0
202	! Define west boundary (from ii=2 to ii=1+nb_rimwidth):
203	DO ir = 1, nb_rimwidth
204	DO ij = 3, jpjglo-2
205	icount=icount+1
206	nbidta(icount,:) = ir + 1 + (jpizoom-1)
207	nbjdta(icount,:) = ij + (jpjzoom-1)
208	nbrdta(icount,:) = ir
209	END DO
210	END DO
211
212	! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth):
213	DO ir=1,nb_rimwidth
214	DO ij=3,jpjglo-2
215	icount=icount+1
216	nbidta(icount,:) = jpiglo-ir + (jpizoom-1)
217	nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points
218	nbjdta(icount,:) = ij + (jpjzoom-1)
219	nbrdta(icount,:) = ir
220	END DO
221	END DO
222
223	ELSE ! Read indices and distances in unstructured boundary data files
224
225	IF( ln_bdy_tides ) THEN
226	! Read tides input files for preference in case there are
227	! no bdydata files.
228	clfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc'
229	clfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc'
230	clfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc'
231	ELSE
232	clfile(1) = filbdy_data_T
233	clfile(2) = filbdy_data_U
234	clfile(3) = filbdy_data_V
235	ENDIF
236
237	! how many files are we to read in?
238	igrd_start = 1
239	igrd_end = 3
240	IF(.NOT. ln_bdy_tides ) THEN
241	IF(.NOT. (ln_bdy_dyn_fla) .AND..NOT. (ln_bdy_tra_frs)) THEN
242	! No T-grid file.
243	igrd_start = 2
244	ELSEIF ( .NOT. ln_bdy_dyn_frs .AND..NOT. ln_bdy_dyn_fla ) THEN
245	! No U-grid or V-grid file.
246	igrd_end = 1
247	ENDIF
248	ENDIF
249
250	DO igrd = igrd_start, igrd_end
251	CALL iom_open( clfile(igrd), inum )
252	id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz )
253	WRITE(numout,*) 'kdimsz : ',kdimsz
254	ib_len = kdimsz(1)
255	IF( ib_len > jpbdta) CALL ctl_stop( &
256	'Boundary data array in file too long.', &
257	'File :', TRIM(clfile(igrd)), &
258	'increase parameter jpbdta.' )
259
260	CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) )
261	DO ii = 1,ib_len
262	nbidta(ii,igrd) = INT( zdta(ii,1) )
263	END DO
264	CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) )
265	DO ii = 1,ib_len
266	nbjdta(ii,igrd) = INT( zdta(ii,1) )
267	END DO
268	CALL iom_get ( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) )
269	DO ii = 1,ib_len
270	nbrdta(ii,igrd) = INT( zdta(ii,1) )
271	END DO
272	CALL iom_close( inum )
273
274	! Check that rimwidth in file is big enough:
275	ibr_max = MAXVAL( nbrdta(:,igrd) )
276	IF(lwp) WRITE(numout,*)
277	IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
278	IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth
279	IF (ibr_max < nb_rimwidth) CALL ctl_stop( &
280	'nb_rimwidth is larger than maximum rimwidth in file' )
281	!
282	END DO
283	!
284	ENDIF
285
286	! Dispatch mapping indices and discrete distances on each processor
287	! *****************************************************************
288
289	iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
290	ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
291	is = mjg(1) + 1 ! if monotasking and no zoom, is=2
292	in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
293
294	DO igrd = igrd_start, igrd_end
295	icount = 0
296	icountr = 0
297	nblen(igrd) = 0
298	nblenrim(igrd) = 0
299	nblendta(igrd) = 0
300	DO ir=1, nb_rimwidth
301	DO ib = 1, jpbdta
302	! check if point is in local domain and equals ir
303	IF( nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND. &
304	& nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND. &
305	& nbrdta(ib,igrd) == ir ) THEN
306	!
307	icount = icount + 1
308	!
309	IF( ir == 1 ) icountr = icountr+1
310	IF (icount > jpbdim) THEN
311	IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'
312	nstop = nstop + 1
313	ELSE
314	nbi(icount, igrd) = nbidta(ib,igrd)- mig(1)+1
315	nbj(icount, igrd) = nbjdta(ib,igrd)- mjg(1)+1
316	nbr(icount, igrd) = nbrdta(ib,igrd)
317	nbmap(icount,igrd) = ib
318	ENDIF
319	ENDIF
320	END DO
321	END DO
322	nblenrim(igrd) = icountr !: length of rim boundary data on each proc
323	nblen (igrd) = icount !: length of boundary data on each proc
324	END DO
325
326	! Compute rim weights
327	! -------------------
328	DO igrd = igrd_start, igrd_end
329	DO ib = 1, nblen(igrd)
330	! tanh formulation
331	nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 )
332	! quadratic
333	! nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2
334	! linear
335	! nbw(ib,igrd) = FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth)
336	END DO
337	END DO
338
339	! Mask corrections
340	! ----------------
341	DO ik = 1, jpkm1
342	DO ij = 1, jpj
343	DO ii = 1, jpi
344	tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
345	umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
346	vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
347	bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij)
348	END DO
349	END DO
350	END DO
351
352	DO ik = 1, jpkm1
353	DO ij = 2, jpjm1
354	DO ii = 2, jpim1
355	fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) &
356	& * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
357	END DO
358	END DO
359	END DO
360
361	tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
362	bdytmask(:,:) = tmask(:,:,1)
363
364	! bdy masks and bmask are now set to zero on boundary points:
365	igrd = 1 ! In the free surface case, bmask is at T-points
366	DO ib = 1, nblenrim(igrd)
367	bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
368	END DO
369	!
370	igrd = 1
371	DO ib = 1, nblenrim(igrd)
372	bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
373	END DO
374	!
375	igrd = 2
376	DO ib = 1, nblenrim(igrd)
377	bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
378	END DO
379	!
380	igrd = 3
381	DO ib = 1, nblenrim(igrd)
382	bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
383	END DO
384
385	! Lateral boundary conditions
386	CALL lbc_lnk( fmask , 'F', 1. )
387	CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
388	CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
389	CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
390
391	IF( ln_bdy_vol .OR. ln_bdy_dyn_fla ) THEN ! Indices and directions of rim velocity components
392	!
393	!flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
394	!flagu = 0 : u is tangential
395	!flagu = 1 : u is normal to the boundary and is direction is inward
396	icount = 0
397	flagu(:) = 0.e0
398
399	igrd = 2 ! u-component
400	DO ib = 1, nblenrim(igrd)
401	zefl=bdytmask(nbi(ib,igrd) , nbj(ib,igrd))
402	zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd))
403	IF( zefl + zwfl ==2 ) THEN
404	icount = icount +1
405	ELSE
406	flagu(ib)=-zefl+zwfl
407	ENDIF
408	END DO
409
410	!flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
411	!flagv = 0 : u is tangential
412	!flagv = 1 : u is normal to the boundary and is direction is inward
413	flagv(:) = 0.e0
414
415	igrd = 3 ! v-component
416	DO ib = 1, nblenrim(igrd)
417	znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd))
418	zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1)
419	IF( znfl + zsfl ==2 ) THEN
420	icount = icount + 1
421	ELSE
422	flagv(ib) = -znfl + zsfl
423	END IF
424	END DO
425
426	IF( icount /= 0 ) THEN
427	IF(lwp) WRITE(numout,*)
428	IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', &
429	' are not boundary points. Check nbi, nbj, indices.'
430	IF(lwp) WRITE(numout,*) ' ========== '
431	IF(lwp) WRITE(numout,*)
432	nstop = nstop + 1
433	ENDIF
434
435	ENDIF
436
437	! Compute total lateral surface for volume correction:
438	! ----------------------------------------------------
439
440	bdysurftot = 0.e0
441	IF( ln_bdy_vol ) THEN
442	igrd = 2 ! Lateral surface at U-points
443	DO ib = 1, nblenrim(igrd)
444	bdysurftot = bdysurftot + hu (nbi(ib,igrd) ,nbj(ib,igrd)) &
445	& * e2u (nbi(ib,igrd) ,nbj(ib,igrd)) * ABS( flagu(ib) ) &
446	& * tmask_i(nbi(ib,igrd) ,nbj(ib,igrd)) &
447	& * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd))
448	END DO
449
450	igrd=3 ! Add lateral surface at V-points
451	DO ib = 1, nblenrim(igrd)
452	bdysurftot = bdysurftot + hv (nbi(ib,igrd),nbj(ib,igrd) ) &
453	& * e1v (nbi(ib,igrd),nbj(ib,igrd) ) * ABS( flagv(ib) ) &
454	& * tmask_i(nbi(ib,igrd),nbj(ib,igrd) ) &
455	& * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1)
456	END DO
457
458	IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
459	END IF
460
461	! Initialise bdy data arrays
462	! --------------------------
463	tbdy(:,:) = 0.e0
464	sbdy(:,:) = 0.e0
465	ubdy(:,:) = 0.e0
466	vbdy(:,:) = 0.e0
467	sshbdy(:) = 0.e0
468	ubtbdy(:) = 0.e0
469	vbtbdy(:) = 0.e0
470
471	! Read in tidal constituents and adjust for model start time
472	! ----------------------------------------------------------
473	IF( ln_bdy_tides ) CALL tide_data
474	!
475	END SUBROUTINE bdy_init
476
477	#else
478	!!---------------------------------------------------------------------------------
479	!! Dummy module NO unstructured open boundaries
480	!!---------------------------------------------------------------------------------
481	CONTAINS
482	SUBROUTINE bdy_init ! Dummy routine
483	END SUBROUTINE bdy_init
484	#endif
485
486	!!=================================================================================
487	END MODULE bdyini

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

Download in other formats: