Context Navigation

traadv_fct.F90 @ 12396

Last change on this file since 12396 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: 32.5 KB

Rev	Line
[5770]	1	MODULE traadv_fct
[3]	2	!!==============================================================================
[5770]	3	!! * MODULE traadv_fct *
	4	!! Ocean tracers: horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
[3]	5	!!==============================================================================
[5770]	6	!! History : 3.7 ! 2015-09 (L. Debreu, G. Madec) original code (inspired from traadv_tvd.F90)
[503]	7	!!----------------------------------------------------------------------
[3]	8
	9	!!----------------------------------------------------------------------
[5770]	10	!! tra_adv_fct : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
	11	!! with sub-time-stepping in the vertical direction
	12	!! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm
	13	!! interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
[3]	14	!!----------------------------------------------------------------------
[3625]	15	USE oce ! ocean dynamics and active tracers
	16	USE dom_oce ! ocean space and time domain
[4990]	17	USE trc_oce ! share passive tracers/Ocean variables
	18	USE trd_oce ! trends: ocean variables
[3625]	19	USE trdtra ! tracers trends
[4990]	20	USE diaptr ! poleward transport diagnostics
[7646]	21	USE diaar5 ! AR5 diagnostics
[9019]	22	USE phycst , ONLY : rau0_rcp
[11407]	23	USE zdf_oce , ONLY : ln_zad_Aimp
[4990]	24	!
[5770]	25	USE in_out_manager ! I/O manager
[9019]	26	USE iom !
[3625]	27	USE lib_mpp ! MPP library
	28	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[5770]	29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[9019]	34	PUBLIC tra_adv_fct ! called by traadv.F90
	35	PUBLIC interp_4th_cpt ! called by traadv_cen.F90
[3]	36
[5770]	37	LOGICAL :: l_trd ! flag to compute trends
[7646]	38	LOGICAL :: l_ptr ! flag to compute poleward transport
	39	LOGICAL :: l_hst ! flag to compute heat/salt transport
[5770]	40	REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
[2528]	41
[7646]	42	! ! tridiag solver associated indices:
	43	INTEGER, PARAMETER :: np_NH = 0 ! Neumann homogeneous boundary condition
	44	INTEGER, PARAMETER :: np_CEN2 = 1 ! 2nd order centered boundary condition
	45
[3]	46	!! * Substitutions
[12377]	47	# include "do_loop_substitute.h90"
[3]	48	!!----------------------------------------------------------------------
[9598]	49	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	50	!! $Id$
[10068]	51	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	52	!!----------------------------------------------------------------------
	53	CONTAINS
	54
[12377]	55	SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pU, pV, pW, &
	56	& Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
[3]	57	!!----------------------------------------------------------------------
[5770]	58	!! * ROUTINE tra_adv_fct *
[3]	59	!!
[6140]	60	!! ** Purpose : Compute the now trend due to total advection of tracers
	61	!! and add it to the general trend of tracer equations
[3]	62	!!
[5770]	63	!! ** Method : - 2nd or 4th FCT scheme on the horizontal direction
	64	!! (choice through the value of kn_fct)
[6140]	65	!! - on the vertical the 4th order is a compact scheme
[5770]	66	!! - corrected flux (monotonic correction)
[3]	67	!!
[12377]	68	!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
[9019]	69	!! - send trends to trdtra module for further diagnostics (l_trdtra=T)
[12377]	70	!! - poleward advective heat and salt transport (ln_diaptr=T)
[503]	71	!!----------------------------------------------------------------------
[12377]	72	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	73	INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
	74	INTEGER , INTENT(in ) :: kit000 ! first time step index
	75	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	76	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	77	INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4)
	78	INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4)
	79	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
	80	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
	81	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
[2715]	82	!
[5770]	83	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[6140]	84	REAL(wp) :: ztra ! local scalar
[5770]	85	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - -
	86	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - -
[9019]	87	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
	88	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry
[11407]	89	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup
	90	LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection
[3]	91	!!----------------------------------------------------------------------
[3294]	92	!
	93	IF( kt == kit000 ) THEN
[2528]	94	IF(lwp) WRITE(numout,*)
[5770]	95	IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
[2528]	96	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3]	97	ENDIF
[2528]	98	!
[9019]	99	l_trd = .FALSE. ! set local switches
[7646]	100	l_hst = .FALSE.
	101	l_ptr = .FALSE.
[11407]	102	ll_zAimp = .FALSE.
[12377]	103	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
	104	IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
	105	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
[9019]	106	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
[5770]	107	!
[7646]	108	IF( l_trd .OR. l_hst ) THEN
[9019]	109	ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
[7753]	110	ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp
[3294]	111	ENDIF
[2528]	112	!
[7646]	113	IF( l_ptr ) THEN
[9019]	114	ALLOCATE( zptry(jpi,jpj,jpk) )
[7753]	115	zptry(:,:,:) = 0._wp
[7646]	116	ENDIF
[6140]	117	! ! surface & bottom value : flux set to zero one for all
[7753]	118	zwz(:,:, 1 ) = 0._wp
	119	zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp
[2528]	120	!
[7753]	121	zwi(:,:,:) = 0._wp
	122	!
[11407]	123	! If adaptive vertical advection, check if it is needed on this PE at this time
	124	IF( ln_zad_Aimp ) THEN
	125	IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
	126	END IF
	127	! If active adaptive vertical advection, build tridiagonal matrix
	128	IF( ll_zAimp ) THEN
	129	ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
[12377]	130	DO_3D_00_00( 1, jpkm1 )
	131	zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs)
	132	zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
	133	zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
	134	END_3D
[11407]	135	END IF
	136	!
[6140]	137	DO jn = 1, kjpt !== loop over the tracers ==!
[5770]	138	!
	139	! !== upstream advection with initial mass fluxes & intermediate update ==!
	140	! !* upstream tracer flux in the i and j direction
[12377]	141	DO_3D_10_10( 1, jpkm1 )
	142	! upstream scheme
	143	zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
	144	zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
	145	zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
	146	zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
	147	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj ,jk,jn,Kbb) )
	148	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) )
	149	END_3D
[5770]	150	! !* upstream tracer flux in the k direction *!
[12377]	151	DO_3D_11_11( 2, jpkm1 )
	152	zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
	153	zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
	154	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
	155	END_3D
[6140]	156	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)
[5770]	157	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
[12377]	158	DO_2D_11_11
	159	zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface
	160	END_2D
[5770]	161	ELSE ! no cavities: only at the ocean surface
[12377]	162	zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
[5770]	163	ENDIF
[5120]	164	ENDIF
[5770]	165	!
[12377]	166	DO_3D_00_00( 1, jpkm1 )
	167	! ! total intermediate advective trends
	168	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	169	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	170	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
	171	! ! update and guess with monotonic sheme
	172	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
	173	zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
	174	END_3D
[11407]	175
	176	IF ( ll_zAimp ) THEN
	177	CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
	178	!
[11411]	179	ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
[12377]	180	DO_3D_00_00( 2, jpkm1 )
	181	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
	182	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
	183	ztw(ji,jj,jk) = 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
	184	zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
	185	END_3D
	186	DO_3D_00_00( 1, jpkm1 )
	187	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) &
	188	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
	189	END_3D
[11407]	190	!
	191	END IF
[5770]	192	!
[7646]	193	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
[7753]	194	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
[2528]	195	END IF
[5770]	196	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[9019]	197	IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:)
[5770]	198	!
	199	! !== anti-diffusive flux : high order minus low order ==!
	200	!
[6140]	201	SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes
[5770]	202	!
[6140]	203	CASE( 2 ) !- 2nd order centered
[12377]	204	DO_3D_10_10( 1, jpkm1 )
	205	zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
	206	zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
	207	END_3D
[5770]	208	!
[6140]	209	CASE( 4 ) !- 4th order centered
[7753]	210	zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
	211	zltv(:,:,jpk) = 0._wp
[6140]	212	DO jk = 1, jpkm1 ! Laplacian
[12377]	213	DO_2D_10_10
	214	ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
	215	ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
	216	END_2D
	217	DO_2D_00_00
	218	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6
	219	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6
	220	END_2D
[503]	221	END DO
[10425]	222	CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
[5770]	223	!
[12377]	224	DO_3D_10_10( 1, jpkm1 )
	225	zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points
	226	zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
	227	! ! C4 minus upstream advective fluxes
	228	zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
	229	zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
	230	END_3D
[5770]	231	!
[6140]	232	CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
[7753]	233	ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
	234	ztv(:,:,jpk) = 0._wp
[12377]	235	DO_3D_10_10( 1, jpkm1 )
	236	ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
	237	ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
	238	END_3D
[10425]	239	CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn)
[5770]	240	!
[12377]	241	DO_3D_00_00( 1, jpkm1 )
	242	zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points (x2)
	243	zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
	244	! ! C4 interpolation of T at u- & v-points (x2)
	245	zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) )
	246	zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) )
	247	! ! C4 minus upstream advective fluxes
	248	zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
	249	zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
	250	END_3D
[5770]	251	!
	252	END SELECT
[6140]	253	!
	254	SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values)
[5770]	255	!
[6140]	256	CASE( 2 ) !- 2nd order centered
[12377]	257	DO_3D_00_00( 2, jpkm1 )
	258	zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) &
	259	& - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
	260	END_3D
[5770]	261	!
[6140]	262	CASE( 4 ) !- 4th order COMPACT
[12377]	263	CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point
	264	DO_3D_00_00( 2, jpkm1 )
	265	zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
	266	END_3D
[5770]	267	!
	268	END SELECT
[6140]	269	IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0
[7753]	270	zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked
[6140]	271	ENDIF
[11407]	272	!
	273	IF ( ll_zAimp ) THEN
[12377]	274	DO_3D_00_00( 1, jpkm1 )
	275	! ! total intermediate advective trends
	276	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	277	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	278	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
	279	ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
	280	END_3D
[11407]	281	!
[11411]	282	CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
[11407]	283	!
[12377]	284	DO_3D_00_00( 2, jpkm1 )
	285	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
	286	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
	287	zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk)
	288	END_3D
[11407]	289	END IF
[11411]	290	!
	291	CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. )
	292	!
[5770]	293	! !== monotonicity algorithm ==!
	294	!
[12377]	295	CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt )
[6140]	296	!
[5770]	297	! !== final trend with corrected fluxes ==!
	298	!
[12377]	299	DO_3D_00_00( 1, jpkm1 )
	300	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	301	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	302	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
	303	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
	304	zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
	305	END_3D
[5770]	306	!
[11407]	307	IF ( ll_zAimp ) THEN
	308	!
[11411]	309	ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
[12377]	310	DO_3D_00_00( 2, jpkm1 )
	311	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
	312	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
	313	ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
	314	zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
	315	END_3D
	316	DO_3D_00_00( 1, jpkm1 )
	317	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) &
	318	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
	319	END_3D
[11407]	320	END IF
	321	!
[9019]	322	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport
	323	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes
	324	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes
	325	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) !
[5770]	326	!
[9019]	327	IF( l_trd ) THEN ! trend diagnostics
[12377]	328	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
	329	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
	330	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
[9019]	331	ENDIF
	332	! ! heat/salt transport
	333	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
[5770]	334	!
[9019]	335	ENDIF
	336	IF( l_ptr ) THEN ! "Poleward" transports
	337	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes
[7646]	338	CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
[2528]	339	ENDIF
[503]	340	!
[6140]	341	END DO ! end of tracer loop
[503]	342	!
[11407]	343	IF ( ll_zAimp ) THEN
	344	DEALLOCATE( zwdia, zwinf, zwsup )
	345	ENDIF
[10024]	346	IF( l_trd .OR. l_hst ) THEN
	347	DEALLOCATE( ztrdx, ztrdy, ztrdz )
	348	ENDIF
	349	IF( l_ptr ) THEN
	350	DEALLOCATE( zptry )
	351	ENDIF
	352	!
[5770]	353	END SUBROUTINE tra_adv_fct
[3]	354
[5737]	355
[12377]	356	SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt )
[3]	357	!!---------------------------------------------------------------------
	358	!! * ROUTINE nonosc *
	359	!!
	360	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	361	!! scheme and the before field by a nonoscillatory algorithm
	362	!!
	363	!! ** Method : ... ???
	364	!! warning : pbef and paft must be masked, but the boundaries
	365	!! conditions on the fluxes are not necessary zalezak (1979)
	366	!! drange (1995) multi-dimensional forward-in-time and upstream-
	367	!! in-space based differencing for fluid
	368	!!----------------------------------------------------------------------
[12377]	369	INTEGER , INTENT(in ) :: Kmm ! time level index
[6140]	370	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[2528]	371	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
	372	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
[2715]	373	!
[4990]	374	INTEGER :: ji, jj, jk ! dummy loop indices
	375	INTEGER :: ikm1 ! local integer
[6140]	376	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
[2715]	377	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
[9019]	378	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
[3]	379	!!----------------------------------------------------------------------
[3294]	380	!
[2715]	381	zbig = 1.e+40_wp
	382	zrtrn = 1.e-15_wp
[7753]	383	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
[785]	384
[3]	385	! Search local extrema
	386	! --------------------
[785]	387	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
[4990]	388	zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), &
	389	& paft * tmask - zbig * ( 1._wp - tmask ) )
	390	zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), &
	391	& paft * tmask + zbig * ( 1._wp - tmask ) )
[785]	392
[5120]	393	DO jk = 1, jpkm1
	394	ikm1 = MAX(jk-1,1)
[12377]	395	DO_2D_00_00
[5120]	396
[12377]	397	! search maximum in neighbourhood
	398	zup = MAX( zbup(ji ,jj ,jk ), &
	399	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
	400	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
	401	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
[3]	402
[12377]	403	! search minimum in neighbourhood
	404	zdo = MIN( zbdo(ji ,jj ,jk ), &
	405	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
	406	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
	407	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
[3]	408
[12377]	409	! positive part of the flux
	410	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
	411	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
	412	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
[785]	413
[12377]	414	! negative part of the flux
	415	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
	416	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
	417	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
[785]	418
[12377]	419	! up & down beta terms
	420	zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
	421	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
	422	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
	423	END_2D
[3]	424	END DO
[10425]	425	CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
[3]	426
[237]	427	! 3. monotonic flux in the i & j direction (paa & pbb)
	428	! ----------------------------------------
[12377]	429	DO_3D_00_00( 1, jpkm1 )
	430	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
	431	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
	432	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
	433	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
[3]	434
[12377]	435	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
	436	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
	437	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
	438	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
[3]	439
[12377]	440	! monotonic flux in the k direction, i.e. pcc
	441	! -------------------------------------------
	442	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
	443	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
	444	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
	445	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
	446	END_3D
[10425]	447	CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
[503]	448	!
[3]	449	END SUBROUTINE nonosc
	450
[5770]	451
[7646]	452	SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
[5770]	453	!!----------------------------------------------------------------------
[7646]	454	!! * ROUTINE interp_4th_cpt_org *
[5770]	455	!!
	456	!! ** Purpose : Compute the interpolation of tracer at w-point
	457	!!
	458	!! ** Method : 4th order compact interpolation
	459	!!----------------------------------------------------------------------
	460	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields
	461	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts
	462	!
	463	INTEGER :: ji, jj, jk ! dummy loop integers
	464	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
	465	!!----------------------------------------------------------------------
	466
[12377]	467	DO_3D_11_11( 3, jpkm1 )
	468	zwd (ji,jj,jk) = 4._wp
	469	zwi (ji,jj,jk) = 1._wp
	470	zws (ji,jj,jk) = 1._wp
	471	zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	472	!
	473	IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom
[5770]	474	zwd (ji,jj,jk) = 1._wp
	475	zwi (ji,jj,jk) = 0._wp
	476	zws (ji,jj,jk) = 0._wp
[12377]	477	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	478	ENDIF
	479	END_3D
[5770]	480	!
[12377]	481	jk = 2 ! Switch to second order centered at top
	482	DO_2D_11_11
	483	zwd (ji,jj,jk) = 1._wp
	484	zwi (ji,jj,jk) = 0._wp
	485	zws (ji,jj,jk) = 0._wp
	486	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	487	END_2D
	488	!
[5770]	489	! !== tridiagonal solve ==!
[12377]	490	DO_2D_11_11
	491	zwt(ji,jj,2) = zwd(ji,jj,2)
	492	END_2D
	493	DO_3D_11_11( 3, jpkm1 )
	494	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	495	END_3D
[5770]	496	!
[12377]	497	DO_2D_11_11
	498	pt_out(ji,jj,2) = zwrm(ji,jj,2)
	499	END_2D
	500	DO_3D_11_11( 3, jpkm1 )
	501	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
	502	END_3D
[5770]	503
[12377]	504	DO_2D_11_11
	505	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
	506	END_2D
	507	DO_3DS_11_11( jpk-2, 2, -1 )
	508	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
	509	END_3D
[5770]	510	!
[7646]	511	END SUBROUTINE interp_4th_cpt_org
	512
	513
	514	SUBROUTINE interp_4th_cpt( pt_in, pt_out )
	515	!!----------------------------------------------------------------------
	516	!! * ROUTINE interp_4th_cpt *
	517	!!
	518	!! ** Purpose : Compute the interpolation of tracer at w-point
	519	!!
	520	!! ** Method : 4th order compact interpolation
	521	!!----------------------------------------------------------------------
	522	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point
	523	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point
	524	!
	525	INTEGER :: ji, jj, jk ! dummy loop integers
	526	INTEGER :: ikt, ikb ! local integers
	527	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
	528	!!----------------------------------------------------------------------
	529	!
	530	! !== build the three diagonal matrix & the RHS ==!
	531	!
[12377]	532	DO_3D_00_00( 3, jpkm1 )
	533	zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal
	534	zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal
	535	zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal
	536	zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS
	537	& * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
	538	END_3D
[7646]	539	!
	540	!!gm
	541	! SELECT CASE( kbc ) !* boundary condition
	542	! CASE( np_NH ) ! Neumann homogeneous at top & bottom
	543	! CASE( np_CEN2 ) ! 2nd order centered at top & bottom
	544	! END SELECT
	545	!!gm
	546	!
[9901]	547	IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case
	548	zwd(:,:,2) = 1._wp ; zwi(:,:,2) = 0._wp ; zws(:,:,2) = 0._wp ; zwrm(:,:,2) = 0._wp
	549	END IF
	550	!
[12377]	551	DO_2D_00_00
	552	ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point
	553	ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point
	554	!
	555	zwd (ji,jj,ikt) = 1._wp ! top
	556	zwi (ji,jj,ikt) = 0._wp
	557	zws (ji,jj,ikt) = 0._wp
	558	zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) )
	559	!
	560	zwd (ji,jj,ikb) = 1._wp ! bottom
	561	zwi (ji,jj,ikb) = 0._wp
	562	zws (ji,jj,ikb) = 0._wp
	563	zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )
	564	END_2D
[7646]	565	!
	566	! !== tridiagonal solver ==!
	567	!
[12377]	568	DO_2D_00_00
	569	zwt(ji,jj,2) = zwd(ji,jj,2)
	570	END_2D
	571	DO_3D_00_00( 3, jpkm1 )
	572	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	573	END_3D
[7646]	574	!
[12377]	575	DO_2D_00_00
	576	pt_out(ji,jj,2) = zwrm(ji,jj,2)
	577	END_2D
	578	DO_3D_00_00( 3, jpkm1 )
	579	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
	580	END_3D
[7646]	581
[12377]	582	DO_2D_00_00
	583	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
	584	END_2D
	585	DO_3DS_00_00( jpk-2, 2, -1 )
	586	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
	587	END_3D
[7646]	588	!
[5770]	589	END SUBROUTINE interp_4th_cpt
[7646]	590
	591
	592	SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
	593	!!----------------------------------------------------------------------
	594	!! * ROUTINE tridia_solver *
	595	!!
	596	!! ** Purpose : solve a symmetric 3diagonal system
	597	!!
	598	!! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk )
	599	!!
	600	!! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 )
	601	!! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 )
	602	!! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 )
	603	!! ( ... )( ... ) ( ... )
	604	!! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k )
	605	!!
	606	!! M is decomposed in the product of an upper and lower triangular matrix.
	607	!! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL
	608	!! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
	609	!! The solution is pta.
	610	!! The 3d array zwt is used as a work space array.
	611	!!----------------------------------------------------------------------
	612	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix
	613	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side
	614	REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)
	615	INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level
	616	! ! =0 pt at t-level
	617	INTEGER :: ji, jj, jk ! dummy loop integers
	618	INTEGER :: kstart ! local indices
	619	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array
	620	!!----------------------------------------------------------------------
	621	!
	622	kstart = 1 + klev
	623	!
[12377]	624	DO_2D_00_00
	625	zwt(ji,jj,kstart) = pD(ji,jj,kstart)
	626	END_2D
	627	DO_3D_00_00( kstart+1, jpkm1 )
	628	zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	629	END_3D
[7646]	630	!
[12377]	631	DO_2D_00_00
	632	pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
	633	END_2D
	634	DO_3D_00_00( kstart+1, jpkm1 )
	635	pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
	636	END_3D
[7646]	637
[12377]	638	DO_2D_00_00
	639	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
	640	END_2D
	641	DO_3DS_00_00( jpk-2, kstart, -1 )
	642	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
	643	END_3D
[7646]	644	!
	645	END SUBROUTINE tridia_solver
	646
[3]	647	!!======================================================================
[5770]	648	END MODULE traadv_fct

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/TRA/traadv_fct.F90 @ 12396

Download in other formats: