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 @ 3294

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

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

Download in other formats: