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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2392

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

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

Download in other formats: