Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynldf_iso.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 20.6 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
[3294]	30	USE timing ! Timing
[3]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[2715]	35	PUBLIC dyn_ldf_iso ! called by step.F90
	36	PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90
[3]	37
[2715]	38	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso)
	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - -
	40
[3]	41	!! * Substitutions
	42	# include "vectopt_loop_substitute.h90"
	43	!!----------------------------------------------------------------------
[2715]	44	!! NEMO/OPA 3.3 , NEMO Consortium (2011)
[1156]	45	!! $Id$
[2715]	46	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	47	!!----------------------------------------------------------------------
	48	CONTAINS
	49
[2715]	50	INTEGER FUNCTION dyn_ldf_iso_alloc()
	51	!!----------------------------------------------------------------------
	52	!! * ROUTINE dyn_ldf_iso_alloc *
	53	!!----------------------------------------------------------------------
	54	ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , &
	55	& zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
	56	!
	57	IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
	58	END FUNCTION dyn_ldf_iso_alloc
	59
	60
[3]	61	SUBROUTINE dyn_ldf_iso( kt )
	62	!!----------------------------------------------------------------------
	63	!! * ROUTINE dyn_ldf_iso *
	64	!!
[455]	65	!! ** Purpose : Compute the before trend of the rotated laplacian
	66	!! operator of lateral momentum diffusion except the diagonal
	67	!! vertical term that will be computed in dynzdf module. Add it
	68	!! to the general trend of momentum equation.
[3]	69	!!
	70	!! ** Method :
[455]	71	!! The momentum lateral diffusive trend is provided by a 2nd
	72	!! order operator rotated along neutral or geopotential surfaces
	73	!! (in s-coordinates).
[3]	74	!! It is computed using before fields (forward in time) and isopyc-
[455]	75	!! nal or geopotential slopes computed in routine ldfslp.
[3]	76	!! Here, u and v components are considered as 2 independent scalar
	77	!! fields. Therefore, the property of splitting divergent and rota-
	78	!! tional part of the flow of the standard, z-coordinate laplacian
	79	!! momentum diffusion is lost.
	80	!! horizontal fluxes associated with the rotated lateral mixing:
	81	!! u-component:
[5836]	82	!! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ ub ]
	83	!! - ahmt e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
	84	!! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ ub ]
	85	!! - ahmf e1f * mi(vslp) dk[ mj(mk(ub)) ]
[3]	86	!! v-component:
[5836]	87	!! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vb ]
	88	!! - ahmf e2t * mj(uslp) dk[ mi(mk(vb)) ]
	89	!! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ ub ]
	90	!! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
[3]	91	!! take the horizontal divergence of the fluxes:
	92	!! diffu = 1/(e1ue2ue3u) { di [ ziut ] + dj-1[ zjuf ] }
	93	!! diffv = 1/(e1ve2ve3v) { di-1[ zivf ] + dj [ zjvt ] }
	94	!! Add this trend to the general trend (ua,va):
	95	!! ua = ua + diffu
[455]	96	!! CAUTION: here the isopycnal part is with a coeff. of aht. This
	97	!! should be modified for applications others than orca_r2 (!!bug)
[3]	98	!!
	99	!! ** Action :
	100	!! Update (ua,va) arrays with the before geopotential biharmonic
	101	!! mixing trend.
[455]	102	!! Update (avmu,avmv) to accompt for the diagonal vertical component
	103	!! of the rotated operator in dynzdf module
[3]	104	!!----------------------------------------------------------------------
[2715]	105	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	106	!
	107	INTEGER :: ji, jj, jk ! dummy loop indices
	108	REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars
[5836]	109	REAL(wp) :: zmskt, zmskf ! - -
[2715]	110	REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - -
	111	REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - -
[3294]	112	!
[7910]	113	REAL(wp), DIMENSION(jpi,jpj) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
[2715]	114	!!----------------------------------------------------------------------
[3294]	115	!
	116	IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_iso')
	117	!
	118	!
[3]	119	IF( kt == nit000 ) THEN
	120	IF(lwp) WRITE(numout,*)
	121	IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
	122	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
[2715]	123	! ! allocate dyn_ldf_bilap arrays
	124	IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
[3]	125	ENDIF
[216]	126
[5836]	127	!!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED
	128	! s-coordinate: Iso-level diffusion on momentum but not on tracer
[455]	129	IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
[2715]	130	!
	131	DO jk = 1, jpk ! set the slopes of iso-level
[455]	132	DO jj = 2, jpjm1
[4488]	133	DO ji = 2, jpim1
[6140]	134	uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
	135	vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
	136	wslpi(ji,jj,jk) = - ( gdepw_b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
	137	wslpj(ji,jj,jk) = - ( gdepw_b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
[455]	138	END DO
	139	END DO
	140	END DO
	141	! Lateral boundary conditions on the slopes
	142	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )
	143	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
	144
	145	!!bug
[2715]	146	IF( kt == nit000 ) then
	147	IF(lwp) WRITE(numout,) ' max slop: u', SQRT( MAXVAL(uslpuslp)), ' v ', SQRT(MAXVAL(vslp)), &
	148	& ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj))
[455]	149	endif
	150	!!end
[216]	151	ENDIF
[455]	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
	162	zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
	163	zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
	164
	165	IF( jk == 1 ) THEN
	166	zdku(:,:) = zdk1u(:,:)
	167	zdkv(:,:) = zdk1v(:,:)
	168	ELSE
	169	zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
	170	zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
	171	ENDIF
	172
	173	! -----f-----
	174	! Horizontal fluxes on U \|
	175	! --------------------=== t u t
	176	! \|
	177	! i-flux at t-point -----f-----
	178
[455]	179	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	180	DO jj = 2, jpjm1
	181	DO ji = fs_2, jpi ! vector opt.
[6140]	182	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u_n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * r1_e1t(ji,jj)
[3]	183
[5836]	184	zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	185	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	186
[5836]	187	zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	188
[5836]	189	ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &
	190	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
[455]	191	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	192	END DO
	193	END DO
	194	ELSE ! other coordinate system (zco or sco) : e3t
	195	DO jj = 2, jpjm1
	196	DO ji = fs_2, jpi ! vector opt.
[6140]	197	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj)
[3]	198
[5836]	199	zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) &
	200	& + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp )
[455]	201
[5836]	202	zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	203
	204	ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &
	205	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	206	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	207	END DO
[3]	208	END DO
[455]	209	ENDIF
[3]	210
	211	! j-flux at f-point
	212	DO jj = 1, jpjm1
	213	DO ji = 1, fs_jpim1 ! vector opt.
[6140]	214	zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj)
[3]	215
[5836]	216	zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	217	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	218
[5836]	219	zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]	220
[455]	221	zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) &
	222	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	223	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
[3]	224	END DO
	225	END DO
	226
	227	! \| t \|
	228	! Horizontal fluxes on V \| \|
	229	! --------------------=== f---v---f
	230	! \| \|
	231	! i-flux at f-point \| t \|
	232
	233	DO jj = 2, jpjm1
	234	DO ji = 1, fs_jpim1 ! vector opt.
[6140]	235	zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj)
[3]	236
[5836]	237	zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	238	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	239
[5836]	240	zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]	241
[5836]	242	zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &
	243	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
	244	& + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
[3]	245	END DO
	246	END DO
	247
	248	! j-flux at t-point
[455]	249	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	250	DO jj = 2, jpj
	251	DO ji = 1, fs_jpim1 ! vector opt.
[6140]	252	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v_n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * r1_e2t(ji,jj)
[3]	253
[5836]	254	zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	255	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	256
[5836]	257	zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]	258
[5836]	259	zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	260	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
[455]	261	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	262	END DO
[3]	263	END DO
[455]	264	ELSE ! other coordinate system (zco or sco) : e3t
	265	DO jj = 2, jpj
	266	DO ji = 1, fs_jpim1 ! vector opt.
[6140]	267	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj)
[3]	268
[455]	269	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	270	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	271
[5836]	272	zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[455]	273
	274	zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	275	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	276	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	277	END DO
	278	END DO
	279	ENDIF
	280
	281
[3]	282	! Second derivative (divergence) and add to the general trend
	283	! -----------------------------------------------------------
	284	DO jj = 2, jpjm1
[5836]	285	DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug
[6140]	286	ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
	287	& + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
	288	va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
	289	& + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[3]	290	END DO
	291	END DO
	292	! ! ===============
	293	END DO ! End of slab
	294	! ! ===============
[216]	295
[455]	296	! print sum trends (used for debugging)
	297	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
	298	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]	299
	300
[455]	301	! ! ===============
	302	DO jj = 2, jpjm1 ! Vertical slab
	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
	315	DO ji = 2, jpi
	316	! i-gradient of u at jj
	317	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) )
	318	! j-gradient of u and v at jj
	319	zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) )
	320	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) )
	321	! j-gradient of u and v at jj+1
	322	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) )
	323	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) )
	324	END DO
	325	END DO
	326	DO jk = 1, jpk
	327	DO ji = 1, jpim1
	328	! i-gradient of v at jj
	329	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) )
	330	END DO
	331	END DO
	332
	333
	334	! I.2 Vertical fluxes
	335	! -------------------
	336
	337	! Surface and bottom vertical fluxes set to zero
	338	DO ji = 1, jpi
	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
	347	DO ji = 2, jpim1
[5836]	348	zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk)
	349	!
[455]	350	zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	351	zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
[5836]	352	!
[455]	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. )
[5836]	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
	358	zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
	359	zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
	360	! vertical flux on u field
	361	zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	362	+zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	363	+ zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	364	+zdj1u(ji,jk ) + zdju (ji ,jk ) )
	365	! update avmu (add isopycnal vertical coefficient to avmu)
[5836]	366	! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
	367	avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
[455]	368	END DO
	369	END DO
	370
	371	! interior (2=<jk=<jpk-1) on V field
	372	DO jk = 2, jpkm1
	373	DO ji = 2, jpim1
[5836]	374	zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk)
[455]	375
	376	zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	377	zvwslpj = zcoef0 * ( 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. )
	383
	384	zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
	385	zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
	386	! vertical flux on v field
	387	zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
[5836]	388	& +zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	389	& + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	390	& +zdjv (ji,jk ) + zdj1v(ji ,jk ) )
[455]	391	! update avmv (add isopycnal vertical coefficient to avmv)
[5836]	392	! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
	393	avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_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
	401	DO ji = 2, jpim1
[6140]	402	ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
	403	va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[455]	404	END DO
	405	END DO
	406	! ! ===============
	407	END DO ! End of slab
	408	! ! ===============
[2715]	409	!
[3294]	410	IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_iso')
	411	!
[3]	412	END SUBROUTINE dyn_ldf_iso
	413
	414	!!======================================================================
	415	END MODULE dynldf_iso

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 7910

Download in other formats: