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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago
Update NEMOGCM from branch nemo_v3_3_beta
Property svn:keywords set to `Id`
File size: 26.7 KB

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

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

Download in other formats: