Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynadv_ubs.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 13.0 KB

Rev	Line
[643]	1	MODULE dynadv_ubs
	2	!!======================================================================
	3	!! * MODULE dynadv_ubs *
	4	!! Ocean dynamics: Update the momentum trend with the flux form advection
	5	!! trend using a 3rd order upstream biased scheme
	6	!!======================================================================
[1566]	7	!! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
	8	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
[643]	9	!!----------------------------------------------------------------------
	10
	11	!!----------------------------------------------------------------------
	12	!! dyn_adv_ubs : flux form momentum advection using (ln_dynadv=T)
	13	!! an 3rd order Upstream Biased Scheme or Quick scheme
	14	!! combined with 2nd or 4th order finite differences
	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and tracers
	17	USE dom_oce ! ocean space and time domain
[5062]	18	USE trd_oce ! trends: ocean variables
	19	USE trddyn ! trend manager: dynamics
	20	!
[2715]	21	USE in_out_manager ! I/O manager
[1129]	22	USE prtctl ! Print control
[2715]	23	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	24	USE lib_mpp ! MPP library
[5062]	25	USE timing ! Timing
[643]	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
[3294]	30	REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
[4153]	31	REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0 2nd order ; =1/32 4th order centred
[643]	32
[1566]	33	PUBLIC dyn_adv_ubs ! routine called by step.F90
[643]	34
	35	!! * Substitutions
	36	# include "vectopt_loop_substitute.h90"
	37	!!----------------------------------------------------------------------
[2715]	38	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]	39	!! $Id$
[2715]	40	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[643]	41	!!----------------------------------------------------------------------
	42	CONTAINS
	43
	44	SUBROUTINE dyn_adv_ubs( kt )
	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE dyn_adv_ubs *
	47	!!
	48	!! ** Purpose : Compute the now momentum advection trend in flux form
[1566]	49	!! and the general trend of the momentum equation.
[643]	50	!!
	51	!! ** Method : The scheme is the one implemeted in ROMS. It depends
	52	!! on two parameter gamma1 and gamma2. The former control the
	53	!! upstream baised part of the scheme and the later the centred
	54	!! part: gamma1 = 0 pure centered (no diffusive part)
	55	!! = 1/4 Quick scheme
	56	!! = 1/3 3rd order Upstream biased scheme
	57	!! gamma2 = 0 2nd order finite differencing
[4153]	58	!! = 1/32 4th order finite differencing
[643]	59	!! For stability reasons, the first term of the fluxes which cor-
	60	!! responds to a second order centered scheme is evaluated using
	61	!! the now velocity (centered in time) while the second term which
	62	!! is the diffusive part of the scheme, is evaluated using the
	63	!! before velocity (forward in time).
	64	!! Default value (hard coded in the begining of the module) are
[4153]	65	!! gamma1=1/3 and gamma2=1/32.
[643]	66	!!
[1566]	67	!! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]	68	!!
	69	!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
	70	!!----------------------------------------------------------------------
[6140]	71	INTEGER, INTENT(in) :: kt ! ocean time-step index
[2715]	72	!
[6140]	73	INTEGER :: ji, jj, jk ! dummy loop indices
	74	REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars
[7910]	75	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfv
	76	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw
	77	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zlu_uu, zlv_vv, zlu_uv, zlv_vu
[643]	78	!!----------------------------------------------------------------------
[3294]	79	!
	80	IF( nn_timing == 1 ) CALL timing_start('dyn_adv_ubs')
	81	!
	82	!
[643]	83	IF( kt == nit000 ) THEN
	84	IF(lwp) WRITE(numout,*)
	85	IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
	86	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	87	ENDIF
[3294]	88	!
[2715]	89	zfu_t(:,:,:) = 0._wp
	90	zfv_t(:,:,:) = 0._wp
	91	zfu_f(:,:,:) = 0._wp
	92	zfv_f(:,:,:) = 0._wp
[1566]	93	!
[2715]	94	zlu_uu(:,:,:,:) = 0._wp
	95	zlv_vv(:,:,:,:) = 0._wp
	96	zlu_uv(:,:,:,:) = 0._wp
	97	zlv_vu(:,:,:,:) = 0._wp
[6140]	98	!
	99	IF( l_trddyn ) THEN ! trends: store the input trends
[1129]	100	zfu_uw(:,:,:) = ua(:,:,:)
	101	zfv_vw(:,:,:) = va(:,:,:)
	102	ENDIF
[1566]	103	! ! =========================== !
	104	DO jk = 1, jpkm1 ! Laplacian of the velocity !
	105	! ! =========================== !
	106	! ! horizontal volume fluxes
[6140]	107	zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
	108	zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
[1566]	109	!
	110	DO jj = 2, jpjm1 ! laplacian
[643]	111	DO ji = fs_2, fs_jpim1 ! vector opt.
[5069]	112	zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj ,jk) - 2.ub (ji,jj,jk) + ub (ji-1,jj ,jk) ) umask(ji,jj,jk)
	113	zlv_vv(ji,jj,jk,1) = ( vb (ji ,jj+1,jk) - 2.vb (ji,jj,jk) + vb (ji ,jj-1,jk) ) vmask(ji,jj,jk)
	114	zlu_uv(ji,jj,jk,1) = ( ub (ji ,jj+1,jk) - ub (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
	115	& - ( ub (ji ,jj ,jk) - ub (ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk)
	116	zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj ,jk) - vb (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
	117	& - ( vb (ji ,jj ,jk) - vb (ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk)
[2715]	118	!
[5069]	119	zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj ,jk) - 2.zfu(ji,jj,jk) + zfu(ji-1,jj ,jk) ) umask(ji,jj,jk)
	120	zlv_vv(ji,jj,jk,2) = ( zfv(ji ,jj+1,jk) - 2.zfv(ji,jj,jk) + zfv(ji ,jj-1,jk) ) vmask(ji,jj,jk)
	121	zlu_uv(ji,jj,jk,2) = ( zfu(ji ,jj+1,jk) - zfu(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
	122	& - ( zfu(ji ,jj ,jk) - zfu(ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk)
	123	zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj ,jk) - zfv(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
	124	& - ( zfv(ji ,jj ,jk) - zfv(ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk)
[643]	125	END DO
	126	END DO
[1566]	127	END DO
	128	CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
	129	CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
	130	CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
	131	CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. )
[6140]	132	!
[1566]	133	! ! ====================== !
	134	! ! Horizontal advection !
	135	DO jk = 1, jpkm1 ! ====================== !
	136	! ! horizontal volume fluxes
[6140]	137	zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
	138	zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
[1566]	139	!
	140	DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point
[643]	141	DO ji = 1, fs_jpim1 ! vector opt.
	142	zui = ( un(ji,jj,jk) + un(ji+1,jj ,jk) )
	143	zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) )
[1566]	144	!
[6140]	145	IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk,1)
	146	ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1)
[643]	147	ENDIF
[6140]	148	IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,jk,1)
	149	ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1)
[643]	150	ENDIF
[1566]	151	!
[643]	152	zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) &
	153	& - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj ,jk,2) ) ) &
	154	& * ( zui - gamma1 * zl_u)
	155	zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji ,jj+1,jk) &
	156	& - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji ,jj+1,jk,2) ) ) &
	157	& * ( zvj - gamma1 * zl_v)
[1566]	158	!
[643]	159	zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) )
	160	zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) )
[6140]	161	IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1)
	162	ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1)
[643]	163	ENDIF
[6140]	164	IF( zfvi > 0 ) THEN ; zl_u = zlu_uv( ji,jj ,jk,1)
	165	ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1)
[643]	166	ENDIF
[1566]	167	!
[643]	168	zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) &
	169	& * ( un(ji,jj,jk) + un(ji ,jj+1,jk) - gamma1 * zl_u )
	170	zfu_f(ji ,jj ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji ,jj+1,jk,2) ) ) &
	171	& * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) - gamma1 * zl_v )
	172	END DO
	173	END DO
[1566]	174	DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes
[643]	175	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	176	ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
	177	& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
	178	va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
	179	& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[643]	180	END DO
	181	END DO
[1566]	182	END DO
[6140]	183	IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic
[1129]	184	zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
	185	zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
[5062]	186	CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
[1129]	187	zfu_t(:,:,:) = ua(:,:,:)
	188	zfv_t(:,:,:) = va(:,:,:)
	189	ENDIF
[1566]	190	! ! ==================== !
	191	! ! Vertical advection !
[6140]	192	! ! ==================== !
	193	DO jj = 2, jpjm1 ! surface/bottom advective fluxes set to zero
	194	DO ji = fs_2, fs_jpim1
	195	zfu_uw(ji,jj,jpk) = 0._wp
	196	zfv_vw(ji,jj,jpk) = 0._wp
	197	zfu_uw(ji,jj, 1 ) = 0._wp
	198	zfv_vw(ji,jj, 1 ) = 0._wp
	199	END DO
	200	END DO
	201	IF( ln_linssh ) THEN ! constant volume : advection through the surface
	202	DO jj = 2, jpjm1
	203	DO ji = fs_2, fs_jpim1
	204	zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1)
	205	zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1)
[643]	206	END DO
[6140]	207	END DO
	208	ENDIF
	209	DO jk = 2, jpkm1 ! interior fluxes
[6750]	210	DO jj = 2, jpj
	211	DO ji = 2, jpi
[6140]	212	zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
	213	END DO
	214	END DO
	215	DO jj = 2, jpjm1
	216	DO ji = fs_2, fs_jpim1 ! vector opt.
	217	zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
	218	zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
	219	END DO
	220	END DO
[643]	221	END DO
[6140]	222	DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence
	223	DO jj = 2, jpjm1
[643]	224	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	225	ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
	226	va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[643]	227	END DO
	228	END DO
	229	END DO
[1566]	230	!
[6140]	231	IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
[1129]	232	zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
	233	zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
[5062]	234	CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
[1129]	235	ENDIF
[6140]	236	! ! Control print
[1129]	237	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
	238	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]	239	!
[2715]	240	!
[3294]	241	IF( nn_timing == 1 ) CALL timing_stop('dyn_adv_ubs')
	242	!
[643]	243	END SUBROUTINE dyn_adv_ubs
	244
	245	!!==============================================================================
	246	END MODULE dynadv_ubs

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 7910

Download in other formats: