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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_iso.F90 @ 10946

Last change on this file since 10946 was 10928, checked in by davestorkey, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Finish converting DYN module, including some updates to previously processed modules, but excluding dynnxt.F90 (which needs to be completely rewritten) and wet_dry.F90 - I need to talk to Enda about this.
Property svn:keywords set to `Id`
File size: 20.9 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
	43	# include "vectopt_loop_substitute.h90"
	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
[10928]	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:
[10928]	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:
[10928]	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 ] }
[10928]	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 :
[10928]	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	!!----------------------------------------------------------------------
[10928]	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	!
	129	DO jk = 1, jpk ! set the slopes of iso-level
[455]	130	DO jj = 2, jpjm1
[4488]	131	DO ji = 2, jpim1
[10928]	132	uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
	133	vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
	134	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
	135	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
[455]	136	END DO
	137	END DO
	138	END DO
	139	! Lateral boundary conditions on the slopes
[10425]	140	CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1., vslp , 'V', -1., wslpi, 'W', -1., wslpj, 'W', -1. )
[9019]	141	!
	142	ENDIF
[9490]	143
	144	zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max
	145
[3]	146	! ! ===============
	147	DO jk = 1, jpkm1 ! Horizontal slab
	148	! ! ===============
	149
	150	! Vertical u- and v-shears at level jk and jk+1
	151	! ---------------------------------------------
	152	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	153	! zdkv(jk=1)=zdkv(jk=2)
	154
[10928]	155	zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
	156	zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
[3]	157
	158	IF( jk == 1 ) THEN
	159	zdku(:,:) = zdk1u(:,:)
	160	zdkv(:,:) = zdk1v(:,:)
	161	ELSE
[10928]	162	zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
	163	zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
[3]	164	ENDIF
	165
	166	! -----f-----
	167	! Horizontal fluxes on U \|
	168	! --------------------=== t u t
	169	! \|
	170	! i-flux at t-point -----f-----
	171
[455]	172	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	173	DO jj = 2, jpjm1
	174	DO ji = fs_2, jpi ! vector opt.
[10928]	175	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]	176
[5836]	177	zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	178	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	179
[9490]	180	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	181
[10928]	182	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
[5836]	183	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
[455]	184	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	185	END DO
	186	END DO
	187	ELSE ! other coordinate system (zco or sco) : e3t
	188	DO jj = 2, jpjm1
	189	DO ji = fs_2, jpi ! vector opt.
[10928]	190	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
[3]	191
[5836]	192	zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) &
	193	& + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp )
[455]	194
[9490]	195	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	196
[10928]	197	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
[455]	198	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	199	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	200	END DO
[3]	201	END DO
[455]	202	ENDIF
[3]	203
	204	! j-flux at f-point
	205	DO jj = 1, jpjm1
	206	DO ji = 1, fs_jpim1 ! vector opt.
[10928]	207	zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
[3]	208
[5836]	209	zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	210	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	211
[9490]	212	zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]	213
[10928]	214	zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) &
[455]	215	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	216	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
[3]	217	END DO
	218	END DO
	219
	220	! \| t \|
	221	! Horizontal fluxes on V \| \|
	222	! --------------------=== f---v---f
	223	! \| \|
	224	! i-flux at f-point \| t \|
	225
	226	DO jj = 2, jpjm1
	227	DO ji = 1, fs_jpim1 ! vector opt.
[10928]	228	zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
[3]	229
[5836]	230	zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	231	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	232
[9490]	233	zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]	234
[10928]	235	zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) &
[5836]	236	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
	237	& + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
[3]	238	END DO
	239	END DO
	240
	241	! j-flux at t-point
[455]	242	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	243	DO jj = 2, jpj
	244	DO ji = 1, fs_jpim1 ! vector opt.
[10928]	245	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]	246
[5836]	247	zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	248	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	249
[9490]	250	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]	251
[10928]	252	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
[5836]	253	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
[455]	254	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	255	END DO
[3]	256	END DO
[455]	257	ELSE ! other coordinate system (zco or sco) : e3t
	258	DO jj = 2, jpj
	259	DO ji = 1, fs_jpim1 ! vector opt.
[10928]	260	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
[3]	261
[455]	262	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	263	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	264
[9490]	265	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[455]	266
[10928]	267	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
[455]	268	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	269	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	270	END DO
	271	END DO
	272	ENDIF
	273
	274
[3]	275	! Second derivative (divergence) and add to the general trend
	276	! -----------------------------------------------------------
	277	DO jj = 2, jpjm1
[5836]	278	DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug
[10928]	279	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
	280	& + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
	281	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
	282	& + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
[3]	283	END DO
	284	END DO
	285	! ! ===============
	286	END DO ! End of slab
	287	! ! ===============
[216]	288
[455]	289	! print sum trends (used for debugging)
[10928]	290	IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
	291	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]	292
	293
[455]	294	! ! ===============
	295	DO jj = 2, jpjm1 ! Vertical slab
	296	! ! ===============
	297
	298
	299	! I. vertical trends associated with the lateral mixing
	300	! =====================================================
	301	! (excluding the vertical flux proportional to dk[t]
	302
	303
	304	! I.1 horizontal momentum gradient
	305	! --------------------------------
	306
	307	DO jk = 1, jpk
	308	DO ji = 2, jpi
	309	! i-gradient of u at jj
[10928]	310	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) )
[455]	311	! j-gradient of u and v at jj
[10928]	312	zdju (ji,jk) = fmask(ji,jj ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) )
	313	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pvv(ji,jj ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) )
[455]	314	! j-gradient of u and v at jj+1
[10928]	315	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) )
	316	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	317	END DO
	318	END DO
	319	DO jk = 1, jpk
	320	DO ji = 1, jpim1
	321	! i-gradient of v at jj
[10928]	322	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	323	END DO
	324	END DO
	325
	326
	327	! I.2 Vertical fluxes
	328	! -------------------
	329
	330	! Surface and bottom vertical fluxes set to zero
	331	DO ji = 1, jpi
	332	zfuw(ji, 1 ) = 0.e0
	333	zfvw(ji, 1 ) = 0.e0
	334	zfuw(ji,jpk) = 0.e0
	335	zfvw(ji,jpk) = 0.e0
	336	END DO
	337
	338	! interior (2=<jk=<jpk-1) on U field
	339	DO jk = 2, jpkm1
	340	DO ji = 2, jpim1
[9490]	341	zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
[5836]	342	!
[9019]	343	zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	344	zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
[5836]	345	!
[9019]	346	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	347	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. )
	348	zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) &
	349	+ fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. )
[455]	350
[9019]	351	zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
	352	zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
[455]	353	! vertical flux on u field
[9019]	354	zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	355	& + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	356	& + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	357	& + zdj1u(ji,jk ) + zdju (ji ,jk ) )
	358	! vertical mixing coefficient (akzu)
[9490]	359	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	360	akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
[455]	361	END DO
	362	END DO
	363
	364	! interior (2=<jk=<jpk-1) on V field
	365	DO jk = 2, jpkm1
	366	DO ji = 2, jpim1
[9490]	367	zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
[9019]	368	!
	369	zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	370	zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	371	!
	372	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	373	& + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ) , 1. )
	374	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	375	& + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ) , 1. )
[455]	376
[9019]	377	zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
	378	zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
[455]	379	! vertical flux on v field
[9019]	380	zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	381	& + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	382	& + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	383	& + zdjv (ji,jk ) + zdj1v(ji ,jk ) )
	384	! vertical mixing coefficient (akzv)
[9490]	385	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	386	akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
[455]	387	END DO
	388	END DO
	389
	390
	391	! I.3 Divergence of vertical fluxes added to the general tracer trend
	392	! -------------------------------------------------------------------
	393	DO jk = 1, jpkm1
	394	DO ji = 2, jpim1
[10928]	395	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)
	396	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]	397	END DO
	398	END DO
	399	! ! ===============
	400	END DO ! End of slab
	401	! ! ===============
[3]	402	END SUBROUTINE dyn_ldf_iso
	403
	404	!!======================================================================
	405	END MODULE dynldf_iso

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

Download in other formats: