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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2633

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

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

Download in other formats: