Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynldf_iso.F90 @ 8576

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

Download in other formats: