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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 8280

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

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

Download in other formats: