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.

dommsk.F90 in branches/2011/dev_LOCEAN_CMCC_2011/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2011/dev_LOCEAN_CMCC_2011/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 3084

Last change on this file since 3084 was 3084, checked in by cetlod, 12 years ago
dev_LOCEAN_CMCC_2011: add in changes CMCC devlopments into the new branch
Property svn:keywords set to `Id`
File size: 28.3 KB

Line
1	MODULE dommsk
2	!!======================================================================
3	!! * MODULE dommsk *
4	!! Ocean initialization : domain land/sea mask
5	!!======================================================================
6	!! History : OPA ! 1987-07 (G. Madec) Original code
7	!! 6.0 ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon)
8	!! 7.0 ! 1996-01 (G. Madec) suppression of common work arrays
9	!! - ! 1996-05 (G. Madec) mask computed from tmask and sup-
10	!! ! pression of the double computation of bmask
11	!! 8.0 ! 1997-02 (G. Madec) mesh information put in domhgr.F
12	!! 8.1 ! 1997-07 (G. Madec) modification of mbathy and fmask
13	!! - ! 1998-05 (G. Roullet) free surface
14	!! 8.2 ! 2000-03 (G. Madec) no slip accurate
15	!! - ! 2001-09 (J.-M. Molines) Open boundaries
16	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
17	!! - ! 2005-11 (V. Garnier) Surface pressure gradient organization
18	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
19	!!----------------------------------------------------------------------
20
21	!!----------------------------------------------------------------------
22	!! dom_msk : compute land/ocean mask
23	!! dom_msk_nsa : update land/ocean mask when no-slip accurate option is used.
24	!!----------------------------------------------------------------------
25	USE oce ! ocean dynamics and tracers
26	USE dom_oce ! ocean space and time domain
27	USE obc_oce ! ocean open boundary conditions
28	USE dynvor ! vorticity term
29	USE in_out_manager ! I/O manager
30	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
31	USE lib_mpp
32	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
33
34	IMPLICIT NONE
35	PRIVATE
36
37	PUBLIC dom_msk ! routine called by inidom.F90
38	PUBLIC dom_msk_alloc ! routine called by nemogcm.F90
39
40	! !!* Namelist namlbc : lateral boundary condition *
41	REAL(wp) :: rn_shlat = 2. ! type of lateral boundary condition on velocity
42
43	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icoord ! Workspace for dom_msk_nsa()
44
45	!! * Substitutions
46	# include "vectopt_loop_substitute.h90"
47	!!----------------------------------------------------------------------
48	!! NEMO/OPA 3.2 , LODYC-IPSL (2009)
49	!! $Id$
50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
51	!!----------------------------------------------------------------------
52	CONTAINS
53
54	INTEGER FUNCTION dom_msk_alloc()
55	!!---------------------------------------------------------------------
56	!! * FUNCTION dom_msk_alloc *
57	!!---------------------------------------------------------------------
58	dom_msk_alloc = 0
59	#if defined key_noslip_accurate
60	ALLOCATE(icoord(jpijpjjpk,3), STAT=dom_msk_alloc)
61	#endif
62	IF( dom_msk_alloc /= 0 ) CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
63	!
64	END FUNCTION dom_msk_alloc
65
66
67	SUBROUTINE dom_msk
68	!!---------------------------------------------------------------------
69	!! * ROUTINE dom_msk *
70	!!
71	!! ** Purpose : Compute land/ocean mask arrays at tracer points, hori-
72	!! zontal velocity points (u & v), vorticity points (f) and baro-
73	!! tropic stream function points (b).
74	!!
75	!! ** Method : The ocean/land mask is computed from the basin bathy-
76	!! metry in level (mbathy) which is defined or read in dommba.
77	!! mbathy equals 0 over continental T-point
78	!! and the number of ocean level over the ocean.
79	!!
80	!! At a given position (ji,jj,jk) the ocean/land mask is given by:
81	!! t-point : 0. IF mbathy( ji ,jj) =< 0
82	!! 1. IF mbathy( ji ,jj) >= jk
83	!! u-point : 0. IF mbathy( ji ,jj) or mbathy(ji+1, jj ) =< 0
84	!! 1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
85	!! v-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) =< 0
86	!! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
87	!! f-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1)
88	!! or mbathy(ji+1,jj) or mbathy(ji+1,jj+1) =< 0
89	!! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
90	!! and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
91	!! b-point : the same definition as for f-point of the first ocean
92	!! level (surface level) but with 0 along coastlines.
93	!! tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
94	!! rows/lines due to cyclic or North Fold boundaries as well
95	!! as MPP halos.
96	!!
97	!! The lateral friction is set through the value of fmask along
98	!! the coast and topography. This value is defined by rn_shlat, a
99	!! namelist parameter:
100	!! rn_shlat = 0, free slip (no shear along the coast)
101	!! rn_shlat = 2, no slip (specified zero velocity at the coast)
102	!! 0 < rn_shlat < 2, partial slip \| non-linear velocity profile
103	!! 2 < rn_shlat, strong slip \| in the lateral boundary layer
104	!!
105	!! N.B. If nperio not equal to 0, the land/ocean mask arrays
106	!! are defined with the proper value at lateral domain boundaries,
107	!! but bmask. indeed, bmask defined the domain over which the
108	!! barotropic stream function is computed. this domain cannot
109	!! contain identical columns because the matrix associated with
110	!! the barotropic stream function equation is then no more inverti-
111	!! ble. therefore bmask is set to 0 along lateral domain boundaries
112	!! even IF nperio is not zero.
113	!!
114	!! In case of open boundaries (lk_obc=T):
115	!! - tmask is set to 1 on the points to be computed bay the open
116	!! boundaries routines.
117	!! - bmask is set to 0 on the open boundaries.
118	!!
119	!! ** Action : tmask : land/ocean mask at t-point (=0. or 1.)
120	!! umask : land/ocean mask at u-point (=0. or 1.)
121	!! vmask : land/ocean mask at v-point (=0. or 1.)
122	!! fmask : land/ocean mask at f-point (=0. or 1.)
123	!! =rn_shlat along lateral boundaries
124	!! bmask : land/ocean mask at barotropic stream
125	!! function point (=0. or 1.) and set to 0 along lateral boundaries
126	!! tmask_i : interior ocean mask
127	!!----------------------------------------------------------------------
128	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released
129	USE wrk_nemo, ONLY: zwf => wrk_2d_1 ! 2D real workspace
130	USE wrk_nemo, ONLY: imsk => iwrk_2d_1 ! 2D integer workspace
131	!
132	INTEGER :: ji, jj, jk ! dummy loop indices
133	INTEGER :: iif, iil, ii0, ii1, ii ! local integers
134	INTEGER :: ijf, ijl, ij0, ij1 ! - -
135	!!
136	NAMELIST/namlbc/ rn_shlat
137	!!---------------------------------------------------------------------
138
139	IF( wrk_in_use(2, 1) .OR. iwrk_in_use(2, 1) ) THEN
140	CALL ctl_stop('dom_msk: requested workspace arrays unavailable') ; RETURN
141	ENDIF
142
143	REWIND( numnam ) ! Namelist namlbc : lateral momentum boundary condition
144	READ ( numnam, namlbc )
145
146	IF(lwp) THEN ! control print
147	WRITE(numout,*)
148	WRITE(numout,*) 'dommsk : ocean mask '
149	WRITE(numout,*) '~~~~~~'
150	WRITE(numout,*) ' Namelist namlbc'
151	WRITE(numout,*) ' lateral momentum boundary cond. rn_shlat = ',rn_shlat
152	ENDIF
153
154	IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral free-slip '
155	ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral no-slip '
156	ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip '
157	ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip '
158	ELSE
159	WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
160	CALL ctl_stop( ctmp1 )
161	ENDIF
162
163	! 1. Ocean/land mask at t-point (computed from mbathy)
164	! -----------------------------
165	! N.B. tmask has already the right boundary conditions since mbathy is ok
166	!
167	tmask(:,:,:) = 0._wp
168	DO jk = 1, jpk
169	DO jj = 1, jpj
170	DO ji = 1, jpi
171	IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) tmask(ji,jj,jk) = 1._wp
172	END DO
173	END DO
174	END DO
175
176	!!gm ????
177	#if defined key_zdfkpp
178	IF( cp_cfg == 'orca' ) THEN
179	IF( jp_cfg == 2 ) THEN ! land point on Bab el Mandeb zonal section
180	ij0 = 87 ; ij1 = 88
181	ii0 = 160 ; ii1 = 161
182	tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
183	ELSE
184	IF(lwp) WRITE(numout,*)
185	IF(lwp) WRITE(numout,cform_war)
186	IF(lwp) WRITE(numout,*)
187	IF(lwp) WRITE(numout,*)' A mask must be applied on Bab el Mandeb strait'
188	IF(lwp) WRITE(numout,*)' in case of ORCAs configurations'
189	IF(lwp) WRITE(numout,*)' This is a problem which is not yet solved'
190	IF(lwp) WRITE(numout,*)
191	ENDIF
192	ENDIF
193	#endif
194	!!gm end
195
196	! Interior domain mask (used for global sum)
197	! --------------------
198	tmask_i(:,:) = tmask(:,:,1)
199	iif = jpreci ! ???
200	iil = nlci - jpreci + 1
201	ijf = jprecj ! ???
202	ijl = nlcj - jprecj + 1
203
204	tmask_i( 1 :iif, : ) = 0._wp ! first columns
205	tmask_i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)
206	tmask_i( : , 1 :ijf) = 0._wp ! first rows
207	tmask_i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)
208
209	! north fold mask
210	! ---------------
211	tpol(1:jpiglo) = 1._wp
212	fpol(1:jpiglo) = 1._wp
213	IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot
214	tpol(jpiglo/2+1:jpiglo) = 0._wp
215	fpol( 1 :jpiglo) = 0._wp
216	IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row
217	DO ji = iif+1, iil-1
218	tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
219	END DO
220	ENDIF
221	ENDIF
222	IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot
223	tpol( 1 :jpiglo) = 0._wp
224	fpol(jpiglo/2+1:jpiglo) = 0._wp
225	ENDIF
226
227	! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask)
228	! -------------------------------------------
229	DO jk = 1, jpk
230	DO jj = 1, jpjm1
231	DO ji = 1, fs_jpim1 ! vector loop
232	umask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk)
233	vmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji ,jj+1,jk)
234	END DO
235	DO ji = 1, jpim1 ! NO vector opt.
236	fmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk) &
237	& * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
238	END DO
239	END DO
240	END DO
241	CALL lbc_lnk( umask, 'U', 1._wp ) ! Lateral boundary conditions
242	CALL lbc_lnk( vmask, 'V', 1._wp )
243	CALL lbc_lnk( fmask, 'F', 1._wp )
244
245
246	! 4. ocean/land mask for the elliptic equation
247	! --------------------------------------------
248	bmask(:,:) = tmask(:,:,1) ! elliptic equation is written at t-point
249	!
250	! ! Boundary conditions
251	! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
252	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
253	bmask( 1 ,:) = 0._wp
254	bmask(jpi,:) = 0._wp
255	ENDIF
256	IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1
257	bmask(:, 1 ) = 0._wp
258	ENDIF
259	! ! north fold :
260	IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
261	DO ji = 1, jpi
262	ii = ji + nimpp - 1
263	bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
264	bmask(ji,jpj ) = 0._wp
265	END DO
266	ENDIF
267	IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
268	bmask(:,jpj) = 0._wp
269	ENDIF
270	!
271	IF( lk_mpp ) THEN ! mpp specificities
272	! ! bmask is set to zero on the overlap region
273	IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp
274	IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp
275	IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp
276	IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp
277	!
278	IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
279	DO ji = 1, nlci
280	ii = ji + nimpp - 1
281	bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
282	bmask(ji,nlcj ) = 0._wp
283	END DO
284	ENDIF
285	IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
286	DO ji = 1, nlci
287	bmask(ji,nlcj ) = 0._wp
288	END DO
289	ENDIF
290	ENDIF
291
292
293	! mask for second order calculation of vorticity
294	! ----------------------------------------------
295	CALL dom_msk_nsa
296
297
298	! Lateral boundary conditions on velocity (modify fmask)
299	! ---------------------------------------
300	DO jk = 1, jpk
301	zwf(:,:) = fmask(:,:,jk)
302	DO jj = 2, jpjm1
303	DO ji = fs_2, fs_jpim1 ! vector opt.
304	IF( fmask(ji,jj,jk) == 0._wp ) THEN
305	fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
306	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
307	ENDIF
308	END DO
309	END DO
310	DO jj = 2, jpjm1
311	IF( fmask(1,jj,jk) == 0._wp ) THEN
312	fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
313	ENDIF
314	IF( fmask(jpi,jj,jk) == 0._wp ) THEN
315	fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
316	ENDIF
317	END DO
318	DO ji = 2, jpim1
319	IF( fmask(ji,1,jk) == 0._wp ) THEN
320	fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
321	ENDIF
322	IF( fmask(ji,jpj,jk) == 0._wp ) THEN
323	fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
324	ENDIF
325	END DO
326	END DO
327	!
328	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration
329	! ! Increased lateral friction near of some straits
330	IF( nn_cla == 0 ) THEN
331	! ! Gibraltar strait : partial slip (fmask=0.5)
332	ij0 = 101 ; ij1 = 101
333	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
334	ij0 = 102 ; ij1 = 102
335	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
336	!
337	! ! Bab el Mandeb : partial slip (fmask=1)
338	ij0 = 87 ; ij1 = 88
339	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
340	ij0 = 88 ; ij1 = 88
341	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
342	!
343	ENDIF
344	! ! Danish straits : strong slip (fmask > 2)
345	! We keep this as an example but it is instable in this case
346	! ij0 = 115 ; ij1 = 115
347	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
348	! ij0 = 116 ; ij1 = 116
349	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
350	!
351	ENDIF
352	!
353	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
354	! ! Increased lateral friction near of some straits
355	IF(lwp) WRITE(numout,*)
356	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
357	IF(lwp) WRITE(numout,*) ' Gibraltar '
358	ii0 = 283 ; ii1 = 284 ! Gibraltar Strait
359	ij0 = 200 ; ij1 = 200 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
360
361	IF(lwp) WRITE(numout,*) ' Bhosporus '
362	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
363	ij0 = 208 ; ij1 = 208 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
364
365	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
366	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
367	ij0 = 149 ; ij1 = 150 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
368
369	IF(lwp) WRITE(numout,*) ' Lombok '
370	ii0 = 44 ; ii1 = 44 ! Lombok Strait
371	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
372
373	IF(lwp) WRITE(numout,*) ' Ombai '
374	ii0 = 53 ; ii1 = 53 ! Ombai Strait
375	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
376
377	IF(lwp) WRITE(numout,*) ' Timor Passage '
378	ii0 = 56 ; ii1 = 56 ! Timor Passage
379	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
380
381	IF(lwp) WRITE(numout,*) ' West Halmahera '
382	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
383	ij0 = 141 ; ij1 = 142 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
384
385	IF(lwp) WRITE(numout,*) ' East Halmahera '
386	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
387	ij0 = 141 ; ij1 = 142 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
388	!
389	ENDIF
390	!
391	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
392
393
394	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
395	! at angles with three ocean points and one land point
396	IF( ln_dynvor_con .OR. ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix) THEN
397	DO jj = 2, jpjm1
398	DO ji = 2, jpim1
399	IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) THEN
400	fmask(ji,jj,jk) = 1._wp
401	ENDIF
402	END DO
403	END DO
404	ENDIF
405
406
407	IF( nprint == 1 .AND. lwp ) THEN ! Control print
408	imsk(:,:) = INT( tmask_i(:,:) )
409	WRITE(numout,*) ' tmask_i : '
410	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
411	& 1, jpj, 1, 1, numout)
412	WRITE (numout,*)
413	WRITE (numout,*) ' dommsk: tmask for each level'
414	WRITE (numout,*) ' ----------------------------'
415	DO jk = 1, jpk
416	imsk(:,:) = INT( tmask(:,:,jk) )
417
418	WRITE(numout,*)
419	WRITE(numout,*) ' level = ',jk
420	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
421	& 1, jpj, 1, 1, numout)
422	END DO
423	WRITE(numout,*)
424	WRITE(numout,*) ' dom_msk: vmask for each level'
425	WRITE(numout,*) ' -----------------------------'
426	DO jk = 1, jpk
427	imsk(:,:) = INT( vmask(:,:,jk) )
428	WRITE(numout,*)
429	WRITE(numout,*) ' level = ',jk
430	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
431	& 1, jpj, 1, 1, numout)
432	END DO
433	WRITE(numout,*)
434	WRITE(numout,*) ' dom_msk: fmask for each level'
435	WRITE(numout,*) ' -----------------------------'
436	DO jk = 1, jpk
437	imsk(:,:) = INT( fmask(:,:,jk) )
438	WRITE(numout,*)
439	WRITE(numout,*) ' level = ',jk
440	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
441	& 1, jpj, 1, 1, numout )
442	END DO
443	WRITE(numout,*)
444	WRITE(numout,*) ' dom_msk: bmask '
445	WRITE(numout,*) ' ---------------'
446	WRITE(numout,*)
447	imsk(:,:) = INT( bmask(:,:) )
448	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
449	& 1, jpj, 1, 1, numout )
450	ENDIF
451	!
452	IF( wrk_not_released(2, 1) .OR. &
453	iwrk_not_released(2, 1) ) CALL ctl_stop('dom_msk: failed to release workspace arrays')
454	!
455	END SUBROUTINE dom_msk
456
457	#if defined key_noslip_accurate
458	!!----------------------------------------------------------------------
459	!! 'key_noslip_accurate' : accurate no-slip boundary condition
460	!!----------------------------------------------------------------------
461
462	SUBROUTINE dom_msk_nsa
463	!!---------------------------------------------------------------------
464	!! * ROUTINE dom_msk_nsa *
465	!!
466	!! ** Purpose :
467	!!
468	!! ** Method :
469	!!
470	!! ** Action :
471	!!----------------------------------------------------------------------
472	INTEGER :: ji, jj, jk, jl ! dummy loop indices
473	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
474	REAL(wp) :: zaa
475	!!---------------------------------------------------------------------
476
477	IF(lwp) WRITE(numout,*)
478	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
479	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
480	IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
481
482	! mask for second order calculation of vorticity
483	! ----------------------------------------------
484	! noslip boundary condition: fmask=1 at convex corner, store
485	! index of straight coast meshes ( 'west', refering to a coast,
486	! means west of the ocean, aso)
487
488	DO jk = 1, jpk
489	DO jl = 1, 4
490	npcoa(jl,jk) = 0
491	DO ji = 1, 2*(jpi+jpj)
492	nicoa(ji,jl,jk) = 0
493	njcoa(ji,jl,jk) = 0
494	END DO
495	END DO
496	END DO
497
498	IF( jperio == 2 ) THEN
499	WRITE(numout,*) ' '
500	WRITE(numout,*) ' symetric boundary conditions need special'
501	WRITE(numout,*) ' treatment not implemented. we stop.'
502	STOP
503	ENDIF
504
505	! convex corners
506
507	DO jk = 1, jpkm1
508	DO jj = 1, jpjm1
509	DO ji = 1, jpim1
510	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
511	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
512	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
513	END DO
514	END DO
515	END DO
516
517	! north-south straight coast
518
519	DO jk = 1, jpkm1
520	inw = 0
521	ine = 0
522	DO jj = 2, jpjm1
523	DO ji = 2, jpim1
524	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
525	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
526	inw = inw + 1
527	nicoa(inw,1,jk) = ji
528	njcoa(inw,1,jk) = jj
529	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
530	ENDIF
531	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
532	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
533	ine = ine + 1
534	nicoa(ine,2,jk) = ji
535	njcoa(ine,2,jk) = jj
536	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
537	ENDIF
538	END DO
539	END DO
540	npcoa(1,jk) = inw
541	npcoa(2,jk) = ine
542	END DO
543
544	! west-east straight coast
545
546	DO jk = 1, jpkm1
547	ins = 0
548	inn = 0
549	DO jj = 2, jpjm1
550	DO ji =2, jpim1
551	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
552	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
553	ins = ins + 1
554	nicoa(ins,3,jk) = ji
555	njcoa(ins,3,jk) = jj
556	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
557	ENDIF
558	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
559	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
560	inn = inn + 1
561	nicoa(inn,4,jk) = ji
562	njcoa(inn,4,jk) = jj
563	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
564	ENDIF
565	END DO
566	END DO
567	npcoa(3,jk) = ins
568	npcoa(4,jk) = inn
569	END DO
570
571	itest = 2 * ( jpi + jpj )
572	DO jk = 1, jpk
573	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
574	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
575
576	WRITE(ctmp1,*) ' level jk = ',jk
577	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
578	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
579	& npcoa(3,jk), npcoa(4,jk)
580	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
581	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
582	ENDIF
583	END DO
584
585	ierror = 0
586	iind = 0
587	ijnd = 0
588	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
589	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
590	DO jk = 1, jpk
591	DO jl = 1, npcoa(1,jk)
592	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
593	ierror = ierror+1
594	icoord(ierror,1) = nicoa(jl,1,jk)
595	icoord(ierror,2) = njcoa(jl,1,jk)
596	icoord(ierror,3) = jk
597	ENDIF
598	END DO
599	DO jl = 1, npcoa(2,jk)
600	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
601	ierror = ierror + 1
602	icoord(ierror,1) = nicoa(jl,2,jk)
603	icoord(ierror,2) = njcoa(jl,2,jk)
604	icoord(ierror,3) = jk
605	ENDIF
606	END DO
607	DO jl = 1, npcoa(3,jk)
608	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
609	ierror = ierror + 1
610	icoord(ierror,1) = nicoa(jl,3,jk)
611	icoord(ierror,2) = njcoa(jl,3,jk)
612	icoord(ierror,3) = jk
613	ENDIF
614	END DO
615	DO jl = 1, npcoa(4,jk)
616	IF( njcoa(jl,4,jk)-2 < 1) THEN
617	ierror=ierror + 1
618	icoord(ierror,1) = nicoa(jl,4,jk)
619	icoord(ierror,2) = njcoa(jl,4,jk)
620	icoord(ierror,3) = jk
621	ENDIF
622	END DO
623	END DO
624
625	IF( ierror > 0 ) THEN
626	IF(lwp) WRITE(numout,*)
627	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
628	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
629	DO jl = 1, ierror
630	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
631	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
632	END DO
633	CALL ctl_stop( 'We stop...' )
634	ENDIF
635	!
636	END SUBROUTINE dom_msk_nsa
637
638	#else
639	!!----------------------------------------------------------------------
640	!! Default option : Empty routine
641	!!----------------------------------------------------------------------
642	SUBROUTINE dom_msk_nsa
643	END SUBROUTINE dom_msk_nsa
644	#endif
645
646	!!======================================================================
647	END MODULE dommsk

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

Download in other formats: