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 @ 13286

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

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

Download in other formats: