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_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5600

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

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

Download in other formats: