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.

traadv_cen2.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 785

Last change on this file since 785 was 719, checked in by ctlod, 17 years ago
get back to the nemo_v2_3 version for trunk
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.6 KB

Rev	Line
[3]	1	MODULE traadv_cen2
[719]	2	!!==============================================================================
	3	!! * MODULE traadv_cen2 *
[3]	4	!! Ocean active tracers: horizontal & vertical advective trend
[719]	5	!!==============================================================================
	6	!! History : 8.2 ! 01-08 (G. Madec, E. Durand) trahad+trazad = traadv
	7	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
	8	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
	9	!! " " ! 06-04 (R. Benshila, G. Madec) Step reorganization
[3]	10	!!----------------------------------------------------------------------
[503]	11
	12	!!----------------------------------------------------------------------
[457]	13	!! tra_adv_cen2 : update the tracer trend with the horizontal and
	14	!! vertical advection trends using a seconder order
	15	!! centered scheme. (k-j-i loops)
[3]	16	!!----------------------------------------------------------------------
	17	USE oce ! ocean dynamics and active tracers
	18	USE dom_oce ! ocean space and time domain
[719]	19	USE trdmod ! ocean active tracers trends
[708]	20	USE trdmod_oce ! ocean variables trends
[719]	21	USE flxrnf !
[3]	22	USE trabbl ! advective term in the BBL
[74]	23	USE ocfzpt !
[3]	24	USE lib_mpp
[74]	25	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[719]	26	USE in_out_manager ! I/O manager
[132]	27	USE diaptr ! poleward transport diagnostics
[719]	28	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[258]	29	USE prtctl ! Print control
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[719]	34	PUBLIC tra_adv_cen2 ! routine called by step.F90
[3]	35
[503]	36	REAL(wp), DIMENSION(jpi,jpj) :: btr2 ! inverse of T-point surface [1/(e1t*e2t)]
	37
[3]	38	!! * Substitutions
	39	# include "domzgr_substitute.h90"
	40	# include "vectopt_loop_substitute.h90"
	41	!!----------------------------------------------------------------------
[719]	42	!! OPA 9.0 , LOCEAN-IPSL (2005)
	43	!! $Header$
[503]	44	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	45	!!----------------------------------------------------------------------
	46
	47	CONTAINS
	48
[457]	49	SUBROUTINE tra_adv_cen2( kt, pun, pvn, pwn )
[3]	50	!!----------------------------------------------------------------------
	51	!! * ROUTINE tra_adv_cen2 *
	52	!!
	53	!! ** Purpose : Compute the now trend due to the advection of tracers
	54	!! and add it to the general trend of passive tracer equations.
	55	!!
	56	!! ** Method : The advection is evaluated by a second order centered
	57	!! scheme using now fields (leap-frog scheme). In specific areas
	58	!! (vicinity of major river mouths, some straits, or where tn is
[457]	59	!! approaching the freezing point) it is mixed with an upstream
[3]	60	!! scheme for stability reasons.
[457]	61	!! Part 0 : compute the upstream / centered flag
	62	!! (3D array, zind, defined at T-point (0<zind<1))
	63	!! Part I : horizontal advection
	64	!! * centered flux:
[3]	65	!! zcenu = e2u*e3u un mi(tn)
	66	!! zcenv = e1v*e3v vn mj(tn)
[457]	67	!! * upstream flux:
[3]	68	!! zupsu = e2u*e3u un (tb(i) or tb(i-1) ) [un>0 or <0]
	69	!! zupsv = e1v*e3v vn (tb(j) or tb(j-1) ) [vn>0 or <0]
[457]	70	!! * mixed upstream / centered horizontal advection scheme
[3]	71	!! zcofi = max(zind(i+1), zind(i))
	72	!! zcofj = max(zind(j+1), zind(j))
	73	!! zwx = zcofi * zupsu + (1-zcofi) * zcenu
	74	!! zwy = zcofj * zupsv + (1-zcofj) * zcenv
[457]	75	!! * horizontal advective trend (divergence of the fluxes)
[3]	76	!! zta = 1/(e1te2te3t) { di-1[zwx] + dj-1[zwy] }
[457]	77	!! * Add this trend now to the general trend of tracer (ta,sa):
[3]	78	!! (ta,sa) = (ta,sa) + ( zta , zsa )
[457]	79	!! * trend diagnostic ('key_trdtra' defined): the trend is
	80	!! saved for diagnostics. The trends saved is expressed as
	81	!! Uh.gradh(T), i.e.
	82	!! save trend = zta + tn divn
[3]	83	!! In addition, the advective trend in the two horizontal direc-
	84	!! tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is
	85	!! equal to (in s-coordinates, and similarly in z-coord.):
	86	!! zta+tndivn=1/(e1te2te3t) { mi-1( e2ue3u un di[tn] )
	87	!! +mj-1( e1v*e3v vn mj[tn] ) }
[457]	88	!! NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so
	89	!! they vanish from the expression of the flux and divergence.
[3]	90	!!
	91	!! Part II : vertical advection
	92	!! For temperature (idem for salinity) the advective trend is com-
	93	!! puted as follows :
	94	!! zta = 1/e3t dk+1[ zwz ]
	95	!! where the vertical advective flux, zwz, is given by :
	96	!! zwz = zcofk * zupst + (1-zcofk) * zcent
[457]	97	!! with
[3]	98	!! zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
	99	!! zcenu = centered flux = wn * mk(tn)
[457]	100	!! The surface boundary condition is :
	101	!! rigid-lid (lk_dynspg_frd = T) : zero advective flux
	102	!! free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1)
[3]	103	!! Add this trend now to the general trend of tracer (ta,sa):
	104	!! (ta,sa) = (ta,sa) + ( zta , zsa )
[457]	105	!! Trend diagnostic ('key_trdtra' defined): the trend is
	106	!! saved for diagnostics. The trends saved is expressed as :
[3]	107	!! save trend = w.gradz(T) = zta - tn divn.
	108	!!
[457]	109	!! ** Action : - update (ta,sa) with the now advective tracer trends
[503]	110	!! - save trends in (ztrdt,ztrds) ('key_trdtra')
	111	!!----------------------------------------------------------------------
	112	USE oce, ONLY : zwx => ua ! use ua as workspace
	113	USE oce, ONLY : zwy => va ! use va as workspace
[3]	114	!!
[503]	115	INTEGER , INTENT(in) :: kt ! ocean time-step index
	116	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component
	117	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component
	118	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component
	119	!!
	120	INTEGER :: ji, jj, jk ! dummy loop indices
[719]	121	REAL(wp) :: &
	122	zbtr, zta, zsa, zfui, zfvj, & ! temporary scalars
	123	zhw, ze3tr, zcofi, zcofj, & ! " "
	124	zupsut, zupsvt, zupsus, zupsvs, & ! " "
	125	zfp_ui, zfp_vj, zfm_ui, zfm_vj, & ! " "
	126	zcofk, zupst, zupss, zcent, & ! " "
	127	zcens, zfp_w, zfm_w, & ! " "
	128	zcenut, zcenvt, zcenus, zcenvs, & ! " "
	129	z_hdivn_x, z_hdivn_y, z_hdivn
	130	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz, ztrdt, zind ! 3D workspace
[503]	131	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zww, ztrds ! " "
[3]	132	!!----------------------------------------------------------------------
	133
	134	IF( kt == nit000 ) THEN
	135	IF(lwp) WRITE(numout,*)
	136	IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
	137	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ Vector optimization case'
	138	IF(lwp) WRITE(numout,*)
[719]	139	!
	140	btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
[3]	141	ENDIF
	142
	143	! Upstream / centered scheme indicator
	144	! ------------------------------------
	145	DO jk = 1, jpk
	146	DO jj = 1, jpj
	147	DO ji = 1, jpi
[719]	148	zind(ji,jj,jk) = MAX ( upsrnfh(ji,jj) * upsrnfz(jk), & ! changing advection scheme near runoff
	149	& upsadv(ji,jj) & ! in the vicinity of some straits
[3]	150	#if defined key_ice_lim
[719]	151	& , tmask(ji,jj,jk) & ! half upstream tracer fluxes
	152	& * MAX( 0., SIGN( 1., fzptn(ji,jj) & ! if tn < ("freezing"+0.1 )
	153	& +0.1-tn(ji,jj,jk) ) ) &
[3]	154	#endif
	155	& )
	156	END DO
	157	END DO
	158	END DO
	159
[719]	160
	161	! Horizontal advective fluxes
	162	! -----------------------------
[3]	163	! ! ===============
	164	DO jk = 1, jpkm1 ! Horizontal slab
	165	! ! ===============
	166	DO jj = 1, jpjm1
	167	DO ji = 1, fs_jpim1 ! vector opt.
	168	! upstream indicator
	169	zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
	170	zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
	171	! volume fluxes * 1/2
[457]	172	#if defined key_zco
	173	zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
	174	zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
[3]	175	#else
[457]	176	zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
	177	zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
[3]	178	#endif
	179	! upstream scheme
	180	zfp_ui = zfui + ABS( zfui )
	181	zfp_vj = zfvj + ABS( zfvj )
	182	zfm_ui = zfui - ABS( zfui )
	183	zfm_vj = zfvj - ABS( zfvj )
	184	zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj ,jk)
	185	zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji ,jj+1,jk)
	186	zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj ,jk)
	187	zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji ,jj+1,jk)
	188	! centered scheme
	189	zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj ,jk) )
	190	zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji ,jj+1,jk) )
	191	zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj ,jk) )
	192	zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji ,jj+1,jk) )
	193	! mixed centered / upstream scheme
	194	zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
	195	zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
	196	zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
	197	zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
	198	END DO
	199	END DO
	200
[503]	201	! Tracer flux divergence at t-point added to the general trend
	202	! --------------------------------------------------------------
[3]	203	DO jj = 2, jpjm1
	204	DO ji = fs_2, fs_jpim1 ! vector opt.
[457]	205	#if defined key_zco
[503]	206	zbtr = btr2(ji,jj)
[457]	207	#else
[503]	208	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
[3]	209	#endif
[719]	210	! horizontal advective trends
	211	zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) &
[3]	212	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) )
	213	zsa = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) &
	214	& + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) )
[719]	215	! add it to the general tracer trends
	216	ta(ji,jj,jk) = ta(ji,jj,jk) + zta
[3]	217	sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
	218	END DO
	219	END DO
	220	! ! ===============
	221	END DO ! End of slab
	222	! ! ===============
	223
[503]	224	! Save the horizontal advective trends for diagnostic
	225	! -----------------------------------------------------
	226	IF( l_trdtra ) THEN
	227	! T/S ZONAL advection trends
	228	ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0
	229	!
	230	DO jk = 1, jpkm1
	231	DO jj = 2, jpjm1
	232	DO ji = fs_2, fs_jpim1 ! vector opt.
	233	!-- Compute zonal divergence by splitting hdivn (see divcur.F90)
	234	! N.B. This computation is not valid along OBCs (if any)
	235	#if defined key_zco
	236	zbtr = btr2(ji,jj)
	237	z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) &
	238	& - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
	239	#else
	240	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	241	z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) &
	242	& - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
	243	#endif
	244	ztrdt(ji,jj,jk) = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x
	245	ztrds(ji,jj,jk) = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x
	246	END DO
	247	END DO
	248	END DO
	249	CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
	250	!
	251	! T/S MERIDIONAL advection trends
	252	DO jk = 1, jpkm1
	253	DO jj = 2, jpjm1
	254	DO ji = fs_2, fs_jpim1 ! vector opt.
	255	!-- Compute merid. divergence by splitting hdivn (see divcur.F90)
	256	! N.B. This computation is not valid along OBCs (if any)
	257	#if defined key_zco
	258	zbtr = btr2(ji,jj)
	259	z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) &
	260	& - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
	261	#else
	262	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	263	z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) &
	264	& - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
	265	#endif
	266	ztrdt(ji,jj,jk) = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y
	267	ztrds(ji,jj,jk) = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y
	268	END DO
	269	END DO
	270	END DO
	271	CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
	272	!
	273	! Save the horizontal up-to-date ta/sa trends
	274	ztrdt(:,:,:) = ta(:,:,:)
	275	ztrds(:,:,:) = sa(:,:,:)
[216]	276	ENDIF
	277
[457]	278	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 had - Ta: ', mask1=tmask, &
	279	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
[3]	280
[719]	281	! 4. "zonal" mean advective heat and salt transport
	282	! -------------------------------------------------
	283
[132]	284	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
[457]	285	IF( lk_zco ) THEN
	286	DO jk = 1, jpkm1
	287	DO jj = 2, jpjm1
	288	DO ji = fs_2, fs_jpim1 ! vector opt.
	289	zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
	290	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
	291	END DO
[3]	292	END DO
	293	END DO
[457]	294	ENDIF
[132]	295	pht_adv(:) = ptr_vj( zwy(:,:,:) )
	296	pst_adv(:) = ptr_vj( zwz(:,:,:) )
[3]	297	ENDIF
	298
	299	! II. Vertical advection
	300	! ----------------------
	301
	302	! Bottom value : flux set to zero
	303	zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0
	304
	305	! Surface value
[592]	306	IF( lk_dynspg_rl .OR. lk_vvl ) THEN
	307	! rigid lid or variable volume: flux set to zero
[359]	308	zwx(:,:, 1 ) = 0.e0 ; zwy(:,:, 1 ) = 0.e0
	309	ELSE
	310	! free surface
[457]	311	zwx(:,:, 1 ) = pwn(:,:,1) * tn(:,:,1)
	312	zwy(:,:, 1 ) = pwn(:,:,1) * sn(:,:,1)
[200]	313	ENDIF
[3]	314
[719]	315	! 1. Vertical advective fluxes
[3]	316	! ----------------------------
[719]	317	! Second order centered tracer flux at w-point
[3]	318	DO jk = 2, jpk
	319	DO jj = 2, jpjm1
	320	DO ji = fs_2, fs_jpim1 ! vector opt.
[719]	321	! upstream indicator
	322	zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
	323	! velocity * 1/2
	324	zhw = 0.5 * pwn(ji,jj,jk)
	325	! upstream scheme
	326	zfp_w = zhw + ABS( zhw )
[3]	327	zfm_w = zhw - ABS( zhw )
	328	zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
	329	zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
[719]	330	! centered scheme
	331	zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
[3]	332	zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
[719]	333	! mixed centered / upstream scheme
	334	zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
[3]	335	zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
	336	END DO
	337	END DO
	338	END DO
	339
	340	! 2. Tracer flux divergence at t-point added to the general trend
	341	! -------------------------
	342	DO jk = 1, jpkm1
	343	DO jj = 2, jpjm1
	344	DO ji = fs_2, fs_jpim1 ! vector opt.
	345	ze3tr = 1. / fse3t(ji,jj,jk)
[719]	346	! vertical advective trends
	347	zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
[3]	348	zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) )
[719]	349	! add it to the general tracer trends
	350	ta(ji,jj,jk) = ta(ji,jj,jk) + zta
[3]	351	sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
	352	END DO
	353	END DO
	354	END DO
	355
[216]	356	! 3. Save the vertical advective trends for diagnostic
	357	! ----------------------------------------------------
	358	IF( l_trdtra ) THEN
	359	! Recompute the vertical advection zta & zsa trends computed
	360	! at the step 2. above in making the difference between the new
[503]	361	! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
[216]	362	! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
	363
[503]	364	DO jk = 1, jpkm1
	365	DO jj = 2, jpjm1
	366	DO ji = fs_2, fs_jpim1 ! vector opt.
	367	#if defined key_zco
	368	zbtr = btr2(ji,jj)
	369	z_hdivn_x = e2u(ji,jj)pun(ji,jj,jk) - e2u(ji-1,jj)pun(ji-1,jj,jk)
	370	z_hdivn_y = e1v(ji,jj)pvn(ji,jj,jk) - e1v(ji,jj-1)pvn(ji,jj-1,jk)
	371	#else
	372	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	373	z_hdivn_x = e2u(ji,jj)fse3u(ji,jj,jk)pun(ji,jj,jk) - e2u(ji-1,jj)fse3u(ji-1,jj,jk)pun(ji-1,jj,jk)
	374	z_hdivn_y = e1v(ji,jj)fse3v(ji,jj,jk)pvn(ji,jj,jk) - e1v(ji,jj-1)fse3v(ji,jj-1,jk)pvn(ji,jj-1,jk)
	375	#endif
	376	z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr
	377	ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
	378	ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
	379	END DO
	380	END DO
	381	END DO
	382	CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
[216]	383	ENDIF
	384
[457]	385	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 zad - Ta: ', mask1=tmask, &
	386	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
[503]	387	!
[3]	388	END SUBROUTINE tra_adv_cen2
	389
	390	!!======================================================================
	391	END MODULE traadv_cen2

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

Download in other formats: