Context Navigation

dommsk.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 33.4 KB

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

Download in other formats: