Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

sshwzv.F90

Last change on this file was 15150, checked in by smasson, 3 years ago
trunk: remove some of the changes done in [15145]
Property svn:keywords set to `Id`
File size: 22.2 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
[6140]	19	!! wzv : compute now vertical velocity
[1438]	20	!!----------------------------------------------------------------------
[6140]	21	USE oce ! ocean dynamics and tracers variables
[12377]	22	USE isf_oce ! ice shelf
[6140]	23	USE dom_oce ! ocean space and time domain variables
	24	USE sbc_oce ! surface boundary condition: ocean
	25	USE domvvl ! Variable volume
	26	USE divhor ! horizontal divergence
	27	USE phycst ! physical constants
[9019]	28	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
[6140]	29	USE bdydyn2d ! bdy_ssh routine
[14139]	30	USE wet_dry ! Wetting/Drying flux limiting
[2528]	31	#if defined key_agrif
[13286]	32	USE agrif_oce
[9570]	33	USE agrif_oce_interp
[2528]	34	#endif
[6140]	35	!
[10364]	36	USE iom
[6140]	37	USE in_out_manager ! I/O manager
	38	USE restart ! only for lrst_oce
	39	USE prtctl ! Print control
	40	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	41	USE lib_mpp ! MPP library
	42	USE timing ! Timing
[14053]	43
[3]	44	IMPLICIT NONE
	45	PRIVATE
	46
[14053]	47	PUBLIC ssh_nxt ! called by step.F90
	48	PUBLIC wzv ! called by step.F90
	49	PUBLIC wAimp ! called by step.F90
	50	PUBLIC ssh_atf ! called by step.F90
[3]	51
	52	!! * Substitutions
[12377]	53	# include "do_loop_substitute.h90"
[13237]	54	# include "domzgr_substitute.h90"
[3]	55	!!----------------------------------------------------------------------
[9598]	56	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	57	!! $Id$
[10068]	58	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	59	!!----------------------------------------------------------------------
[3]	60	CONTAINS
	61
[12377]	62	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
[3]	63	!!----------------------------------------------------------------------
[4292]	64	!! * ROUTINE ssh_nxt *
[1438]	65	!!
[12377]	66	!! ** Purpose : compute the after ssh (ssh(Kaa))
[3]	67	!!
[4292]	68	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
	69	!! is computed by integrating the horizontal divergence and multiply by
	70	!! by the time step.
[3]	71	!!
[12377]	72	!! ** action : ssh(:,:,Kaa), after sea surface height
[2528]	73	!!
	74	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	75	!!----------------------------------------------------------------------
[12377]	76	INTEGER , INTENT(in ) :: kt ! time step
	77	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
	78	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
[4292]	79	!
[14834]	80	INTEGER :: ji, jj, jk ! dummy loop index
[12489]	81	REAL(wp) :: zcoef ! local scalar
[9019]	82	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[3]	83	!!----------------------------------------------------------------------
[3294]	84	!
[9019]	85	IF( ln_timing ) CALL timing_start('ssh_nxt')
[3294]	86	!
[3]	87	IF( kt == nit000 ) THEN
	88	IF(lwp) WRITE(numout,*)
[4292]	89	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	90	IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]	91	ENDIF
[2528]	92	!
[12489]	93	zcoef = 0.5_wp * r1_rho0
[3]	94
[1438]	95	! !------------------------------!
	96	! ! After Sea Surface Height !
	97	! !------------------------------!
[9023]	98	IF(ln_wd_il) THEN
[12489]	99	CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
[7753]	100	ENDIF
[7646]	101
[12377]	102	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
[7646]	103	!
[7753]	104	zhdiv(:,:) = 0._wp
[15150]	105	DO_3D( 1, nn_hls, 1, nn_hls, 1, jpkm1 ) ! Horizontal divergence of barotropic transports
[14834]	106	zhdiv(ji,jj) = zhdiv(ji,jj) + e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)
	107	END_3D
[1438]	108	! ! Sea surface elevation time stepping
[4338]	109	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	110	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
	111	!
[15150]	112	DO_2D_OVR( 1, nn_hls, 1, nn_hls ) ! Loop bounds limited by hdiv definition in div_hor
[14834]	113	pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
	114	END_2D
	115	! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp)
	116	IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )
[9023]	117	!
	118	#if defined key_agrif
[13237]	119	Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
	120	CALL agrif_ssh( kt )
[9023]	121	#endif
	122	!
[5930]	123	IF ( .NOT.ln_dynspg_ts ) THEN
[7646]	124	IF( ln_bdy ) THEN
[14834]	125	IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
[12377]	126	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
[5930]	127	ENDIF
[4292]	128	ENDIF
	129	! !------------------------------!
	130	! ! outputs !
	131	! !------------------------------!
	132	!
[12377]	133	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
[4292]	134	!
[9019]	135	IF( ln_timing ) CALL timing_stop('ssh_nxt')
[4292]	136	!
	137	END SUBROUTINE ssh_nxt
	138
	139
[13237]	140	SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
[4292]	141	!!----------------------------------------------------------------------
	142	!! * ROUTINE wzv *
	143	!!
	144	!! ** Purpose : compute the now vertical velocity
	145	!!
	146	!! ** Method : - Using the incompressibility hypothesis, the vertical
	147	!! velocity is computed by integrating the horizontal divergence
	148	!! from the bottom to the surface minus the scale factor evolution.
	149	!! The boundary conditions are w=0 at the bottom (no flux) and.
	150	!!
[12377]	151	!! ** action : pww : now vertical velocity
[4292]	152	!!
	153	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	154	!!----------------------------------------------------------------------
[12377]	155	INTEGER , INTENT(in) :: kt ! time step
	156	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
[13237]	157	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
[4292]	158	!
[5836]	159	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	160	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
[4292]	161	!!----------------------------------------------------------------------
	162	!
[9019]	163	IF( ln_timing ) CALL timing_start('wzv')
[5836]	164	!
[4292]	165	IF( kt == nit000 ) THEN
	166	IF(lwp) WRITE(numout,*)
	167	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
	168	IF(lwp) WRITE(numout,*) '~~~~~ '
	169	!
[12377]	170	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4292]	171	ENDIF
	172	! !------------------------------!
	173	! ! Now Vertical Velocity !
	174	! !------------------------------!
	175	!
[13237]	176	! !===============================!
	177	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
	178	! !===============================!
[9019]	179	ALLOCATE( zhdiv(jpi,jpj,jpk) )
[4292]	180	!
	181	DO jk = 1, jpkm1
	182	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	183	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[15058]	184	DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
[12377]	185	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) )
	186	END_2D
[592]	187	END DO
[15058]	188	IF( nn_hls == 1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
[4292]	189	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
[12377]	190	! ! Same question holds for hdiv. Perhaps just for security
[15055]	191	! ! clem: yes it is a problem because ww is used in many other places where we need the halos
	192	!
[15058]	193	DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
[4292]	194	! computation of w
[14834]	195	pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
	196	& + zhdiv(ji,jj,jk) &
	197	& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
	198	& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
	199	END_3D
[12377]	200	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
[9019]	201	DEALLOCATE( zhdiv )
[13237]	202	! !=================================!
	203	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
	204	! !=================================!
[15058]	205	DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
[14834]	206	pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) ) * tmask(ji,jj,jk)
	207	END_3D
[13237]	208	! !==========================================!
	209	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
	210	! !==========================================!
[15058]	211	DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
[14834]	212	pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
	213	& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
	214	& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
	215	END_3D
[1438]	216	ENDIF
[592]	217
[7646]	218	IF( ln_bdy ) THEN
[4327]	219	DO jk = 1, jpkm1
[15058]	220	DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
[15055]	221	pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
	222	END_2D
[4327]	223	END DO
	224	ENDIF
[4292]	225	!
[13286]	226	#if defined key_agrif
	227	IF( .NOT. AGRIF_Root() ) THEN
	228	!
[12965]	229	! Mask vertical velocity at first/last columns/row
	230	! inside computational domain (cosmetic)
[13286]	231	DO jk = 1, jpkm1
	232	IF( lk_west ) THEN ! --- West --- !
	233	DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
	234	DO jj = 1, jpj
	235	pww(ji,jj,jk) = 0._wp
	236	END DO
	237	END DO
	238	ENDIF
	239	IF( lk_east ) THEN ! --- East --- !
	240	DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
	241	DO jj = 1, jpj
	242	pww(ji,jj,jk) = 0._wp
	243	END DO
	244	END DO
	245	ENDIF
	246	IF( lk_south ) THEN ! --- South --- !
	247	DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
	248	DO ji = 1, jpi
	249	pww(ji,jj,jk) = 0._wp
	250	END DO
	251	END DO
	252	ENDIF
	253	IF( lk_north ) THEN ! --- North --- !
	254	DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
	255	DO ji = 1, jpi
	256	pww(ji,jj,jk) = 0._wp
	257	END DO
	258	END DO
	259	ENDIF
	260	!
	261	END DO
[12965]	262	!
[9023]	263	ENDIF
[13286]	264	#endif
[5836]	265	!
[9124]	266	IF( ln_timing ) CALL timing_stop('wzv')
[9023]	267	!
[5836]	268	END SUBROUTINE wzv
[592]	269
	270
[14205]	271	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
[1438]	272	!!----------------------------------------------------------------------
[12377]	273	!! * ROUTINE ssh_atf *
[1438]	274	!!
[12377]	275	!! ** Purpose : Apply Asselin time filter to now SSH.
[1438]	276	!!
[2528]	277	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	278	!! from the filter, see Leclair and Madec 2010) and swap :
[12489]	279	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
	280	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
[1438]	281	!!
[12377]	282	!! ** action : - pssh(:,:,Kmm) time filtered
[2528]	283	!!
	284	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	285	!!----------------------------------------------------------------------
[12377]	286	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	287	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
[14205]	288	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field
[6140]	289	!
	290	REAL(wp) :: zcoef ! local scalar
[1438]	291	!!----------------------------------------------------------------------
[3294]	292	!
[12377]	293	IF( ln_timing ) CALL timing_start('ssh_atf')
[3294]	294	!
[1438]	295	IF( kt == nit000 ) THEN
	296	IF(lwp) WRITE(numout,*)
[12377]	297	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
[1438]	298	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	299	ENDIF
[14205]	300	!
	301	IF( .NOT.l_1st_euler ) THEN ! Apply Asselin time filter on Kmm field (not on euler 1st)
[14053]	302	!
[14205]	303	IF( ln_linssh ) THEN ! filtered "now" field
	304	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
	305	!
	306	ELSE ! filtered "now" field with forcing removed
[12489]	307	zcoef = rn_atfp * rn_Dt * r1_rho0
[14205]	308	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) &
	309	& - zcoef * ( emp_b(:,:) - emp(:,:) &
	310	& - rnf_b(:,:) + rnf(:,:) &
[15004]	311	& - fwfisf_cav_b(:,:) + fwfisf_cav(:,:) &
	312	& - fwfisf_par_b(:,:) + fwfisf_par(:,:) ) * ssmask(:,:)
[12377]	313
	314	! ice sheet coupling
[14053]	315	IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 ) &
[14205]	316	& pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0._wp ) * ssmask(:,:)
[12377]	317
[6140]	318	ENDIF
[1438]	319	ENDIF
	320	!
[14053]	321	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf - pssh(:,:,Kmm): ', mask1=tmask )
[2528]	322	!
[12377]	323	IF( ln_timing ) CALL timing_stop('ssh_atf')
[3294]	324	!
[12377]	325	END SUBROUTINE ssh_atf
[3]	326
[13237]	327
[12377]	328	SUBROUTINE wAimp( kt, Kmm )
[10364]	329	!!----------------------------------------------------------------------
	330	!! * ROUTINE wAimp *
	331	!!
	332	!! ** Purpose : compute the Courant number and partition vertical velocity
	333	!! if a proportion needs to be treated implicitly
	334	!!
	335	!! ** Method : -
	336	!!
[12377]	337	!! ** action : ww : now vertical velocity (to be handled explicitly)
[10364]	338	!! : wi : now vertical velocity (for implicit treatment)
	339	!!
[11414]	340	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
	341	!! implicit scheme for vertical advection in oceanic modeling.
	342	!! Ocean Modelling, 91, 38-69.
[10364]	343	!!----------------------------------------------------------------------
	344	INTEGER, INTENT(in) :: kt ! time step
[12377]	345	INTEGER, INTENT(in) :: Kmm ! time level index
[10364]	346	!
	347	INTEGER :: ji, jj, jk ! dummy loop indices
[13237]	348	REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars
[10364]	349	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
[11407]	350	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
[10364]	351	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
	352	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
	353	!!----------------------------------------------------------------------
	354	!
	355	IF( ln_timing ) CALL timing_start('wAimp')
	356	!
	357	IF( kt == nit000 ) THEN
	358	IF(lwp) WRITE(numout,*)
	359	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
	360	IF(lwp) WRITE(numout,*) '~~~~~ '
	361	ENDIF
	362	!
[11414]	363	! Calculate Courant numbers
[13237]	364	zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability)
[11414]	365	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[15102]	366	DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
[12377]	367	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[13237]	368	Cu_adv(ji,jj,jk) = zdt * &
	369	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
	370	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
	371	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
	372	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
	373	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
[12377]	374	& * r1_e1e2t(ji,jj) &
[13237]	375	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
	376	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
	377	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
	378	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
[12377]	379	& * r1_e1e2t(ji,jj) &
	380	& ) * z1_e3t
	381	END_3D
[11414]	382	ELSE
[15102]	383	DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
[12377]	384	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[13237]	385	Cu_adv(ji,jj,jk) = zdt * &
	386	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
[12377]	387	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
	388	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
	389	& * r1_e1e2t(ji,jj) &
	390	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
	391	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
	392	& * r1_e1e2t(ji,jj) &
	393	& ) * z1_e3t
	394	END_3D
[11414]	395	ENDIF
[10364]	396	CALL iom_put("Courant",Cu_adv)
	397	!
	398	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
[15102]	399	DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary
[12377]	400	!
	401	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
[11293]	402	! alt:
[12377]	403	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
[11293]	404	! zCu = Cu_adv(ji,jj,jk)
	405	! ELSE
	406	! zCu = Cu_adv(ji,jj,jk-1)
	407	! ENDIF
[12377]	408	!
	409	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
	410	zcff = 0._wp
	411	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
	412	zcff = ( zCu - Cu_min )**2
	413	zcff = zcff / ( Fcu + zcff )
	414	ELSE !<-- Mostly implicit
	415	zcff = ( zCu - Cu_max )/ zCu
	416	ENDIF
	417	zcff = MIN(1._wp, zcff)
	418	!
	419	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
	420	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
	421	!
	422	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
	423	END_3D
[11293]	424	Cu_adv(:,:,1) = 0._wp
[10364]	425	ELSE
	426	! Fully explicit everywhere
[15102]	427	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
[10907]	428	wi (:,:,:) = 0._wp
[10364]	429	ENDIF
	430	CALL iom_put("wimp",wi)
	431	CALL iom_put("wi_cff",Cu_adv)
[12377]	432	CALL iom_put("wexp",ww)
[10364]	433	!
	434	IF( ln_timing ) CALL timing_stop('wAimp')
	435	!
	436	END SUBROUTINE wAimp
[14053]	437
[3]	438	!!======================================================================
[1565]	439	END MODULE sshwzv

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: