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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_diag_bitcomp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 9608

Last change on this file since 9608 was 9608, checked in by frrh, 6 years ago
Introduce and apply a 3D interior mask for T, U, V and W grid diagnostics.
File size: 32.7 KB

Rev	Line
[3]	1	MODULE dommsk
[1566]	2	!!======================================================================
[3]	3	!! * MODULE dommsk *
	4	!! Ocean initialization : domain land/sea mask
[1566]	5	!!======================================================================
	6	!! History : OPA ! 1987-07 (G. Madec) Original code
[2528]	7	!! 6.0 ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon)
	8	!! 7.0 ! 1996-01 (G. Madec) suppression of common work arrays
[1566]	9	!! - ! 1996-05 (G. Madec) mask computed from tmask and sup-
	10	!! ! pression of the double computation of bmask
[2528]	11	!! 8.0 ! 1997-02 (G. Madec) mesh information put in domhgr.F
	12	!! 8.1 ! 1997-07 (G. Madec) modification of mbathy and fmask
[1566]	13	!! - ! 1998-05 (G. Roullet) free surface
[2528]	14	!! 8.2 ! 2000-03 (G. Madec) no slip accurate
[1566]	15	!! - ! 2001-09 (J.-M. Molines) Open boundaries
	16	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	17	!! - ! 2005-11 (V. Garnier) Surface pressure gradient organization
	18	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
	19	!!----------------------------------------------------------------------
[3]	20
	21	!!----------------------------------------------------------------------
	22	!! dom_msk : compute land/ocean mask
[1566]	23	!! dom_msk_nsa : update land/ocean mask when no-slip accurate option is used.
[3]	24	!!----------------------------------------------------------------------
	25	USE oce ! ocean dynamics and tracers
	26	USE dom_oce ! ocean space and time domain
	27	USE in_out_manager ! I/O manager
	28	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[32]	29	USE lib_mpp
[367]	30	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[3294]	31	USE wrk_nemo ! Memory allocation
	32	USE timing ! Timing
[3]	33
	34	IMPLICIT NONE
	35	PRIVATE
	36
[2715]	37	PUBLIC dom_msk ! routine called by inidom.F90
	38	PUBLIC dom_msk_alloc ! routine called by nemogcm.F90
[3]	39
[1601]	40	! !!* Namelist namlbc : lateral boundary condition *
[4147]	41	REAL(wp) :: rn_shlat ! type of lateral boundary condition on velocity
	42	LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition
[3294]	43	! with analytical eqs.
[2715]	44
[3294]	45
[2715]	46	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icoord ! Workspace for dom_msk_nsa()
	47
[3]	48	!! * Substitutions
	49	# include "vectopt_loop_substitute.h90"
[1566]	50	!!----------------------------------------------------------------------
	51	!! NEMO/OPA 3.2 , LODYC-IPSL (2009)
[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
[9608]	235	tmask_i_3d(:,:,:) = tmask(:,:,:) ! Initialise 3D interior tmask with standard t mask
	236	! Now mask out any wrap columns
	237	tmask_i_3d( 1 :iif,:,:) = 0._wp ! first columns
	238	tmask_i_3d(iil:jpi,:,:) = 0._wp ! last columns (including mpp extra columns)
	239	! Now mask out any extra rows
	240	tmask_i_3d(:,1:ijf,:) = 0._wp ! first rows
	241	tmask_i_3d(:,ijl:jpj,:) = 0._wp ! last rows (including mpp extra rows)
	242
	243
[3]	244	! north fold mask
[1566]	245	! ---------------
[2528]	246	tpol(1:jpiglo) = 1._wp
	247	fpol(1:jpiglo) = 1._wp
[3]	248	IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot
[2528]	249	tpol(jpiglo/2+1:jpiglo) = 0._wp
	250	fpol( 1 :jpiglo) = 0._wp
[1566]	251	IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row
[291]	252	DO ji = iif+1, iil-1
	253	tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
[9608]	254	tmask_i_3d(ji,nlej-1,:) = tmask_i_3d(ji,nlej-1,:) * tpol(mig(ji))
[291]	255	END DO
	256	ENDIF
[3]	257	ENDIF
	258	IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot
[2528]	259	tpol( 1 :jpiglo) = 0._wp
	260	fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]	261	ENDIF
	262
	263	! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask)
	264	! -------------------------------------------
	265	DO jk = 1, jpk
	266	DO jj = 1, jpjm1
	267	DO ji = 1, fs_jpim1 ! vector loop
	268	umask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk)
	269	vmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji ,jj+1,jk)
[1271]	270	END DO
[1694]	271	DO ji = 1, jpim1 ! NO vector opt.
[3]	272	fmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk) &
[62]	273	& * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]	274	END DO
	275	END DO
	276	END DO
[4990]	277	! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
	278	DO jj = 1, jpjm1
	279	DO ji = 1, fs_jpim1 ! vector loop
	280	umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))
	281	vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))
	282	END DO
	283	DO ji = 1, jpim1 ! NO vector opt.
	284	fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &
	285	& * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
	286	END DO
	287	END DO
[2528]	288	CALL lbc_lnk( umask, 'U', 1._wp ) ! Lateral boundary conditions
	289	CALL lbc_lnk( vmask, 'V', 1._wp )
	290	CALL lbc_lnk( fmask, 'F', 1._wp )
[4990]	291	CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions
	292	CALL lbc_lnk( vmask_i, 'V', 1._wp )
	293	CALL lbc_lnk( fmask_i, 'F', 1._wp )
[3]	294
[5120]	295	! 3. Ocean/land mask at wu-, wv- and w points
	296	!----------------------------------------------
	297	wmask (:,:,1) = tmask(:,:,1) ! ????????
	298	wumask(:,:,1) = umask(:,:,1) ! ????????
	299	wvmask(:,:,1) = vmask(:,:,1) ! ????????
	300	DO jk=2,jpk
	301	wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
	302	wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)
	303	wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
	304	END DO
[3]	305
	306	! 4. ocean/land mask for the elliptic equation
	307	! --------------------------------------------
[4990]	308	bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point
[1566]	309	!
	310	! ! Boundary conditions
	311	! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
[3]	312	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
[2528]	313	bmask( 1 ,:) = 0._wp
	314	bmask(jpi,:) = 0._wp
[3]	315	ENDIF
[1566]	316	IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1
[2528]	317	bmask(:, 1 ) = 0._wp
[3]	318	ENDIF
[1566]	319	! ! north fold :
	320	IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
	321	DO ji = 1, jpi
[1528]	322	ii = ji + nimpp - 1
	323	bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
[2528]	324	bmask(ji,jpj ) = 0._wp
[1528]	325	END DO
[3]	326	ENDIF
[1566]	327	IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[2528]	328	bmask(:,jpj) = 0._wp
[3]	329	ENDIF
[1566]	330	!
	331	IF( lk_mpp ) THEN ! mpp specificities
	332	! ! bmask is set to zero on the overlap region
[2528]	333	IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp
	334	IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp
	335	IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp
	336	IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp
[1566]	337	!
	338	IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
[1528]	339	DO ji = 1, nlci
	340	ii = ji + nimpp - 1
	341	bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
[2528]	342	bmask(ji,nlcj ) = 0._wp
[1528]	343	END DO
[32]	344	ENDIF
[1566]	345	IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[1528]	346	DO ji = 1, nlci
[2528]	347	bmask(ji,nlcj ) = 0._wp
[1528]	348	END DO
[32]	349	ENDIF
[3]	350	ENDIF
	351
	352
	353	! mask for second order calculation of vorticity
	354	! ----------------------------------------------
	355	CALL dom_msk_nsa
	356
	357
	358	! Lateral boundary conditions on velocity (modify fmask)
[1566]	359	! ---------------------------------------
[3]	360	DO jk = 1, jpk
[1566]	361	zwf(:,:) = fmask(:,:,jk)
[3]	362	DO jj = 2, jpjm1
	363	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	364	IF( fmask(ji,jj,jk) == 0._wp ) THEN
	365	fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	366	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
[3]	367	ENDIF
	368	END DO
	369	END DO
	370	DO jj = 2, jpjm1
[2528]	371	IF( fmask(1,jj,jk) == 0._wp ) THEN
	372	fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]	373	ENDIF
[2528]	374	IF( fmask(jpi,jj,jk) == 0._wp ) THEN
	375	fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]	376	ENDIF
[1566]	377	END DO
[3]	378	DO ji = 2, jpim1
[2528]	379	IF( fmask(ji,1,jk) == 0._wp ) THEN
	380	fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]	381	ENDIF
[2528]	382	IF( fmask(ji,jpj,jk) == 0._wp ) THEN
	383	fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]	384	ENDIF
	385	END DO
	386	END DO
[1566]	387	!
	388	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration
	389	! ! Increased lateral friction near of some straits
[2528]	390	IF( nn_cla == 0 ) THEN
[1273]	391	! ! Gibraltar strait : partial slip (fmask=0.5)
	392	ij0 = 101 ; ij1 = 101
[2528]	393	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[32]	394	ij0 = 102 ; ij1 = 102
[2528]	395	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[1273]	396	!
	397	! ! Bab el Mandeb : partial slip (fmask=1)
	398	ij0 = 87 ; ij1 = 88
[2528]	399	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	400	ij0 = 88 ; ij1 = 88
[2528]	401	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	402	!
[3]	403	ENDIF
[1707]	404	! ! Danish straits : strong slip (fmask > 2)
	405	! We keep this as an example but it is instable in this case
	406	! ij0 = 115 ; ij1 = 115
[2528]	407	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	408	! ij0 = 116 ; ij1 = 116
[2528]	409	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	410	!
	411	ENDIF
[1566]	412	!
[2528]	413	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	414	! ! Increased lateral friction near of some straits
[5506]	415	! This dirty section will be suppressed by simplification process:
	416	! all this will come back in input files
	417	! Currently these hard-wired indices relate to configuration with
	418	! extend grid (jpjglo=332)
	419	!
	420	isrow = 332 - jpjglo
	421	!
[2528]	422	IF(lwp) WRITE(numout,*)
	423	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	424	IF(lwp) WRITE(numout,*) ' Gibraltar '
[5385]	425	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait
[6487]	426	ij0 = 241 - isrow ; ij1 = 241 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[1566]	427
[2528]	428	IF(lwp) WRITE(numout,*) ' Bhosporus '
[5385]	429	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
[6487]	430	ij0 = 248 - isrow ; ij1 = 248 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	431
	432	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
[5385]	433	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
[6487]	434	ij0 = 189 - isrow ; ij1 = 190 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	435
	436	IF(lwp) WRITE(numout,*) ' Lombok '
[5385]	437	ii0 = 44 ; ii1 = 44 ! Lombok Strait
[6487]	438	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	439
	440	IF(lwp) WRITE(numout,*) ' Ombai '
[5385]	441	ii0 = 53 ; ii1 = 53 ! Ombai Strait
[6487]	442	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	443
	444	IF(lwp) WRITE(numout,*) ' Timor Passage '
[5385]	445	ii0 = 56 ; ii1 = 56 ! Timor Passage
[6487]	446	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	447
	448	IF(lwp) WRITE(numout,*) ' West Halmahera '
[5385]	449	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
[6487]	450	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	451
	452	IF(lwp) WRITE(numout,*) ' East Halmahera '
[5385]	453	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
[6487]	454	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	455	!
	456	ENDIF
	457	!
[6491]	458	IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN
	459	! ! ORCA_R025 configuration
	460	! ! Increased lateral friction on parts of Antarctic coastline
	461	! ! for increased stability
	462	! ! NB. This only works to do this here if we have free slip
	463	! ! generally, so fmask is zero at coast points.
	464	IF(lwp) WRITE(numout,*)
	465	IF(lwp) WRITE(numout,*) ' orca_r025: increase friction in following regions : '
	466	IF(lwp) WRITE(numout,*) ' whole Antarctic coastline: partial slip shlat=1 '
	467
	468	zphi_drake_passage = -58.0_wp
	469	zshlat_antarc = 1.0_wp
	470	zwf(:,:) = fmask(:,:,1)
	471	DO jj = 2, jpjm1
	472	DO ji = fs_2, fs_jpim1 ! vector opt.
	473	IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN
	474	fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	475	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
	476	ENDIF
	477	END DO
	478	END DO
	479	END IF
	480	!
[2528]	481	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	482
[3294]	483	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	484
[1566]	485	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	486	imsk(:,:) = INT( tmask_i(:,:) )
	487	WRITE(numout,*) ' tmask_i : '
	488	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	489	& 1, jpj, 1, 1, numout)
	490	WRITE (numout,*)
	491	WRITE (numout,*) ' dommsk: tmask for each level'
	492	WRITE (numout,*) ' ----------------------------'
	493	DO jk = 1, jpk
	494	imsk(:,:) = INT( tmask(:,:,jk) )
	495
	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: vmask for each level'
	503	WRITE(numout,*) ' -----------------------------'
	504	DO jk = 1, jpk
	505	imsk(:,:) = INT( vmask(:,:,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: fmask for each level'
	513	WRITE(numout,*) ' -----------------------------'
	514	DO jk = 1, jpk
	515	imsk(:,:) = INT( fmask(:,:,jk) )
	516	WRITE(numout,*)
	517	WRITE(numout,*) ' level = ',jk
	518	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	519	& 1, jpj, 1, 1, numout )
	520	END DO
	521	WRITE(numout,*)
	522	WRITE(numout,*) ' dom_msk: bmask '
	523	WRITE(numout,*) ' ---------------'
	524	WRITE(numout,*)
	525	imsk(:,:) = INT( bmask(:,:) )
	526	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	527	& 1, jpj, 1, 1, numout )
[3]	528	ENDIF
[1566]	529	!
[3294]	530	CALL wrk_dealloc( jpi, jpj, imsk )
	531	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	532	!
[3294]	533	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	534	!
[3]	535	END SUBROUTINE dom_msk
	536
	537	#if defined key_noslip_accurate
	538	!!----------------------------------------------------------------------
	539	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	540	!!----------------------------------------------------------------------
	541
	542	SUBROUTINE dom_msk_nsa
	543	!!---------------------------------------------------------------------
	544	!! * ROUTINE dom_msk_nsa *
	545	!!
	546	!! ** Purpose :
	547	!!
	548	!! ** Method :
	549	!!
	550	!! ** Action :
	551	!!----------------------------------------------------------------------
[2715]	552	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	553	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	554	REAL(wp) :: zaa
[3]	555	!!---------------------------------------------------------------------
[3294]	556	!
	557	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	558	!
[2715]	559	IF(lwp) WRITE(numout,*)
	560	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	561	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
[8280]	562	IF( lk_mpp ) CALL ctl_stop('STOP', ' mpp version is not yet implemented' )
[3]	563
	564	! mask for second order calculation of vorticity
	565	! ----------------------------------------------
	566	! noslip boundary condition: fmask=1 at convex corner, store
	567	! index of straight coast meshes ( 'west', refering to a coast,
	568	! means west of the ocean, aso)
	569
	570	DO jk = 1, jpk
	571	DO jl = 1, 4
	572	npcoa(jl,jk) = 0
	573	DO ji = 1, 2*(jpi+jpj)
	574	nicoa(ji,jl,jk) = 0
	575	njcoa(ji,jl,jk) = 0
	576	END DO
	577	END DO
	578	END DO
	579
	580	IF( jperio == 2 ) THEN
	581	WRITE(numout,*) ' '
	582	WRITE(numout,*) ' symetric boundary conditions need special'
	583	WRITE(numout,*) ' treatment not implemented. we stop.'
[8280]	584	CALL ctl_stop('STOP', 'NEMO abort from dom_msk_nsa')
[3]	585	ENDIF
	586
	587	! convex corners
	588
	589	DO jk = 1, jpkm1
	590	DO jj = 1, jpjm1
	591	DO ji = 1, jpim1
	592	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	593	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	594	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	595	END DO
	596	END DO
	597	END DO
	598
	599	! north-south straight coast
	600
	601	DO jk = 1, jpkm1
	602	inw = 0
	603	ine = 0
	604	DO jj = 2, jpjm1
	605	DO ji = 2, jpim1
	606	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	607	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	608	inw = inw + 1
	609	nicoa(inw,1,jk) = ji
	610	njcoa(inw,1,jk) = jj
	611	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	612	ENDIF
	613	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	614	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	615	ine = ine + 1
	616	nicoa(ine,2,jk) = ji
	617	njcoa(ine,2,jk) = jj
	618	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	619	ENDIF
	620	END DO
	621	END DO
	622	npcoa(1,jk) = inw
	623	npcoa(2,jk) = ine
	624	END DO
	625
	626	! west-east straight coast
	627
	628	DO jk = 1, jpkm1
	629	ins = 0
	630	inn = 0
	631	DO jj = 2, jpjm1
	632	DO ji =2, jpim1
	633	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	634	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	635	ins = ins + 1
	636	nicoa(ins,3,jk) = ji
	637	njcoa(ins,3,jk) = jj
	638	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	639	ENDIF
	640	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	641	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	642	inn = inn + 1
	643	nicoa(inn,4,jk) = ji
	644	njcoa(inn,4,jk) = jj
	645	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	646	ENDIF
	647	END DO
	648	END DO
	649	npcoa(3,jk) = ins
	650	npcoa(4,jk) = inn
	651	END DO
	652
	653	itest = 2 * ( jpi + jpj )
	654	DO jk = 1, jpk
	655	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	656	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	657
	658	WRITE(ctmp1,*) ' level jk = ',jk
	659	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	660	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	661	& npcoa(3,jk), npcoa(4,jk)
[474]	662	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	663	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	664	ENDIF
	665	END DO
	666
	667	ierror = 0
	668	iind = 0
	669	ijnd = 0
[2528]	670	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	671	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	672	DO jk = 1, jpk
	673	DO jl = 1, npcoa(1,jk)
	674	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	675	ierror = ierror+1
	676	icoord(ierror,1) = nicoa(jl,1,jk)
	677	icoord(ierror,2) = njcoa(jl,1,jk)
	678	icoord(ierror,3) = jk
	679	ENDIF
	680	END DO
	681	DO jl = 1, npcoa(2,jk)
	682	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	683	ierror = ierror + 1
	684	icoord(ierror,1) = nicoa(jl,2,jk)
	685	icoord(ierror,2) = njcoa(jl,2,jk)
	686	icoord(ierror,3) = jk
	687	ENDIF
	688	END DO
	689	DO jl = 1, npcoa(3,jk)
	690	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	691	ierror = ierror + 1
	692	icoord(ierror,1) = nicoa(jl,3,jk)
	693	icoord(ierror,2) = njcoa(jl,3,jk)
	694	icoord(ierror,3) = jk
	695	ENDIF
	696	END DO
[2528]	697	DO jl = 1, npcoa(4,jk)
[3]	698	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	699	ierror=ierror + 1
	700	icoord(ierror,1) = nicoa(jl,4,jk)
	701	icoord(ierror,2) = njcoa(jl,4,jk)
	702	icoord(ierror,3) = jk
[3]	703	ENDIF
	704	END DO
	705	END DO
	706
	707	IF( ierror > 0 ) THEN
	708	IF(lwp) WRITE(numout,*)
	709	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	710	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	711	DO jl = 1, ierror
	712	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	713	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	714	END DO
[474]	715	CALL ctl_stop( 'We stop...' )
[3]	716	ENDIF
[2715]	717	!
[3294]	718	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	719	!
[3]	720	END SUBROUTINE dom_msk_nsa
	721
	722	#else
	723	!!----------------------------------------------------------------------
	724	!! Default option : Empty routine
	725	!!----------------------------------------------------------------------
	726	SUBROUTINE dom_msk_nsa
	727	END SUBROUTINE dom_msk_nsa
	728	#endif
	729
	730	!!======================================================================
	731	END MODULE dommsk

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

Download in other formats: