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

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

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

Download in other formats: