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

Context Navigation

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5802

Last change on this file since 5802 was 5802, checked in by mathiot, 9 years ago
ice sheet coupling: reproducibility fixed in MISOMIP configuration + contrain bathy to be constant during all the run
Property svn:keywords set to `Id`
File size: 30.9 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
[5385]	136	INTEGER :: isrow ! index for ORCA1 starting row
[3294]	137	INTEGER , POINTER, DIMENSION(:,:) :: imsk
	138	REAL(wp), POINTER, DIMENSION(:,:) :: zwf
[1601]	139	!!
[3294]	140	NAMELIST/namlbc/ rn_shlat, ln_vorlat
[3]	141	!!---------------------------------------------------------------------
[3294]	142	!
	143	IF( nn_timing == 1 ) CALL timing_start('dom_msk')
	144	!
	145	CALL wrk_alloc( jpi, jpj, imsk )
	146	CALL wrk_alloc( jpi, jpj, zwf )
	147	!
[4147]	148	REWIND( numnam_ref ) ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
	149	READ ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
	150	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
	151
	152	REWIND( numnam_cfg ) ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
	153	READ ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
	154	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
[4624]	155	IF(lwm) WRITE ( numond, namlbc )
[1566]	156
	157	IF(lwp) THEN ! control print
[3]	158	WRITE(numout,*)
	159	WRITE(numout,*) 'dommsk : ocean mask '
	160	WRITE(numout,*) '~~~~~~'
[1566]	161	WRITE(numout,*) ' Namelist namlbc'
[3294]	162	WRITE(numout,*) ' lateral momentum boundary cond. rn_shlat = ',rn_shlat
	163	WRITE(numout,*) ' consistency with analytical form ln_vorlat = ',ln_vorlat
[3]	164	ENDIF
	165
[2528]	166	IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral free-slip '
[1601]	167	ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral no-slip '
	168	ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip '
	169	ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip '
	170	ELSE
	171	WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
	172	CALL ctl_stop( ctmp1 )
[3]	173	ENDIF
	174
	175	! 1. Ocean/land mask at t-point (computed from mbathy)
	176	! -----------------------------
[1566]	177	! N.B. tmask has already the right boundary conditions since mbathy is ok
	178	!
[2528]	179	tmask(:,:,:) = 0._wp
[3]	180	DO jk = 1, jpk
	181	DO jj = 1, jpj
	182	DO ji = 1, jpi
[2528]	183	IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) tmask(ji,jj,jk) = 1._wp
[3]	184	END DO
	185	END DO
	186	END DO
[4990]	187
	188	! (ISF) define barotropic mask and mask the ice shelf point
	189	ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
	190
	191	DO jk = 1, jpk
	192	DO jj = 1, jpj
	193	DO ji = 1, jpi
	194	IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp ) THEN
	195	tmask(ji,jj,jk) = 0._wp
	196	END IF
	197	END DO
	198	END DO
	199	END DO
[3]	200
[1566]	201	!!gm ????
[255]	202	#if defined key_zdfkpp
[1601]	203	IF( cp_cfg == 'orca' ) THEN
	204	IF( jp_cfg == 2 ) THEN ! land point on Bab el Mandeb zonal section
[291]	205	ij0 = 87 ; ij1 = 88
	206	ii0 = 160 ; ii1 = 161
[2528]	207	tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
[291]	208	ELSE
	209	IF(lwp) WRITE(numout,*)
	210	IF(lwp) WRITE(numout,cform_war)
	211	IF(lwp) WRITE(numout,*)
	212	IF(lwp) WRITE(numout,*)' A mask must be applied on Bab el Mandeb strait'
	213	IF(lwp) WRITE(numout,*)' in case of ORCAs configurations'
	214	IF(lwp) WRITE(numout,*)' This is a problem which is not yet solved'
	215	IF(lwp) WRITE(numout,*)
	216	ENDIF
	217	ENDIF
[255]	218	#endif
[1566]	219	!!gm end
[291]	220
[3]	221	! Interior domain mask (used for global sum)
	222	! --------------------
[4990]	223	tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf
[5802]	224
	225	tmask_h(:,:) = 1._wp ! 0 on the halo and 1 elsewhere
[3]	226	iif = jpreci ! ???
	227	iil = nlci - jpreci + 1
	228	ijf = jprecj ! ???
	229	ijl = nlcj - jprecj + 1
	230
[5802]	231	tmask_h( 1 :iif, : ) = 0._wp ! first columns
	232	tmask_h(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)
	233	tmask_h( : , 1 :ijf) = 0._wp ! first rows
	234	tmask_h( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)
[3]	235
	236	! north fold mask
[1566]	237	! ---------------
[2528]	238	tpol(1:jpiglo) = 1._wp
	239	fpol(1:jpiglo) = 1._wp
[3]	240	IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot
[2528]	241	tpol(jpiglo/2+1:jpiglo) = 0._wp
	242	fpol( 1 :jpiglo) = 0._wp
[1566]	243	IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row
[291]	244	DO ji = iif+1, iil-1
[5802]	245	tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
[291]	246	END DO
	247	ENDIF
[3]	248	ENDIF
[5802]	249
	250	tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:)
	251
[3]	252	IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot
[2528]	253	tpol( 1 :jpiglo) = 0._wp
	254	fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]	255	ENDIF
	256
	257	! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask)
	258	! -------------------------------------------
	259	DO jk = 1, jpk
	260	DO jj = 1, jpjm1
	261	DO ji = 1, fs_jpim1 ! vector loop
	262	umask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk)
	263	vmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji ,jj+1,jk)
[1271]	264	END DO
[1694]	265	DO ji = 1, jpim1 ! NO vector opt.
[3]	266	fmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk) &
[62]	267	& * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]	268	END DO
	269	END DO
	270	END DO
[5619]	271	! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
[4990]	272	DO jj = 1, jpjm1
	273	DO ji = 1, fs_jpim1 ! vector loop
[5619]	274	ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))
	275	ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))
[4990]	276	END DO
	277	DO ji = 1, jpim1 ! NO vector opt.
[5619]	278	ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &
[4990]	279	& * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
	280	END DO
	281	END DO
[5619]	282	CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions
	283	CALL lbc_lnk( vmask , 'V', 1._wp )
	284	CALL lbc_lnk( fmask , 'F', 1._wp )
	285	CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions
	286	CALL lbc_lnk( ssvmask, 'V', 1._wp )
	287	CALL lbc_lnk( ssfmask, 'F', 1._wp )
[3]	288
[5120]	289	! 3. Ocean/land mask at wu-, wv- and w points
	290	!----------------------------------------------
	291	wmask (:,:,1) = tmask(:,:,1) ! ????????
	292	wumask(:,:,1) = umask(:,:,1) ! ????????
	293	wvmask(:,:,1) = vmask(:,:,1) ! ????????
	294	DO jk=2,jpk
	295	wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
	296	wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)
	297	wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
	298	END DO
[3]	299
	300	! 4. ocean/land mask for the elliptic equation
	301	! --------------------------------------------
[4990]	302	bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point
[1566]	303	!
	304	! ! Boundary conditions
	305	! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
[3]	306	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
[2528]	307	bmask( 1 ,:) = 0._wp
	308	bmask(jpi,:) = 0._wp
[3]	309	ENDIF
[1566]	310	IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1
[2528]	311	bmask(:, 1 ) = 0._wp
[3]	312	ENDIF
[1566]	313	! ! north fold :
	314	IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
	315	DO ji = 1, jpi
[1528]	316	ii = ji + nimpp - 1
	317	bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
[2528]	318	bmask(ji,jpj ) = 0._wp
[1528]	319	END DO
[3]	320	ENDIF
[1566]	321	IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[2528]	322	bmask(:,jpj) = 0._wp
[3]	323	ENDIF
[1566]	324	!
	325	IF( lk_mpp ) THEN ! mpp specificities
	326	! ! bmask is set to zero on the overlap region
[2528]	327	IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp
	328	IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp
	329	IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp
	330	IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp
[1566]	331	!
	332	IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
[1528]	333	DO ji = 1, nlci
	334	ii = ji + nimpp - 1
	335	bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
[2528]	336	bmask(ji,nlcj ) = 0._wp
[1528]	337	END DO
[32]	338	ENDIF
[1566]	339	IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[1528]	340	DO ji = 1, nlci
[2528]	341	bmask(ji,nlcj ) = 0._wp
[1528]	342	END DO
[32]	343	ENDIF
[3]	344	ENDIF
	345
	346
	347	! mask for second order calculation of vorticity
	348	! ----------------------------------------------
	349	CALL dom_msk_nsa
	350
	351
	352	! Lateral boundary conditions on velocity (modify fmask)
[1566]	353	! ---------------------------------------
[3]	354	DO jk = 1, jpk
[1566]	355	zwf(:,:) = fmask(:,:,jk)
[3]	356	DO jj = 2, jpjm1
	357	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	358	IF( fmask(ji,jj,jk) == 0._wp ) THEN
	359	fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
	360	& zwf(ji-1,jj), zwf(ji,jj-1) ) )
[3]	361	ENDIF
	362	END DO
	363	END DO
	364	DO jj = 2, jpjm1
[2528]	365	IF( fmask(1,jj,jk) == 0._wp ) THEN
	366	fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]	367	ENDIF
[2528]	368	IF( fmask(jpi,jj,jk) == 0._wp ) THEN
	369	fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]	370	ENDIF
[1566]	371	END DO
[3]	372	DO ji = 2, jpim1
[2528]	373	IF( fmask(ji,1,jk) == 0._wp ) THEN
	374	fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]	375	ENDIF
[2528]	376	IF( fmask(ji,jpj,jk) == 0._wp ) THEN
	377	fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]	378	ENDIF
	379	END DO
	380	END DO
[1566]	381	!
	382	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration
	383	! ! Increased lateral friction near of some straits
[2528]	384	IF( nn_cla == 0 ) THEN
[1273]	385	! ! Gibraltar strait : partial slip (fmask=0.5)
	386	ij0 = 101 ; ij1 = 101
[2528]	387	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[32]	388	ij0 = 102 ; ij1 = 102
[2528]	389	ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
[1273]	390	!
	391	! ! Bab el Mandeb : partial slip (fmask=1)
	392	ij0 = 87 ; ij1 = 88
[2528]	393	ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	394	ij0 = 88 ; ij1 = 88
[2528]	395	ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
[1273]	396	!
[3]	397	ENDIF
[1707]	398	! ! Danish straits : strong slip (fmask > 2)
	399	! We keep this as an example but it is instable in this case
	400	! ij0 = 115 ; ij1 = 115
[2528]	401	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]	402	! ij0 = 116 ; ij1 = 116
[2528]	403	! ii0 = 145 ; ii1 = 146 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]	404	!
	405	ENDIF
[1566]	406	!
[2528]	407	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	408	! ! Increased lateral friction near of some straits
[5506]	409	! This dirty section will be suppressed by simplification process:
	410	! all this will come back in input files
	411	! Currently these hard-wired indices relate to configuration with
	412	! extend grid (jpjglo=332)
	413	!
	414	isrow = 332 - jpjglo
	415	!
[2528]	416	IF(lwp) WRITE(numout,*)
	417	IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
	418	IF(lwp) WRITE(numout,*) ' Gibraltar '
[5385]	419	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait
[5552]	420	ij0 = 241 - isrow ; ij1 = 241 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[1566]	421
[2528]	422	IF(lwp) WRITE(numout,*) ' Bhosporus '
[5385]	423	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
[5552]	424	ij0 = 248 - isrow ; ij1 = 248 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	425
	426	IF(lwp) WRITE(numout,*) ' Makassar (Top) '
[5385]	427	ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
[5552]	428	ij0 = 189 - isrow ; ij1 = 190 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	429
	430	IF(lwp) WRITE(numout,*) ' Lombok '
[5385]	431	ii0 = 44 ; ii1 = 44 ! Lombok Strait
[5552]	432	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	433
	434	IF(lwp) WRITE(numout,*) ' Ombai '
[5385]	435	ii0 = 53 ; ii1 = 53 ! Ombai Strait
[5552]	436	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	437
	438	IF(lwp) WRITE(numout,*) ' Timor Passage '
[5385]	439	ii0 = 56 ; ii1 = 56 ! Timor Passage
[5552]	440	ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
[2528]	441
	442	IF(lwp) WRITE(numout,*) ' West Halmahera '
[5385]	443	ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
[5552]	444	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	445
	446	IF(lwp) WRITE(numout,*) ' East Halmahera '
[5385]	447	ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
[5552]	448	ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
[2528]	449	!
	450	ENDIF
	451	!
	452	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	453
[3294]	454	! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]	455
[1566]	456	IF( nprint == 1 .AND. lwp ) THEN ! Control print
[3]	457	imsk(:,:) = INT( tmask_i(:,:) )
	458	WRITE(numout,*) ' tmask_i : '
	459	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	460	& 1, jpj, 1, 1, numout)
	461	WRITE (numout,*)
	462	WRITE (numout,*) ' dommsk: tmask for each level'
	463	WRITE (numout,*) ' ----------------------------'
	464	DO jk = 1, jpk
	465	imsk(:,:) = INT( tmask(:,:,jk) )
	466
	467	WRITE(numout,*)
	468	WRITE(numout,*) ' level = ',jk
	469	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	470	& 1, jpj, 1, 1, numout)
	471	END DO
	472	WRITE(numout,*)
	473	WRITE(numout,*) ' dom_msk: vmask for each level'
	474	WRITE(numout,*) ' -----------------------------'
	475	DO jk = 1, jpk
	476	imsk(:,:) = INT( vmask(:,:,jk) )
	477	WRITE(numout,*)
	478	WRITE(numout,*) ' level = ',jk
	479	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	480	& 1, jpj, 1, 1, numout)
	481	END DO
	482	WRITE(numout,*)
	483	WRITE(numout,*) ' dom_msk: fmask for each level'
	484	WRITE(numout,*) ' -----------------------------'
	485	DO jk = 1, jpk
	486	imsk(:,:) = INT( fmask(:,:,jk) )
	487	WRITE(numout,*)
	488	WRITE(numout,*) ' level = ',jk
	489	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
	490	& 1, jpj, 1, 1, numout )
	491	END DO
	492	WRITE(numout,*)
	493	WRITE(numout,*) ' dom_msk: bmask '
	494	WRITE(numout,*) ' ---------------'
	495	WRITE(numout,*)
	496	imsk(:,:) = INT( bmask(:,:) )
	497	CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &
[2528]	498	& 1, jpj, 1, 1, numout )
[3]	499	ENDIF
[1566]	500	!
[3294]	501	CALL wrk_dealloc( jpi, jpj, imsk )
	502	CALL wrk_dealloc( jpi, jpj, zwf )
[2715]	503	!
[3294]	504	IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
	505	!
[3]	506	END SUBROUTINE dom_msk
	507
	508	#if defined key_noslip_accurate
	509	!!----------------------------------------------------------------------
	510	!! 'key_noslip_accurate' : accurate no-slip boundary condition
	511	!!----------------------------------------------------------------------
	512
	513	SUBROUTINE dom_msk_nsa
	514	!!---------------------------------------------------------------------
	515	!! * ROUTINE dom_msk_nsa *
	516	!!
	517	!! ** Purpose :
	518	!!
	519	!! ** Method :
	520	!!
	521	!! ** Action :
	522	!!----------------------------------------------------------------------
[2715]	523	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[1566]	524	INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd
	525	REAL(wp) :: zaa
[3]	526	!!---------------------------------------------------------------------
[3294]	527	!
	528	IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')
	529	!
[2715]	530	IF(lwp) WRITE(numout,*)
	531	IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
	532	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'
	533	IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
[3]	534
	535	! mask for second order calculation of vorticity
	536	! ----------------------------------------------
	537	! noslip boundary condition: fmask=1 at convex corner, store
	538	! index of straight coast meshes ( 'west', refering to a coast,
	539	! means west of the ocean, aso)
	540
	541	DO jk = 1, jpk
	542	DO jl = 1, 4
	543	npcoa(jl,jk) = 0
	544	DO ji = 1, 2*(jpi+jpj)
	545	nicoa(ji,jl,jk) = 0
	546	njcoa(ji,jl,jk) = 0
	547	END DO
	548	END DO
	549	END DO
	550
	551	IF( jperio == 2 ) THEN
	552	WRITE(numout,*) ' '
	553	WRITE(numout,*) ' symetric boundary conditions need special'
	554	WRITE(numout,*) ' treatment not implemented. we stop.'
	555	STOP
	556	ENDIF
	557
	558	! convex corners
	559
	560	DO jk = 1, jpkm1
	561	DO jj = 1, jpjm1
	562	DO ji = 1, jpim1
	563	zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
[32]	564	&+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	565	IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp
[3]	566	END DO
	567	END DO
	568	END DO
	569
	570	! north-south straight coast
	571
	572	DO jk = 1, jpkm1
	573	inw = 0
	574	ine = 0
	575	DO jj = 2, jpjm1
	576	DO ji = 2, jpim1
	577	zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]	578	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	579	inw = inw + 1
	580	nicoa(inw,1,jk) = ji
	581	njcoa(inw,1,jk) = jj
	582	IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj
	583	ENDIF
	584	zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
[2528]	585	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	586	ine = ine + 1
	587	nicoa(ine,2,jk) = ji
	588	njcoa(ine,2,jk) = jj
	589	IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj
	590	ENDIF
	591	END DO
	592	END DO
	593	npcoa(1,jk) = inw
	594	npcoa(2,jk) = ine
	595	END DO
	596
	597	! west-east straight coast
	598
	599	DO jk = 1, jpkm1
	600	ins = 0
	601	inn = 0
	602	DO jj = 2, jpjm1
	603	DO ji =2, jpim1
	604	zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
[2528]	605	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	606	ins = ins + 1
	607	nicoa(ins,3,jk) = ji
	608	njcoa(ins,3,jk) = jj
	609	IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
	610	ENDIF
	611	zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
[2528]	612	IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]	613	inn = inn + 1
	614	nicoa(inn,4,jk) = ji
	615	njcoa(inn,4,jk) = jj
	616	IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
	617	ENDIF
	618	END DO
	619	END DO
	620	npcoa(3,jk) = ins
	621	npcoa(4,jk) = inn
	622	END DO
	623
	624	itest = 2 * ( jpi + jpj )
	625	DO jk = 1, jpk
	626	IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &
	627	npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
[474]	628
	629	WRITE(ctmp1,*) ' level jk = ',jk
	630	WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
	631	WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &
[32]	632	& npcoa(3,jk), npcoa(4,jk)
[474]	633	WRITE(ctmp4,) ' 2(jpi+jpj) = ',itest,'. we stop.'
	634	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]	635	ENDIF
	636	END DO
	637
	638	ierror = 0
	639	iind = 0
	640	ijnd = 0
[2528]	641	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
	642	IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
[3]	643	DO jk = 1, jpk
	644	DO jl = 1, npcoa(1,jk)
	645	IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
	646	ierror = ierror+1
	647	icoord(ierror,1) = nicoa(jl,1,jk)
	648	icoord(ierror,2) = njcoa(jl,1,jk)
	649	icoord(ierror,3) = jk
	650	ENDIF
	651	END DO
	652	DO jl = 1, npcoa(2,jk)
	653	IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
	654	ierror = ierror + 1
	655	icoord(ierror,1) = nicoa(jl,2,jk)
	656	icoord(ierror,2) = njcoa(jl,2,jk)
	657	icoord(ierror,3) = jk
	658	ENDIF
	659	END DO
	660	DO jl = 1, npcoa(3,jk)
	661	IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
	662	ierror = ierror + 1
	663	icoord(ierror,1) = nicoa(jl,3,jk)
	664	icoord(ierror,2) = njcoa(jl,3,jk)
	665	icoord(ierror,3) = jk
	666	ENDIF
	667	END DO
[2528]	668	DO jl = 1, npcoa(4,jk)
[3]	669	IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]	670	ierror=ierror + 1
	671	icoord(ierror,1) = nicoa(jl,4,jk)
	672	icoord(ierror,2) = njcoa(jl,4,jk)
	673	icoord(ierror,3) = jk
[3]	674	ENDIF
	675	END DO
	676	END DO
	677
	678	IF( ierror > 0 ) THEN
	679	IF(lwp) WRITE(numout,*)
	680	IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'
	681	IF(lwp) WRITE(numout,*) ' Bad marking off at points:'
	682	DO jl = 1, ierror
	683	IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &
[32]	684	& ' Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]	685	END DO
[474]	686	CALL ctl_stop( 'We stop...' )
[3]	687	ENDIF
[2715]	688	!
[3294]	689	IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')
	690	!
[3]	691	END SUBROUTINE dom_msk_nsa
	692
	693	#else
	694	!!----------------------------------------------------------------------
	695	!! Default option : Empty routine
	696	!!----------------------------------------------------------------------
	697	SUBROUTINE dom_msk_nsa
	698	END SUBROUTINE dom_msk_nsa
	699	#endif
	700
	701	!!======================================================================
	702	END MODULE dommsk

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

Download in other formats: