Context Navigation

dommsk.F90 @ 11248

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

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

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

Context Navigation

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 11248

Download in other formats: