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

Context Navigation

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

Initial restructuring.
Use fldread to read open boundary data.

Property svn:keywords set to Id

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

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

Download in other formats: