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/2021/dev_r14318_RK3_stage1/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90 @ 14618

Last change on this file since 14618 was 14618, checked in by techene, 3 years ago
#2506 RK3 work in progess : for 2 configurag
Property svn:keywords set to `Id`
File size: 28.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
[14053]	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
	10	!! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work
	11	!! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection
	12	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering.
	13	!! - ! 2020-08 (S. Techene, G. Madec) add here ssh initiatlisation
[3]	14	!!----------------------------------------------------------------------
[1438]	15
[3]	16	!!----------------------------------------------------------------------
[6140]	17	!! ssh_nxt : after ssh
[12377]	18	!! ssh_atf : time filter the ssh arrays
[14547]	19	!! wzv : generic interface of vertical velocity calculation
	20	!! wzv_MLF : MLF: compute NOW vertical velocity
	21	!! wzv_RK3 : RK3: Compute a vertical velocity
[1438]	22	!!----------------------------------------------------------------------
[6140]	23	USE oce ! ocean dynamics and tracers variables
[12377]	24	USE isf_oce ! ice shelf
[6140]	25	USE dom_oce ! ocean space and time domain variables
	26	USE sbc_oce ! surface boundary condition: ocean
	27	USE domvvl ! Variable volume
	28	USE divhor ! horizontal divergence
	29	USE phycst ! physical constants
[9019]	30	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
[6140]	31	USE bdydyn2d ! bdy_ssh routine
[14139]	32	USE wet_dry ! Wetting/Drying flux limiting
[2528]	33	#if defined key_agrif
[13286]	34	USE agrif_oce
[9570]	35	USE agrif_oce_interp
[2528]	36	#endif
[6140]	37	!
[10364]	38	USE iom
[6140]	39	USE in_out_manager ! I/O manager
	40	USE restart ! only for lrst_oce
	41	USE prtctl ! Print control
	42	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	43	USE lib_mpp ! MPP library
	44	USE timing ! Timing
[14053]	45
[3]	46	IMPLICIT NONE
	47	PRIVATE
[14547]	48	! !! * Interface
	49	INTERFACE wzv
	50	MODULE PROCEDURE wzv_MLF, wzv_RK3
	51	END INTERFACE
[3]	52
[14053]	53	PUBLIC ssh_nxt ! called by step.F90
	54	PUBLIC wzv ! called by step.F90
	55	PUBLIC wAimp ! called by step.F90
	56	PUBLIC ssh_atf ! called by step.F90
[3]	57
	58	!! * Substitutions
[12377]	59	# include "do_loop_substitute.h90"
[13237]	60	# include "domzgr_substitute.h90"
[3]	61	!!----------------------------------------------------------------------
[9598]	62	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	63	!! $Id$
[10068]	64	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	65	!!----------------------------------------------------------------------
[3]	66	CONTAINS
	67
[12377]	68	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
[3]	69	!!----------------------------------------------------------------------
[4292]	70	!! * ROUTINE ssh_nxt *
[1438]	71	!!
[12377]	72	!! ** Purpose : compute the after ssh (ssh(Kaa))
[3]	73	!!
[4292]	74	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
	75	!! is computed by integrating the horizontal divergence and multiply by
	76	!! by the time step.
[3]	77	!!
[12377]	78	!! ** action : ssh(:,:,Kaa), after sea surface height
[2528]	79	!!
	80	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	81	!!----------------------------------------------------------------------
[12377]	82	INTEGER , INTENT(in ) :: kt ! time step
	83	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
	84	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
[14547]	85	!
[12489]	86	INTEGER :: jk ! dummy loop index
	87	REAL(wp) :: zcoef ! local scalar
[9019]	88	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[14547]	89	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[3]	90	!!----------------------------------------------------------------------
[3294]	91	!
[9019]	92	IF( ln_timing ) CALL timing_start('ssh_nxt')
[3294]	93	!
[3]	94	IF( kt == nit000 ) THEN
	95	IF(lwp) WRITE(numout,*)
[4292]	96	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	97	IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]	98	ENDIF
[2528]	99	!
[12489]	100	zcoef = 0.5_wp * r1_rho0
[1438]	101	! !------------------------------!
	102	! ! After Sea Surface Height !
	103	! !------------------------------!
[9023]	104	IF(ln_wd_il) THEN
[14547]	105	CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
[7753]	106	ENDIF
[7646]	107
[12377]	108	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
[7646]	109	!
[7753]	110	zhdiv(:,:) = 0._wp
[1438]	111	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[12377]	112	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
[1438]	113	END DO
	114	! ! Sea surface elevation time stepping
[4338]	115	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	116	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
	117	!
[12489]	118	pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
[9023]	119	!
	120	#if defined key_agrif
[13237]	121	Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
	122	CALL agrif_ssh( kt )
[9023]	123	#endif
	124	!
[5930]	125	IF ( .NOT.ln_dynspg_ts ) THEN
[7646]	126	IF( ln_bdy ) THEN
[13226]	127	CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
[12377]	128	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
[5930]	129	ENDIF
[4292]	130	ENDIF
	131	! !------------------------------!
	132	! ! outputs !
	133	! !------------------------------!
	134	!
[12377]	135	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
[4292]	136	!
[9019]	137	IF( ln_timing ) CALL timing_stop('ssh_nxt')
[4292]	138	!
	139	END SUBROUTINE ssh_nxt
	140
[14547]	141
	142	SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww )
[4292]	143	!!----------------------------------------------------------------------
[14547]	144	!! * ROUTINE wzv_MLF *
[4292]	145	!!
	146	!! ** Purpose : compute the now vertical velocity
	147	!!
	148	!! ** Method : - Using the incompressibility hypothesis, the vertical
	149	!! velocity is computed by integrating the horizontal divergence
	150	!! from the bottom to the surface minus the scale factor evolution.
	151	!! The boundary conditions are w=0 at the bottom (no flux) and.
	152	!!
[12377]	153	!! ** action : pww : now vertical velocity
[4292]	154	!!
	155	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	156	!!----------------------------------------------------------------------
[12377]	157	INTEGER , INTENT(in) :: kt ! time step
	158	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
[13237]	159	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
[4292]	160	!
[5836]	161	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	162	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
[4292]	163	!!----------------------------------------------------------------------
	164	!
[14547]	165	IF( ln_timing ) CALL timing_start('wzv_MLF')
[5836]	166	!
[4292]	167	IF( kt == nit000 ) THEN
	168	IF(lwp) WRITE(numout,*)
[14547]	169	IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
	170	IF(lwp) WRITE(numout,*) '~~~~~~~'
[4292]	171	!
[12377]	172	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4292]	173	ENDIF
	174	! !------------------------------!
	175	! ! Now Vertical Velocity !
	176	! !------------------------------!
	177	!
[13237]	178	! !===============================!
	179	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
	180	! !===============================!
[9019]	181	ALLOCATE( zhdiv(jpi,jpj,jpk) )
[4292]	182	!
	183	DO jk = 1, jpkm1
	184	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	185	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[13295]	186	DO_2D( 0, 0, 0, 0 )
[12377]	187	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) )
	188	END_2D
[592]	189	END DO
[13226]	190	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
[4292]	191	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
[12377]	192	! ! Same question holds for hdiv. Perhaps just for security
[4292]	193	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	194	! computation of w
[13237]	195	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	196	& + zhdiv(:,:,jk) &
	197	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
	198	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
[4292]	199	END DO
[12377]	200	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
[14547]	201	DEALLOCATE( zhdiv )
[13237]	202	! !=================================!
	203	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
	204	! !=================================!
	205	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	206	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk)
	207	END DO
	208	! !==========================================!
	209	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
	210	! !==========================================!
[13497]	211	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
[14547]	212	! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
	213	!!gm old
	214	! pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	215	! & + r1_Dt * ( e3t(:,:,jk,Kaa) &
	216	! & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
	217	!!gm slightly faster :
[14618]	218	#if defined key_qco
[14547]	219	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	220	& + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) ) ) * tmask(:,:,jk)
[14618]	221	#else
	222	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	223	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
	224	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
	225	#endif
	226	!!gm end !!st need to be there when not using qco
[14547]	227
	228
[4292]	229	END DO
[1438]	230	ENDIF
[592]	231
[7646]	232	IF( ln_bdy ) THEN
[4327]	233	DO jk = 1, jpkm1
[12377]	234	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
[4327]	235	END DO
	236	ENDIF
[4292]	237	!
[13286]	238	#if defined key_agrif
	239	IF( .NOT. AGRIF_Root() ) THEN
	240	!
[12965]	241	! Mask vertical velocity at first/last columns/row
	242	! inside computational domain (cosmetic)
[13286]	243	DO jk = 1, jpkm1
	244	IF( lk_west ) THEN ! --- West --- !
	245	DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
	246	DO jj = 1, jpj
	247	pww(ji,jj,jk) = 0._wp
	248	END DO
	249	END DO
	250	ENDIF
	251	IF( lk_east ) THEN ! --- East --- !
	252	DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
	253	DO jj = 1, jpj
	254	pww(ji,jj,jk) = 0._wp
	255	END DO
	256	END DO
	257	ENDIF
	258	IF( lk_south ) THEN ! --- South --- !
	259	DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
	260	DO ji = 1, jpi
	261	pww(ji,jj,jk) = 0._wp
	262	END DO
	263	END DO
	264	ENDIF
	265	IF( lk_north ) THEN ! --- North --- !
	266	DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
	267	DO ji = 1, jpi
	268	pww(ji,jj,jk) = 0._wp
	269	END DO
	270	END DO
	271	ENDIF
	272	!
	273	END DO
[12965]	274	!
[9023]	275	ENDIF
[13286]	276	#endif
[5836]	277	!
[14547]	278	IF( ln_timing ) CALL timing_stop('wzv_MLF')
[9023]	279	!
[14547]	280	END SUBROUTINE wzv_MLF
[592]	281
	282
[14547]	283	SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
	284	!!----------------------------------------------------------------------
	285	!! * ROUTINE wzv_RK3 *
	286	!!
	287	!! ** Purpose : compute the now vertical velocity
	288	!!
	289	!! ** Method : - Using the incompressibility hypothesis, the vertical
	290	!! velocity is computed by integrating the horizontal divergence
	291	!! from the bottom to the surface minus the scale factor evolution.
	292	!! The boundary conditions are w=0 at the bottom (no flux) and.
	293	!!
	294	!! ** action : pww : now vertical velocity
	295	!!
	296	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	297	!!----------------------------------------------------------------------
	298	INTEGER , INTENT(in ) :: kt ! time step
	299	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices
	300	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: puu, pvv ! horizontal velocity at Kmm
	301	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
	302	!
	303	INTEGER :: ji, jj, jk ! dummy loop indices
	304	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
	305	REAL(wp) , DIMENSION(jpi,jpj,jpk) :: ze3div
	306	!!----------------------------------------------------------------------
	307	!
	308	IF( ln_timing ) CALL timing_start('wzv_RK3')
	309	!
	310	IF( kt == nit000 ) THEN
	311	IF(lwp) WRITE(numout,*)
	312	IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity '
	313	IF(lwp) WRITE(numout,*) '~~~~~ '
	314	!
	315	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
	316	ENDIF
	317	!
	318	CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div )
	319	! !------------------------------!
	320	! ! Now Vertical Velocity !
	321	! !------------------------------!
	322	!
	323	! !===============================!
	324	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
	325	! !===============================!
	326	ALLOCATE( zhdiv(jpi,jpj,jpk) )
	327	!
	328	DO jk = 1, jpkm1
	329	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
	330	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
	331	DO_2D( 0, 0, 0, 0 )
	332	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) )
	333	END_2D
	334	END DO
	335	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
	336	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
	337	! ! Same question holds for hdiv. Perhaps just for security
	338	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	339	! computation of w
	340	pww(:,:,jk) = pww(:,:,jk+1) - ( ze3div(:,:,jk) + zhdiv(:,:,jk) &
	341	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
	342	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
	343	END DO
	344	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
	345	DEALLOCATE( zhdiv )
	346	! !=================================!
	347	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
	348	! !=================================!
	349	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	350	pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk)
	351	END DO
	352	! !==========================================!
	353	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
	354	! !==========================================!
	355	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	356	! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
	357	pww(:,:,jk) = pww(:,:,jk+1) - ( ze3div(:,:,jk) &
	358	& + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) ) ) * tmask(:,:,jk)
	359	END DO
	360	ENDIF
	361
	362	IF( ln_bdy ) THEN
	363	DO jk = 1, jpkm1
	364	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
	365	END DO
	366	ENDIF
	367	!
	368	#if defined key_agrif
	369	IF( .NOT. AGRIF_Root() ) THEN
	370	!
	371	! Mask vertical velocity at first/last columns/row
	372	! inside computational domain (cosmetic)
	373	DO jk = 1, jpkm1
	374	IF( lk_west ) THEN ! --- West --- !
	375	DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
	376	DO jj = 1, jpj
	377	pww(ji,jj,jk) = 0._wp
	378	END DO
	379	END DO
	380	ENDIF
	381	IF( lk_east ) THEN ! --- East --- !
	382	DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
	383	DO jj = 1, jpj
	384	pww(ji,jj,jk) = 0._wp
	385	END DO
	386	END DO
	387	ENDIF
	388	IF( lk_south ) THEN ! --- South --- !
	389	DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
	390	DO ji = 1, jpi
	391	pww(ji,jj,jk) = 0._wp
	392	END DO
	393	END DO
	394	ENDIF
	395	IF( lk_north ) THEN ! --- North --- !
	396	DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
	397	DO ji = 1, jpi
	398	pww(ji,jj,jk) = 0._wp
	399	END DO
	400	END DO
	401	ENDIF
	402	!
	403	END DO
	404	!
	405	ENDIF
	406	#endif
	407	!
	408	IF( ln_timing ) CALL timing_stop('wzv_RK3')
	409	!
	410	END SUBROUTINE wzv_RK3
	411
	412
[14205]	413	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
[1438]	414	!!----------------------------------------------------------------------
[12377]	415	!! * ROUTINE ssh_atf *
[1438]	416	!!
[12377]	417	!! ** Purpose : Apply Asselin time filter to now SSH.
[1438]	418	!!
[2528]	419	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	420	!! from the filter, see Leclair and Madec 2010) and swap :
[12489]	421	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
	422	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
[1438]	423	!!
[12377]	424	!! ** action : - pssh(:,:,Kmm) time filtered
[2528]	425	!!
	426	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	427	!!----------------------------------------------------------------------
[12377]	428	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	429	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
[14205]	430	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field
[6140]	431	!
	432	REAL(wp) :: zcoef ! local scalar
[1438]	433	!!----------------------------------------------------------------------
[3294]	434	!
[12377]	435	IF( ln_timing ) CALL timing_start('ssh_atf')
[3294]	436	!
[1438]	437	IF( kt == nit000 ) THEN
	438	IF(lwp) WRITE(numout,*)
[12377]	439	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
[1438]	440	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	441	ENDIF
[14205]	442	!
	443	IF( .NOT.l_1st_euler ) THEN ! Apply Asselin time filter on Kmm field (not on euler 1st)
[14053]	444	!
[14205]	445	IF( ln_linssh ) THEN ! filtered "now" field
	446	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
	447	!
	448	ELSE ! filtered "now" field with forcing removed
[12489]	449	zcoef = rn_atfp * rn_Dt * r1_rho0
[14205]	450	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) &
	451	& - zcoef * ( emp_b(:,:) - emp(:,:) &
	452	& - rnf_b(:,:) + rnf(:,:) &
	453	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
	454	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
[12377]	455
	456	! ice sheet coupling
[14053]	457	IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 ) &
[14205]	458	& pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0._wp ) * ssmask(:,:)
[12377]	459
[6140]	460	ENDIF
[1438]	461	ENDIF
	462	!
[14053]	463	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf - pssh(:,:,Kmm): ', mask1=tmask )
[2528]	464	!
[12377]	465	IF( ln_timing ) CALL timing_stop('ssh_atf')
[3294]	466	!
[12377]	467	END SUBROUTINE ssh_atf
[3]	468
[13237]	469
[12377]	470	SUBROUTINE wAimp( kt, Kmm )
[10364]	471	!!----------------------------------------------------------------------
	472	!! * ROUTINE wAimp *
	473	!!
	474	!! ** Purpose : compute the Courant number and partition vertical velocity
	475	!! if a proportion needs to be treated implicitly
	476	!!
	477	!! ** Method : -
	478	!!
[12377]	479	!! ** action : ww : now vertical velocity (to be handled explicitly)
[10364]	480	!! : wi : now vertical velocity (for implicit treatment)
	481	!!
[11414]	482	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
	483	!! implicit scheme for vertical advection in oceanic modeling.
	484	!! Ocean Modelling, 91, 38-69.
[10364]	485	!!----------------------------------------------------------------------
	486	INTEGER, INTENT(in) :: kt ! time step
[12377]	487	INTEGER, INTENT(in) :: Kmm ! time level index
[10364]	488	!
	489	INTEGER :: ji, jj, jk ! dummy loop indices
[13237]	490	REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars
[10364]	491	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
[11407]	492	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
[10364]	493	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
	494	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
	495	!!----------------------------------------------------------------------
	496	!
	497	IF( ln_timing ) CALL timing_start('wAimp')
	498	!
	499	IF( kt == nit000 ) THEN
	500	IF(lwp) WRITE(numout,*)
	501	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
	502	IF(lwp) WRITE(numout,*) '~~~~~ '
[11293]	503	wi(:,:,:) = 0._wp
[10364]	504	ENDIF
	505	!
[11414]	506	! Calculate Courant numbers
[13237]	507	zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability)
[11414]	508	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[13295]	509	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	510	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[13237]	511	Cu_adv(ji,jj,jk) = zdt * &
	512	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
	513	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
	514	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
	515	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
	516	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
[12377]	517	& * r1_e1e2t(ji,jj) &
[13237]	518	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
	519	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
	520	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
	521	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
[12377]	522	& * r1_e1e2t(ji,jj) &
	523	& ) * z1_e3t
	524	END_3D
[11414]	525	ELSE
[13295]	526	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	527	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[13237]	528	Cu_adv(ji,jj,jk) = zdt * &
	529	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
[12377]	530	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
	531	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
	532	& * r1_e1e2t(ji,jj) &
	533	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
	534	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
	535	& * r1_e1e2t(ji,jj) &
	536	& ) * z1_e3t
	537	END_3D
[11414]	538	ENDIF
[13226]	539	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
[10364]	540	!
	541	CALL iom_put("Courant",Cu_adv)
	542	!
	543	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
[13497]	544	DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary
[12377]	545	!
	546	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
[11293]	547	! alt:
[12377]	548	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
[11293]	549	! zCu = Cu_adv(ji,jj,jk)
	550	! ELSE
	551	! zCu = Cu_adv(ji,jj,jk-1)
	552	! ENDIF
[12377]	553	!
	554	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
	555	zcff = 0._wp
	556	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
	557	zcff = ( zCu - Cu_min )**2
	558	zcff = zcff / ( Fcu + zcff )
	559	ELSE !<-- Mostly implicit
	560	zcff = ( zCu - Cu_max )/ zCu
	561	ENDIF
	562	zcff = MIN(1._wp, zcff)
	563	!
	564	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
	565	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
	566	!
	567	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
	568	END_3D
[11293]	569	Cu_adv(:,:,1) = 0._wp
[10364]	570	ELSE
	571	! Fully explicit everywhere
[11407]	572	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
[10907]	573	wi (:,:,:) = 0._wp
[10364]	574	ENDIF
	575	CALL iom_put("wimp",wi)
	576	CALL iom_put("wi_cff",Cu_adv)
[12377]	577	CALL iom_put("wexp",ww)
[10364]	578	!
	579	IF( ln_timing ) CALL timing_stop('wAimp')
	580	!
	581	END SUBROUTINE wAimp
[14053]	582
[3]	583	!!======================================================================
[1565]	584	END MODULE sshwzv

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

Download in other formats: