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

Context Navigation

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5770

Last change on this file since 5770 was 5770, checked in by gm, 9 years ago
#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)
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
[5770]	19	!! 3.6 ! 2015-05 (P. Mathiot) ISF: add wmask,wumask and wvmask
[1566]	20	!!----------------------------------------------------------------------
[3]	21
	22	!!----------------------------------------------------------------------
	23	!! dom_msk : compute land/ocean mask
[1566]	24	!! dom_msk_nsa : update land/ocean mask when no-slip accurate option is used.
[3]	25	!!----------------------------------------------------------------------
	26	USE oce ! ocean dynamics and tracers
	27	USE dom_oce ! ocean space and time domain
	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
[3294]	32	USE wrk_nemo ! Memory allocation
	33	USE timing ! Timing
[3]	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
[2715]	38	PUBLIC dom_msk ! routine called by inidom.F90
	39	PUBLIC dom_msk_alloc ! routine called by nemogcm.F90
[3]	40
[1601]	41	! !!* Namelist namlbc : lateral boundary condition *
[4147]	42	REAL(wp) :: rn_shlat ! type of lateral boundary condition on velocity
	43	LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition
[3294]	44	! with analytical eqs.
[2715]	45
[3294]	46
[2715]	47	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icoord ! Workspace for dom_msk_nsa()
	48
[3]	49	!! * Substitutions
	50	# include "vectopt_loop_substitute.h90"
[1566]	51	!!----------------------------------------------------------------------
	52	!! NEMO/OPA 3.2 , LODYC-IPSL (2009)
[1152]	53	!! $Id$
[2528]	54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1566]	55	!!----------------------------------------------------------------------
[3]	56	CONTAINS
	57
[2715]	58	INTEGER FUNCTION dom_msk_alloc()
	59	!!---------------------------------------------------------------------
	60	!! * FUNCTION dom_msk_alloc *
	61	!!---------------------------------------------------------------------
	62	dom_msk_alloc = 0
	63	#if defined key_noslip_accurate
	64	ALLOCATE(icoord(jpijpjjpk,3), STAT=dom_msk_alloc)
	65	#endif
	66	IF( dom_msk_alloc /= 0 ) CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
	67	!
	68	END FUNCTION dom_msk_alloc
	69
	70
[3]	71	SUBROUTINE dom_msk
	72	!!---------------------------------------------------------------------
	73	!! * ROUTINE dom_msk *
	74	!!
	75	!! ** Purpose : Compute land/ocean mask arrays at tracer points, hori-
	76	!! zontal velocity points (u & v), vorticity points (f) and baro-
	77	!! tropic stream function points (b).
	78	!!
	79	!! ** Method : The ocean/land mask is computed from the basin bathy-
	80	!! metry in level (mbathy) which is defined or read in dommba.
[1528]	81	!! mbathy equals 0 over continental T-point
	82	!! and the number of ocean level over the ocean.
[3]	83	!!
	84	!! At a given position (ji,jj,jk) the ocean/land mask is given by:
	85	!! t-point : 0. IF mbathy( ji ,jj) =< 0
	86	!! 1. IF mbathy( ji ,jj) >= jk
	87	!! u-point : 0. IF mbathy( ji ,jj) or mbathy(ji+1, jj ) =< 0
	88	!! 1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
	89	!! v-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) =< 0
	90	!! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
	91	!! f-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1)
	92	!! or mbathy(ji+1,jj) or mbathy(ji+1,jj+1) =< 0
	93	!! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
[2528]	94	!! and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
[3]	95	!! b-point : the same definition as for f-point of the first ocean
	96	!! level (surface level) but with 0 along coastlines.
[2528]	97	!! tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
	98	!! rows/lines due to cyclic or North Fold boundaries as well
	99	!! as MPP halos.
[3]	100	!!
	101	!! The lateral friction is set through the value of fmask along
[1601]	102	!! the coast and topography. This value is defined by rn_shlat, a
[3]	103	!! namelist parameter:
[1601]	104	!! rn_shlat = 0, free slip (no shear along the coast)
	105	!! rn_shlat = 2, no slip (specified zero velocity at the coast)
	106	!! 0 < rn_shlat < 2, partial slip \| non-linear velocity profile
	107	!! 2 < rn_shlat, strong slip \| in the lateral boundary layer
[3]	108	!!
	109	!! N.B. If nperio not equal to 0, the land/ocean mask arrays
	110	!! are defined with the proper value at lateral domain boundaries,
	111	!! but bmask. indeed, bmask defined the domain over which the
	112	!! barotropic stream function is computed. this domain cannot
	113	!! contain identical columns because the matrix associated with
	114	!! the barotropic stream function equation is then no more inverti-
	115	!! ble. therefore bmask is set to 0 along lateral domain boundaries
	116	!! even IF nperio is not zero.
	117	!!
[4328]	118	!! In case of open boundaries (lk_bdy=T):
[3]	119	!! - tmask is set to 1 on the points to be computed bay the open
	120	!! boundaries routines.
	121	!! - bmask is set to 0 on the open boundaries.
	122	!!
[1566]	123	!! ** Action : tmask : land/ocean mask at t-point (=0. or 1.)
	124	!! umask : land/ocean mask at u-point (=0. or 1.)
	125	!! vmask : land/ocean mask at v-point (=0. or 1.)
	126	!! fmask : land/ocean mask at f-point (=0. or 1.)
[1601]	127	!! =rn_shlat along lateral boundaries
[1566]	128	!! bmask : land/ocean mask at barotropic stream
	129	!! function point (=0. or 1.) and set to 0 along lateral boundaries
[2528]	130	!! tmask_i : interior ocean mask
[3]	131	!!----------------------------------------------------------------------
[5770]	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
	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
[4990]	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	! --------------------
[4990]	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
[4990]	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 )
[4990]	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
[5120]	284	! 3. Ocean/land mask at wu-, wv- and w points
	285	!----------------------------------------------
[5770]	286	wmask (:,:,1) = tmask(:,:,1) ! surface
	287	wumask(:,:,1) = umask(:,:,1)
	288	wvmask(:,:,1) = vmask(:,:,1)
	289	DO jk = 2, jpk ! interior values
	290	wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
	291	wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)
	292	wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
[5120]	293	END DO
[3]	294
	295	! 4. ocean/land mask for the elliptic equation
	296	! --------------------------------------------
[4990]	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
[5770]	379	! ! Gibraltar strait : partial slip (fmask=0.5)
	380	ij0 = 101 ; ij1 = 101
	381	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
	382	ij0 = 102 ; ij1 = 102
	383	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
	384	!
	385	! ! Bab el Mandeb : partial slip (fmask=1)
	386	ij0 = 87 ; ij1 = 88
	387	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
	388	ij0 = 88 ; ij1 = 88
	389	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
	390	!
[1707]	391	! ! Danish straits : strong slip (fmask > 2)
	392	! We keep this as an example but it is instable in this case
	393	! ij0 = 115 ; ij1 = 115
[2528]	394	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	395	! ij0 = 116 ; ij1 = 116
[2528]	396	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	397	!
	398	ENDIF
[1566]	399	!
[2528]	400	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	401	! ! Increased lateral friction near of some straits
[5506]	402	! This dirty section will be suppressed by simplification process:
	403	! all this will come back in input files
	404	! Currently these hard-wired indices relate to configuration with
	405	! extend grid (jpjglo=332)
	406	!
	407	isrow = 332 - jpjglo
	408	!
[2528]	409	IF(lwp) WRITE(numout,*)
	410	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	411	IF(lwp) WRITE(numout,*) ' Gibraltar '
[5385]	412	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait
[5552]	413	ij0 = 241 - isrow ; ij1 = 241 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[1566]	414
[2528]	415	IF(lwp) WRITE(numout,*) ' Bhosporus '
[5385]	416	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
[5552]	417	ij0 = 248 - isrow ; ij1 = 248 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	418
	419	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
[5385]	420	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
[5552]	421	ij0 = 189 - isrow ; ij1 = 190 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	422
	423	IF(lwp) WRITE(numout,*) ' Lombok '
[5385]	424	ii0 = 44 ; ii1 = 44 ! Lombok Strait
[5552]	425	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	426
	427	IF(lwp) WRITE(numout,*) ' Ombai '
[5385]	428	ii0 = 53 ; ii1 = 53 ! Ombai Strait
[5552]	429	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	430
	431	IF(lwp) WRITE(numout,*) ' Timor Passage '
[5385]	432	ii0 = 56 ; ii1 = 56 ! Timor Passage
[5552]	433	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	434
	435	IF(lwp) WRITE(numout,*) ' West Halmahera '
[5385]	436	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
[5552]	437	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	438
	439	IF(lwp) WRITE(numout,*) ' East Halmahera '
[5385]	440	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
[5552]	441	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	442	!
	443	ENDIF
	444	!
	445	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	446
[3294]	447	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	448
[1566]	449	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	450	imsk(:,:) = INT( tmask_i(:,:) )
	451	WRITE(numout,*) ' tmask_i : '
	452	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	453	& 1, jpj, 1, 1, numout)
	454	WRITE (numout,*)
	455	WRITE (numout,*) ' dommsk: tmask for each level'
	456	WRITE (numout,*) ' ----------------------------'
	457	DO jk = 1, jpk
	458	imsk(:,:) = INT( tmask(:,:,jk) )
	459
	460	WRITE(numout,*)
	461	WRITE(numout,*) ' level = ',jk
	462	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	463	& 1, jpj, 1, 1, numout)
	464	END DO
	465	WRITE(numout,*)
	466	WRITE(numout,*) ' dom_msk: vmask for each level'
	467	WRITE(numout,*) ' -----------------------------'
	468	DO jk = 1, jpk
	469	imsk(:,:) = INT( vmask(:,:,jk) )
	470	WRITE(numout,*)
	471	WRITE(numout,*) ' level = ',jk
	472	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	473	& 1, jpj, 1, 1, numout)
	474	END DO
	475	WRITE(numout,*)
	476	WRITE(numout,*) ' dom_msk: fmask for each level'
	477	WRITE(numout,*) ' -----------------------------'
	478	DO jk = 1, jpk
	479	imsk(:,:) = INT( fmask(:,:,jk) )
	480	WRITE(numout,*)
	481	WRITE(numout,*) ' level = ',jk
	482	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	483	& 1, jpj, 1, 1, numout )
	484	END DO
	485	WRITE(numout,*)
	486	WRITE(numout,*) ' dom_msk: bmask '
	487	WRITE(numout,*) ' ---------------'
	488	WRITE(numout,*)
	489	imsk(:,:) = INT( bmask(:,:) )
	490	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	491	& 1, jpj, 1, 1, numout )
[3]	492	ENDIF
[1566]	493	!
[3294]	494	CALL wrk_dealloc( jpi, jpj, imsk )
	495	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	496	!
[3294]	497	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	498	!
[3]	499	END SUBROUTINE dom_msk
	500
	501	#if defined key_noslip_accurate
	502	!!----------------------------------------------------------------------
	503	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	504	!!----------------------------------------------------------------------
	505
	506	SUBROUTINE dom_msk_nsa
	507	!!---------------------------------------------------------------------
	508	!! * ROUTINE dom_msk_nsa *
	509	!!
	510	!! ** Purpose :
	511	!!
	512	!! ** Method :
	513	!!
	514	!! ** Action :
	515	!!----------------------------------------------------------------------
[2715]	516	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	517	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	518	REAL(wp) :: zaa
[3]	519	!!---------------------------------------------------------------------
[3294]	520	!
	521	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	522	!
[2715]	523	IF(lwp) WRITE(numout,*)
	524	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	525	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
	526	IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
[3]	527
	528	! mask for second order calculation of vorticity
	529	! ----------------------------------------------
	530	! noslip boundary condition: fmask=1 at convex corner, store
	531	! index of straight coast meshes ( 'west', refering to a coast,
	532	! means west of the ocean, aso)
	533
	534	DO jk = 1, jpk
	535	DO jl = 1, 4
	536	npcoa(jl,jk) = 0
	537	DO ji = 1, 2*(jpi+jpj)
	538	nicoa(ji,jl,jk) = 0
	539	njcoa(ji,jl,jk) = 0
	540	END DO
	541	END DO
	542	END DO
	543
	544	IF( jperio == 2 ) THEN
	545	WRITE(numout,*) ' '
	546	WRITE(numout,*) ' symetric boundary conditions need special'
	547	WRITE(numout,*) ' treatment not implemented. we stop.'
	548	STOP
	549	ENDIF
	550
	551	! convex corners
	552
	553	DO jk = 1, jpkm1
	554	DO jj = 1, jpjm1
	555	DO ji = 1, jpim1
	556	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	557	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	558	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	559	END DO
	560	END DO
	561	END DO
	562
	563	! north-south straight coast
	564
	565	DO jk = 1, jpkm1
	566	inw = 0
	567	ine = 0
	568	DO jj = 2, jpjm1
	569	DO ji = 2, jpim1
	570	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	571	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	572	inw = inw + 1
	573	nicoa(inw,1,jk) = ji
	574	njcoa(inw,1,jk) = jj
	575	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	576	ENDIF
	577	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	578	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	579	ine = ine + 1
	580	nicoa(ine,2,jk) = ji
	581	njcoa(ine,2,jk) = jj
	582	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	583	ENDIF
	584	END DO
	585	END DO
	586	npcoa(1,jk) = inw
	587	npcoa(2,jk) = ine
	588	END DO
	589
	590	! west-east straight coast
	591
	592	DO jk = 1, jpkm1
	593	ins = 0
	594	inn = 0
	595	DO jj = 2, jpjm1
	596	DO ji =2, jpim1
	597	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	598	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	599	ins = ins + 1
	600	nicoa(ins,3,jk) = ji
	601	njcoa(ins,3,jk) = jj
	602	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	603	ENDIF
	604	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	605	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	606	inn = inn + 1
	607	nicoa(inn,4,jk) = ji
	608	njcoa(inn,4,jk) = jj
	609	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	610	ENDIF
	611	END DO
	612	END DO
	613	npcoa(3,jk) = ins
	614	npcoa(4,jk) = inn
	615	END DO
	616
	617	itest = 2 * ( jpi + jpj )
	618	DO jk = 1, jpk
	619	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	620	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	621
	622	WRITE(ctmp1,*) ' level jk = ',jk
	623	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	624	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	625	& npcoa(3,jk), npcoa(4,jk)
[474]	626	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	627	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	628	ENDIF
	629	END DO
	630
	631	ierror = 0
	632	iind = 0
	633	ijnd = 0
[2528]	634	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	635	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	636	DO jk = 1, jpk
	637	DO jl = 1, npcoa(1,jk)
	638	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	639	ierror = ierror+1
	640	icoord(ierror,1) = nicoa(jl,1,jk)
	641	icoord(ierror,2) = njcoa(jl,1,jk)
	642	icoord(ierror,3) = jk
	643	ENDIF
	644	END DO
	645	DO jl = 1, npcoa(2,jk)
	646	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	647	ierror = ierror + 1
	648	icoord(ierror,1) = nicoa(jl,2,jk)
	649	icoord(ierror,2) = njcoa(jl,2,jk)
	650	icoord(ierror,3) = jk
	651	ENDIF
	652	END DO
	653	DO jl = 1, npcoa(3,jk)
	654	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	655	ierror = ierror + 1
	656	icoord(ierror,1) = nicoa(jl,3,jk)
	657	icoord(ierror,2) = njcoa(jl,3,jk)
	658	icoord(ierror,3) = jk
	659	ENDIF
	660	END DO
[2528]	661	DO jl = 1, npcoa(4,jk)
[3]	662	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	663	ierror=ierror + 1
	664	icoord(ierror,1) = nicoa(jl,4,jk)
	665	icoord(ierror,2) = njcoa(jl,4,jk)
	666	icoord(ierror,3) = jk
[3]	667	ENDIF
	668	END DO
	669	END DO
	670
	671	IF( ierror > 0 ) THEN
	672	IF(lwp) WRITE(numout,*)
	673	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	674	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	675	DO jl = 1, ierror
	676	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	677	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	678	END DO
[474]	679	CALL ctl_stop( 'We stop...' )
[3]	680	ENDIF
[2715]	681	!
[3294]	682	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	683	!
[3]	684	END SUBROUTINE dom_msk_nsa
	685
	686	#else
	687	!!----------------------------------------------------------------------
	688	!! Default option : Empty routine
	689	!!----------------------------------------------------------------------
	690	SUBROUTINE dom_msk_nsa
	691	END SUBROUTINE dom_msk_nsa
	692	#endif
	693
	694	!!======================================================================
	695	END MODULE dommsk

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

Download in other formats: