New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dommsk.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5120

Last change on this file since 5120 was 5120, checked in by acc, 10 years ago
#1473. Re-organisation and optimisation of ice shelf cavity option. This commit merges changes from the dev_r5094_UKMO_ISFCLEAN branch onto the trunk. Results will change, even with ln_isfcav=F, due to a return to original definitions of the vertical metrics. All changes have been reviewed and SETTE tested.
Property svn:keywords set to `Id`
File size: 30.2 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
[4990]	186
	187	! (ISF) define barotropic mask and mask the ice shelf point
	188	ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
	189
	190	DO jk = 1, jpk
	191	DO jj = 1, jpj
	192	DO ji = 1, jpi
	193	IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp ) THEN
	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	! --------------------
[4990]	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
[4990]	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
	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,:)))
	270	END DO
	271	DO ji = 1, jpim1 ! NO vector opt.
	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,:)))
	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 )
[4990]	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
[5120]	283	! 3. Ocean/land mask at wu-, wv- and w points
	284	!----------------------------------------------
	285	wmask (:,:,1) = tmask(:,:,1) ! ????????
	286	wumask(:,:,1) = umask(:,:,1) ! ????????
	287	wvmask(:,:,1) = vmask(:,:,1) ! ????????
	288	DO jk=2,jpk
	289	wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
	290	wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)
	291	wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
	292	END DO
[3]	293
	294	! 4. ocean/land mask for the elliptic equation
	295	! --------------------------------------------
[4990]	296	bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point
[1566]	297	!
	298	! ! Boundary conditions
	299	! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
[3]	300	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
[2528]	301	bmask( 1 ,:) = 0._wp
	302	bmask(jpi,:) = 0._wp
[3]	303	ENDIF
[1566]	304	IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1
[2528]	305	bmask(:, 1 ) = 0._wp
[3]	306	ENDIF
[1566]	307	! ! north fold :
	308	IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
	309	DO ji = 1, jpi
[1528]	310	ii = ji + nimpp - 1
	311	bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
[2528]	312	bmask(ji,jpj ) = 0._wp
[1528]	313	END DO
[3]	314	ENDIF
[1566]	315	IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[2528]	316	bmask(:,jpj) = 0._wp
[3]	317	ENDIF
[1566]	318	!
	319	IF( lk_mpp ) THEN ! mpp specificities
	320	! ! bmask is set to zero on the overlap region
[2528]	321	IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp
	322	IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp
	323	IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp
	324	IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp
[1566]	325	!
	326	IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
[1528]	327	DO ji = 1, nlci
	328	ii = ji + nimpp - 1
	329	bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
[2528]	330	bmask(ji,nlcj ) = 0._wp
[1528]	331	END DO
[32]	332	ENDIF
[1566]	333	IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[1528]	334	DO ji = 1, nlci
[2528]	335	bmask(ji,nlcj ) = 0._wp
[1528]	336	END DO
[32]	337	ENDIF
[3]	338	ENDIF
	339
	340
	341	! mask for second order calculation of vorticity
	342	! ----------------------------------------------
	343	CALL dom_msk_nsa
	344
	345
	346	! Lateral boundary conditions on velocity (modify fmask)
[1566]	347	! ---------------------------------------
[3]	348	DO jk = 1, jpk
[1566]	349	zwf(:,:) = fmask(:,:,jk)
[3]	350	DO jj = 2, jpjm1
	351	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	352	IF( fmask(ji,jj,jk) == 0._wp ) THEN
	353	fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	354	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
[3]	355	ENDIF
	356	END DO
	357	END DO
	358	DO jj = 2, jpjm1
[2528]	359	IF( fmask(1,jj,jk) == 0._wp ) THEN
	360	fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]	361	ENDIF
[2528]	362	IF( fmask(jpi,jj,jk) == 0._wp ) THEN
	363	fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]	364	ENDIF
[1566]	365	END DO
[3]	366	DO ji = 2, jpim1
[2528]	367	IF( fmask(ji,1,jk) == 0._wp ) THEN
	368	fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]	369	ENDIF
[2528]	370	IF( fmask(ji,jpj,jk) == 0._wp ) THEN
	371	fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]	372	ENDIF
	373	END DO
	374	END DO
[1566]	375	!
	376	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration
	377	! ! Increased lateral friction near of some straits
[2528]	378	IF( nn_cla == 0 ) THEN
[1273]	379	! ! Gibraltar strait : partial slip (fmask=0.5)
	380	ij0 = 101 ; ij1 = 101
[2528]	381	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[32]	382	ij0 = 102 ; ij1 = 102
[2528]	383	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[1273]	384	!
	385	! ! Bab el Mandeb : partial slip (fmask=1)
	386	ij0 = 87 ; ij1 = 88
[2528]	387	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	388	ij0 = 88 ; ij1 = 88
[2528]	389	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	390	!
[3]	391	ENDIF
[1707]	392	! ! Danish straits : strong slip (fmask > 2)
	393	! We keep this as an example but it is instable in this case
	394	! ij0 = 115 ; ij1 = 115
[2528]	395	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	396	! ij0 = 116 ; ij1 = 116
[2528]	397	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	398	!
	399	ENDIF
[1566]	400	!
[2528]	401	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	402	! ! Increased lateral friction near of some straits
	403	IF(lwp) WRITE(numout,*)
	404	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	405	IF(lwp) WRITE(numout,*) ' Gibraltar '
	406	ii0 = 283 ; ii1 = 284 ! Gibraltar Strait
	407	ij0 = 200 ; ij1 = 200 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
[1566]	408
[2528]	409	IF(lwp) WRITE(numout,*) ' Bhosporus '
	410	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
	411	ij0 = 208 ; ij1 = 208 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	412
	413	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
	414	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
	415	ij0 = 149 ; ij1 = 150 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
	416
	417	IF(lwp) WRITE(numout,*) ' Lombok '
	418	ii0 = 44 ; ii1 = 44 ! Lombok Strait
	419	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	420
	421	IF(lwp) WRITE(numout,*) ' Ombai '
	422	ii0 = 53 ; ii1 = 53 ! Ombai Strait
	423	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	424
	425	IF(lwp) WRITE(numout,*) ' Timor Passage '
	426	ii0 = 56 ; ii1 = 56 ! Timor Passage
	427	ij0 = 124 ; ij1 = 125 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp
	428
	429	IF(lwp) WRITE(numout,*) ' West Halmahera '
	430	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
	431	ij0 = 141 ; ij1 = 142 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
	432
	433	IF(lwp) WRITE(numout,*) ' East Halmahera '
	434	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
	435	ij0 = 141 ; ij1 = 142 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp
	436	!
	437	ENDIF
	438	!
	439	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	440
[3294]	441	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	442
[1566]	443	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	444	imsk(:,:) = INT( tmask_i(:,:) )
	445	WRITE(numout,*) ' tmask_i : '
	446	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	447	& 1, jpj, 1, 1, numout)
	448	WRITE (numout,*)
	449	WRITE (numout,*) ' dommsk: tmask for each level'
	450	WRITE (numout,*) ' ----------------------------'
	451	DO jk = 1, jpk
	452	imsk(:,:) = INT( tmask(:,:,jk) )
	453
	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: vmask for each level'
	461	WRITE(numout,*) ' -----------------------------'
	462	DO jk = 1, jpk
	463	imsk(:,:) = INT( vmask(:,:,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: fmask for each level'
	471	WRITE(numout,*) ' -----------------------------'
	472	DO jk = 1, jpk
	473	imsk(:,:) = INT( fmask(:,:,jk) )
	474	WRITE(numout,*)
	475	WRITE(numout,*) ' level = ',jk
	476	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	477	& 1, jpj, 1, 1, numout )
	478	END DO
	479	WRITE(numout,*)
	480	WRITE(numout,*) ' dom_msk: bmask '
	481	WRITE(numout,*) ' ---------------'
	482	WRITE(numout,*)
	483	imsk(:,:) = INT( bmask(:,:) )
	484	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	485	& 1, jpj, 1, 1, numout )
[3]	486	ENDIF
[1566]	487	!
[3294]	488	CALL wrk_dealloc( jpi, jpj, imsk )
	489	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	490	!
[3294]	491	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	492	!
[3]	493	END SUBROUTINE dom_msk
	494
	495	#if defined key_noslip_accurate
	496	!!----------------------------------------------------------------------
	497	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	498	!!----------------------------------------------------------------------
	499
	500	SUBROUTINE dom_msk_nsa
	501	!!---------------------------------------------------------------------
	502	!! * ROUTINE dom_msk_nsa *
	503	!!
	504	!! ** Purpose :
	505	!!
	506	!! ** Method :
	507	!!
	508	!! ** Action :
	509	!!----------------------------------------------------------------------
[2715]	510	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	511	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	512	REAL(wp) :: zaa
[3]	513	!!---------------------------------------------------------------------
[3294]	514	!
	515	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	516	!
[2715]	517	IF(lwp) WRITE(numout,*)
	518	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	519	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
	520	IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
[3]	521
	522	! mask for second order calculation of vorticity
	523	! ----------------------------------------------
	524	! noslip boundary condition: fmask=1 at convex corner, store
	525	! index of straight coast meshes ( 'west', refering to a coast,
	526	! means west of the ocean, aso)
	527
	528	DO jk = 1, jpk
	529	DO jl = 1, 4
	530	npcoa(jl,jk) = 0
	531	DO ji = 1, 2*(jpi+jpj)
	532	nicoa(ji,jl,jk) = 0
	533	njcoa(ji,jl,jk) = 0
	534	END DO
	535	END DO
	536	END DO
	537
	538	IF( jperio == 2 ) THEN
	539	WRITE(numout,*) ' '
	540	WRITE(numout,*) ' symetric boundary conditions need special'
	541	WRITE(numout,*) ' treatment not implemented. we stop.'
	542	STOP
	543	ENDIF
	544
	545	! convex corners
	546
	547	DO jk = 1, jpkm1
	548	DO jj = 1, jpjm1
	549	DO ji = 1, jpim1
	550	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	551	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	552	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	553	END DO
	554	END DO
	555	END DO
	556
	557	! north-south straight coast
	558
	559	DO jk = 1, jpkm1
	560	inw = 0
	561	ine = 0
	562	DO jj = 2, jpjm1
	563	DO ji = 2, jpim1
	564	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	565	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	566	inw = inw + 1
	567	nicoa(inw,1,jk) = ji
	568	njcoa(inw,1,jk) = jj
	569	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	570	ENDIF
	571	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	572	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	573	ine = ine + 1
	574	nicoa(ine,2,jk) = ji
	575	njcoa(ine,2,jk) = jj
	576	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	577	ENDIF
	578	END DO
	579	END DO
	580	npcoa(1,jk) = inw
	581	npcoa(2,jk) = ine
	582	END DO
	583
	584	! west-east straight coast
	585
	586	DO jk = 1, jpkm1
	587	ins = 0
	588	inn = 0
	589	DO jj = 2, jpjm1
	590	DO ji =2, jpim1
	591	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	592	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	593	ins = ins + 1
	594	nicoa(ins,3,jk) = ji
	595	njcoa(ins,3,jk) = jj
	596	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	597	ENDIF
	598	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	599	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	600	inn = inn + 1
	601	nicoa(inn,4,jk) = ji
	602	njcoa(inn,4,jk) = jj
	603	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	604	ENDIF
	605	END DO
	606	END DO
	607	npcoa(3,jk) = ins
	608	npcoa(4,jk) = inn
	609	END DO
	610
	611	itest = 2 * ( jpi + jpj )
	612	DO jk = 1, jpk
	613	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	614	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	615
	616	WRITE(ctmp1,*) ' level jk = ',jk
	617	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	618	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	619	& npcoa(3,jk), npcoa(4,jk)
[474]	620	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	621	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	622	ENDIF
	623	END DO
	624
	625	ierror = 0
	626	iind = 0
	627	ijnd = 0
[2528]	628	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	629	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	630	DO jk = 1, jpk
	631	DO jl = 1, npcoa(1,jk)
	632	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	633	ierror = ierror+1
	634	icoord(ierror,1) = nicoa(jl,1,jk)
	635	icoord(ierror,2) = njcoa(jl,1,jk)
	636	icoord(ierror,3) = jk
	637	ENDIF
	638	END DO
	639	DO jl = 1, npcoa(2,jk)
	640	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	641	ierror = ierror + 1
	642	icoord(ierror,1) = nicoa(jl,2,jk)
	643	icoord(ierror,2) = njcoa(jl,2,jk)
	644	icoord(ierror,3) = jk
	645	ENDIF
	646	END DO
	647	DO jl = 1, npcoa(3,jk)
	648	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	649	ierror = ierror + 1
	650	icoord(ierror,1) = nicoa(jl,3,jk)
	651	icoord(ierror,2) = njcoa(jl,3,jk)
	652	icoord(ierror,3) = jk
	653	ENDIF
	654	END DO
[2528]	655	DO jl = 1, npcoa(4,jk)
[3]	656	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	657	ierror=ierror + 1
	658	icoord(ierror,1) = nicoa(jl,4,jk)
	659	icoord(ierror,2) = njcoa(jl,4,jk)
	660	icoord(ierror,3) = jk
[3]	661	ENDIF
	662	END DO
	663	END DO
	664
	665	IF( ierror > 0 ) THEN
	666	IF(lwp) WRITE(numout,*)
	667	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	668	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	669	DO jl = 1, ierror
	670	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	671	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	672	END DO
[474]	673	CALL ctl_stop( 'We stop...' )
[3]	674	ENDIF
[2715]	675	!
[3294]	676	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	677	!
[3]	678	END SUBROUTINE dom_msk_nsa
	679
	680	#else
	681	!!----------------------------------------------------------------------
	682	!! Default option : Empty routine
	683	!!----------------------------------------------------------------------
	684	SUBROUTINE dom_msk_nsa
	685	END SUBROUTINE dom_msk_nsa
	686	#endif
	687
	688	!!======================================================================
	689	END MODULE dommsk

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

Download in other formats: