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.

ldfslp.F90 in trunk/NEMO/OPA_SRC/LDF – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 1793

Last change on this file since 1793 was 1515, checked in by ctlod, 15 years ago
wrong sign in the lbc related to pn2 in the computation of the slope just below the mixed layer (step II), see ticket: #435
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 31.9 KB

Rev	Line
[3]	1	MODULE ldfslp
	2	!!======================================================================
	3	!! * MODULE ldfslp *
	4	!! Ocean physics: slopes of neutral surfaces
	5	!!======================================================================
[1515]	6	!! History : OPA ! 1994-12 (G. Madec, M. Imbard) Original code
	7	!! 8.0 ! 1997-06 (G. Madec) optimization, lbc
	8	!! 8.1 ! 1999-10 (A. Jouzeau) NEW profile in the mixed layer
	9	!! NEMO 0.5 ! 2002-10 (G. Madec) Free form, F90
	10	!! 1.0 ! 2005-10 (A. Beckmann) correction for s-coordinates
	11	!!----------------------------------------------------------------------
[3]	12	#if defined key_ldfslp \|\| defined key_esopa
	13	!!----------------------------------------------------------------------
	14	!! 'key_ldfslp' Rotation of lateral mixing tensor
	15	!!----------------------------------------------------------------------
	16	!! ldf_slp : compute the slopes of neutral surface
	17	!! ldf_slp_mxl : compute the slopes of iso-neutral surface
	18	!! ldf_slp_init : initialization of the slopes computation
	19	!!----------------------------------------------------------------------
	20	USE oce ! ocean dynamics and tracers
	21	USE dom_oce ! ocean space and time domain
	22	USE ldftra_oce
	23	USE ldfdyn_oce
	24	USE phycst ! physical constants
	25	USE zdfmxl ! mixed layer depth
	26	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	27	USE in_out_manager ! I/O manager
[258]	28	USE prtctl ! Print control
[3]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[1515]	33	PUBLIC ldf_slp ! routine called by step.F90
[3]	34
[1515]	35	LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag
	36	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: uslp, wslpi !: i_slope at U- and W-points
	37	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: vslp, wslpj !: j-slope at V- and W-points
[3]	38
[1515]	39	REAL(wp), DIMENSION(jpi,jpj,jpk) :: omlmask ! mask of the surface mixed layer at T-pt
	40	REAL(wp), DIMENSION(jpi,jpj) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer
	41	REAL(wp), DIMENSION(jpi,jpj) :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer
[3]	42
	43	!! * Substitutions
	44	# include "domzgr_substitute.h90"
	45	# include "vectopt_loop_substitute.h90"
	46	!!----------------------------------------------------------------------
[1515]	47	!! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
[1156]	48	!! $Id$
	49	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	50	!!----------------------------------------------------------------------
	51
	52	CONTAINS
	53
	54	SUBROUTINE ldf_slp( kt, prd, pn2 )
	55	!!----------------------------------------------------------------------
	56	!! * ROUTINE ldf_slp *
	57	!!
[1515]	58	!! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal
	59	!! surfaces referenced locally) ('key_traldfiso').
	60	!!
[3]	61	!! ** Method : The slope in the i-direction is computed at U- and
	62	!! W-points (uslp, wslpi) and the slope in the j-direction is
	63	!! computed at V- and W-points (vslp, wslpj).
	64	!! They are bounded by 1/100 over the whole ocean, and within the
	65	!! surface layer they are bounded by the distance to the surface
	66	!! ( slope<= depth/l where l is the length scale of horizontal
	67	!! diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
	68	!! of 10cm/s)
	69	!! A horizontal shapiro filter is applied to the slopes
[461]	70	!! ln_sco=T, s-coordinate, add to the previously computed slopes
[3]	71	!! the slope of the model level surface.
	72	!! macro-tasked on horizontal slab (jk-loop) (2, jpk-1)
	73	!! [slopes already set to zero at level 1, and to zero or the ocean
[461]	74	!! bottom slope (ln_sco=T) at level jpk in inildf]
[3]	75	!!
	76	!! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes
	77	!! of now neutral surfaces at u-, w- and v- w-points, resp.
[1515]	78	!!----------------------------------------------------------------------
	79	USE oce , zgru => ua ! use ua as workspace
	80	USE oce , zgrv => va ! use va as workspace
	81	USE oce , zwy => ta ! use ta as workspace
	82	USE oce , zwz => sa ! use sa as workspace
[3]	83	!!
[1515]	84	INTEGER , INTENT(in) :: kt ! ocean time-step index
	85	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: prd ! in situ density
	86	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pn2 ! Brunt-Vaisala frequency (locally ref.)
	87	!!
	88	INTEGER :: ji , jj , jk ! dummy loop indices
	89	INTEGER :: ii0, ii1, iku ! temporary integer
	90	INTEGER :: ij0, ij1, ikv ! temporary integer
	91	REAL(wp) :: zeps, zmg, zm05g, zalpha ! temporary scalars
	92	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! - -
	93	REAL(wp) :: zcofu , zcofv , zcofw ! - -
	94	REAL(wp) :: zau, zbu, zai, zbi, z1u, z1wu ! - -
	95	REAL(wp) :: zav, zbv, zaj, zbj, z1v, z1wv !
	96	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zww ! 3D workspace
[3]	97	!!----------------------------------------------------------------------
	98
[1515]	99	IF( kt == nit000 ) CALL ldf_slp_init ! initialization (first time-step only)
[3]	100
[1515]	101	zeps = 1.e-20 ! Local constant initialization
[32]	102	zmg = -1.0 / grav
	103	zm05g = -0.5 / grav
[1515]	104	!
[3]	105	zww(:,:,:) = 0.e0
	106	zwz(:,:,:) = 0.e0
[1515]	107	! ! horizontal density gradient computation
[3]	108	DO jk = 1, jpk
	109	DO jj = 1, jpjm1
	110	DO ji = 1, fs_jpim1 ! vector opt.
	111	zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) )
	112	zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) )
	113	END DO
	114	END DO
	115	END DO
[1515]	116	IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level
[789]	117	# if defined key_vectopt_loop
[1515]	118	DO jj = 1, 1
	119	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[3]	120	# else
[461]	121	DO jj = 1, jpjm1
	122	DO ji = 1, jpim1
[3]	123	# endif
[1515]	124	iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1 ! last ocean level
[461]	125	ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
	126	zgru(ji,jj,iku) = gru(ji,jj)
	127	zgrv(ji,jj,ikv) = grv(ji,jj)
	128	END DO
[3]	129	END DO
[461]	130	ENDIF
[1515]	131
	132	CALL ldf_slp_mxl( prd, pn2 ) ! Slopes of isopycnal surfaces just below the mixed layer
[3]	133
	134
[1515]	135	! I. slopes at u and v point \| uslp = d/di( prd ) / d/dz( prd )
	136	! =========================== \| vslp = d/dj( prd ) / d/dz( prd )
	137	!
	138	! !* Local vertical density gradient evaluated from N^2
	139	DO jk = 2, jpkm1 ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
[3]	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
[1515]	142	zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * ( pn2 (ji,jj,jk) + pn2 (ji,jj,jk+1) ) &
	143	& / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. )
[3]	144	END DO
	145	END DO
[1515]	146	END DO
	147	DO jk = 2, jpkm1 !* Slopes at u and v points
[3]	148	DO jj = 2, jpjm1
	149	DO ji = fs_2, fs_jpim1 ! vector opt.
	150	! horizontal and vertical density gradient at u- and v-points
	151	zau = 1. / e1u(ji,jj) * zgru(ji,jj,jk)
	152	zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk)
	153	zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj ,jk) )
	154	zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(ji ,jj+1,jk) )
	155	! bound the slopes: abs(zw.)<= 1/100 and zb..<0
	156	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	157	zbu = MIN( zbu, -100.ABS( zau ), -7.e+3/fse3u(ji,jj,jk)ABS( zau ) )
	158	zbv = MIN( zbv, -100.ABS( zav ), -7.e+3/fse3v(ji,jj,jk)ABS( zav ) )
	159	! uslp and vslp output in zwz and zww, resp.
	160	zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
[1515]	161	zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps ) &
	162	& + zalpha * uslpml(ji,jj) &
	163	& * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) &
	164	& / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk)
[3]	165	zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
[1515]	166	zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps ) &
	167	& + zalpha * vslpml(ji,jj) &
	168	& * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) &
	169	& / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
[3]	170	END DO
	171	END DO
[1515]	172	END DO
	173	CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions
	174	!
	175	zcofu = 1. / 16. !* horizontal Shapiro filter
	176	zcofv = 1. / 16.
	177	DO jk = 2, jpkm1
	178	DO jj = 2, jpjm1, jpj-3 ! rows jj=2 and =jpjm1 only
[3]	179	DO ji = 2, jpim1
[1515]	180	uslp(ji,jj,jk) = zcofu * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	181	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	182	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	183	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	184	& + 4.* zwz(ji ,jj ,jk) )
	185	vslp(ji,jj,jk) = zcofv * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	186	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	187	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	188	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	189	& + 4.* zww(ji,jj ,jk) )
[3]	190	END DO
	191	END DO
[1515]	192	DO jj = 3, jpj-2 ! other rows
[3]	193	DO ji = fs_2, fs_jpim1 ! vector opt.
	194	uslp(ji,jj,jk) = zcofu * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	195	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	196	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	197	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	198	& + 4.* zwz(ji ,jj ,jk) )
	199	vslp(ji,jj,jk) = zcofv * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	200	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	201	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	202	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	203	& + 4.* zww(ji,jj ,jk) )
	204	END DO
	205	END DO
[1515]	206	! !* decrease along coastal boundaries
[3]	207	DO jj = 2, jpjm1
	208	DO ji = fs_2, fs_jpim1 ! vector opt.
	209	z1u = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5
	210	z1v = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5
	211	z1wu = ( umask(ji,jj,jk) + umask(ji,jj,jk+1) )*.5
	212	z1wv = ( vmask(ji,jj,jk) + vmask(ji,jj,jk+1) )*.5
	213	uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu
	214	vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv
	215	END DO
	216	END DO
[1515]	217	END DO
[3]	218
	219
[1515]	220	! II. slopes at w point \| wslpi = mij( d/di( prd ) / d/dz( prd )
	221	! =========================== \| wslpj = mij( d/dj( prd ) / d/dz( prd )
	222	!
	223	! !* Local vertical density gradient evaluated from N^2
	224	DO jk = 2, jpkm1 ! zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
[3]	225	DO jj = 1, jpj
	226	DO ji = 1, jpi
[1515]	227	zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
[3]	228	END DO
	229	END DO
[1515]	230	END DO
	231	DO jk = 2, jpkm1 !* Slopes at w point
[3]	232	DO jj = 2, jpjm1
	233	DO ji = fs_2, fs_jpim1 ! vector opt.
[1515]	234	! ! horizontal density i-gradient at w-points
[3]	235	zcoef1 = MAX( zeps, umask(ji-1,jj,jk )+umask(ji,jj,jk ) &
	236	& +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) )
	237	zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
	238	zai = zcoef1 * ( zgru(ji ,jj,jk ) + zgru(ji ,jj,jk-1) &
	239	& + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk ) ) * tmask (ji,jj,jk)
[1515]	240	! ! horizontal density j-gradient at w-points
[3]	241	zcoef2 = MAX( zeps, vmask(ji,jj-1,jk )+vmask(ji,jj,jk-1) &
	242	& +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk ) )
	243	zcoef2 = 1.0 / ( zcoef2 * e2t (ji,jj) )
	244	zaj = zcoef2 * ( zgrv(ji,jj ,jk ) + zgrv(ji,jj ,jk-1) &
	245	& + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk ) ) * tmask (ji,jj,jk)
[1515]	246	! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
	247	! ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
[3]	248	zbi = MIN( zwy (ji,jj,jk),- 100.ABS(zai), -7.e+3/fse3w(ji,jj,jk)ABS(zai) )
	249	zbj = MIN( zwy (ji,jj,jk), -100.ABS(zaj), -7.e+3/fse3w(ji,jj,jk)ABS(zaj) )
[1515]	250	! ! wslpi and wslpj output in zwz and zww, resp.
[461]	251	zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )
	252	zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
	253	zwz(ji,jj,jk) = ( zai / ( zbi - zeps) * ( 1. - zalpha ) &
	254	& + zcoef3 * wslpiml(ji,jj) * zalpha ) * tmask (ji,jj,jk)
	255	zww(ji,jj,jk) = ( zaj / ( zbj - zeps) * ( 1. - zalpha ) &
	256	& + zcoef3 * wslpjml(ji,jj) * zalpha ) * tmask (ji,jj,jk)
[3]	257	END DO
	258	END DO
[1515]	259	END DO
	260	CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions on zwz and zww
	261	!
	262	! !* horizontal Shapiro filter
	263	DO jk = 2, jpkm1
	264	DO jj = 2, jpjm1, jpj-3 ! rows jj=2 and =jpjm1
[3]	265	DO ji = 2, jpim1
[1515]	266	zcofw = tmask(ji,jj,jk) / 16.
[3]	267	wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	268	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	269	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	270	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	271	& + 4.* zwz(ji ,jj ,jk) ) * zcofw
	272
	273	wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	274	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	275	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	276	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	277	& + 4.* zww(ji ,jj ,jk) ) * zcofw
	278	END DO
[1515]	279	END DO
	280	DO jj = 3, jpj-2 ! other rows
[3]	281	DO ji = fs_2, fs_jpim1 ! vector opt.
[1515]	282	zcofw = tmask(ji,jj,jk) / 16.
[3]	283	wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	284	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	285	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	286	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	287	& + 4.* zwz(ji ,jj ,jk) ) * zcofw
	288
	289	wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	290	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	291	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	292	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	293	& + 4.* zww(ji ,jj ,jk) ) * zcofw
	294	END DO
	295	END DO
[1515]	296	! !* decrease along coastal boundaries
[3]	297	DO jj = 2, jpjm1
	298	DO ji = fs_2, fs_jpim1 ! vector opt.
	299	z1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5
	300	z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5
	301	wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z1u * z1v
	302	wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z1u * z1v
	303	END DO
	304	END DO
[1515]	305	END DO
[3]	306
	307
[1515]	308	! III. Specific grid points
	309	! ===========================
	310	!
	311	IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area
	312	! ! Gibraltar Strait
	313	ij0 = 50 ; ij1 = 53
	314	ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	315	ij0 = 51 ; ij1 = 53
	316	ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	317	ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	318	ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	319	!
	320	! ! Mediterrannean Sea
	321	ij0 = 49 ; ij1 = 56
	322	ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	323	ij0 = 50 ; ij1 = 56
	324	ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	325	ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	326	ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
	327	ENDIF
[3]	328
	329
[1515]	330	! IV. Lateral boundary conditions
	331	! ===============================
[461]	332	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )
	333	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
[3]	334
[1515]	335
[258]	336	IF(ln_ctl) THEN
	337	CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
	338	CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
[49]	339	ENDIF
[1515]	340	!
[3]	341	END SUBROUTINE ldf_slp
	342
	343
	344	SUBROUTINE ldf_slp_mxl( prd, pn2 )
	345	!!----------------------------------------------------------------------
	346	!! * ROUTINE ldf_slp_mxl *
	347	!!
[1515]	348	!! ** Purpose : Compute the slopes of iso-neutral surface just below
	349	!! the mixed layer.
	350	!!
[3]	351	!! ** Method :
	352	!! The slope in the i-direction is computed at u- and w-points
	353	!! (uslp, wslpi) and the slope in the j-direction is computed at
	354	!! v- and w-points (vslp, wslpj).
	355	!! They are bounded by 1/100 over the whole ocean, and within the
	356	!! surface layer they are bounded by the distance to the surface
	357	!! ( slope<= depth/l where l is the length scale of horizontal
	358	!! diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
	359	!! of 10cm/s)
	360	!!
[1515]	361	!! ** Action : Compute uslp, wslpi, and vslp, wslpj, the i- and j-slopes
[3]	362	!! of now neutral surfaces at u-, w- and v- w-points, resp.
[1515]	363	!!----------------------------------------------------------------------
	364	USE oce , zgru => ua ! ua, va used as workspace and set to hor.
	365	USE oce , zgrv => va ! density gradient in ldf_slp
	366	!!
	367	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: prd ! in situ density
	368	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pn2 ! Brunt-Vaisala frequency (locally ref.)
[3]	369	!!
[1515]	370	INTEGER :: ji, jj, jk ! dummy loop indices
	371	INTEGER :: ik, ikm1 ! temporary integers
	372	REAL(wp) :: zeps, zmg, zm05g ! temporary scalars
	373	REAL(wp) :: zcoef1, zcoef2 ! - -
	374	REAL(wp) :: zau, zbu, zai, zbi ! - -
	375	REAL(wp) :: zav, zbv, zaj, zbj ! - -
	376	REAL(wp), DIMENSION(jpi,jpj) :: zwy ! 2D workspace
[3]	377	!!----------------------------------------------------------------------
	378
[1515]	379	zeps = 1.e-20 ! Local constant initialization
[32]	380	zmg = -1.0 / grav
	381	zm05g = -0.5 / grav
[1515]	382	!
[3]	383	uslpml (1,:) = 0.e0 ; uslpml (jpi,:) = 0.e0
	384	vslpml (1,:) = 0.e0 ; vslpml (jpi,:) = 0.e0
	385	wslpiml(1,:) = 0.e0 ; wslpiml(jpi,:) = 0.e0
	386	wslpjml(1,:) = 0.e0 ; wslpjml(jpi,:) = 0.e0
	387
[1515]	388	! ! surface mixed layer mask
	389	DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise
[789]	390	# if defined key_vectopt_loop
[1515]	391	DO jj = 1, 1
	392	DO ji = 1, jpij ! vector opt. (forced unrolling)
[3]	393	# else
	394	DO jj = 1, jpj
	395	DO ji = 1, jpi
	396	# endif
	397	ik = nmln(ji,jj) - 1
[1515]	398	IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1.e0
	399	ELSE ; omlmask(ji,jj,jk) = 0.e0
[3]	400	ENDIF
	401	END DO
	402	END DO
	403	END DO
	404
	405
	406	! Slopes of isopycnal surfaces just before bottom of mixed layer
	407	! --------------------------------------------------------------
[1515]	408	! The slope are computed as in the 3D case.
	409	! A key point here is the definition of the mixed layer at u- and v-points.
	410	! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
	411	! Otherwise, a n2 value inside the mixed layer can be involved in the computation
	412	! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
	413	! induce velocity field near the base of the mixed layer.
[3]	414	!-----------------------------------------------------------------------
[1515]	415	!
	416	zwy(:,jpj) = 0.e0 !* vertical density gradient for u-slope (from N^2)
[3]	417	zwy(jpi,:) = 0.e0
[789]	418	# if defined key_vectopt_loop
[1515]	419	DO jj = 1, 1
	420	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[3]	421	# else
	422	DO jj = 1, jpjm1
	423	DO ji = 1, jpim1
	424	# endif
[1515]	425	ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) ! avoid spurious recirculation
	426	ik = MIN( ik, jpkm1 ) ! if ik = jpk take jpkm1 values
	427	zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * ( pn2 (ji,jj,ik) + pn2 (ji,jj,ik+1) ) &
	428	& / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
[3]	429	END DO
	430	END DO
[1515]	431	CALL lbc_lnk( zwy, 'U', 1. ) ! lateral boundary conditions NO sign change
[3]	432
[1515]	433	! !* Slope at u points
[789]	434	# if defined key_vectopt_loop
[1515]	435	DO jj = 1, 1
	436	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3]	437	# else
	438	DO jj = 2, jpjm1
	439	DO ji = 2, jpim1
	440	# endif
	441	! horizontal and vertical density gradient at u-points
	442	ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
[1515]	443	ik = MIN( ik, jpkm1 )
[3]	444	zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
	445	zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
	446	! bound the slopes: abs(zw.)<= 1/100 and zb..<0
	447	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	448	zbu = MIN( zbu, -100.ABS(zau), -7.e+3/fse3u(ji,jj,ik)ABS(zau) )
	449	! uslpml
	450	uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
	451	END DO
	452	END DO
[1515]	453	CALL lbc_lnk( uslpml, 'U', -1. ) ! lateral boundary conditions (i-gradient => sign change)
[3]	454
[1515]	455	! !* vertical density gradient for v-slope (from N^2)
[789]	456	# if defined key_vectopt_loop
[1515]	457	DO jj = 1, 1
	458	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[3]	459	# else
	460	DO jj = 1, jpjm1
	461	DO ji = 1, jpim1
	462	# endif
	463	ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
[1515]	464	ik = MIN( ik, jpkm1 )
	465	zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * ( pn2 (ji,jj,ik) + pn2 (ji,jj,ik+1) ) &
	466	& / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
[3]	467	END DO
	468	END DO
[1515]	469	CALL lbc_lnk( zwy, 'V', 1. ) ! lateral boundary conditions NO sign change
[3]	470
[1515]	471	! !* Slope at v points
[789]	472	# if defined key_vectopt_loop
[1515]	473	DO jj = 1, 1
	474	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3]	475	# else
	476	DO jj = 2, jpjm1
	477	DO ji = 2, jpim1
	478	# endif
	479	! horizontal and vertical density gradient at v-points
	480	ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
	481	ik = MIN( ik,jpkm1 )
	482	zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
	483	zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
	484	! bound the slopes: abs(zw.)<= 1/100 and zb..<0
	485	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	486	zbv = MIN( zbv, -100.ABS(zav), -7.e+3/fse3v(ji,jj,ik)ABS( zav ) )
	487	! vslpml
	488	vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
	489	END DO
	490	END DO
[1515]	491	CALL lbc_lnk( vslpml, 'V', -1. ) ! lateral boundary conditions (j-gradient => sign change)
[3]	492
	493
[1515]	494	! !* vertical density gradient for w-slope (from N^2)
[789]	495	# if defined key_vectopt_loop
[1515]	496	DO jj = 1, 1
	497	DO ji = 1, jpij ! vector opt. (forced unrolling)
[3]	498	# else
	499	DO jj = 1, jpj
	500	DO ji = 1, jpi
	501	# endif
[1515]	502	ik = nmln(ji,jj) + 1
	503	ik = MIN( ik, jpk )
[3]	504	ikm1 = MAX ( 1, ik-1)
	505	zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) * &
	506	& ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
	507	END DO
	508	END DO
	509
[1515]	510	! !* Slopes at w points
[789]	511	# if defined key_vectopt_loop
[1515]	512	DO jj = 1, 1
	513	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3]	514	# else
	515	DO jj = 2, jpjm1
	516	DO ji = 2, jpim1
	517	# endif
[1515]	518	ik = nmln(ji,jj) + 1
	519	ik = MIN( ik, jpk )
	520	ikm1 = MAX ( 1, ik-1 )
[3]	521	! horizontal density i-gradient at w-points
	522	zcoef1 = MAX( zeps, umask(ji-1,jj,ik )+umask(ji,jj,ik ) &
	523	& +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) )
	524	zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
	525	zai = zcoef1 * ( zgru(ji ,jj,ik ) + zgru(ji ,jj,ikm1) &
	526	& + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik ) ) * tmask (ji,jj,ik)
	527	! horizontal density j-gradient at w-points
	528	zcoef2 = MAX( zeps, vmask(ji,jj-1,ik )+vmask(ji,jj,ikm1) &
	529	& +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik ) )
	530	zcoef2 = 1.0 / ( zcoef2 * e2t (ji,jj) )
	531	zaj = zcoef2 * ( zgrv(ji,jj ,ik ) + zgrv(ji,jj ,ikm1) &
	532	& + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik ) ) * tmask (ji,jj,ik)
	533	! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
	534	! static instability:
	535	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	536	zbi = MIN ( zwy (ji,jj),- 100.ABS(zai), -7.e+3/fse3w(ji,jj,ik)ABS(zai) )
	537	zbj = MIN ( zwy (ji,jj), -100.ABS(zaj), -7.e+3/fse3w(ji,jj,ik)ABS(zaj) )
	538	! wslpiml and wslpjml
	539	wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
	540	wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
	541	END DO
	542	END DO
[1515]	543	CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions
	544	!
[3]	545	END SUBROUTINE ldf_slp_mxl
	546
	547
	548	SUBROUTINE ldf_slp_init
	549	!!----------------------------------------------------------------------
	550	!! * ROUTINE ldf_slp_init *
	551	!!
	552	!! ** Purpose : Initialization for the isopycnal slopes computation
	553	!!
	554	!! ** Method : read the nammbf namelist and check the parameter
	555	!! values called by tra_dmp at the first timestep (nit000)
	556	!!----------------------------------------------------------------------
	557	INTEGER :: ji, jj, jk ! dummy loop indices
	558	!!----------------------------------------------------------------------
	559
[1515]	560	IF(lwp) THEN
[3]	561	WRITE(numout,*)
	562	WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
	563	WRITE(numout,*) '~~~~~~~'
	564	ENDIF
	565
	566	! Direction of lateral diffusion (tracers and/or momentum)
	567	! ------------------------------
	568	! set the slope to zero (even in s-coordinates)
	569
	570	uslp (:,:,:) = 0.e0
	571	vslp (:,:,:) = 0.e0
	572	wslpi(:,:,:) = 0.e0
	573	wslpj(:,:,:) = 0.e0
	574
	575	uslpml (:,:) = 0.e0
	576	vslpml (:,:) = 0.e0
	577	wslpiml(:,:) = 0.e0
	578	wslpjml(:,:) = 0.e0
	579
[592]	580	IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
[461]	581	IF(lwp) THEN
	582	WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
	583	ENDIF
[3]	584
	585	! geopotential diffusion in s-coordinates on tracers and/or momentum
	586	! The slopes of s-surfaces are computed once (no call to ldfslp in step)
	587	! The slopes for momentum diffusion are i- or j- averaged of those on tracers
	588
	589	! set the slope of diffusion to the slope of s-surfaces
	590	! ( c a u t i o n : minus sign as fsdep has positive value )
	591	DO jk = 1, jpk
	592	DO jj = 2, jpjm1
	593	DO ji = fs_2, fs_jpim1 ! vector opt.
[461]	594	uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
	595	vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
	596	wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
	597	wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
[3]	598	END DO
	599	END DO
	600	END DO
	601	! Lateral boundary conditions on the slopes
[461]	602	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )
	603	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
[3]	604	ENDIF
[1515]	605	!
[3]	606	END SUBROUTINE ldf_slp_init
	607
	608	#else
	609	!!------------------------------------------------------------------------
	610	!! Dummy module : NO Rotation of lateral mixing tensor
	611	!!------------------------------------------------------------------------
[32]	612	LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag
[3]	613	CONTAINS
	614	SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine
	615	INTEGER, INTENT(in) :: kt
[1515]	616	REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
[32]	617	WRITE(,) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
[3]	618	END SUBROUTINE ldf_slp
	619	#endif
	620
	621	!!======================================================================
	622	END MODULE ldfslp

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

Download in other formats: