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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 10302

Last change on this file since 10302 was 10302, checked in by dford, 5 years ago
Merge in revisions 8447:10159 of dev_r5518_GO6_package.
File size: 33.3 KB

Rev	Line
[10302]	1
[3]	2	MODULE dommsk
[1566]	3	!!======================================================================
[3]	4	!! * MODULE dommsk *
	5	!! Ocean initialization : domain land/sea mask
[1566]	6	!!======================================================================
	7	!! History : OPA ! 1987-07 (G. Madec) Original code
[2528]	8	!! 6.0 ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon)
	9	!! 7.0 ! 1996-01 (G. Madec) suppression of common work arrays
[1566]	10	!! - ! 1996-05 (G. Madec) mask computed from tmask and sup-
	11	!! ! pression of the double computation of bmask
[2528]	12	!! 8.0 ! 1997-02 (G. Madec) mesh information put in domhgr.F
	13	!! 8.1 ! 1997-07 (G. Madec) modification of mbathy and fmask
[1566]	14	!! - ! 1998-05 (G. Roullet) free surface
[2528]	15	!! 8.2 ! 2000-03 (G. Madec) no slip accurate
[1566]	16	!! - ! 2001-09 (J.-M. Molines) Open boundaries
	17	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	18	!! - ! 2005-11 (V. Garnier) Surface pressure gradient organization
	19	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
	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
[10302]	33	USE domwri
[3294]	34	USE timing ! Timing
[3]	35
	36	IMPLICIT NONE
	37	PRIVATE
	38
[2715]	39	PUBLIC dom_msk ! routine called by inidom.F90
	40	PUBLIC dom_msk_alloc ! routine called by nemogcm.F90
[3]	41
[1601]	42	! !!* Namelist namlbc : lateral boundary condition *
[4147]	43	REAL(wp) :: rn_shlat ! type of lateral boundary condition on velocity
	44	LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition
[3294]	45	! with analytical eqs.
[2715]	46
[3294]	47
[2715]	48	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icoord ! Workspace for dom_msk_nsa()
	49
[3]	50	!! * Substitutions
	51	# include "vectopt_loop_substitute.h90"
[1566]	52	!!----------------------------------------------------------------------
	53	!! NEMO/OPA 3.2 , LODYC-IPSL (2009)
[6486]	54	!! $Id$
[2528]	55	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1566]	56	!!----------------------------------------------------------------------
[3]	57	CONTAINS
	58
[2715]	59	INTEGER FUNCTION dom_msk_alloc()
	60	!!---------------------------------------------------------------------
	61	!! * FUNCTION dom_msk_alloc *
	62	!!---------------------------------------------------------------------
	63	dom_msk_alloc = 0
	64	#if defined key_noslip_accurate
	65	ALLOCATE(icoord(jpijpjjpk,3), STAT=dom_msk_alloc)
	66	#endif
	67	IF( dom_msk_alloc /= 0 ) CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
	68	!
	69	END FUNCTION dom_msk_alloc
	70
	71
[3]	72	SUBROUTINE dom_msk
	73	!!---------------------------------------------------------------------
	74	!! * ROUTINE dom_msk *
	75	!!
	76	!! ** Purpose : Compute land/ocean mask arrays at tracer points, hori-
	77	!! zontal velocity points (u & v), vorticity points (f) and baro-
	78	!! tropic stream function points (b).
	79	!!
	80	!! ** Method : The ocean/land mask is computed from the basin bathy-
	81	!! metry in level (mbathy) which is defined or read in dommba.
[1528]	82	!! mbathy equals 0 over continental T-point
	83	!! and the number of ocean level over the ocean.
[3]	84	!!
	85	!! At a given position (ji,jj,jk) the ocean/land mask is given by:
	86	!! t-point : 0. IF mbathy( ji ,jj) =< 0
	87	!! 1. IF mbathy( ji ,jj) >= jk
	88	!! u-point : 0. IF mbathy( ji ,jj) or mbathy(ji+1, jj ) =< 0
	89	!! 1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
	90	!! v-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) =< 0
	91	!! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
	92	!! f-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1)
	93	!! or mbathy(ji+1,jj) or mbathy(ji+1,jj+1) =< 0
	94	!! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
[2528]	95	!! and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
[3]	96	!! b-point : the same definition as for f-point of the first ocean
	97	!! level (surface level) but with 0 along coastlines.
[2528]	98	!! tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
	99	!! rows/lines due to cyclic or North Fold boundaries as well
	100	!! as MPP halos.
[3]	101	!!
	102	!! The lateral friction is set through the value of fmask along
[1601]	103	!! the coast and topography. This value is defined by rn_shlat, a
[3]	104	!! namelist parameter:
[1601]	105	!! rn_shlat = 0, free slip (no shear along the coast)
	106	!! rn_shlat = 2, no slip (specified zero velocity at the coast)
	107	!! 0 < rn_shlat < 2, partial slip \| non-linear velocity profile
	108	!! 2 < rn_shlat, strong slip \| in the lateral boundary layer
[3]	109	!!
	110	!! N.B. If nperio not equal to 0, the land/ocean mask arrays
	111	!! are defined with the proper value at lateral domain boundaries,
	112	!! but bmask. indeed, bmask defined the domain over which the
	113	!! barotropic stream function is computed. this domain cannot
	114	!! contain identical columns because the matrix associated with
	115	!! the barotropic stream function equation is then no more inverti-
	116	!! ble. therefore bmask is set to 0 along lateral domain boundaries
	117	!! even IF nperio is not zero.
	118	!!
[4328]	119	!! In case of open boundaries (lk_bdy=T):
[3]	120	!! - tmask is set to 1 on the points to be computed bay the open
	121	!! boundaries routines.
	122	!! - bmask is set to 0 on the open boundaries.
	123	!!
[1566]	124	!! ** Action : tmask : land/ocean mask at t-point (=0. or 1.)
	125	!! umask : land/ocean mask at u-point (=0. or 1.)
	126	!! vmask : land/ocean mask at v-point (=0. or 1.)
	127	!! fmask : land/ocean mask at f-point (=0. or 1.)
[1601]	128	!! =rn_shlat along lateral boundaries
[1566]	129	!! bmask : land/ocean mask at barotropic stream
	130	!! function point (=0. or 1.) and set to 0 along lateral boundaries
[2528]	131	!! tmask_i : interior ocean mask
[3]	132	!!----------------------------------------------------------------------
[2715]	133	!
[454]	134	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	135	INTEGER :: iif, iil, ii0, ii1, ii ! local integers
	136	INTEGER :: ijf, ijl, ij0, ij1 ! - -
[4147]	137	INTEGER :: ios
[5385]	138	INTEGER :: isrow ! index for ORCA1 starting row
[3294]	139	INTEGER , POINTER, DIMENSION(:,:) :: imsk
[6491]	140	REAL(wp) :: zphi_drake_passage, zshlat_antarc
[3294]	141	REAL(wp), POINTER, DIMENSION(:,:) :: zwf
[10302]	142	REAL(wp) :: uvt(jpi,jpj) ! dummy array for masking purposes.
[1601]	143	!!
[3294]	144	NAMELIST/namlbc/ rn_shlat, ln_vorlat
[3]	145	!!---------------------------------------------------------------------
[3294]	146	!
	147	IF( nn_timing == 1 ) CALL timing_start('dom_msk')
	148	!
	149	CALL wrk_alloc( jpi, jpj, imsk )
	150	CALL wrk_alloc( jpi, jpj, zwf )
	151	!
[4147]	152	REWIND( numnam_ref ) ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
	153	READ ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
	154	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
	155
	156	REWIND( numnam_cfg ) ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
	157	READ ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
	158	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
[4624]	159	IF(lwm) WRITE ( numond, namlbc )
[1566]	160
	161	IF(lwp) THEN ! control print
[3]	162	WRITE(numout,*)
	163	WRITE(numout,*) 'dommsk : ocean mask '
	164	WRITE(numout,*) '~~~~~~'
[1566]	165	WRITE(numout,*) ' Namelist namlbc'
[3294]	166	WRITE(numout,*) ' lateral momentum boundary cond. rn_shlat = ',rn_shlat
	167	WRITE(numout,*) ' consistency with analytical form ln_vorlat = ',ln_vorlat
[3]	168	ENDIF
	169
[2528]	170	IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral free-slip '
[1601]	171	ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral no-slip '
	172	ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip '
	173	ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip '
	174	ELSE
	175	WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
	176	CALL ctl_stop( ctmp1 )
[3]	177	ENDIF
	178
	179	! 1. Ocean/land mask at t-point (computed from mbathy)
	180	! -----------------------------
[1566]	181	! N.B. tmask has already the right boundary conditions since mbathy is ok
	182	!
[2528]	183	tmask(:,:,:) = 0._wp
[3]	184	DO jk = 1, jpk
	185	DO jj = 1, jpj
	186	DO ji = 1, jpi
[2528]	187	IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) tmask(ji,jj,jk) = 1._wp
[3]	188	END DO
	189	END DO
	190	END DO
[4990]	191
	192	! (ISF) define barotropic mask and mask the ice shelf point
	193	ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
	194
	195	DO jk = 1, jpk
	196	DO jj = 1, jpj
	197	DO ji = 1, jpi
	198	IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp ) THEN
	199	tmask(ji,jj,jk) = 0._wp
	200	END IF
	201	END DO
	202	END DO
	203	END DO
[3]	204
[1566]	205	!!gm ????
[255]	206	#if defined key_zdfkpp
[1601]	207	IF( cp_cfg == 'orca' ) THEN
	208	IF( jp_cfg == 2 ) THEN ! land point on Bab el Mandeb zonal section
[291]	209	ij0 = 87 ; ij1 = 88
	210	ii0 = 160 ; ii1 = 161
[2528]	211	tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
[291]	212	ELSE
	213	IF(lwp) WRITE(numout,*)
	214	IF(lwp) WRITE(numout,cform_war)
	215	IF(lwp) WRITE(numout,*)
	216	IF(lwp) WRITE(numout,*)' A mask must be applied on Bab el Mandeb strait'
	217	IF(lwp) WRITE(numout,*)' in case of ORCAs configurations'
	218	IF(lwp) WRITE(numout,*)' This is a problem which is not yet solved'
	219	IF(lwp) WRITE(numout,*)
	220	ENDIF
	221	ENDIF
[255]	222	#endif
[1566]	223	!!gm end
[291]	224
[3]	225	! Interior domain mask (used for global sum)
	226	! --------------------
[4990]	227	tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf
[10302]	228
[3]	229	iif = jpreci ! ???
	230	iil = nlci - jpreci + 1
	231	ijf = jprecj ! ???
	232	ijl = nlcj - jprecj + 1
	233
[2528]	234	tmask_i( 1 :iif, : ) = 0._wp ! first columns
	235	tmask_i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)
	236	tmask_i( : , 1 :ijf) = 0._wp ! first rows
	237	tmask_i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)
[3]	238
	239	! north fold mask
[1566]	240	! ---------------
[2528]	241	tpol(1:jpiglo) = 1._wp
	242	fpol(1:jpiglo) = 1._wp
[3]	243	IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot
[2528]	244	tpol(jpiglo/2+1:jpiglo) = 0._wp
	245	fpol( 1 :jpiglo) = 0._wp
[1566]	246	IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row
[291]	247	DO ji = iif+1, iil-1
	248	tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
	249	END DO
	250	ENDIF
[3]	251	ENDIF
[10302]	252
	253
[3]	254	IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot
[2528]	255	tpol( 1 :jpiglo) = 0._wp
	256	fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]	257	ENDIF
	258
	259	! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask)
	260	! -------------------------------------------
	261	DO jk = 1, jpk
	262	DO jj = 1, jpjm1
	263	DO ji = 1, fs_jpim1 ! vector loop
	264	umask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk)
	265	vmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji ,jj+1,jk)
[1271]	266	END DO
[1694]	267	DO ji = 1, jpim1 ! NO vector opt.
[3]	268	fmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk) &
[62]	269	& * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]	270	END DO
[10302]	271	END DO
[3]	272	END DO
[4990]	273	! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
	274	DO jj = 1, jpjm1
	275	DO ji = 1, fs_jpim1 ! vector loop
	276	umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))
	277	vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))
	278	END DO
	279	DO ji = 1, jpim1 ! NO vector opt.
	280	fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &
	281	& * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
	282	END DO
	283	END DO
[2528]	284	CALL lbc_lnk( umask, 'U', 1._wp ) ! Lateral boundary conditions
	285	CALL lbc_lnk( vmask, 'V', 1._wp )
	286	CALL lbc_lnk( fmask, 'F', 1._wp )
[4990]	287	CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions
	288	CALL lbc_lnk( vmask_i, 'V', 1._wp )
	289	CALL lbc_lnk( fmask_i, 'F', 1._wp )
[3]	290
[10302]	291
	292	! Set up mask for diagnostics on T points, to exclude duplicate
	293	! data points in wrap and N-fold regions.
	294	CALL dom_uniq( uvt, 'T' )
	295	DO jk = 1, jpk
	296	tmask_i_diag(:,:,jk) = tmask(:,:,jk) * uvt(:,:)
	297	END DO
	298
	299	! Set up mask for diagnostics on U points, to exclude duplicate
	300	! data points in wrap and N-fold regions.
	301	umask_i_diag(:,:,:) = 1.0
	302	umask_i_diag(2:jpim1,:,:) = tmask(2:jpim1,:,:) + tmask(3:jpi,:,:)
	303	CALL lbc_lnk( umask_i_diag, 'U', 1. )
	304
	305	! Now mask out any duplicate points
	306	CALL dom_uniq( uvt, 'U' )
	307	DO jk = 1, jpk
	308	umask_i_diag(:,:,jk) = umask_i_diag(:,:,jk) * uvt(:,:)
	309	END DO
	310
	311
	312	! Set up mask for diagnostics on V points, to exclude duplicate
	313	! data points in wrap and N-fold regions.
	314	vmask_i_diag(:,:,:) = 1.0
	315	vmask_i_diag(:,2:jpjm1,:) = tmask(:,2:jpjm1,:) + tmask(:,3:jpj,:)
	316	CALL lbc_lnk( vmask_i_diag, 'V', 1. )
	317
	318	! Now mask out any duplicate points
	319	CALL dom_uniq( uvt, 'V' )
	320	DO jk = 1, jpk
	321	vmask_i_diag(:,:,jk) = vmask_i_diag(:,:,jk) * uvt(:,:)
	322	END DO
	323
	324
	325
[5120]	326	! 3. Ocean/land mask at wu-, wv- and w points
	327	!----------------------------------------------
	328	wmask (:,:,1) = tmask(:,:,1) ! ????????
	329	wumask(:,:,1) = umask(:,:,1) ! ????????
	330	wvmask(:,:,1) = vmask(:,:,1) ! ????????
	331	DO jk=2,jpk
	332	wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
	333	wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)
	334	wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
	335	END DO
[3]	336
	337	! 4. ocean/land mask for the elliptic equation
	338	! --------------------------------------------
[4990]	339	bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point
[1566]	340	!
	341	! ! Boundary conditions
	342	! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
[3]	343	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
[2528]	344	bmask( 1 ,:) = 0._wp
	345	bmask(jpi,:) = 0._wp
[3]	346	ENDIF
[1566]	347	IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1
[2528]	348	bmask(:, 1 ) = 0._wp
[3]	349	ENDIF
[1566]	350	! ! north fold :
	351	IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
	352	DO ji = 1, jpi
[1528]	353	ii = ji + nimpp - 1
	354	bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
[2528]	355	bmask(ji,jpj ) = 0._wp
[1528]	356	END DO
[3]	357	ENDIF
[1566]	358	IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[2528]	359	bmask(:,jpj) = 0._wp
[3]	360	ENDIF
[1566]	361	!
	362	IF( lk_mpp ) THEN ! mpp specificities
	363	! ! bmask is set to zero on the overlap region
[2528]	364	IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp
	365	IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp
	366	IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp
	367	IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp
[1566]	368	!
	369	IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
[1528]	370	DO ji = 1, nlci
	371	ii = ji + nimpp - 1
	372	bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
[2528]	373	bmask(ji,nlcj ) = 0._wp
[1528]	374	END DO
[32]	375	ENDIF
[1566]	376	IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[1528]	377	DO ji = 1, nlci
[2528]	378	bmask(ji,nlcj ) = 0._wp
[1528]	379	END DO
[32]	380	ENDIF
[3]	381	ENDIF
	382
	383
	384	! mask for second order calculation of vorticity
	385	! ----------------------------------------------
	386	CALL dom_msk_nsa
	387
	388
	389	! Lateral boundary conditions on velocity (modify fmask)
[1566]	390	! ---------------------------------------
[3]	391	DO jk = 1, jpk
[1566]	392	zwf(:,:) = fmask(:,:,jk)
[3]	393	DO jj = 2, jpjm1
	394	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	395	IF( fmask(ji,jj,jk) == 0._wp ) THEN
	396	fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	397	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
[3]	398	ENDIF
	399	END DO
	400	END DO
	401	DO jj = 2, jpjm1
[2528]	402	IF( fmask(1,jj,jk) == 0._wp ) THEN
	403	fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]	404	ENDIF
[2528]	405	IF( fmask(jpi,jj,jk) == 0._wp ) THEN
	406	fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]	407	ENDIF
[1566]	408	END DO
[3]	409	DO ji = 2, jpim1
[2528]	410	IF( fmask(ji,1,jk) == 0._wp ) THEN
	411	fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]	412	ENDIF
[2528]	413	IF( fmask(ji,jpj,jk) == 0._wp ) THEN
	414	fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]	415	ENDIF
	416	END DO
	417	END DO
[1566]	418	!
	419	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration
	420	! ! Increased lateral friction near of some straits
[2528]	421	IF( nn_cla == 0 ) THEN
[1273]	422	! ! Gibraltar strait : partial slip (fmask=0.5)
	423	ij0 = 101 ; ij1 = 101
[2528]	424	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[32]	425	ij0 = 102 ; ij1 = 102
[2528]	426	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[1273]	427	!
	428	! ! Bab el Mandeb : partial slip (fmask=1)
	429	ij0 = 87 ; ij1 = 88
[2528]	430	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	431	ij0 = 88 ; ij1 = 88
[2528]	432	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	433	!
[3]	434	ENDIF
[1707]	435	! ! Danish straits : strong slip (fmask > 2)
	436	! We keep this as an example but it is instable in this case
	437	! ij0 = 115 ; ij1 = 115
[2528]	438	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	439	! ij0 = 116 ; ij1 = 116
[2528]	440	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	441	!
	442	ENDIF
[1566]	443	!
[2528]	444	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	445	! ! Increased lateral friction near of some straits
[5506]	446	! This dirty section will be suppressed by simplification process:
	447	! all this will come back in input files
	448	! Currently these hard-wired indices relate to configuration with
	449	! extend grid (jpjglo=332)
	450	!
	451	isrow = 332 - jpjglo
	452	!
[2528]	453	IF(lwp) WRITE(numout,*)
	454	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	455	IF(lwp) WRITE(numout,*) ' Gibraltar '
[5385]	456	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait
[6487]	457	ij0 = 241 - isrow ; ij1 = 241 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[1566]	458
[2528]	459	IF(lwp) WRITE(numout,*) ' Bhosporus '
[5385]	460	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
[6487]	461	ij0 = 248 - isrow ; ij1 = 248 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	462
	463	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
[5385]	464	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
[6487]	465	ij0 = 189 - isrow ; ij1 = 190 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	466
	467	IF(lwp) WRITE(numout,*) ' Lombok '
[5385]	468	ii0 = 44 ; ii1 = 44 ! Lombok Strait
[6487]	469	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	470
	471	IF(lwp) WRITE(numout,*) ' Ombai '
[5385]	472	ii0 = 53 ; ii1 = 53 ! Ombai Strait
[6487]	473	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	474
	475	IF(lwp) WRITE(numout,*) ' Timor Passage '
[5385]	476	ii0 = 56 ; ii1 = 56 ! Timor Passage
[6487]	477	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	478
	479	IF(lwp) WRITE(numout,*) ' West Halmahera '
[5385]	480	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
[6487]	481	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	482
	483	IF(lwp) WRITE(numout,*) ' East Halmahera '
[5385]	484	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
[6487]	485	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	486	!
	487	ENDIF
	488	!
[6491]	489	IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN
	490	! ! ORCA_R025 configuration
	491	! ! Increased lateral friction on parts of Antarctic coastline
	492	! ! for increased stability
	493	! ! NB. This only works to do this here if we have free slip
	494	! ! generally, so fmask is zero at coast points.
	495	IF(lwp) WRITE(numout,*)
	496	IF(lwp) WRITE(numout,*) ' orca_r025: increase friction in following regions : '
	497	IF(lwp) WRITE(numout,*) ' whole Antarctic coastline: partial slip shlat=1 '
	498
	499	zphi_drake_passage = -58.0_wp
	500	zshlat_antarc = 1.0_wp
	501	zwf(:,:) = fmask(:,:,1)
	502	DO jj = 2, jpjm1
	503	DO ji = fs_2, fs_jpim1 ! vector opt.
	504	IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN
	505	fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	506	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
	507	ENDIF
	508	END DO
	509	END DO
	510	END IF
	511	!
[2528]	512	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	513
[3294]	514	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	515
[1566]	516	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	517	imsk(:,:) = INT( tmask_i(:,:) )
	518	WRITE(numout,*) ' tmask_i : '
	519	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	520	& 1, jpj, 1, 1, numout)
	521	WRITE (numout,*)
	522	WRITE (numout,*) ' dommsk: tmask for each level'
	523	WRITE (numout,*) ' ----------------------------'
	524	DO jk = 1, jpk
	525	imsk(:,:) = INT( tmask(:,:,jk) )
	526
	527	WRITE(numout,*)
	528	WRITE(numout,*) ' level = ',jk
	529	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	530	& 1, jpj, 1, 1, numout)
	531	END DO
	532	WRITE(numout,*)
	533	WRITE(numout,*) ' dom_msk: vmask for each level'
	534	WRITE(numout,*) ' -----------------------------'
	535	DO jk = 1, jpk
	536	imsk(:,:) = INT( vmask(:,:,jk) )
	537	WRITE(numout,*)
	538	WRITE(numout,*) ' level = ',jk
	539	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	540	& 1, jpj, 1, 1, numout)
	541	END DO
	542	WRITE(numout,*)
	543	WRITE(numout,*) ' dom_msk: fmask for each level'
	544	WRITE(numout,*) ' -----------------------------'
	545	DO jk = 1, jpk
	546	imsk(:,:) = INT( fmask(:,:,jk) )
	547	WRITE(numout,*)
	548	WRITE(numout,*) ' level = ',jk
	549	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	550	& 1, jpj, 1, 1, numout )
	551	END DO
	552	WRITE(numout,*)
	553	WRITE(numout,*) ' dom_msk: bmask '
	554	WRITE(numout,*) ' ---------------'
	555	WRITE(numout,*)
	556	imsk(:,:) = INT( bmask(:,:) )
	557	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	558	& 1, jpj, 1, 1, numout )
[3]	559	ENDIF
[1566]	560	!
[3294]	561	CALL wrk_dealloc( jpi, jpj, imsk )
	562	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	563	!
[3294]	564	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	565	!
[3]	566	END SUBROUTINE dom_msk
	567
	568	#if defined key_noslip_accurate
	569	!!----------------------------------------------------------------------
	570	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	571	!!----------------------------------------------------------------------
	572
	573	SUBROUTINE dom_msk_nsa
	574	!!---------------------------------------------------------------------
	575	!! * ROUTINE dom_msk_nsa *
	576	!!
	577	!! ** Purpose :
	578	!!
	579	!! ** Method :
	580	!!
	581	!! ** Action :
	582	!!----------------------------------------------------------------------
[2715]	583	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	584	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	585	REAL(wp) :: zaa
[3]	586	!!---------------------------------------------------------------------
[3294]	587	!
	588	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	589	!
[2715]	590	IF(lwp) WRITE(numout,*)
	591	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	592	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
[8280]	593	IF( lk_mpp ) CALL ctl_stop('STOP', ' mpp version is not yet implemented' )
[3]	594
	595	! mask for second order calculation of vorticity
	596	! ----------------------------------------------
	597	! noslip boundary condition: fmask=1 at convex corner, store
	598	! index of straight coast meshes ( 'west', refering to a coast,
	599	! means west of the ocean, aso)
	600
	601	DO jk = 1, jpk
	602	DO jl = 1, 4
	603	npcoa(jl,jk) = 0
	604	DO ji = 1, 2*(jpi+jpj)
	605	nicoa(ji,jl,jk) = 0
	606	njcoa(ji,jl,jk) = 0
	607	END DO
	608	END DO
	609	END DO
	610
	611	IF( jperio == 2 ) THEN
	612	WRITE(numout,*) ' '
	613	WRITE(numout,*) ' symetric boundary conditions need special'
	614	WRITE(numout,*) ' treatment not implemented. we stop.'
[8280]	615	CALL ctl_stop('STOP', 'NEMO abort from dom_msk_nsa')
[3]	616	ENDIF
	617
	618	! convex corners
	619
	620	DO jk = 1, jpkm1
	621	DO jj = 1, jpjm1
	622	DO ji = 1, jpim1
	623	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	624	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	625	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	626	END DO
	627	END DO
	628	END DO
	629
	630	! north-south straight coast
	631
	632	DO jk = 1, jpkm1
	633	inw = 0
	634	ine = 0
	635	DO jj = 2, jpjm1
	636	DO ji = 2, jpim1
	637	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	638	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	639	inw = inw + 1
	640	nicoa(inw,1,jk) = ji
	641	njcoa(inw,1,jk) = jj
	642	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	643	ENDIF
	644	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	645	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	646	ine = ine + 1
	647	nicoa(ine,2,jk) = ji
	648	njcoa(ine,2,jk) = jj
	649	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	650	ENDIF
	651	END DO
	652	END DO
	653	npcoa(1,jk) = inw
	654	npcoa(2,jk) = ine
	655	END DO
	656
	657	! west-east straight coast
	658
	659	DO jk = 1, jpkm1
	660	ins = 0
	661	inn = 0
	662	DO jj = 2, jpjm1
	663	DO ji =2, jpim1
	664	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	665	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	666	ins = ins + 1
	667	nicoa(ins,3,jk) = ji
	668	njcoa(ins,3,jk) = jj
	669	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	670	ENDIF
	671	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	672	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	673	inn = inn + 1
	674	nicoa(inn,4,jk) = ji
	675	njcoa(inn,4,jk) = jj
	676	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	677	ENDIF
	678	END DO
	679	END DO
	680	npcoa(3,jk) = ins
	681	npcoa(4,jk) = inn
	682	END DO
	683
	684	itest = 2 * ( jpi + jpj )
	685	DO jk = 1, jpk
	686	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	687	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	688
	689	WRITE(ctmp1,*) ' level jk = ',jk
	690	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	691	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	692	& npcoa(3,jk), npcoa(4,jk)
[474]	693	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	694	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	695	ENDIF
	696	END DO
	697
	698	ierror = 0
	699	iind = 0
	700	ijnd = 0
[2528]	701	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	702	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	703	DO jk = 1, jpk
	704	DO jl = 1, npcoa(1,jk)
	705	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	706	ierror = ierror+1
	707	icoord(ierror,1) = nicoa(jl,1,jk)
	708	icoord(ierror,2) = njcoa(jl,1,jk)
	709	icoord(ierror,3) = jk
	710	ENDIF
	711	END DO
	712	DO jl = 1, npcoa(2,jk)
	713	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	714	ierror = ierror + 1
	715	icoord(ierror,1) = nicoa(jl,2,jk)
	716	icoord(ierror,2) = njcoa(jl,2,jk)
	717	icoord(ierror,3) = jk
	718	ENDIF
	719	END DO
	720	DO jl = 1, npcoa(3,jk)
	721	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	722	ierror = ierror + 1
	723	icoord(ierror,1) = nicoa(jl,3,jk)
	724	icoord(ierror,2) = njcoa(jl,3,jk)
	725	icoord(ierror,3) = jk
	726	ENDIF
	727	END DO
[2528]	728	DO jl = 1, npcoa(4,jk)
[3]	729	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	730	ierror=ierror + 1
	731	icoord(ierror,1) = nicoa(jl,4,jk)
	732	icoord(ierror,2) = njcoa(jl,4,jk)
	733	icoord(ierror,3) = jk
[3]	734	ENDIF
	735	END DO
	736	END DO
	737
	738	IF( ierror > 0 ) THEN
	739	IF(lwp) WRITE(numout,*)
	740	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	741	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	742	DO jl = 1, ierror
	743	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	744	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	745	END DO
[474]	746	CALL ctl_stop( 'We stop...' )
[3]	747	ENDIF
[2715]	748	!
[3294]	749	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	750	!
[3]	751	END SUBROUTINE dom_msk_nsa
	752
	753	#else
	754	!!----------------------------------------------------------------------
	755	!! Default option : Empty routine
	756	!!----------------------------------------------------------------------
	757	SUBROUTINE dom_msk_nsa
	758	END SUBROUTINE dom_msk_nsa
	759	#endif
	760
	761	!!======================================================================
	762	END MODULE dommsk

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

Download in other formats: