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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90 @ 10969

Last change on this file since 10969 was 10969, checked in by davestorkey, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : ASM and knock-on changes. NB. Only tested for compilation - no SETTE tests for ASM.
Property svn:keywords set to `Id`
File size: 17.0 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
[10969]	56	SUBROUTINE ssh_nxt( kt, Kbb, Kmm )
[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	!!----------------------------------------------------------------------
[10969]	70	INTEGER, INTENT(in) :: kt ! time step
	71	INTEGER, INTENT(in) :: Kbb, Kmm ! time level index
[4292]	72	!
[7753]	73	INTEGER :: jk ! dummy loop indice
[5836]	74	REAL(wp) :: z2dt, zcoef ! local scalars
[9019]	75	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[3]	76	!!----------------------------------------------------------------------
[3294]	77	!
[9019]	78	IF( ln_timing ) CALL timing_start('ssh_nxt')
[3294]	79	!
[3]	80	IF( kt == nit000 ) THEN
	81	IF(lwp) WRITE(numout,*)
[4292]	82	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	83	IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]	84	ENDIF
[2528]	85	!
[7646]	86	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
[2715]	87	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
[7646]	88	zcoef = 0.5_wp * r1_rau0
[3]	89
[1438]	90	! !------------------------------!
	91	! ! After Sea Surface Height !
	92	! !------------------------------!
[9023]	93	IF(ln_wd_il) THEN
[7753]	94	CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
	95	ENDIF
[7646]	96
[10969]	97	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
[7646]	98	!
[7753]	99	zhdiv(:,:) = 0._wp
[1438]	100	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[7753]	101	zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
[1438]	102	END DO
	103	! ! Sea surface elevation time stepping
[4338]	104	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	105	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
	106	!
[7753]	107	ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
[9023]	108	!
	109	#if defined key_agrif
	110	CALL agrif_ssh( kt )
	111	#endif
	112	!
[5930]	113	IF ( .NOT.ln_dynspg_ts ) THEN
[7646]	114	IF( ln_bdy ) THEN
[10425]	115	CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. ) ! Not sure that's necessary
[5930]	116	CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
	117	ENDIF
[4292]	118	ENDIF
	119	! !------------------------------!
	120	! ! outputs !
	121	! !------------------------------!
	122	!
[9440]	123	IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask )
[4292]	124	!
[9019]	125	IF( ln_timing ) CALL timing_stop('ssh_nxt')
[4292]	126	!
	127	END SUBROUTINE ssh_nxt
	128
	129
	130	SUBROUTINE wzv( kt )
	131	!!----------------------------------------------------------------------
	132	!! * ROUTINE wzv *
	133	!!
	134	!! ** Purpose : compute the now vertical velocity
	135	!!
	136	!! ** Method : - Using the incompressibility hypothesis, the vertical
	137	!! velocity is computed by integrating the horizontal divergence
	138	!! from the bottom to the surface minus the scale factor evolution.
	139	!! The boundary conditions are w=0 at the bottom (no flux) and.
	140	!!
	141	!! ** action : wn : now vertical velocity
	142	!!
	143	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	144	!!----------------------------------------------------------------------
[5836]	145	INTEGER, INTENT(in) :: kt ! time step
[4292]	146	!
[5836]	147	INTEGER :: ji, jj, jk ! dummy loop indices
	148	REAL(wp) :: z1_2dt ! local scalars
[9019]	149	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
[4292]	150	!!----------------------------------------------------------------------
	151	!
[9019]	152	IF( ln_timing ) CALL timing_start('wzv')
[5836]	153	!
[4292]	154	IF( kt == nit000 ) THEN
	155	IF(lwp) WRITE(numout,*)
	156	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
	157	IF(lwp) WRITE(numout,*) '~~~~~ '
	158	!
[7753]	159	wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4292]	160	ENDIF
	161	! !------------------------------!
	162	! ! Now Vertical Velocity !
	163	! !------------------------------!
	164	z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)
	165	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt
	166	!
	167	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases
[9019]	168	ALLOCATE( zhdiv(jpi,jpj,jpk) )
[4292]	169	!
	170	DO jk = 1, jpkm1
	171	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	172	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[4292]	173	DO jj = 2, jpjm1
	174	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	175	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]	176	END DO
[592]	177	END DO
	178	END DO
[10425]	179	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
[4292]	180	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
	181	! ! Same question holds for hdivn. Perhaps just for security
	182	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	183	! computation of w
[7753]	184	wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) &
	185	& + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)
[4292]	186	END DO
	187	! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
[9019]	188	DEALLOCATE( zhdiv )
[4292]	189	ELSE ! z_star and linear free surface cases
	190	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
[7753]	191	! computation of w
	192	wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) &
	193	& + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)
[4292]	194	END DO
[1438]	195	ENDIF
[592]	196
[7646]	197	IF( ln_bdy ) THEN
[4327]	198	DO jk = 1, jpkm1
[7753]	199	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
[4327]	200	END DO
	201	ENDIF
[4292]	202	!
[9023]	203	#if defined key_agrif
	204	IF( .NOT. AGRIF_Root() ) THEN
	205	IF ((nbondi == 1).OR.(nbondi == 2)) wn(nlci-1 , : ,:) = 0.e0 ! east
	206	IF ((nbondi == -1).OR.(nbondi == 2)) wn(2 , : ,:) = 0.e0 ! west
	207	IF ((nbondj == 1).OR.(nbondj == 2)) wn(: ,nlcj-1 ,:) = 0.e0 ! north
	208	IF ((nbondj == -1).OR.(nbondj == 2)) wn(: ,2 ,:) = 0.e0 ! south
	209	ENDIF
	210	#endif
[5836]	211	!
[9124]	212	IF( ln_timing ) CALL timing_stop('wzv')
[9023]	213	!
[5836]	214	END SUBROUTINE wzv
[592]	215
	216
[4292]	217	SUBROUTINE ssh_swp( kt )
[1438]	218	!!----------------------------------------------------------------------
	219	!! * ROUTINE ssh_nxt *
	220	!!
	221	!! ** Purpose : achieve the sea surface height time stepping by
	222	!! applying Asselin time filter and swapping the arrays
[4292]	223	!! ssha already computed in ssh_nxt
[1438]	224	!!
[2528]	225	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	226	!! from the filter, see Leclair and Madec 2010) and swap :
	227	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
	228	!! - atfp * rdt * ( emp_b - emp ) / rau0
	229	!! sshn = ssha
[1438]	230	!!
	231	!! ** action : - sshb, sshn : before & now sea surface height
	232	!! ready for the next time step
[2528]	233	!!
	234	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	235	!!----------------------------------------------------------------------
[2528]	236	INTEGER, INTENT(in) :: kt ! ocean time-step index
[6140]	237	!
	238	REAL(wp) :: zcoef ! local scalar
[1438]	239	!!----------------------------------------------------------------------
[3294]	240	!
[9124]	241	IF( ln_timing ) CALL timing_start('ssh_swp')
[3294]	242	!
[1438]	243	IF( kt == nit000 ) THEN
	244	IF(lwp) WRITE(numout,*)
[4292]	245	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
[1438]	246	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	247	ENDIF
[6140]	248	! !== Euler time-stepping: no filter, just swap ==!
[9239]	249	IF ( neuler == 0 .AND. kt == nit000 ) THEN
[7753]	250	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
[5836]	251	!
[6140]	252	ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==!
	253	! ! before <-- now filtered
[7753]	254	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
[6140]	255	IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed
	256	zcoef = atfp * rdt * r1_rau0
[7753]	257	sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) &
	258	& - rnf_b(:,:) + rnf (:,:) &
	259	& + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
[6140]	260	ENDIF
[7753]	261	sshn(:,:) = ssha(:,:) ! now <-- after
[1438]	262	ENDIF
	263	!
[9440]	264	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask )
[2528]	265	!
[9019]	266	IF( ln_timing ) CALL timing_stop('ssh_swp')
[3294]	267	!
[4292]	268	END SUBROUTINE ssh_swp
[3]	269
[10364]	270	SUBROUTINE wAimp( kt )
	271	!!----------------------------------------------------------------------
	272	!! * ROUTINE wAimp *
	273	!!
	274	!! ** Purpose : compute the Courant number and partition vertical velocity
	275	!! if a proportion needs to be treated implicitly
	276	!!
	277	!! ** Method : -
	278	!!
	279	!! ** action : wn : now vertical velocity (to be handled explicitly)
	280	!! : wi : now vertical velocity (for implicit treatment)
	281	!!
	282	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	283	!!----------------------------------------------------------------------
	284	INTEGER, INTENT(in) :: kt ! time step
	285	!
	286	INTEGER :: ji, jj, jk ! dummy loop indices
	287	REAL(wp) :: zCu, zcff, z1_e3w ! local scalars
	288	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
	289	REAL(wp) , PARAMETER :: Cu_max = 0.27 ! local parameters
	290	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
	291	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
	292	!!----------------------------------------------------------------------
	293	!
	294	IF( ln_timing ) CALL timing_start('wAimp')
	295	!
	296	IF( kt == nit000 ) THEN
	297	IF(lwp) WRITE(numout,*)
	298	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
	299	IF(lwp) WRITE(numout,*) '~~~~~ '
	300	!
	301	Cu_adv(:,:,jpk) = 0._wp ! bottom value : Cu_adv=0 (set once for all)
	302	ENDIF
	303	!
	304	DO jk = 1, jpkm1 ! calculate Courant numbers
	305	DO jj = 2, jpjm1
	306	DO ji = 2, fs_jpim1 ! vector opt.
	307	z1_e3w = 1._wp / e3w_n(ji,jj,jk)
	308	Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) &
	309	& + ( MAX( e2u(ji ,jj)e3uw_n(ji ,jj,jk)un(ji ,jj,jk), 0._wp ) - &
	310	& MIN( e2u(ji-1,jj)e3uw_n(ji-1,jj,jk)un(ji-1,jj,jk), 0._wp ) ) &
	311	& * r1_e1e2t(ji,jj) &
	312	& + ( MAX( e1v(ji,jj )e3vw_n(ji,jj ,jk)vn(ji,jj ,jk), 0._wp ) - &
	313	& MIN( e1v(ji,jj-1)e3vw_n(ji,jj-1,jk)vn(ji,jj-1,jk), 0._wp ) ) &
	314	& * r1_e1e2t(ji,jj) &
	315	& ) * z1_e3w
	316	END DO
	317	END DO
	318	END DO
	319	!
	320	CALL iom_put("Courant",Cu_adv)
	321	!
	322	wi(:,:,:) = 0._wp ! Includes top and bottom values set to zero
	323	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
	324	DO jk = 1, jpkm1 ! or scan Courant criterion and partition
	325	DO jj = 2, jpjm1 ! w where necessary
	326	DO ji = 2, fs_jpim1 ! vector opt.
	327	!
	328	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
	329	!
	330	IF( zCu < Cu_min ) THEN !<-- Fully explicit
	331	zcff = 0._wp
	332	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
	333	zcff = ( zCu - Cu_min )**2
	334	zcff = zcff / ( Fcu + zcff )
	335	ELSE !<-- Mostly implicit
	336	zcff = ( zCu - Cu_max )/ zCu
	337	ENDIF
	338	zcff = MIN(1._wp, zcff)
	339	!
	340	wi(ji,jj,jk) = zcff * wn(ji,jj,jk)
	341	wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
	342	!
	343	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient
	344	END DO
	345	END DO
	346	END DO
	347	ELSE
	348	! Fully explicit everywhere
	349	Cu_adv = 0.0_wp ! Reuse array to output coefficient
	350	ENDIF
	351	CALL iom_put("wimp",wi)
	352	CALL iom_put("wi_cff",Cu_adv)
	353	CALL iom_put("wexp",wn)
	354	!
	355	IF( ln_timing ) CALL timing_stop('wAimp')
	356	!
	357	END SUBROUTINE wAimp
[3]	358	!!======================================================================
[1565]	359	END MODULE sshwzv

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

Download in other formats: