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.

traldf_triad.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90 @ 10980

Last change on this file since 10980 was 10980, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : interface changes to traldf routines for design compliance and consistency. SETTE tested (GYRE_PISCES, only)
Property svn:keywords set to `Id`
File size: 23.7 KB

Rev	Line
[5758]	1	MODULE traldf_triad
[2371]	2	!!======================================================================
[5758]	3	!! * MODULE traldf_triad *
[2371]	4	!! Ocean tracers: horizontal component of the lateral tracer mixing trend
	5	!!======================================================================
[5758]	6	!! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) Griffies operator (original code)
	7	!! 3.7 ! 2013-12 (F. Lemarie, G. Madec) triad operator (Griffies) + Method of Stabilizing Correction
[2205]	8	!!----------------------------------------------------------------------
[5758]	9
[2205]	10	!!----------------------------------------------------------------------
[5758]	11	!! tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator
[2205]	12	!!----------------------------------------------------------------------
[6140]	13	USE oce ! ocean dynamics and active tracers
	14	USE dom_oce ! ocean space and time domain
	15	USE phycst ! physical constants
	16	USE trc_oce ! share passive tracers/Ocean variables
	17	USE zdf_oce ! ocean vertical physics
	18	USE ldftra ! lateral physics: eddy diffusivity
	19	USE ldfslp ! lateral physics: iso-neutral slopes
	20	USE traldf_iso ! lateral diffusion (Madec operator) (tra_ldf_iso routine)
	21	USE diaptr ! poleward transport diagnostics
[7646]	22	USE diaar5 ! AR5 diagnostics
[6140]	23	USE zpshde ! partial step: hor. derivative (zps_hde routine)
[5758]	24	!
[6140]	25	USE in_out_manager ! I/O manager
	26	USE iom ! I/O library
	27	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	28	USE lib_mpp ! MPP library
[2205]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[5758]	33	PUBLIC tra_ldf_triad ! routine called by traldf.F90
[2205]	34
[5758]	35	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels
[2371]	36
[7646]	37	LOGICAL :: l_ptr ! flag to compute poleward transport
	38	LOGICAL :: l_hst ! flag to compute heat transport
	39
	40
[2205]	41	!! * Substitutions
[2371]	42	# include "vectopt_loop_substitute.h90"
[2205]	43	!!----------------------------------------------------------------------
[9598]	44	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2287]	45	!! $Id$
[10068]	46	!! Software governed by the CeCILL license (see ./LICENSE)
[2205]	47	!!----------------------------------------------------------------------
	48	CONTAINS
	49
[10980]	50	SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, &
	51	& pgu , pgv , pgui, pgvi , &
	52	& pt , pt2, pt_rhs, kjpt, kpass )
[2450]	53	!!----------------------------------------------------------------------
[5758]	54	!! * ROUTINE tra_ldf_triad *
[2450]	55	!!
[3294]	56	!! ** Purpose : Compute the before horizontal tracer (t & s) diffusive
	57	!! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
[2450]	58	!! add it to the general trend of tracer equation.
	59	!!
[3294]	60	!! ** Method : The horizontal component of the lateral diffusive trends
[2450]	61	!! is provided by a 2nd order operator rotated along neural or geopo-
	62	!! tential surfaces to which an eddy induced advection can be added
	63	!! It is computed using before fields (forward in time) and isopyc-
	64	!! nal or geopotential slopes computed in routine ldfslp.
	65	!!
[5758]	66	!! see documentation for the desciption
[2450]	67	!!
[10980]	68	!! ** Action : pt_rhs updated with the before rotated diffusion
[5758]	69	!! ah_wslp2 ....
	70	!! akz stabilizing vertical diffusivity coefficient (used in trazdf_imp)
[2450]	71	!!----------------------------------------------------------------------
	72	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	73	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2450]	74	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	75	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[5758]	76	INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage
[10922]	77	INTEGER , INTENT(in) :: Kmm ! ocean time level indices
[5758]	78	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s]
[5777]	79	REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels
[5758]	80	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels
[10980]	81	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2)
	82	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt2 ! tracer (only used in kpass=2)
	83	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[2715]	84	!
[5758]	85	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	86	INTEGER :: ip,jp,kp ! dummy loop indices
	87	INTEGER :: ierr ! local integer
	88	REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars
	89	REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - -
	90	REAL(wp) :: zcoef0, ze3w_2, zsign, z2dt, z1_2dt ! - -
[2371]	91	!
[2454]	92	REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv
[5758]	93	REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
[2454]	94	REAL(wp) :: zah, zah_slp, zaei_slp
[9019]	95	REAL(wp), DIMENSION(jpi,jpj ) :: z2d ! 2D workspace
	96	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw ! 3D -
[2205]	97	!!----------------------------------------------------------------------
[3294]	98	!
[5758]	99	IF( .NOT.ALLOCATED(zdkt3d) ) THEN
	100	ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr )
[10425]	101	CALL mpp_sum ( 'traldf_triad', ierr )
[5758]	102	IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays')
[2450]	103	ENDIF
[5758]	104	!
	105	IF( kpass == 1 .AND. kt == kit000 ) THEN
	106	IF(lwp) WRITE(numout,*)
	107	IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype
	108	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
	109	ENDIF
[7646]	110	!
	111	l_hst = .FALSE.
	112	l_ptr = .FALSE.
	113	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
	114	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
	115	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
	116	!
	117	! ! set time step size (Euler/Leapfrog)
[6140]	118	IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt ! at nit000 (Euler)
	119	ELSE ; z2dt = 2.* rdt ! (Leapfrog)
[5758]	120	ENDIF
	121	z1_2dt = 1._wp / z2dt
	122	!
	123	IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0)
	124	ELSE ; zsign = -1._wp
	125	ENDIF
[6140]	126	!
[2205]	127	!!----------------------------------------------------------------------
[5758]	128	!! 0 - calculate ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw
[2371]	129	!!----------------------------------------------------------------------
[5758]	130	!
	131	IF( kpass == 1 ) THEN !== first pass only and whatever the tracer is ==!
	132	!
	133	akz (:,:,:) = 0._wp
	134	ah_wslp2(:,:,:) = 0._wp
	135	IF( ln_ldfeiv_dia ) THEN
	136	zpsi_uw(:,:,:) = 0._wp
	137	zpsi_vw(:,:,:) = 0._wp
	138	ENDIF
	139	!
	140	DO ip = 0, 1 ! i-k triads
	141	DO kp = 0, 1
	142	DO jk = 1, jpkm1
	143	DO jj = 1, jpjm1
	144	DO ji = 1, fs_jpim1
[10954]	145	ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
	146	zbu = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
[5758]	147	zah = 0.25_wp * pahu(ji,jj,jk)
	148	zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
	149	! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by adding gradient of depth)
[10954]	150	zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
[5758]	151	zslope2 = zslope2 *zslope2
	152	ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
	153	akz (ji+ip,jj,jk+kp) = akz (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj) &
	154	& * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
[6140]	155	!
[5758]	156	IF( ln_ldfeiv_dia ) zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &
	157	& + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew
	158	END DO
[2450]	159	END DO
	160	END DO
	161	END DO
	162	END DO
[5758]	163	!
	164	DO jp = 0, 1 ! j-k triads
	165	DO kp = 0, 1
	166	DO jk = 1, jpkm1
	167	DO jj = 1, jpjm1
	168	DO ji = 1, fs_jpim1
[10954]	169	ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
	170	zbv = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
[5758]	171	zah = 0.25_wp * pahv(ji,jj,jk)
	172	zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
	173	! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
	174	! (do this by adding gradient of depth)
[10954]	175	zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
[5758]	176	zslope2 = zslope2 * zslope2
	177	ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
	178	akz (ji,jj+jp,jk+kp) = akz (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj) &
	179	& * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
	180	!
	181	IF( ln_ldfeiv_dia ) zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) &
	182	& + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew
	183	END DO
[2450]	184	END DO
	185	END DO
	186	END DO
	187	END DO
[5147]	188	!
[5758]	189	IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient
	190	!
	191	IF( ln_traldf_blp ) THEN ! bilaplacian operator
	192	DO jk = 2, jpkm1
	193	DO jj = 1, jpjm1
	194	DO ji = 1, fs_jpim1
	195	akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) &
[10954]	196	& * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) ) )
[5758]	197	END DO
[5147]	198	END DO
[3294]	199	END DO
[5758]	200	ELSEIF( ln_traldf_lap ) THEN ! laplacian operator
	201	DO jk = 2, jpkm1
	202	DO jj = 1, jpjm1
	203	DO ji = 1, fs_jpim1
[10954]	204	ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
[5758]	205	zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )
	206	akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
	207	END DO
	208	END DO
	209	END DO
	210	ENDIF
	211	!
	212	ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
	213	akz(:,:,:) = ah_wslp2(:,:,:)
[5147]	214	ENDIF
	215	!
[10922]	216	IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
[5758]	217	!
	218	ENDIF !== end 1st pass only ==!
	219	!
	220	! ! ===========
	221	DO jn = 1, kjpt ! tracer loop
	222	! ! ===========
[2371]	223	! Zero fluxes for each tracer
[5758]	224	!!gm this should probably be done outside the jn loop
[2371]	225	ztfw(:,:,:) = 0._wp
	226	zftu(:,:,:) = 0._wp
	227	zftv(:,:,:) = 0._wp
[3294]	228	!
[5758]	229	DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==!
[2371]	230	DO jj = 1, jpjm1
	231	DO ji = 1, fs_jpim1 ! vector opt.
[10980]	232	zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
	233	zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
[2371]	234	END DO
[2205]	235	END DO
	236	END DO
[5758]	237	IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction at top/bottom ocean level
	238	DO jj = 1, jpjm1 ! bottom level
	239	DO ji = 1, fs_jpim1 ! vector opt.
[3294]	240	zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
	241	zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
[2371]	242	END DO
	243	END DO
[5758]	244	IF( ln_isfcav ) THEN ! top level (ocean cavities only)
	245	DO jj = 1, jpjm1
	246	DO ji = 1, fs_jpim1 ! vector opt.
	247	IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)
	248	IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)
	249	END DO
	250	END DO
	251	ENDIF
[2371]	252	ENDIF
[6140]	253	!
[2371]	254	!!----------------------------------------------------------------------
	255	!! II - horizontal trend (full)
	256	!!----------------------------------------------------------------------
	257	!
	258	DO jk = 1, jpkm1
	259	! !== Vertical tracer gradient at level jk and jk+1
[10980]	260	zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
[2371]	261	!
[3294]	262	! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
	263	IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1)
[10980]	264	ELSE ; zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
[2371]	265	ENDIF
[5758]	266	!
	267	zaei_slp = 0._wp
	268	!
	269	IF( ln_botmix_triad ) THEN
[3294]	270	DO ip = 0, 1 !== Horizontal & vertical fluxes
	271	DO kp = 0, 1
	272	DO jj = 1, jpjm1
	273	DO ji = 1, fs_jpim1
[5758]	274	ze1ur = r1_e1u(ji,jj)
[3294]	275	zdxt = zdit(ji,jj,jk) * ze1ur
[10954]	276	ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
[3294]	277	zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr
	278	zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
[5758]	279	zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp)
[6140]	280	!
[10954]	281	zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
[5758]	282	! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked....
	283	zah = pahu(ji,jj,jk)
[3294]	284	zah_slp = zah * zslope_iso
[5758]	285	IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew
	286	zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
	287	ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt * zbu * ze3wr
[3294]	288	END DO
[2371]	289	END DO
	290	END DO
	291	END DO
[6140]	292	!
[3294]	293	DO jp = 0, 1
	294	DO kp = 0, 1
	295	DO jj = 1, jpjm1
	296	DO ji = 1, fs_jpim1
[5758]	297	ze2vr = r1_e2v(ji,jj)
[3294]	298	zdyt = zdjt(ji,jj,jk) * ze2vr
[10954]	299	ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
[3294]	300	zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr
	301	zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
	302	zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp)
[10954]	303	zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
[5758]	304	! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked...
	305	zah = pahv(ji,jj,jk)
[3294]	306	zah_slp = zah * zslope_iso
[5758]	307	IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew
	308	zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
	309	ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt * zbv * ze3wr
[3294]	310	END DO
[2371]	311	END DO
	312	END DO
	313	END DO
[6140]	314	!
[3294]	315	ELSE
[6140]	316	!
[5758]	317	DO ip = 0, 1 !== Horizontal & vertical fluxes
[3294]	318	DO kp = 0, 1
	319	DO jj = 1, jpjm1
	320	DO ji = 1, fs_jpim1
[5758]	321	ze1ur = r1_e1u(ji,jj)
[3294]	322	zdxt = zdit(ji,jj,jk) * ze1ur
[10954]	323	ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
[3294]	324	zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr
	325	zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
	326	zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp)
[6140]	327	!
[10954]	328	zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
[5758]	329	! ln_botmix_triad is .F. mask zah for bottom half cells
	330	zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ????
[3294]	331	zah_slp = zah * zslope_iso
[6140]	332	IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! aeit(ji+ip,jj,jk)*zslope_skew
[5758]	333	zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
[3294]	334	ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
	335	END DO
	336	END DO
	337	END DO
	338	END DO
[6140]	339	!
[3294]	340	DO jp = 0, 1
	341	DO kp = 0, 1
	342	DO jj = 1, jpjm1
	343	DO ji = 1, fs_jpim1
[5758]	344	ze2vr = r1_e2v(ji,jj)
[3294]	345	zdyt = zdjt(ji,jj,jk) * ze2vr
[10954]	346	ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
[3294]	347	zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr
	348	zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
	349	zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp)
[10954]	350	zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
[5758]	351	! ln_botmix_triad is .F. mask zah for bottom half cells
	352	zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ????
[3294]	353	zah_slp = zah * zslope_iso
[6140]	354	IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! aeit(ji,jj+jp,jk)*zslope_skew
[3294]	355	zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
	356	ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
	357	END DO
	358	END DO
	359	END DO
	360	END DO
[5758]	361	ENDIF
	362	! !== horizontal divergence and add to the general trend ==!
[2450]	363	DO jj = 2 , jpjm1
[3294]	364	DO ji = fs_2, fs_jpim1 ! vector opt.
[10980]	365	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) &
[5758]	366	& + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) &
[10954]	367	& / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
[2450]	368	END DO
	369	END DO
	370	!
	371	END DO
	372	!
[5758]	373	! !== add the vertical 33 flux ==!
	374	IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz
	375	DO jk = 2, jpkm1
	376	DO jj = 1, jpjm1
	377	DO ji = fs_2, fs_jpim1
[10954]	378	ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
[5758]	379	& * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) &
[10980]	380	& * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
[5758]	381	END DO
	382	END DO
	383	END DO
	384	ELSE ! bilaplacian
	385	SELECT CASE( kpass )
	386	CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2
	387	DO jk = 2, jpkm1
	388	DO jj = 1, jpjm1
	389	DO ji = fs_2, fs_jpim1
[10954]	390	ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
[10980]	391	& * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
[5758]	392	END DO
	393	END DO
	394	END DO
[10980]	395	CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp.
[5758]	396	DO jk = 2, jpkm1
	397	DO jj = 1, jpjm1
	398	DO ji = fs_2, fs_jpim1
[10954]	399	ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
[10980]	400	& * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) &
	401	& + akz (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) )
[5758]	402	END DO
	403	END DO
	404	END DO
	405	END SELECT
	406	ENDIF
	407	!
[10980]	408	DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pt_rhs ==!
[2450]	409	DO jj = 2, jpjm1
[3294]	410	DO ji = fs_2, fs_jpim1 ! vector opt.
[10980]	411	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &
	412	& / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
[2450]	413	END DO
	414	END DO
	415	END DO
	416	!
[5758]	417	IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==!
	418	( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass (bilaplacian) ==!
	419	!
	420	! ! "Poleward" diffusive heat or salt transports (T-S case only)
[7646]	421	IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:) )
	422	! ! Diffusive heat transports
	423	IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
[5758]	424	!
	425	ENDIF !== end pass selection ==!
[2450]	426	!
[5758]	427	! ! ===============
	428	END DO ! end tracer loop
	429	! ! ===============
	430	END SUBROUTINE tra_ldf_triad
[2371]	431
[2205]	432	!!==============================================================================
[5758]	433	END MODULE traldf_triad

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

Download in other formats: