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

Context Navigation

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 4726

Last change on this file since 4726 was 4726, checked in by mathiot, 10 years ago
ISF branch: change name of 2 variables (icedep => risfdep and lmask => ssmask), cosmetic changes and add ldfslp key
Property svn:keywords set to `Id`
File size: 29.8 KB

Rev	Line
[3]	1	MODULE dommsk
[1566]	2	!!======================================================================
[3]	3	!! * MODULE dommsk *
	4	!! Ocean initialization : domain land/sea mask
[1566]	5	!!======================================================================
	6	!! History : OPA ! 1987-07 (G. Madec) Original code
[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 *
[4147]	41	REAL(wp) :: rn_shlat ! type of lateral boundary condition on velocity
	42	LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition
[3294]	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	!!
[4328]	117	!! In case of open boundaries (lk_bdy=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 ! - -
[4147]	135	INTEGER :: ios
[3294]	136	INTEGER , POINTER, DIMENSION(:,:) :: imsk
	137	REAL(wp), POINTER, DIMENSION(:,:) :: zwf
[1601]	138	!!
[3294]	139	NAMELIST/namlbc/ rn_shlat, ln_vorlat
[3]	140	!!---------------------------------------------------------------------
[3294]	141	!
	142	IF( nn_timing == 1 ) CALL timing_start('dom_msk')
	143	!
	144	CALL wrk_alloc( jpi, jpj, imsk )
	145	CALL wrk_alloc( jpi, jpj, zwf )
	146	!
[4147]	147	REWIND( numnam_ref ) ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
	148	READ ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
	149	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
	150
	151	REWIND( numnam_cfg ) ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
	152	READ ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
	153	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
[4624]	154	IF(lwm) WRITE ( numond, namlbc )
[1566]	155
	156	IF(lwp) THEN ! control print
[3]	157	WRITE(numout,*)
	158	WRITE(numout,*) 'dommsk : ocean mask '
	159	WRITE(numout,*) '~~~~~~'
[1566]	160	WRITE(numout,*) ' Namelist namlbc'
[3294]	161	WRITE(numout,*) ' lateral momentum boundary cond. rn_shlat = ',rn_shlat
	162	WRITE(numout,*) ' consistency with analytical form ln_vorlat = ',ln_vorlat
[3]	163	ENDIF
	164
[2528]	165	IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral free-slip '
[1601]	166	ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral no-slip '
	167	ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip '
	168	ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip '
	169	ELSE
	170	WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
	171	CALL ctl_stop( ctmp1 )
[3]	172	ENDIF
	173
	174	! 1. Ocean/land mask at t-point (computed from mbathy)
	175	! -----------------------------
[1566]	176	! N.B. tmask has already the right boundary conditions since mbathy is ok
	177	!
[2528]	178	tmask(:,:,:) = 0._wp
[3]	179	DO jk = 1, jpk
	180	DO jj = 1, jpj
	181	DO ji = 1, jpi
[2528]	182	IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) tmask(ji,jj,jk) = 1._wp
[3]	183	END DO
	184	END DO
	185	END DO
[4666]	186
	187	! (ISF) define barotropic mask and mask the ice shelf point
[4726]	188	ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
[4666]	189
	190	DO jk = 1, jpk
	191	DO jj = 1, jpj
	192	DO ji = 1, jpi
[4726]	193	IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp ) THEN
[4666]	194	tmask(ji,jj,jk) = 0._wp
	195	END IF
	196	END DO
	197	END DO
	198	END DO
[3]	199
[1566]	200	!!gm ????
[255]	201	#if defined key_zdfkpp
[1601]	202	IF( cp_cfg == 'orca' ) THEN
	203	IF( jp_cfg == 2 ) THEN ! land point on Bab el Mandeb zonal section
[291]	204	ij0 = 87 ; ij1 = 88
	205	ii0 = 160 ; ii1 = 161
[2528]	206	tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
[291]	207	ELSE
	208	IF(lwp) WRITE(numout,*)
	209	IF(lwp) WRITE(numout,cform_war)
	210	IF(lwp) WRITE(numout,*)
	211	IF(lwp) WRITE(numout,*)' A mask must be applied on Bab el Mandeb strait'
	212	IF(lwp) WRITE(numout,*)' in case of ORCAs configurations'
	213	IF(lwp) WRITE(numout,*)' This is a problem which is not yet solved'
	214	IF(lwp) WRITE(numout,*)
	215	ENDIF
	216	ENDIF
[255]	217	#endif
[1566]	218	!!gm end
[291]	219
[3]	220	! Interior domain mask (used for global sum)
	221	! --------------------
[4726]	222	tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf
[3]	223	iif = jpreci ! ???
	224	iil = nlci - jpreci + 1
	225	ijf = jprecj ! ???
	226	ijl = nlcj - jprecj + 1
	227
[2528]	228	tmask_i( 1 :iif, : ) = 0._wp ! first columns
	229	tmask_i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)
	230	tmask_i( : , 1 :ijf) = 0._wp ! first rows
	231	tmask_i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)
[3]	232
	233	! north fold mask
[1566]	234	! ---------------
[2528]	235	tpol(1:jpiglo) = 1._wp
	236	fpol(1:jpiglo) = 1._wp
[3]	237	IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot
[2528]	238	tpol(jpiglo/2+1:jpiglo) = 0._wp
	239	fpol( 1 :jpiglo) = 0._wp
[1566]	240	IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row
[291]	241	DO ji = iif+1, iil-1
	242	tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
	243	END DO
	244	ENDIF
[3]	245	ENDIF
	246	IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot
[2528]	247	tpol( 1 :jpiglo) = 0._wp
	248	fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]	249	ENDIF
	250
	251	! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask)
	252	! -------------------------------------------
	253	DO jk = 1, jpk
	254	DO jj = 1, jpjm1
	255	DO ji = 1, fs_jpim1 ! vector loop
	256	umask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk)
	257	vmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji ,jj+1,jk)
[1271]	258	END DO
[1694]	259	DO ji = 1, jpim1 ! NO vector opt.
[3]	260	fmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk) &
[62]	261	& * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]	262	END DO
	263	END DO
	264	END DO
[4666]	265	! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
	266	DO jj = 1, jpjm1
	267	DO ji = 1, fs_jpim1 ! vector loop
[4726]	268	umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))
	269	vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))
[4666]	270	END DO
	271	DO ji = 1, jpim1 ! NO vector opt.
[4726]	272	fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &
	273	& * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
[4666]	274	END DO
	275	END DO
[2528]	276	CALL lbc_lnk( umask, 'U', 1._wp ) ! Lateral boundary conditions
	277	CALL lbc_lnk( vmask, 'V', 1._wp )
	278	CALL lbc_lnk( fmask, 'F', 1._wp )
[4666]	279	CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions
	280	CALL lbc_lnk( vmask_i, 'V', 1._wp )
	281	CALL lbc_lnk( fmask_i, 'F', 1._wp )
[3]	282
	283
	284	! 4. ocean/land mask for the elliptic equation
	285	! --------------------------------------------
[4726]	286	bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point
[1566]	287	!
	288	! ! Boundary conditions
	289	! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
[3]	290	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
[2528]	291	bmask( 1 ,:) = 0._wp
	292	bmask(jpi,:) = 0._wp
[3]	293	ENDIF
[1566]	294	IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1
[2528]	295	bmask(:, 1 ) = 0._wp
[3]	296	ENDIF
[1566]	297	! ! north fold :
	298	IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
	299	DO ji = 1, jpi
[1528]	300	ii = ji + nimpp - 1
	301	bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
[2528]	302	bmask(ji,jpj ) = 0._wp
[1528]	303	END DO
[3]	304	ENDIF
[1566]	305	IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[2528]	306	bmask(:,jpj) = 0._wp
[3]	307	ENDIF
[1566]	308	!
	309	IF( lk_mpp ) THEN ! mpp specificities
	310	! ! bmask is set to zero on the overlap region
[2528]	311	IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp
	312	IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp
	313	IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp
	314	IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp
[1566]	315	!
	316	IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
[1528]	317	DO ji = 1, nlci
	318	ii = ji + nimpp - 1
	319	bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
[2528]	320	bmask(ji,nlcj ) = 0._wp
[1528]	321	END DO
[32]	322	ENDIF
[1566]	323	IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[1528]	324	DO ji = 1, nlci
[2528]	325	bmask(ji,nlcj ) = 0._wp
[1528]	326	END DO
[32]	327	ENDIF
[3]	328	ENDIF
	329
	330
	331	! mask for second order calculation of vorticity
	332	! ----------------------------------------------
	333	CALL dom_msk_nsa
	334
	335
	336	! Lateral boundary conditions on velocity (modify fmask)
[1566]	337	! ---------------------------------------
[3]	338	DO jk = 1, jpk
[1566]	339	zwf(:,:) = fmask(:,:,jk)
[3]	340	DO jj = 2, jpjm1
	341	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	342	IF( fmask(ji,jj,jk) == 0._wp ) THEN
	343	fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	344	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
[3]	345	ENDIF
	346	END DO
	347	END DO
	348	DO jj = 2, jpjm1
[2528]	349	IF( fmask(1,jj,jk) == 0._wp ) THEN
	350	fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]	351	ENDIF
[2528]	352	IF( fmask(jpi,jj,jk) == 0._wp ) THEN
	353	fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]	354	ENDIF
[1566]	355	END DO
[3]	356	DO ji = 2, jpim1
[2528]	357	IF( fmask(ji,1,jk) == 0._wp ) THEN
	358	fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]	359	ENDIF
[2528]	360	IF( fmask(ji,jpj,jk) == 0._wp ) THEN
	361	fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]	362	ENDIF
	363	END DO
	364	END DO
[1566]	365	!
	366	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration
	367	! ! Increased lateral friction near of some straits
[2528]	368	IF( nn_cla == 0 ) THEN
[1273]	369	! ! Gibraltar strait : partial slip (fmask=0.5)
	370	ij0 = 101 ; ij1 = 101
[2528]	371	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[32]	372	ij0 = 102 ; ij1 = 102
[2528]	373	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[1273]	374	!
	375	! ! Bab el Mandeb : partial slip (fmask=1)
	376	ij0 = 87 ; ij1 = 88
[2528]	377	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	378	ij0 = 88 ; ij1 = 88
[2528]	379	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	380	!
[3]	381	ENDIF
[1707]	382	! ! Danish straits : strong slip (fmask > 2)
	383	! We keep this as an example but it is instable in this case
	384	! ij0 = 115 ; ij1 = 115
[2528]	385	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	386	! ij0 = 116 ; ij1 = 116
[2528]	387	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	388	!
	389	ENDIF
[1566]	390	!
[2528]	391	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	392	! ! Increased lateral friction near of some straits
	393	IF(lwp) WRITE(numout,*)
	394	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	395	IF(lwp) WRITE(numout,*) ' Gibraltar '
	396	ii0 = 283 ; ii1 = 284 ! Gibraltar Strait
	397	ij0 = 200 ; ij1 = 200 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
[1566]	398
[2528]	399	IF(lwp) WRITE(numout,*) ' Bhosporus '
	400	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
	401	ij0 = 208 ; ij1 = 208 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	402
	403	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
	404	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
	405	ij0 = 149 ; ij1 = 150 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
	406
	407	IF(lwp) WRITE(numout,*) ' Lombok '
	408	ii0 = 44 ; ii1 = 44 ! Lombok Strait
	409	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	410
	411	IF(lwp) WRITE(numout,*) ' Ombai '
	412	ii0 = 53 ; ii1 = 53 ! Ombai Strait
	413	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	414
	415	IF(lwp) WRITE(numout,*) ' Timor Passage '
	416	ii0 = 56 ; ii1 = 56 ! Timor Passage
	417	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	418
	419	IF(lwp) WRITE(numout,*) ' West Halmahera '
	420	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
	421	ij0 = 141 ; ij1 = 142 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
	422
	423	IF(lwp) WRITE(numout,*) ' East Halmahera '
	424	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
	425	ij0 = 141 ; ij1 = 142 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
	426	!
	427	ENDIF
	428	!
	429	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	430
[3294]	431	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	432
[1566]	433	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	434	imsk(:,:) = INT( tmask_i(:,:) )
	435	WRITE(numout,*) ' tmask_i : '
	436	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	437	& 1, jpj, 1, 1, numout)
	438	WRITE (numout,*)
	439	WRITE (numout,*) ' dommsk: tmask for each level'
	440	WRITE (numout,*) ' ----------------------------'
	441	DO jk = 1, jpk
	442	imsk(:,:) = INT( tmask(:,:,jk) )
	443
	444	WRITE(numout,*)
	445	WRITE(numout,*) ' level = ',jk
	446	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	447	& 1, jpj, 1, 1, numout)
	448	END DO
	449	WRITE(numout,*)
	450	WRITE(numout,*) ' dom_msk: vmask for each level'
	451	WRITE(numout,*) ' -----------------------------'
	452	DO jk = 1, jpk
	453	imsk(:,:) = INT( vmask(:,:,jk) )
	454	WRITE(numout,*)
	455	WRITE(numout,*) ' level = ',jk
	456	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	457	& 1, jpj, 1, 1, numout)
	458	END DO
	459	WRITE(numout,*)
	460	WRITE(numout,*) ' dom_msk: fmask for each level'
	461	WRITE(numout,*) ' -----------------------------'
	462	DO jk = 1, jpk
	463	imsk(:,:) = INT( fmask(:,:,jk) )
	464	WRITE(numout,*)
	465	WRITE(numout,*) ' level = ',jk
	466	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	467	& 1, jpj, 1, 1, numout )
	468	END DO
	469	WRITE(numout,*)
	470	WRITE(numout,*) ' dom_msk: bmask '
	471	WRITE(numout,*) ' ---------------'
	472	WRITE(numout,*)
	473	imsk(:,:) = INT( bmask(:,:) )
	474	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	475	& 1, jpj, 1, 1, numout )
[3]	476	ENDIF
[1566]	477	!
[3294]	478	CALL wrk_dealloc( jpi, jpj, imsk )
	479	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	480	!
[3294]	481	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	482	!
[3]	483	END SUBROUTINE dom_msk
	484
	485	#if defined key_noslip_accurate
	486	!!----------------------------------------------------------------------
	487	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	488	!!----------------------------------------------------------------------
	489
	490	SUBROUTINE dom_msk_nsa
	491	!!---------------------------------------------------------------------
	492	!! * ROUTINE dom_msk_nsa *
	493	!!
	494	!! ** Purpose :
	495	!!
	496	!! ** Method :
	497	!!
	498	!! ** Action :
	499	!!----------------------------------------------------------------------
[2715]	500	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	501	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	502	REAL(wp) :: zaa
[3]	503	!!---------------------------------------------------------------------
[3294]	504	!
	505	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	506	!
[2715]	507	IF(lwp) WRITE(numout,*)
	508	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	509	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
	510	IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
[3]	511
	512	! mask for second order calculation of vorticity
	513	! ----------------------------------------------
	514	! noslip boundary condition: fmask=1 at convex corner, store
	515	! index of straight coast meshes ( 'west', refering to a coast,
	516	! means west of the ocean, aso)
	517
	518	DO jk = 1, jpk
	519	DO jl = 1, 4
	520	npcoa(jl,jk) = 0
	521	DO ji = 1, 2*(jpi+jpj)
	522	nicoa(ji,jl,jk) = 0
	523	njcoa(ji,jl,jk) = 0
	524	END DO
	525	END DO
	526	END DO
	527
	528	IF( jperio == 2 ) THEN
	529	WRITE(numout,*) ' '
	530	WRITE(numout,*) ' symetric boundary conditions need special'
	531	WRITE(numout,*) ' treatment not implemented. we stop.'
	532	STOP
	533	ENDIF
	534
	535	! convex corners
	536
	537	DO jk = 1, jpkm1
	538	DO jj = 1, jpjm1
	539	DO ji = 1, jpim1
	540	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	541	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	542	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	543	END DO
	544	END DO
	545	END DO
	546
	547	! north-south straight coast
	548
	549	DO jk = 1, jpkm1
	550	inw = 0
	551	ine = 0
	552	DO jj = 2, jpjm1
	553	DO ji = 2, jpim1
	554	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	555	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	556	inw = inw + 1
	557	nicoa(inw,1,jk) = ji
	558	njcoa(inw,1,jk) = jj
	559	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	560	ENDIF
	561	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	562	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	563	ine = ine + 1
	564	nicoa(ine,2,jk) = ji
	565	njcoa(ine,2,jk) = jj
	566	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	567	ENDIF
	568	END DO
	569	END DO
	570	npcoa(1,jk) = inw
	571	npcoa(2,jk) = ine
	572	END DO
	573
	574	! west-east straight coast
	575
	576	DO jk = 1, jpkm1
	577	ins = 0
	578	inn = 0
	579	DO jj = 2, jpjm1
	580	DO ji =2, jpim1
	581	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	582	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	583	ins = ins + 1
	584	nicoa(ins,3,jk) = ji
	585	njcoa(ins,3,jk) = jj
	586	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	587	ENDIF
	588	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	589	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	590	inn = inn + 1
	591	nicoa(inn,4,jk) = ji
	592	njcoa(inn,4,jk) = jj
	593	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	594	ENDIF
	595	END DO
	596	END DO
	597	npcoa(3,jk) = ins
	598	npcoa(4,jk) = inn
	599	END DO
	600
	601	itest = 2 * ( jpi + jpj )
	602	DO jk = 1, jpk
	603	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	604	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	605
	606	WRITE(ctmp1,*) ' level jk = ',jk
	607	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	608	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	609	& npcoa(3,jk), npcoa(4,jk)
[474]	610	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	611	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	612	ENDIF
	613	END DO
	614
	615	ierror = 0
	616	iind = 0
	617	ijnd = 0
[2528]	618	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	619	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	620	DO jk = 1, jpk
	621	DO jl = 1, npcoa(1,jk)
	622	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	623	ierror = ierror+1
	624	icoord(ierror,1) = nicoa(jl,1,jk)
	625	icoord(ierror,2) = njcoa(jl,1,jk)
	626	icoord(ierror,3) = jk
	627	ENDIF
	628	END DO
	629	DO jl = 1, npcoa(2,jk)
	630	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	631	ierror = ierror + 1
	632	icoord(ierror,1) = nicoa(jl,2,jk)
	633	icoord(ierror,2) = njcoa(jl,2,jk)
	634	icoord(ierror,3) = jk
	635	ENDIF
	636	END DO
	637	DO jl = 1, npcoa(3,jk)
	638	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	639	ierror = ierror + 1
	640	icoord(ierror,1) = nicoa(jl,3,jk)
	641	icoord(ierror,2) = njcoa(jl,3,jk)
	642	icoord(ierror,3) = jk
	643	ENDIF
	644	END DO
[2528]	645	DO jl = 1, npcoa(4,jk)
[3]	646	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	647	ierror=ierror + 1
	648	icoord(ierror,1) = nicoa(jl,4,jk)
	649	icoord(ierror,2) = njcoa(jl,4,jk)
	650	icoord(ierror,3) = jk
[3]	651	ENDIF
	652	END DO
	653	END DO
	654
	655	IF( ierror > 0 ) THEN
	656	IF(lwp) WRITE(numout,*)
	657	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	658	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	659	DO jl = 1, ierror
	660	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	661	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	662	END DO
[474]	663	CALL ctl_stop( 'We stop...' )
[3]	664	ENDIF
[2715]	665	!
[3294]	666	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	667	!
[3]	668	END SUBROUTINE dom_msk_nsa
	669
	670	#else
	671	!!----------------------------------------------------------------------
	672	!! Default option : Empty routine
	673	!!----------------------------------------------------------------------
	674	SUBROUTINE dom_msk_nsa
	675	END SUBROUTINE dom_msk_nsa
	676	#endif
	677
	678	!!======================================================================
	679	END MODULE dommsk

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

Download in other formats: