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/trunk/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynldf_iso.F90 @ 14325

Last change on this file since 14325 was 14215, checked in by acc, 3 years ago
trunk changes to swap the order of arguments to the DO LOOP macros. These changes result in a more natural i-j-k ordering as explained in #2595. SETTE is passed before and after these changes and results are unchanged. This fixes #2595
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
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[2715]	34	PUBLIC dyn_ldf_iso ! called by step.F90
	35	PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90
[3]	36
[9019]	37	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity
	38
[2715]	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso)
	40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - -
	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	!!----------------------------------------------------------------------
[9019]	56	ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , &
	57	& akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
[2715]	58	!
	59	IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
	60	END FUNCTION dyn_ldf_iso_alloc
	61
	62
[12377]	63	SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
[3]	64	!!----------------------------------------------------------------------
	65	!! * ROUTINE dyn_ldf_iso *
	66	!!
[455]	67	!! ** Purpose : Compute the before trend of the rotated laplacian
	68	!! operator of lateral momentum diffusion except the diagonal
	69	!! vertical term that will be computed in dynzdf module. Add it
	70	!! to the general trend of momentum equation.
[3]	71	!!
	72	!! ** Method :
[455]	73	!! The momentum lateral diffusive trend is provided by a 2nd
	74	!! order operator rotated along neutral or geopotential surfaces
	75	!! (in s-coordinates).
[3]	76	!! It is computed using before fields (forward in time) and isopyc-
[455]	77	!! nal or geopotential slopes computed in routine ldfslp.
[3]	78	!! Here, u and v components are considered as 2 independent scalar
	79	!! fields. Therefore, the property of splitting divergent and rota-
	80	!! tional part of the flow of the standard, z-coordinate laplacian
	81	!! momentum diffusion is lost.
	82	!! horizontal fluxes associated with the rotated lateral mixing:
	83	!! u-component:
[12377]	84	!! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ uu ]
	85	!! - ahmt e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
	86	!! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ uu ]
	87	!! - ahmf e1f * mi(vslp) dk[ mj(mk(uu)) ]
[3]	88	!! v-component:
[12377]	89	!! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vv ]
	90	!! - ahmf e2t * mj(uslp) dk[ mi(mk(vv)) ]
	91	!! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ vv ]
	92	!! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
[3]	93	!! take the horizontal divergence of the fluxes:
	94	!! diffu = 1/(e1ue2ue3u) { di [ ziut ] + dj-1[ zjuf ] }
	95	!! diffv = 1/(e1ve2ve3v) { di-1[ zivf ] + dj [ zjvt ] }
[12377]	96	!! Add this trend to the general trend (uu(rhs),vv(rhs)):
	97	!! uu(rhs) = uu(rhs) + diffu
[455]	98	!! CAUTION: here the isopycnal part is with a coeff. of aht. This
	99	!! should be modified for applications others than orca_r2 (!!bug)
[3]	100	!!
	101	!! ** Action :
[12377]	102	!! -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
[9019]	103	!! -(akzu,akzv) to accompt for the diagonal vertical component
	104	!! of the rotated operator in dynzdf module
[3]	105	!!----------------------------------------------------------------------
[12377]	106	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	107	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
	108	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
[2715]	109	!
	110	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	111	REAL(wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj ! local scalars
	112	REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - -
[9490]	113	REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0 ! - -
[9019]	114	REAL(wp), DIMENSION(jpi,jpj) :: ziut, zivf, zdku, zdk1u ! 2D workspace
	115	REAL(wp), DIMENSION(jpi,jpj) :: zjuf, zjvt, zdkv, zdk1v ! - -
[2715]	116	!!----------------------------------------------------------------------
[3294]	117	!
[3]	118	IF( kt == nit000 ) THEN
	119	IF(lwp) WRITE(numout,*)
	120	IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
	121	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
[2715]	122	! ! allocate dyn_ldf_bilap arrays
	123	IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
[3]	124	ENDIF
[216]	125
[5836]	126	!!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED
	127	! s-coordinate: Iso-level diffusion on momentum but not on tracer
[455]	128	IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
[2715]	129	!
[13497]	130	DO_3D( 0, 0, 0, 0, 1, jpk ) ! set the slopes of iso-level
[12377]	131	uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
	132	vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
	133	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
	134	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
	135	END_3D
[455]	136	! Lateral boundary conditions on the slopes
[13226]	137	CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
[9019]	138	!
	139	ENDIF
[9490]	140
	141	zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max
	142
[3]	143	! ! ===============
	144	DO jk = 1, jpkm1 ! Horizontal slab
	145	! ! ===============
	146
	147	! Vertical u- and v-shears at level jk and jk+1
	148	! ---------------------------------------------
	149	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	150	! zdkv(jk=1)=zdkv(jk=2)
	151
[12377]	152	zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
	153	zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
[3]	154
	155	IF( jk == 1 ) THEN
	156	zdku(:,:) = zdk1u(:,:)
	157	zdkv(:,:) = zdk1v(:,:)
	158	ELSE
[12377]	159	zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
	160	zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
[3]	161	ENDIF
	162
	163	! -----f-----
	164	! Horizontal fluxes on U \|
	165	! --------------------=== t u t
	166	! \|
	167	! i-flux at t-point -----f-----
	168
[455]	169	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
[14215]	170	DO_2D( 0, 1, 0, 0 )
[13237]	171	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) &
	172	& * MIN( e3u(ji ,jj,jk,Kmm), &
	173	& e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
[3]	174
[12377]	175	zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	176	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	177
[12377]	178	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	179
	180	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
	181	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	182	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	183	END_2D
[455]	184	ELSE ! other coordinate system (zco or sco) : e3t
[14215]	185	DO_2D( 0, 1, 0, 0 )
[13237]	186	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
	187	& * e2t(ji,jj) * e3t(ji,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 )
[455]	191
[12377]	192	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	193
[12377]	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	ENDIF
[3]	199
	200	! j-flux at f-point
[13295]	201	DO_2D( 1, 0, 1, 0 )
[13237]	202	zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
	203	& * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
[3]	204
[12377]	205	zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	206	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	207
[12377]	208	zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]	209
[12377]	210	zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) &
	211	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	212	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
	213	END_2D
[3]	214
	215	! \| t \|
	216	! Horizontal fluxes on V \| \|
	217	! --------------------=== f---v---f
	218	! \| \|
	219	! i-flux at f-point \| t \|
	220
[14215]	221	DO_2D( 1, 0, 0, 0 )
[13237]	222	zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
	223	& * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
[3]	224
[12377]	225	zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	226	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	227
[12377]	228	zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]	229
[12377]	230	zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) &
	231	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
	232	& + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
	233	END_2D
[3]	234
	235	! j-flux at t-point
[455]	236	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
[14215]	237	DO_2D( 1, 0, 0, 1 )
[13237]	238	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) &
	239	& * MIN( e3v(ji,jj ,jk,Kmm), &
	240	& e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
[3]	241
[12377]	242	zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	243	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	244
[12377]	245	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]	246
[12377]	247	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
	248	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	249	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	250	END_2D
[455]	251	ELSE ! other coordinate system (zco or sco) : e3t
[14215]	252	DO_2D( 1, 0, 0, 1 )
[13237]	253	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
	254	& * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
[3]	255
[12377]	256	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	257	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	258
[12377]	259	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[455]	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	ENDIF
	266
	267
[3]	268	! Second derivative (divergence) and add to the general trend
	269	! -----------------------------------------------------------
[13497]	270	DO_2D( 0, 0, 0, 0 ) !!gm Question vectop possible??? !!bug
[12377]	271	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
[13237]	272	& + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) &
	273	& / e3u(ji,jj,jk,Kmm)
[12377]	274	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
[13237]	275	& + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) &
	276	& / e3v(ji,jj,jk,Kmm)
[12377]	277	END_2D
[3]	278	! ! ===============
	279	END DO ! End of slab
	280	! ! ===============
[216]	281
[455]	282	! print sum trends (used for debugging)
[12377]	283	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
	284	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]	285
	286
[455]	287	! ! ===============
	288	DO jj = 2, jpjm1 ! Vertical slab
	289	! ! ===============
	290
	291
	292	! I. vertical trends associated with the lateral mixing
	293	! =====================================================
	294	! (excluding the vertical flux proportional to dk[t]
	295
	296
	297	! I.1 horizontal momentum gradient
	298	! --------------------------------
	299
	300	DO jk = 1, jpk
	301	DO ji = 2, jpi
	302	! i-gradient of u at jj
[12377]	303	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) )
[455]	304	! j-gradient of u and v at jj
[12377]	305	zdju (ji,jk) = fmask(ji,jj ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) )
	306	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pvv(ji,jj ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) )
[455]	307	! j-gradient of u and v at jj+1
[12377]	308	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) )
	309	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	310	END DO
	311	END DO
	312	DO jk = 1, jpk
	313	DO ji = 1, jpim1
	314	! i-gradient of v at jj
[12377]	315	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	316	END DO
	317	END DO
	318
	319
	320	! I.2 Vertical fluxes
	321	! -------------------
	322
	323	! Surface and bottom vertical fluxes set to zero
	324	DO ji = 1, jpi
	325	zfuw(ji, 1 ) = 0.e0
	326	zfvw(ji, 1 ) = 0.e0
	327	zfuw(ji,jpk) = 0.e0
	328	zfvw(ji,jpk) = 0.e0
	329	END DO
	330
	331	! interior (2=<jk=<jpk-1) on U field
	332	DO jk = 2, jpkm1
	333	DO ji = 2, jpim1
[9490]	334	zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
[5836]	335	!
[9019]	336	zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	337	zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
[5836]	338	!
[9019]	339	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	340	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. )
	341	zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) &
	342	+ fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. )
[455]	343
[9019]	344	zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
	345	zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
[455]	346	! vertical flux on u field
[9019]	347	zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	348	& + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	349	& + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	350	& + zdj1u(ji,jk ) + zdju (ji ,jk ) )
	351	! vertical mixing coefficient (akzu)
[9490]	352	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	353	akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
[455]	354	END DO
	355	END DO
	356
	357	! interior (2=<jk=<jpk-1) on V field
	358	DO jk = 2, jpkm1
	359	DO ji = 2, jpim1
[9490]	360	zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
[9019]	361	!
	362	zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	363	zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	364	!
	365	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	366	& + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ) , 1. )
	367	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	368	& + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ) , 1. )
[455]	369
[9019]	370	zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
	371	zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
[455]	372	! vertical flux on v field
[9019]	373	zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	374	& + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	375	& + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	376	& + zdjv (ji,jk ) + zdj1v(ji ,jk ) )
	377	! vertical mixing coefficient (akzv)
[9490]	378	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	379	akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
[455]	380	END DO
	381	END DO
	382
	383
	384	! I.3 Divergence of vertical fluxes added to the general tracer trend
	385	! -------------------------------------------------------------------
	386	DO jk = 1, jpkm1
	387	DO ji = 2, jpim1
[13237]	388	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) &
	389	& / e3u(ji,jj,jk,Kmm)
	390	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) &
	391	& / e3v(ji,jj,jk,Kmm)
[455]	392	END DO
	393	END DO
	394	! ! ===============
	395	END DO ! End of slab
	396	! ! ===============
[3]	397	END SUBROUTINE dyn_ldf_iso
	398
	399	!!======================================================================
	400	END MODULE dynldf_iso

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

Download in other formats: