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/ENHANCE-02_ISF_nemo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90 @ 11931

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

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

Download in other formats: