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.

dynzad.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 7710

Last change on this file since 7710 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 13.7 KB

Rev	Line
[3]	1	MODULE dynzad
	2	!!======================================================================
	3	!! * MODULE dynzad *
	4	!! Ocean dynamics : vertical advection trend
	5	!!======================================================================
[2715]	6	!! History : OPA ! 1991-01 (G. Madec) Original code
	7	!! 7.0 ! 1991-11 (G. Madec)
	8	!! 7.5 ! 1996-01 (G. Madec) statement function for e3
	9	!! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90
[503]	10	!!----------------------------------------------------------------------
[3]	11
	12	!!----------------------------------------------------------------------
[503]	13	!! dyn_zad : vertical advection momentum trend
[3]	14	!!----------------------------------------------------------------------
[503]	15	USE oce ! ocean dynamics and tracers
	16	USE dom_oce ! ocean space and time domain
[888]	17	USE sbc_oce ! surface boundary condition: ocean
[4990]	18	USE trd_oce ! trends: ocean variables
	19	USE trddyn ! trend manager: dynamics
	20	!
[719]	21	USE in_out_manager ! I/O manager
[4990]	22	USE lib_mpp ! MPP library
[503]	23	USE prtctl ! Print control
[4990]	24	USE wrk_nemo ! Memory Allocation
	25	USE timing ! Timing
[3]	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
[4990]	30	PUBLIC dyn_zad ! routine called by dynadv.F90
	31	PUBLIC dyn_zad_zts ! routine called by dynadv.F90
[3]	32
	33	!! * Substitutions
	34	# include "vectopt_loop_substitute.h90"
	35	!!----------------------------------------------------------------------
[2528]	36	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	37	!! $Id$
[2715]	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	39	!!----------------------------------------------------------------------
	40	CONTAINS
	41
	42	SUBROUTINE dyn_zad ( kt )
	43	!!----------------------------------------------------------------------
	44	!! * ROUTINE dynzad *
	45	!!
	46	!! ** Purpose : Compute the now vertical momentum advection trend and
	47	!! add it to the general trend of momentum equation.
	48	!!
	49	!! ** Method : The now vertical advection of momentum is given by:
[5836]	50	!! w dz(u) = ua + 1/(e1e2ue3u) mk+1[ mi(e1e2twn) dk(un) ]
	51	!! w dz(v) = va + 1/(e1e2ve3v) mk+1[ mj(e1e2twn) dk(vn) ]
[3]	52	!! Add this trend to the general trend (ua,va):
	53	!! (ua,va) = (ua,va) + w dz(u,v)
	54	!!
	55	!! ** Action : - Update (ua,va) with the vert. momentum adv. trends
[4990]	56	!! - Send the trends to trddyn for diagnostics (l_trddyn=T)
[3294]	57	!!----------------------------------------------------------------------
[503]	58	INTEGER, INTENT(in) :: kt ! ocean time-step inedx
[2715]	59	!
[503]	60	INTEGER :: ji, jj, jk ! dummy loop indices
	61	REAL(wp) :: zua, zva ! temporary scalars
[3294]	62	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw
	63	REAL(wp), POINTER, DIMENSION(:,: ) :: zww
	64	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
[3]	65	!!----------------------------------------------------------------------
[3294]	66	!
	67	IF( nn_timing == 1 ) CALL timing_start('dyn_zad')
	68	!
	69	CALL wrk_alloc( jpi,jpj, zww )
	70	CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw )
	71	!
[3]	72	IF( kt == nit000 ) THEN
	73	IF(lwp)WRITE(numout,*)
	74	IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
	75	ENDIF
[216]	76
[503]	77	IF( l_trddyn ) THEN ! Save ua and va trends
[3294]	78	CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
[7698]	79	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
	80	DO jk = 1, jpk
	81	DO jj = 1, jpj
	82	DO ji = 1, jpi
	83	ztrdu(ji,jj,jk) = ua(ji,jj,jk)
	84	ztrdv(ji,jj,jk) = va(ji,jj,jk)
	85	END DO
	86	END DO
	87	END DO
[216]	88	ENDIF
[3]	89
[7698]	90	!$OMP PARALLEL
[503]	91	DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical
[7698]	92	!$OMP DO schedule(static) private(jj, ji)
[503]	93	DO jj = 2, jpj ! vertical fluxes
	94	DO ji = fs_2, jpi ! vector opt.
[5836]	95	zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
[3]	96	END DO
	97	END DO
[7698]	98	!$OMP DO schedule(static) private(jj, ji)
[503]	99	DO jj = 2, jpjm1 ! vertical momentum advection at w-point
	100	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	101	zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )
	102	zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )
[3]	103	END DO
	104	END DO
	105	END DO
[7698]	106	!$OMP END PARALLEL
[5120]	107	!
	108	! Surface and bottom advective fluxes set to zero
	109	IF ( ln_isfcav ) THEN
[7698]	110	!$OMP PARALLEL DO schedule(static) private(jj, ji)
[5120]	111	DO jj = 2, jpjm1
	112	DO ji = fs_2, fs_jpim1 ! vector opt.
	113	zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
	114	zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
	115	zwuw(ji,jj,jpk) = 0._wp
	116	zwvw(ji,jj,jpk) = 0._wp
	117	END DO
	118	END DO
	119	ELSE
[7698]	120	!$OMP PARALLEL DO schedule(static) private(jj, ji)
[5120]	121	DO jj = 2, jpjm1
	122	DO ji = fs_2, fs_jpim1 ! vector opt.
	123	zwuw(ji,jj, 1 ) = 0._wp
	124	zwvw(ji,jj, 1 ) = 0._wp
	125	zwuw(ji,jj,jpk) = 0._wp
	126	zwvw(ji,jj,jpk) = 0._wp
	127	END DO
	128	END DO
	129	END IF
[3]	130
[7698]	131	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zua, zva)
[503]	132	DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points
[3]	133	DO jj = 2, jpjm1
[503]	134	DO ji = fs_2, fs_jpim1 ! vector opt.
	135	! ! vertical momentum advective trends
[6140]	136	zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
	137	zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[503]	138	! ! add the trends to the general momentum trends
[3]	139	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	140	va(ji,jj,jk) = va(ji,jj,jk) + zva
	141	END DO
	142	END DO
	143	END DO
	144
[503]	145	IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic
[7698]	146	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
	147	DO jk = 1, jpk
	148	DO jj = 1, jpj
	149	DO ji = 1, jpi
	150	ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
	151	ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
	152	END DO
	153	END DO
	154	END DO
[4990]	155	CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
[3294]	156	CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
[216]	157	ENDIF
[503]	158	! ! Control print
	159	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, &
	160	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	161	!
[3294]	162	CALL wrk_dealloc( jpi,jpj, zww )
	163	CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw )
[2715]	164	!
[3294]	165	IF( nn_timing == 1 ) CALL timing_stop('dyn_zad')
	166	!
[3]	167	END SUBROUTINE dyn_zad
	168
[5836]	169
[4990]	170	SUBROUTINE dyn_zad_zts ( kt )
	171	!!----------------------------------------------------------------------
	172	!! * ROUTINE dynzad_zts *
	173	!!
	174	!! ** Purpose : Compute the now vertical momentum advection trend and
	175	!! add it to the general trend of momentum equation. This version
	176	!! uses sub-timesteps for improved numerical stability with small
	177	!! vertical grid sizes. This is especially relevant when using
	178	!! embedded ice with thin surface boxes.
	179	!!
	180	!! ** Method : The now vertical advection of momentum is given by:
	181	!! w dz(u) = ua + 1/(e1ue2ue3u) mk+1[ mi(e1te2twn) dk(un) ]
	182	!! w dz(v) = va + 1/(e1ve2ve3v) mk+1[ mj(e1te2twn) dk(vn) ]
	183	!! Add this trend to the general trend (ua,va):
	184	!! (ua,va) = (ua,va) + w dz(u,v)
	185	!!
	186	!! ** Action : - Update (ua,va) with the vert. momentum adv. trends
	187	!! - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
	188	!!----------------------------------------------------------------------
	189	INTEGER, INTENT(in) :: kt ! ocean time-step inedx
	190	!
	191	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	192	INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection
	193	INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps
	194	REAL(wp) :: zua, zva ! temporary scalars
	195	REAL(wp) :: zr_rdt ! temporary scalar
	196	REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection
	197	REAL(wp) :: zts ! length of sub-timestep for vertical advection
	198	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww
	199	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
	200	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs
	201	!!----------------------------------------------------------------------
	202	!
	203	IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts')
	204	!
[5836]	205	CALL wrk_alloc( jpi,jpj,jpk, zwuw, zwvw, zww )
	206	CALL wrk_alloc( jpi,jpj,jpk,3, zus , zvs )
[4990]	207	!
	208	IF( kt == nit000 ) THEN
	209	IF(lwp)WRITE(numout,*)
	210	IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
	211	ENDIF
	212
	213	IF( l_trddyn ) THEN ! Save ua and va trends
	214	CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
	215	ztrdu(:,:,:) = ua(:,:,:)
	216	ztrdv(:,:,:) = va(:,:,:)
	217	ENDIF
	218
	219	IF( neuler == 0 .AND. kt == nit000 ) THEN
	220	z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping)
	221	ELSE
	222	z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog)
	223	ENDIF
	224
	225	DO jk = 2, jpkm1 ! Calculate and store vertical fluxes
	226	DO jj = 2, jpj
	227	DO ji = fs_2, jpi ! vector opt.
[5836]	228	zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
[4990]	229	END DO
	230	END DO
	231	END DO
[5836]	232
	233	DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero
[4990]	234	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	235	!!gm missing ISF boundary condition
[5120]	236	zwuw(ji,jj, 1 ) = 0._wp
	237	zwvw(ji,jj, 1 ) = 0._wp
[4990]	238	zwuw(ji,jj,jpk) = 0._wp
	239	zwvw(ji,jj,jpk) = 0._wp
	240	END DO
	241	END DO
	242
	243	! Start with before values and use sub timestepping to reach after values
	244
	245	zus(:,:,:,1) = ub(:,:,:)
	246	zvs(:,:,:,1) = vb(:,:,:)
	247
	248	DO jl = 1, jnzts ! Start of sub timestepping loop
	249
	250	IF( jl == 1 ) THEN ! Euler forward to kick things off
	251	jtb = 1 ; jtn = 1 ; jta = 2
	252	zts = z2dtzts
	253	ELSEIF( jl == 2 ) THEN ! First leapfrog step
	254	jtb = 1 ; jtn = 2 ; jta = 3
	255	zts = 2._wp * z2dtzts
	256	ELSE ! Shuffle pointers for subsequent leapfrog steps
	257	jtb = MOD(jtb,3) + 1
	258	jtn = MOD(jtn,3) + 1
	259	jta = MOD(jta,3) + 1
	260	ENDIF
	261
	262	DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical
	263	DO jj = 2, jpjm1 ! vertical momentum advection at w-point
	264	DO ji = fs_2, fs_jpim1 ! vector opt.
[5120]	265	zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk)
	266	zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk)
[4990]	267	END DO
	268	END DO
	269	END DO
	270	DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points
	271	DO jj = 2, jpjm1
	272	DO ji = fs_2, fs_jpim1 ! vector opt.
	273	! ! vertical momentum advective trends
[6140]	274	zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
	275	zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[4990]	276	zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
	277	zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
	278	END DO
	279	END DO
	280	END DO
	281
	282	END DO ! End of sub timestepping loop
	283
	284	zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
	285	DO jk = 1, jpkm1 ! Recover trends over the outer timestep
	286	DO jj = 2, jpjm1
	287	DO ji = fs_2, fs_jpim1 ! vector opt.
	288	! ! vertical momentum advective trends
	289	! ! add the trends to the general momentum trends
	290	ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
	291	va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
	292	END DO
	293	END DO
	294	END DO
	295
	296	IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic
	297	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	298	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
	299	CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
	300	CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
	301	ENDIF
	302	! ! Control print
	303	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, &
	304	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	305	!
[5836]	306	CALL wrk_dealloc( jpi,jpj,jpk, zwuw, zwvw, zww )
	307	CALL wrk_dealloc( jpi,jpj,jpk,3, zus , zvs )
[4990]	308	!
	309	IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts')
	310	!
	311	END SUBROUTINE dyn_zad_zts
	312
[503]	313	!!======================================================================
[3]	314	END MODULE dynzad

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

Download in other formats: