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 @ 881

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

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

Download in other formats: