Context Navigation

dynldf_iso.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

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

Download in other formats: