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.

sshwzv.F90 in NEMO/trunk/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/DYN/sshwzv.F90 @ 10364

Last change on this file since 10364 was 10364, checked in by acc, 5 years ago
Introduce Adaptive-Implicit vertical advection option to the trunk. This is code merged from branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP (see ticket #2042). The structure for the option is complete but is currently only successful with the flux-limited advection scheme (ln_traadv_mus). The use of this scheme with flux corrected advection schemes is not recommended until improvements to the nonoscillatory algorithm have been completed (work in progress elsewhere). The scheme is activated via a new namelist switch (ln_zad_Aimp) and is off by default.
Property svn:keywords set to `Id`
File size: 16.9 KB

Rev	Line
[1565]	1	MODULE sshwzv
[3]	2	!!==============================================================================
[1438]	3	!! * MODULE sshwzv *
	4	!! Ocean dynamics : sea surface height and vertical velocity
[3]	5	!!==============================================================================
[1438]	6	!! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code
[2528]	7	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA
	8	!! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
	9	!! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module
[4292]	10	!! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work
[3]	11	!!----------------------------------------------------------------------
[1438]	12
[3]	13	!!----------------------------------------------------------------------
[6140]	14	!! ssh_nxt : after ssh
	15	!! ssh_swp : filter ans swap the ssh arrays
	16	!! wzv : compute now vertical velocity
[1438]	17	!!----------------------------------------------------------------------
[6140]	18	USE oce ! ocean dynamics and tracers variables
	19	USE dom_oce ! ocean space and time domain variables
	20	USE sbc_oce ! surface boundary condition: ocean
	21	USE domvvl ! Variable volume
	22	USE divhor ! horizontal divergence
	23	USE phycst ! physical constants
[9019]	24	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
[6140]	25	USE bdydyn2d ! bdy_ssh routine
[2528]	26	#if defined key_agrif
[9570]	27	USE agrif_oce_interp
[2528]	28	#endif
[6140]	29	!
[10364]	30	USE iom
[6140]	31	USE in_out_manager ! I/O manager
	32	USE restart ! only for lrst_oce
	33	USE prtctl ! Print control
	34	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	35	USE lib_mpp ! MPP library
	36	USE timing ! Timing
[9023]	37	USE wet_dry ! Wetting/Drying flux limiting
[592]	38
[3]	39	IMPLICIT NONE
	40	PRIVATE
	41
[1438]	42	PUBLIC ssh_nxt ! called by step.F90
[4292]	43	PUBLIC wzv ! called by step.F90
[10364]	44	PUBLIC wAimp ! called by step.F90
[4292]	45	PUBLIC ssh_swp ! called by step.F90
[3]	46
	47	!! * Substitutions
[1438]	48	# include "vectopt_loop_substitute.h90"
[3]	49	!!----------------------------------------------------------------------
[9598]	50	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	51	!! $Id$
[10068]	52	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	53	!!----------------------------------------------------------------------
[3]	54	CONTAINS
	55
[4292]	56	SUBROUTINE ssh_nxt( kt )
[3]	57	!!----------------------------------------------------------------------
[4292]	58	!! * ROUTINE ssh_nxt *
[1438]	59	!!
[4292]	60	!! ** Purpose : compute the after ssh (ssha)
[3]	61	!!
[4292]	62	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
	63	!! is computed by integrating the horizontal divergence and multiply by
	64	!! by the time step.
[3]	65	!!
[5836]	66	!! ** action : ssha, after sea surface height
[2528]	67	!!
	68	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	69	!!----------------------------------------------------------------------
[5836]	70	INTEGER, INTENT(in) :: kt ! time step
[4292]	71	!
[7753]	72	INTEGER :: jk ! dummy loop indice
[5836]	73	REAL(wp) :: z2dt, zcoef ! local scalars
[9019]	74	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[3]	75	!!----------------------------------------------------------------------
[3294]	76	!
[9019]	77	IF( ln_timing ) CALL timing_start('ssh_nxt')
[3294]	78	!
[3]	79	IF( kt == nit000 ) THEN
	80	IF(lwp) WRITE(numout,*)
[4292]	81	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	82	IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]	83	ENDIF
[2528]	84	!
[7646]	85	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
[2715]	86	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
[7646]	87	zcoef = 0.5_wp * r1_rau0
[3]	88
[1438]	89	! !------------------------------!
	90	! ! After Sea Surface Height !
	91	! !------------------------------!
[9023]	92	IF(ln_wd_il) THEN
[7753]	93	CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
	94	ENDIF
[7646]	95
[7753]	96	CALL div_hor( kt ) ! Horizontal divergence
[7646]	97	!
[7753]	98	zhdiv(:,:) = 0._wp
[1438]	99	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[7753]	100	zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
[1438]	101	END DO
	102	! ! Sea surface elevation time stepping
[4338]	103	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	104	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
	105	!
[7753]	106	ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
[9023]	107	!
	108	#if defined key_agrif
	109	CALL agrif_ssh( kt )
	110	#endif
	111	!
[5930]	112	IF ( .NOT.ln_dynspg_ts ) THEN
[7646]	113	IF( ln_bdy ) THEN
[5930]	114	CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
	115	CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
	116	ENDIF
[4292]	117	ENDIF
	118	! !------------------------------!
	119	! ! outputs !
	120	! !------------------------------!
	121	!
[9440]	122	IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask )
[4292]	123	!
[9019]	124	IF( ln_timing ) CALL timing_stop('ssh_nxt')
[4292]	125	!
	126	END SUBROUTINE ssh_nxt
	127
	128
	129	SUBROUTINE wzv( kt )
	130	!!----------------------------------------------------------------------
	131	!! * ROUTINE wzv *
	132	!!
	133	!! ** Purpose : compute the now vertical velocity
	134	!!
	135	!! ** Method : - Using the incompressibility hypothesis, the vertical
	136	!! velocity is computed by integrating the horizontal divergence
	137	!! from the bottom to the surface minus the scale factor evolution.
	138	!! The boundary conditions are w=0 at the bottom (no flux) and.
	139	!!
	140	!! ** action : wn : now vertical velocity
	141	!!
	142	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	143	!!----------------------------------------------------------------------
[5836]	144	INTEGER, INTENT(in) :: kt ! time step
[4292]	145	!
[5836]	146	INTEGER :: ji, jj, jk ! dummy loop indices
	147	REAL(wp) :: z1_2dt ! local scalars
[9019]	148	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
[4292]	149	!!----------------------------------------------------------------------
	150	!
[9019]	151	IF( ln_timing ) CALL timing_start('wzv')
[5836]	152	!
[4292]	153	IF( kt == nit000 ) THEN
	154	IF(lwp) WRITE(numout,*)
	155	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
	156	IF(lwp) WRITE(numout,*) '~~~~~ '
	157	!
[7753]	158	wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4292]	159	ENDIF
	160	! !------------------------------!
	161	! ! Now Vertical Velocity !
	162	! !------------------------------!
	163	z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)
	164	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt
	165	!
	166	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases
[9019]	167	ALLOCATE( zhdiv(jpi,jpj,jpk) )
[4292]	168	!
	169	DO jk = 1, jpkm1
	170	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	171	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[4292]	172	DO jj = 2, jpjm1
	173	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	174	zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
[4292]	175	END DO
[592]	176	END DO
	177	END DO
[4292]	178	CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
	179	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
	180	! ! Same question holds for hdivn. Perhaps just for security
	181	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	182	! computation of w
[7753]	183	wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) &
	184	& + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)
[4292]	185	END DO
	186	! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
[9019]	187	DEALLOCATE( zhdiv )
[4292]	188	ELSE ! z_star and linear free surface cases
	189	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
[7753]	190	! computation of w
	191	wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) &
	192	& + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)
[4292]	193	END DO
[1438]	194	ENDIF
[592]	195
[7646]	196	IF( ln_bdy ) THEN
[4327]	197	DO jk = 1, jpkm1
[7753]	198	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
[4327]	199	END DO
	200	ENDIF
[4292]	201	!
[9023]	202	#if defined key_agrif
	203	IF( .NOT. AGRIF_Root() ) THEN
	204	IF ((nbondi == 1).OR.(nbondi == 2)) wn(nlci-1 , : ,:) = 0.e0 ! east
	205	IF ((nbondi == -1).OR.(nbondi == 2)) wn(2 , : ,:) = 0.e0 ! west
	206	IF ((nbondj == 1).OR.(nbondj == 2)) wn(: ,nlcj-1 ,:) = 0.e0 ! north
	207	IF ((nbondj == -1).OR.(nbondj == 2)) wn(: ,2 ,:) = 0.e0 ! south
	208	ENDIF
	209	#endif
[5836]	210	!
[9124]	211	IF( ln_timing ) CALL timing_stop('wzv')
[9023]	212	!
[5836]	213	END SUBROUTINE wzv
[592]	214
	215
[4292]	216	SUBROUTINE ssh_swp( kt )
[1438]	217	!!----------------------------------------------------------------------
	218	!! * ROUTINE ssh_nxt *
	219	!!
	220	!! ** Purpose : achieve the sea surface height time stepping by
	221	!! applying Asselin time filter and swapping the arrays
[4292]	222	!! ssha already computed in ssh_nxt
[1438]	223	!!
[2528]	224	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	225	!! from the filter, see Leclair and Madec 2010) and swap :
	226	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
	227	!! - atfp * rdt * ( emp_b - emp ) / rau0
	228	!! sshn = ssha
[1438]	229	!!
	230	!! ** action : - sshb, sshn : before & now sea surface height
	231	!! ready for the next time step
[2528]	232	!!
	233	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	234	!!----------------------------------------------------------------------
[2528]	235	INTEGER, INTENT(in) :: kt ! ocean time-step index
[6140]	236	!
	237	REAL(wp) :: zcoef ! local scalar
[1438]	238	!!----------------------------------------------------------------------
[3294]	239	!
[9124]	240	IF( ln_timing ) CALL timing_start('ssh_swp')
[3294]	241	!
[1438]	242	IF( kt == nit000 ) THEN
	243	IF(lwp) WRITE(numout,*)
[4292]	244	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
[1438]	245	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	246	ENDIF
[6140]	247	! !== Euler time-stepping: no filter, just swap ==!
[9239]	248	IF ( neuler == 0 .AND. kt == nit000 ) THEN
[7753]	249	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
[5836]	250	!
[6140]	251	ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==!
	252	! ! before <-- now filtered
[7753]	253	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
[6140]	254	IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed
	255	zcoef = atfp * rdt * r1_rau0
[7753]	256	sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) &
	257	& - rnf_b(:,:) + rnf (:,:) &
	258	& + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
[6140]	259	ENDIF
[7753]	260	sshn(:,:) = ssha(:,:) ! now <-- after
[1438]	261	ENDIF
	262	!
[9440]	263	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask )
[2528]	264	!
[9019]	265	IF( ln_timing ) CALL timing_stop('ssh_swp')
[3294]	266	!
[4292]	267	END SUBROUTINE ssh_swp
[3]	268
[10364]	269	SUBROUTINE wAimp( kt )
	270	!!----------------------------------------------------------------------
	271	!! * ROUTINE wAimp *
	272	!!
	273	!! ** Purpose : compute the Courant number and partition vertical velocity
	274	!! if a proportion needs to be treated implicitly
	275	!!
	276	!! ** Method : -
	277	!!
	278	!! ** action : wn : now vertical velocity (to be handled explicitly)
	279	!! : wi : now vertical velocity (for implicit treatment)
	280	!!
	281	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	282	!!----------------------------------------------------------------------
	283	INTEGER, INTENT(in) :: kt ! time step
	284	!
	285	INTEGER :: ji, jj, jk ! dummy loop indices
	286	REAL(wp) :: zCu, zcff, z1_e3w ! local scalars
	287	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
	288	REAL(wp) , PARAMETER :: Cu_max = 0.27 ! local parameters
	289	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
	290	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
	291	!!----------------------------------------------------------------------
	292	!
	293	IF( ln_timing ) CALL timing_start('wAimp')
	294	!
	295	IF( kt == nit000 ) THEN
	296	IF(lwp) WRITE(numout,*)
	297	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
	298	IF(lwp) WRITE(numout,*) '~~~~~ '
	299	!
	300	Cu_adv(:,:,jpk) = 0._wp ! bottom value : Cu_adv=0 (set once for all)
	301	ENDIF
	302	!
	303	DO jk = 1, jpkm1 ! calculate Courant numbers
	304	DO jj = 2, jpjm1
	305	DO ji = 2, fs_jpim1 ! vector opt.
	306	z1_e3w = 1._wp / e3w_n(ji,jj,jk)
	307	Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) &
	308	& + ( MAX( e2u(ji ,jj)e3uw_n(ji ,jj,jk)un(ji ,jj,jk), 0._wp ) - &
	309	& MIN( e2u(ji-1,jj)e3uw_n(ji-1,jj,jk)un(ji-1,jj,jk), 0._wp ) ) &
	310	& * r1_e1e2t(ji,jj) &
	311	& + ( MAX( e1v(ji,jj )e3vw_n(ji,jj ,jk)vn(ji,jj ,jk), 0._wp ) - &
	312	& MIN( e1v(ji,jj-1)e3vw_n(ji,jj-1,jk)vn(ji,jj-1,jk), 0._wp ) ) &
	313	& * r1_e1e2t(ji,jj) &
	314	& ) * z1_e3w
	315	END DO
	316	END DO
	317	END DO
	318	!
	319	CALL iom_put("Courant",Cu_adv)
	320	!
	321	wi(:,:,:) = 0._wp ! Includes top and bottom values set to zero
	322	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
	323	DO jk = 1, jpkm1 ! or scan Courant criterion and partition
	324	DO jj = 2, jpjm1 ! w where necessary
	325	DO ji = 2, fs_jpim1 ! vector opt.
	326	!
	327	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
	328	!
	329	IF( zCu < Cu_min ) THEN !<-- Fully explicit
	330	zcff = 0._wp
	331	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
	332	zcff = ( zCu - Cu_min )**2
	333	zcff = zcff / ( Fcu + zcff )
	334	ELSE !<-- Mostly implicit
	335	zcff = ( zCu - Cu_max )/ zCu
	336	ENDIF
	337	zcff = MIN(1._wp, zcff)
	338	!
	339	wi(ji,jj,jk) = zcff * wn(ji,jj,jk)
	340	wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
	341	!
	342	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient
	343	END DO
	344	END DO
	345	END DO
	346	ELSE
	347	! Fully explicit everywhere
	348	Cu_adv = 0.0_wp ! Reuse array to output coefficient
	349	ENDIF
	350	CALL iom_put("wimp",wi)
	351	CALL iom_put("wi_cff",Cu_adv)
	352	CALL iom_put("wexp",wn)
	353	!
	354	IF( ln_timing ) CALL timing_stop('wAimp')
	355	!
	356	END SUBROUTINE wAimp
[3]	357	!!======================================================================
[1565]	358	END MODULE sshwzv

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

Download in other formats: