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.

dynldf_iso.F90 in NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/DYN/dynldf_iso.F90 @ 14856

Last change on this file since 14856 was 14856, checked in by smueller, 3 years ago
Synchronizing with ^{/NEMO/trunk@14854 (ticket #2353)}
Property svn:keywords set to `Id`
File size: 21.2 KB

Rev	Line
[3]	1	MODULE dynldf_iso
	2	!!======================================================================
	3	!! * MODULE dynldf_iso *
[5836]	4	!! Ocean dynamics: lateral viscosity trend (rotated laplacian operator)
[3]	5	!!======================================================================
[2715]	6	!! History : OPA ! 97-07 (G. Madec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	8	!! - ! 2004-08 (C. Talandier) New trends organization
	9	!! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion
[5836]	10	!! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification,
	11	!! ! add velocity dependent coefficient and optional read in file
[2715]	12	!!----------------------------------------------------------------------
[5836]	13
[3]	14	!!----------------------------------------------------------------------
	15	!! dyn_ldf_iso : update the momentum trend with the horizontal part
	16	!! of the lateral diffusion using isopycnal or horizon-
	17	!! tal s-coordinate laplacian operator.
	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and tracers
	20	USE dom_oce ! ocean space and time domain
[5836]	21	USE ldfdyn ! lateral diffusion: eddy viscosity coef.
	22	USE ldftra ! lateral physics: eddy diffusivity
[3]	23	USE zdf_oce ! ocean vertical physics
	24	USE ldfslp ! iso-neutral slopes
[4990]	25	!
[3]	26	USE in_out_manager ! I/O manager
[2715]	27	USE lib_mpp ! MPP library
[5836]	28	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	29	USE prtctl ! Print control
[14856]	30	#if defined key_loop_fusion
	31	USE dynldf_iso_lf, ONLY: dyn_ldf_iso_lf ! lateral mixing - loop fusion version (dyn_ldf_iso routine )
	32	#endif
[3]	33
	34	IMPLICIT NONE
	35	PRIVATE
	36
[2715]	37	PUBLIC dyn_ldf_iso ! called by step.F90
	38	PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90
[3]	39
[9019]	40	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity
[2715]	41
[3]	42	!! * Substitutions
[12377]	43	# include "do_loop_substitute.h90"
[13237]	44	# include "domzgr_substitute.h90"
[3]	45	!!----------------------------------------------------------------------
[9598]	46	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1156]	47	!! $Id$
[10068]	48	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	49	!!----------------------------------------------------------------------
	50	CONTAINS
	51
[2715]	52	INTEGER FUNCTION dyn_ldf_iso_alloc()
	53	!!----------------------------------------------------------------------
	54	!! * ROUTINE dyn_ldf_iso_alloc *
	55	!!----------------------------------------------------------------------
[14856]	56	dyn_ldf_iso_alloc = 0
	57	IF( .NOT. ALLOCATED( akzu ) ) THEN
	58	ALLOCATE( akzu(jpi,jpj,jpk), akzv(jpi,jpj,jpk), STAT=dyn_ldf_iso_alloc )
	59	!
	60	IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
	61	ENDIF
[2715]	62	END FUNCTION dyn_ldf_iso_alloc
	63
	64
[12377]	65	SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
[3]	66	!!----------------------------------------------------------------------
	67	!! * ROUTINE dyn_ldf_iso *
	68	!!
[455]	69	!! ** Purpose : Compute the before trend of the rotated laplacian
	70	!! operator of lateral momentum diffusion except the diagonal
	71	!! vertical term that will be computed in dynzdf module. Add it
	72	!! to the general trend of momentum equation.
[3]	73	!!
	74	!! ** Method :
[455]	75	!! The momentum lateral diffusive trend is provided by a 2nd
	76	!! order operator rotated along neutral or geopotential surfaces
	77	!! (in s-coordinates).
[3]	78	!! It is computed using before fields (forward in time) and isopyc-
[455]	79	!! nal or geopotential slopes computed in routine ldfslp.
[3]	80	!! Here, u and v components are considered as 2 independent scalar
	81	!! fields. Therefore, the property of splitting divergent and rota-
	82	!! tional part of the flow of the standard, z-coordinate laplacian
	83	!! momentum diffusion is lost.
	84	!! horizontal fluxes associated with the rotated lateral mixing:
	85	!! u-component:
[12377]	86	!! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ uu ]
	87	!! - ahmt e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
	88	!! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ uu ]
	89	!! - ahmf e1f * mi(vslp) dk[ mj(mk(uu)) ]
[3]	90	!! v-component:
[12377]	91	!! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vv ]
	92	!! - ahmf e2t * mj(uslp) dk[ mi(mk(vv)) ]
	93	!! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ vv ]
	94	!! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
[3]	95	!! take the horizontal divergence of the fluxes:
	96	!! diffu = 1/(e1ue2ue3u) { di [ ziut ] + dj-1[ zjuf ] }
	97	!! diffv = 1/(e1ve2ve3v) { di-1[ zivf ] + dj [ zjvt ] }
[12377]	98	!! Add this trend to the general trend (uu(rhs),vv(rhs)):
	99	!! uu(rhs) = uu(rhs) + diffu
[455]	100	!! CAUTION: here the isopycnal part is with a coeff. of aht. This
	101	!! should be modified for applications others than orca_r2 (!!bug)
[3]	102	!!
	103	!! ** Action :
[12377]	104	!! -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
[9019]	105	!! -(akzu,akzv) to accompt for the diagonal vertical component
	106	!! of the rotated operator in dynzdf module
[3]	107	!!----------------------------------------------------------------------
[12377]	108	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	109	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
	110	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
[2715]	111	!
	112	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	113	REAL(wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj ! local scalars
	114	REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - -
[9490]	115	REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0 ! - -
[14856]	116	REAL(wp), DIMENSION(A2D(nn_hls)) :: ziut, zivf, zdku, zdk1u ! 2D workspace
	117	REAL(wp), DIMENSION(A2D(nn_hls)) :: zjuf, zjvt, zdkv, zdk1v ! - -
	118	REAL(wp), DIMENSION(A1Di(nn_hls),jpk) :: zfuw, zdiu, zdju, zdj1u ! - -
	119	REAL(wp), DIMENSION(A1Di(nn_hls),jpk) :: zfvw, zdiv, zdjv, zdj1v ! - -
[2715]	120	!!----------------------------------------------------------------------
[3294]	121	!
[14856]	122	#if defined key_loop_fusion
	123	CALL dyn_ldf_iso_lf( kt, Kbb, Kmm, puu, pvv, Krhs )
	124	#else
	125
	126	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
	127	IF( kt == nit000 ) THEN
	128	IF(lwp) WRITE(numout,*)
	129	IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
	130	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
	131	! ! allocate dyn_ldf_iso arrays
	132	IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
	133	ENDIF
[3]	134	ENDIF
[216]	135
[5836]	136	!!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED
	137	! s-coordinate: Iso-level diffusion on momentum but not on tracer
[455]	138	IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
[2715]	139	!
[14856]	140	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) ! set the slopes of iso-level
[12377]	141	uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
	142	vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
	143	wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
	144	wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
	145	END_3D
[455]	146	! Lateral boundary conditions on the slopes
[14856]	147	IF (nn_hls == 1) CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
[9019]	148	!
[14856]	149	ENDIF
[9490]	150
	151	zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max
	152
[3]	153	! ! ===============
	154	DO jk = 1, jpkm1 ! Horizontal slab
	155	! ! ===============
	156
	157	! Vertical u- and v-shears at level jk and jk+1
	158	! ---------------------------------------------
	159	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	160	! zdkv(jk=1)=zdkv(jk=2)
	161
[14856]	162	DO_2D( 1, 1, 1, 1 )
	163	zdk1u(ji,jj) = ( puu(ji,jj,jk,Kbb) -puu(ji,jj,jk+1,Kbb) ) * umask(ji,jj,jk+1)
	164	zdk1v(ji,jj) = ( pvv(ji,jj,jk,Kbb) -pvv(ji,jj,jk+1,Kbb) ) * vmask(ji,jj,jk+1)
	165	END_2D
[3]	166
	167	IF( jk == 1 ) THEN
	168	zdku(:,:) = zdk1u(:,:)
	169	zdkv(:,:) = zdk1v(:,:)
	170	ELSE
[14856]	171	DO_2D( 1, 1, 1, 1 )
	172	zdku(ji,jj) = ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
	173	zdkv(ji,jj) = ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
	174	END_2D
[3]	175	ENDIF
	176
	177	! -----f-----
	178	! Horizontal fluxes on U \|
	179	! --------------------=== t u t
	180	! \|
	181	! i-flux at t-point -----f-----
	182
[455]	183	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
[14822]	184	DO_2D( 0, 1, 0, 0 )
[13237]	185	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) &
	186	& * MIN( e3u(ji ,jj,jk,Kmm), &
	187	& e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
[3]	188
[12377]	189	zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	190	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	191
[12377]	192	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	193
	194	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
	195	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	196	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	197	END_2D
[455]	198	ELSE ! other coordinate system (zco or sco) : e3t
[14822]	199	DO_2D( 0, 1, 0, 0 )
[13237]	200	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
	201	& * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
[3]	202
[12377]	203	zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) &
	204	& + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp )
[455]	205
[12377]	206	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	207
[12377]	208	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
	209	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	210	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	211	END_2D
[455]	212	ENDIF
[3]	213
	214	! j-flux at f-point
[13295]	215	DO_2D( 1, 0, 1, 0 )
[13237]	216	zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
	217	& * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
[3]	218
[12377]	219	zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	220	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	221
[12377]	222	zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]	223
[12377]	224	zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) &
	225	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	226	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
	227	END_2D
[3]	228
	229	! \| t \|
	230	! Horizontal fluxes on V \| \|
	231	! --------------------=== f---v---f
	232	! \| \|
	233	! i-flux at f-point \| t \|
	234
[14822]	235	DO_2D( 1, 0, 0, 0 )
[13237]	236	zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
	237	& * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
[3]	238
[12377]	239	zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	240	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	241
[12377]	242	zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]	243
[12377]	244	zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) &
	245	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
	246	& + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
	247	END_2D
[3]	248
	249	! j-flux at t-point
[455]	250	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
[14822]	251	DO_2D( 1, 0, 0, 1 )
[13237]	252	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) &
	253	& * MIN( e3v(ji,jj ,jk,Kmm), &
	254	& e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
[3]	255
[12377]	256	zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	257	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	258
[12377]	259	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]	260
[12377]	261	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
	262	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	263	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	264	END_2D
[455]	265	ELSE ! other coordinate system (zco or sco) : e3t
[14822]	266	DO_2D( 1, 0, 0, 1 )
[13237]	267	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
	268	& * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
[3]	269
[12377]	270	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	271	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	272
[12377]	273	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[455]	274
[12377]	275	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
	276	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	277	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	278	END_2D
[455]	279	ENDIF
	280
	281
[3]	282	! Second derivative (divergence) and add to the general trend
	283	! -----------------------------------------------------------
[13497]	284	DO_2D( 0, 0, 0, 0 ) !!gm Question vectop possible??? !!bug
[12377]	285	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
[13237]	286	& + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) &
	287	& / e3u(ji,jj,jk,Kmm)
[12377]	288	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
[13237]	289	& + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) &
	290	& / e3v(ji,jj,jk,Kmm)
[12377]	291	END_2D
[3]	292	! ! ===============
	293	END DO ! End of slab
	294	! ! ===============
[216]	295
[455]	296	! print sum trends (used for debugging)
[12377]	297	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
	298	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]	299
	300
[455]	301	! ! ===============
[14856]	302	DO jj = ntsj, ntej ! Vertical slab
[455]	303	! ! ===============
	304
	305
	306	! I. vertical trends associated with the lateral mixing
	307	! =====================================================
	308	! (excluding the vertical flux proportional to dk[t]
	309
	310
	311	! I.1 horizontal momentum gradient
	312	! --------------------------------
	313
	314	DO jk = 1, jpk
[14856]	315	DO ji = ntsi, ntei + nn_hls
[455]	316	! i-gradient of u at jj
[12377]	317	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) )
[455]	318	! j-gradient of u and v at jj
[12377]	319	zdju (ji,jk) = fmask(ji,jj ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) )
	320	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pvv(ji,jj ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) )
[455]	321	! j-gradient of u and v at jj+1
[12377]	322	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) )
	323	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	324	END DO
	325	END DO
	326	DO jk = 1, jpk
[14856]	327	DO ji = ntsi - nn_hls, ntei
[455]	328	! i-gradient of v at jj
[12377]	329	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	330	END DO
	331	END DO
	332
	333
	334	! I.2 Vertical fluxes
	335	! -------------------
	336
	337	! Surface and bottom vertical fluxes set to zero
[14856]	338	DO ji = ntsi - nn_hls, ntei + nn_hls
[455]	339	zfuw(ji, 1 ) = 0.e0
	340	zfvw(ji, 1 ) = 0.e0
	341	zfuw(ji,jpk) = 0.e0
	342	zfvw(ji,jpk) = 0.e0
	343	END DO
	344
	345	! interior (2=<jk=<jpk-1) on U field
	346	DO jk = 2, jpkm1
[14856]	347	DO ji = ntsi, ntei
[9490]	348	zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
[5836]	349	!
[9019]	350	zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	351	zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
[5836]	352	!
[9019]	353	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	354	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. )
	355	zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) &
	356	+ fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. )
[455]	357
[9019]	358	zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
	359	zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
[455]	360	! vertical flux on u field
[9019]	361	zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	362	& + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	363	& + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	364	& + zdj1u(ji,jk ) + zdju (ji ,jk ) )
	365	! vertical mixing coefficient (akzu)
[9490]	366	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	367	akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
[455]	368	END DO
	369	END DO
	370
	371	! interior (2=<jk=<jpk-1) on V field
	372	DO jk = 2, jpkm1
[14856]	373	DO ji = ntsi, ntei
[9490]	374	zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
[9019]	375	!
	376	zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	377	zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	378	!
	379	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	380	& + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ) , 1. )
	381	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	382	& + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ) , 1. )
[455]	383
[9019]	384	zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
	385	zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
[455]	386	! vertical flux on v field
[9019]	387	zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	388	& + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	389	& + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	390	& + zdjv (ji,jk ) + zdj1v(ji ,jk ) )
	391	! vertical mixing coefficient (akzv)
[9490]	392	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	393	akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
[455]	394	END DO
	395	END DO
	396
	397
	398	! I.3 Divergence of vertical fluxes added to the general tracer trend
	399	! -------------------------------------------------------------------
	400	DO jk = 1, jpkm1
[14856]	401	DO ji = ntsi, ntei
[13237]	402	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) &
	403	& / e3u(ji,jj,jk,Kmm)
	404	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) &
	405	& / e3v(ji,jj,jk,Kmm)
[455]	406	END DO
	407	END DO
	408	! ! ===============
	409	END DO ! End of slab
	410	! ! ===============
[14856]	411	#endif
[3]	412	END SUBROUTINE dyn_ldf_iso
	413
	414	!!======================================================================
	415	END MODULE dynldf_iso

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

Download in other formats: