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

Context Navigation

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 9 years ago
Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.
File size: 31.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)
[5234]	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
[5443]	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
[5443]	284	! 3. Ocean/land mask at wu-, wv- and w points
	285	!----------------------------------------------
	286	wmask (:,:,1) = tmask(:,:,1) ! ????????
	287	wumask(:,:,1) = umask(:,:,1) ! ????????
	288	wvmask(:,:,1) = vmask(:,:,1) ! ????????
	289	DO jk=2,jpk
	290	wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
	291	wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)
	292	wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
	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
[2528]	379	IF( nn_cla == 0 ) THEN
[1273]	380	! ! Gibraltar strait : partial slip (fmask=0.5)
	381	ij0 = 101 ; ij1 = 101
[2528]	382	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[32]	383	ij0 = 102 ; ij1 = 102
[2528]	384	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[1273]	385	!
	386	! ! Bab el Mandeb : partial slip (fmask=1)
	387	ij0 = 87 ; ij1 = 88
[2528]	388	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	389	ij0 = 88 ; ij1 = 88
[2528]	390	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	391	!
[3]	392	ENDIF
[1707]	393	! ! Danish straits : strong slip (fmask > 2)
	394	! We keep this as an example but it is instable in this case
	395	! ij0 = 115 ; ij1 = 115
[2528]	396	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	397	! ij0 = 116 ; ij1 = 116
[2528]	398	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	399	!
	400	ENDIF
[1566]	401	!
[2528]	402	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	403	! ! Increased lateral friction near of some straits
[5443]	404	! This dirty section will be suppressed by simplification process: all this will come back in input files
	405	! Currently these hard-wired indices relate to the original (pre-v3.6) configuration
	406	! which had a grid-size of 362x292.
	407	! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6.
	408	! The original domain can still be used optionally if the isf code is not activated.
	409	! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
	410	!
	411	IF ( jpjglo == 292 ) THEN ; isrow = 0 ! Using pre-v3.6 files or adjusted start row from isf-extended grid
	412	ELSEIF( jpjglo == 332 ) THEN ; isrow = 40 ! Using full isfextended domain.
	413	ENDIF
	414
[2528]	415	IF(lwp) WRITE(numout,*)
	416	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	417	IF(lwp) WRITE(numout,*) ' Gibraltar '
[5443]	418	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait
	419	ij0 = 201 + isrow ; ij1 = 201 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[1566]	420
[2528]	421	IF(lwp) WRITE(numout,*) ' Bhosporus '
[5443]	422	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
	423	ij0 = 208 + isrow ; ij1 = 208 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	424
	425	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
[5443]	426	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
	427	ij0 = 149 + isrow ; ij1 = 150 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	428
	429	IF(lwp) WRITE(numout,*) ' Lombok '
[5443]	430	ii0 = 44 ; ii1 = 44 ! Lombok Strait
	431	ij0 = 124 + isrow ; ij1 = 125 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	432
	433	IF(lwp) WRITE(numout,*) ' Ombai '
[5443]	434	ii0 = 53 ; ii1 = 53 ! Ombai Strait
	435	ij0 = 124 + isrow ; ij1 = 125 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	436
	437	IF(lwp) WRITE(numout,*) ' Timor Passage '
[5443]	438	ii0 = 56 ; ii1 = 56 ! Timor Passage
	439	ij0 = 124 + isrow ; ij1 = 125 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	440
	441	IF(lwp) WRITE(numout,*) ' West Halmahera '
[5443]	442	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
	443	ij0 = 141 + isrow ; ij1 = 142 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	444
	445	IF(lwp) WRITE(numout,*) ' East Halmahera '
[5443]	446	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
	447	ij0 = 141 + isrow ; ij1 = 142 + isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	448	!
	449	ENDIF
	450	!
	451	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	452
[3294]	453	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	454
[1566]	455	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	456	imsk(:,:) = INT( tmask_i(:,:) )
	457	WRITE(numout,*) ' tmask_i : '
	458	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	459	& 1, jpj, 1, 1, numout)
	460	WRITE (numout,*)
	461	WRITE (numout,*) ' dommsk: tmask for each level'
	462	WRITE (numout,*) ' ----------------------------'
	463	DO jk = 1, jpk
	464	imsk(:,:) = INT( tmask(:,:,jk) )
	465
	466	WRITE(numout,*)
	467	WRITE(numout,*) ' level = ',jk
	468	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	469	& 1, jpj, 1, 1, numout)
	470	END DO
	471	WRITE(numout,*)
	472	WRITE(numout,*) ' dom_msk: vmask for each level'
	473	WRITE(numout,*) ' -----------------------------'
	474	DO jk = 1, jpk
	475	imsk(:,:) = INT( vmask(:,:,jk) )
	476	WRITE(numout,*)
	477	WRITE(numout,*) ' level = ',jk
	478	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	479	& 1, jpj, 1, 1, numout)
	480	END DO
	481	WRITE(numout,*)
	482	WRITE(numout,*) ' dom_msk: fmask for each level'
	483	WRITE(numout,*) ' -----------------------------'
	484	DO jk = 1, jpk
	485	imsk(:,:) = INT( fmask(:,:,jk) )
	486	WRITE(numout,*)
	487	WRITE(numout,*) ' level = ',jk
	488	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	489	& 1, jpj, 1, 1, numout )
	490	END DO
	491	WRITE(numout,*)
	492	WRITE(numout,*) ' dom_msk: bmask '
	493	WRITE(numout,*) ' ---------------'
	494	WRITE(numout,*)
	495	imsk(:,:) = INT( bmask(:,:) )
	496	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	497	& 1, jpj, 1, 1, numout )
[3]	498	ENDIF
[1566]	499	!
[3294]	500	CALL wrk_dealloc( jpi, jpj, imsk )
	501	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	502	!
[3294]	503	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	504	!
[3]	505	END SUBROUTINE dom_msk
	506
	507	#if defined key_noslip_accurate
	508	!!----------------------------------------------------------------------
	509	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	510	!!----------------------------------------------------------------------
	511
	512	SUBROUTINE dom_msk_nsa
	513	!!---------------------------------------------------------------------
	514	!! * ROUTINE dom_msk_nsa *
	515	!!
	516	!! ** Purpose :
	517	!!
	518	!! ** Method :
	519	!!
	520	!! ** Action :
	521	!!----------------------------------------------------------------------
[2715]	522	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	523	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	524	REAL(wp) :: zaa
[3]	525	!!---------------------------------------------------------------------
[3294]	526	!
	527	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	528	!
[2715]	529	IF(lwp) WRITE(numout,*)
	530	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	531	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
	532	IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
[3]	533
	534	! mask for second order calculation of vorticity
	535	! ----------------------------------------------
	536	! noslip boundary condition: fmask=1 at convex corner, store
	537	! index of straight coast meshes ( 'west', refering to a coast,
	538	! means west of the ocean, aso)
	539
	540	DO jk = 1, jpk
	541	DO jl = 1, 4
	542	npcoa(jl,jk) = 0
	543	DO ji = 1, 2*(jpi+jpj)
	544	nicoa(ji,jl,jk) = 0
	545	njcoa(ji,jl,jk) = 0
	546	END DO
	547	END DO
	548	END DO
	549
	550	IF( jperio == 2 ) THEN
	551	WRITE(numout,*) ' '
	552	WRITE(numout,*) ' symetric boundary conditions need special'
	553	WRITE(numout,*) ' treatment not implemented. we stop.'
	554	STOP
	555	ENDIF
	556
	557	! convex corners
	558
	559	DO jk = 1, jpkm1
	560	DO jj = 1, jpjm1
	561	DO ji = 1, jpim1
	562	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	563	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	564	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	565	END DO
	566	END DO
	567	END DO
	568
	569	! north-south straight coast
	570
	571	DO jk = 1, jpkm1
	572	inw = 0
	573	ine = 0
	574	DO jj = 2, jpjm1
	575	DO ji = 2, jpim1
	576	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	577	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	578	inw = inw + 1
	579	nicoa(inw,1,jk) = ji
	580	njcoa(inw,1,jk) = jj
	581	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	582	ENDIF
	583	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	584	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	585	ine = ine + 1
	586	nicoa(ine,2,jk) = ji
	587	njcoa(ine,2,jk) = jj
	588	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	589	ENDIF
	590	END DO
	591	END DO
	592	npcoa(1,jk) = inw
	593	npcoa(2,jk) = ine
	594	END DO
	595
	596	! west-east straight coast
	597
	598	DO jk = 1, jpkm1
	599	ins = 0
	600	inn = 0
	601	DO jj = 2, jpjm1
	602	DO ji =2, jpim1
	603	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	604	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	605	ins = ins + 1
	606	nicoa(ins,3,jk) = ji
	607	njcoa(ins,3,jk) = jj
	608	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	609	ENDIF
	610	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	611	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	612	inn = inn + 1
	613	nicoa(inn,4,jk) = ji
	614	njcoa(inn,4,jk) = jj
	615	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	616	ENDIF
	617	END DO
	618	END DO
	619	npcoa(3,jk) = ins
	620	npcoa(4,jk) = inn
	621	END DO
	622
	623	itest = 2 * ( jpi + jpj )
	624	DO jk = 1, jpk
	625	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	626	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	627
	628	WRITE(ctmp1,*) ' level jk = ',jk
	629	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	630	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	631	& npcoa(3,jk), npcoa(4,jk)
[474]	632	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	633	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	634	ENDIF
	635	END DO
	636
	637	ierror = 0
	638	iind = 0
	639	ijnd = 0
[2528]	640	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	641	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	642	DO jk = 1, jpk
	643	DO jl = 1, npcoa(1,jk)
	644	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	645	ierror = ierror+1
	646	icoord(ierror,1) = nicoa(jl,1,jk)
	647	icoord(ierror,2) = njcoa(jl,1,jk)
	648	icoord(ierror,3) = jk
	649	ENDIF
	650	END DO
	651	DO jl = 1, npcoa(2,jk)
	652	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	653	ierror = ierror + 1
	654	icoord(ierror,1) = nicoa(jl,2,jk)
	655	icoord(ierror,2) = njcoa(jl,2,jk)
	656	icoord(ierror,3) = jk
	657	ENDIF
	658	END DO
	659	DO jl = 1, npcoa(3,jk)
	660	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	661	ierror = ierror + 1
	662	icoord(ierror,1) = nicoa(jl,3,jk)
	663	icoord(ierror,2) = njcoa(jl,3,jk)
	664	icoord(ierror,3) = jk
	665	ENDIF
	666	END DO
[2528]	667	DO jl = 1, npcoa(4,jk)
[3]	668	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	669	ierror=ierror + 1
	670	icoord(ierror,1) = nicoa(jl,4,jk)
	671	icoord(ierror,2) = njcoa(jl,4,jk)
	672	icoord(ierror,3) = jk
[3]	673	ENDIF
	674	END DO
	675	END DO
	676
	677	IF( ierror > 0 ) THEN
	678	IF(lwp) WRITE(numout,*)
	679	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	680	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	681	DO jl = 1, ierror
	682	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	683	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	684	END DO
[474]	685	CALL ctl_stop( 'We stop...' )
[3]	686	ENDIF
[2715]	687	!
[3294]	688	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	689	!
[3]	690	END SUBROUTINE dom_msk_nsa
	691
	692	#else
	693	!!----------------------------------------------------------------------
	694	!! Default option : Empty routine
	695	!!----------------------------------------------------------------------
	696	SUBROUTINE dom_msk_nsa
	697	END SUBROUTINE dom_msk_nsa
	698	#endif
	699
	700	!!======================================================================
	701	END MODULE dommsk

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

Download in other formats: