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.

domvvl.F90 in NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/DOM – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/DOM/domvvl.F90 @ 14098

Last change on this file since 14098 was 14098, checked in by nicolasmartin, 3 years ago
#2414 Sync merge with trunk
Property svn:keywords set to `Id`
File size: 53.7 KB

Rev	Line
[592]	1	MODULE domvvl
	2	!!======================================================================
	3	!! * MODULE domvvl *
[14098]	4	!! Ocean :
[592]	5	!!======================================================================
[1438]	6	!! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code
	7	!! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate
[9019]	8	!! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
[5120]	9	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
[12377]	10	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
[14066]	11	!! - ! 2020-02 (G. Madec, S. Techene) introduce ssh to h0 ratio
[592]	12	!!----------------------------------------------------------------------
[5836]	13
[592]	14	USE oce ! ocean dynamics and tracers
[6140]	15	USE phycst ! physical constant
[592]	16	USE dom_oce ! ocean space and time domain
[4292]	17	USE sbc_oce ! ocean surface boundary condition
[6152]	18	USE wet_dry ! wetting and drying
[7646]	19	USE usrdef_istate ! user defined initial state (wad only)
[6140]	20	USE restart ! ocean restart
	21	!
[592]	22	USE in_out_manager ! I/O manager
[4292]	23	USE iom ! I/O manager library
[592]	24	USE lib_mpp ! distributed memory computing library
	25	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3294]	26	USE timing ! Timing
[592]	27
[14098]	28	#if defined key_agrif
	29	USE agrif_oce ! initial state interpolation
	30	USE agrif_oce_interp
	31	#endif
	32
[592]	33	IMPLICIT NONE
	34	PRIVATE
	35
[5836]	36	! !!* Namelist nam_vvl
	37	LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
	38	LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate
	39	LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate
	40	LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate
	41	LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate
	42	LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer
	43	! ! conservation: not used yet
	44	REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient
	45	REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days]
	46	REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days]
	47	REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation
	48	LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints
[592]	49
[5836]	50	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
	51	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence
	52	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
	53	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
	54	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors
	55	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence
[1438]	56
[13237]	57	#if defined key_qco
	58	!!----------------------------------------------------------------------
	59	!! 'key_qco' EMPTY MODULE Quasi-Eulerian vertical coordonate
	60	!!----------------------------------------------------------------------
	61	#else
	62	!!----------------------------------------------------------------------
	63	!! Default key Old management of time varying vertical coordinate
	64	!!----------------------------------------------------------------------
[14098]	65
[13237]	66	!!----------------------------------------------------------------------
	67	!! dom_vvl_init : define initial vertical scale factors, depths and column thickness
	68	!! dom_vvl_sf_nxt : Compute next vertical scale factors
	69	!! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid
	70	!! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
	71	!! dom_vvl_rst : read/write restart file
	72	!! dom_vvl_ctl : Check the vvl options
	73	!!----------------------------------------------------------------------
	74
	75	PUBLIC dom_vvl_init ! called by domain.F90
	76	PUBLIC dom_vvl_zgr ! called by isfcpl.F90
	77	PUBLIC dom_vvl_sf_nxt ! called by step.F90
	78	PUBLIC dom_vvl_sf_update ! called by step.F90
	79	PUBLIC dom_vvl_interpol ! called by dynnxt.F90
[14098]	80
[592]	81	!! * Substitutions
[12377]	82	# include "do_loop_substitute.h90"
[592]	83	!!----------------------------------------------------------------------
[10068]	84	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	85	!! $Id$
[10068]	86	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	87	!!----------------------------------------------------------------------
[4292]	88	CONTAINS
	89
[2715]	90	INTEGER FUNCTION dom_vvl_alloc()
	91	!!----------------------------------------------------------------------
[4292]	92	!! * FUNCTION dom_vvl_alloc *
[2715]	93	!!----------------------------------------------------------------------
[6140]	94	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
[4292]	95	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[4338]	96	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
	97	& dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
	98	& STAT = dom_vvl_alloc )
[10425]	99	CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
	100	IF( dom_vvl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
[6140]	101	un_td = 0._wp
	102	vn_td = 0._wp
[4292]	103	ENDIF
	104	IF( ln_vvl_ztilde ) THEN
	105	ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
[10425]	106	CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
	107	IF( dom_vvl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
[4292]	108	ENDIF
[6140]	109	!
[2715]	110	END FUNCTION dom_vvl_alloc
	111
	112
[12377]	113	SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa )
[592]	114	!!----------------------------------------------------------------------
[4292]	115	!! * ROUTINE dom_vvl_init *
[14098]	116	!!
[4292]	117	!! ** Purpose : Initialization of all scale factors, depths
	118	!! and water column heights
	119	!!
	120	!! ** Method : - use restart file and/or initialize
	121	!! - interpolate scale factors
	122	!!
[6140]	123	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
[12377]	124	!! - Regrid: e3[u/v](:,:,:,Kmm)
[14098]	125	!! e3[u/v](:,:,:,Kmm)
	126	!! e3w(:,:,:,Kmm)
[12377]	127	!! e3[u/v]w_b
[14098]	128	!! e3[u/v]w_n
[12377]	129	!! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
[4292]	130	!! - h(t/u/v)_0
	131	!! - frq_rst_e3t and frq_rst_hdv
	132	!!
	133	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
[592]	134	!!----------------------------------------------------------------------
[12377]	135	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
[6140]	136	!
[4292]	137	IF(lwp) WRITE(numout,*)
	138	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
	139	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[6140]	140	!
	141	CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)
	142	!
	143	! ! Allocate module arrays
[4292]	144	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
[6140]	145	!
	146	! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
[12377]	147	CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )
	148	e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all
[6140]	149	!
[12377]	150	CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
	151	!
	152	END SUBROUTINE dom_vvl_init
[13237]	153
	154
[12377]	155	SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa)
	156	!!----------------------------------------------------------------------
	157	!! * ROUTINE dom_vvl_init *
[14098]	158	!!
	159	!! ** Purpose : Interpolation of all scale factors,
[12377]	160	!! depths and water column heights
	161	!!
	162	!! ** Method : - interpolate scale factors
	163	!!
	164	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
	165	!! - Regrid: e3(u/v)_n
[14098]	166	!! e3(u/v)_b
	167	!! e3w_n
	168	!! e3(u/v)w_b
	169	!! e3(u/v)w_n
[12377]	170	!! gdept_n, gdepw_n and gde3w_n
	171	!! - h(t/u/v)_0
	172	!! - frq_rst_e3t and frq_rst_hdv
	173	!!
	174	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
	175	!!----------------------------------------------------------------------
	176	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
	177	!!----------------------------------------------------------------------
	178	INTEGER :: ji, jj, jk
	179	INTEGER :: ii0, ii1, ij0, ij1
	180	REAL(wp):: zcoef
	181	!!----------------------------------------------------------------------
	182	!
[6140]	183	! !== Set of all other vertical scale factors ==! (now and before)
	184	! ! Horizontal interpolation of e3t
[12377]	185	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U
	186	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
[14098]	187	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V
[12377]	188	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
	189	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F
[14098]	190	! ! Vertical interpolation of e3t,u,v
[12377]	191	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W
	192	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' )
	193	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW
	194	CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
	195	CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW
	196	CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
[9449]	197
	198	! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
[12377]	199	e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
	200	e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
	201	e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
[6140]	202	!
	203	! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep)
[12377]	204	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration)
	205	gdepw(:,:,1,Kmm) = 0.0_wp
	206	gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg
	207	gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
	208	gdepw(:,:,1,Kbb) = 0.0_wp
[13497]	209	DO_3D( 1, 1, 1, 1, 2, jpk ) ! vertical sum
[12377]	210	! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	211	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
[14098]	212	! ! 0.5 where jk = mikt
[12377]	213	!!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ??
	214	zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
	215	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
	216	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) &
[14098]	217	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
[12377]	218	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
	219	gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
	220	gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) &
[14098]	221	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb))
[12377]	222	END_3D
[6140]	223	!
	224	! !== thickness of the water column !! (ocean portion only)
[12377]	225	ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....
	226	hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
	227	hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
	228	hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
	229	hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
[6140]	230	DO jk = 2, jpkm1
[12377]	231	ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
	232	hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
	233	hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
	234	hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
	235	hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
[4370]	236	END DO
[6140]	237	!
	238	! !== inverse of water column thickness ==! (u- and v- points)
[12377]	239	r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF
	240	r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
	241	r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
	242	r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
[7753]	243
[6140]	244	! !== z_tilde coordinate case ==! (Restoring frequencies)
[4292]	245	IF( ln_vvl_ztilde ) THEN
[6140]	246	!!gm : idea: add here a READ in a file of custumized restoring frequency
	247	! ! Values in days provided via the namelist
	248	! ! use rsmall to avoid possible division by zero errors with faulty settings
[7753]	249	frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp )
	250	frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
[6140]	251	!
	252	IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile
[7753]	253	frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings
[12489]	254	frq_rst_hdv(:,:) = 1._wp / rn_Dt
[4292]	255	ENDIF
[6140]	256	IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator
[13295]	257	DO_2D( 1, 1, 1, 1 )
[6140]	258	!!gm case \|gphi\| >= 6 degrees is useless initialized just above by default
[12377]	259	IF( ABS(gphit(ji,jj)) >= 6.) THEN
	260	! values outside the equatorial band and transition zone (ztilde)
	261	frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )
	262	frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
	263	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star
	264	! values inside the equatorial band (ztilde as zstar)
	265	frq_rst_e3t(ji,jj) = 0.0_wp
[12489]	266	frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt
[12377]	267	ELSE ! transition band (2.5 to 6 degrees N/S)
	268	! ! (linearly transition from z-tilde to z-star)
	269	frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp &
	270	& * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
	271	& * 180._wp / 3.5_wp ) )
[12489]	272	frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) &
	273	& + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp &
[12377]	274	& * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
	275	& * 180._wp / 3.5_wp ) )
	276	ENDIF
	277	END_2D
[10213]	278	IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
	279	IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
[14098]	280	ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1
[13286]	281	ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls
[10213]	282	frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
[12489]	283	frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt
[10213]	284	ENDIF
[4292]	285	ENDIF
	286	ENDIF
	287	ENDIF
[6140]	288	!
[12377]	289	END SUBROUTINE dom_vvl_zgr
[4292]	290
	291
[14098]	292	SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
[4292]	293	!!----------------------------------------------------------------------
	294	!! * ROUTINE dom_vvl_sf_nxt *
[14098]	295	!!
[4292]	296	!! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,
	297	!! tranxt and dynspg routines
	298	!!
	299	!! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.
[14098]	300	!! - z_tilde_case: after scale factor increment =
[4292]	301	!! high frequency part of horizontal divergence
	302	!! + retsoring towards the background grid
	303	!! + thickness difusion
	304	!! Then repartition of ssh INCREMENT proportionnaly
	305	!! to the "baroclinic" level thickness.
	306	!!
	307	!! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case
[14098]	308	!! - tilde_e3t_a: after increment of vertical scale factor
[4292]	309	!! in z_tilde case
[6140]	310	!! - e3(t/u/v)_a
[4292]	311	!!
	312	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
	313	!!----------------------------------------------------------------------
[12377]	314	INTEGER, INTENT( in ) :: kt ! time step
	315	INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step
	316	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
[6140]	317	!
	318	INTEGER :: ji, jj, jk ! dummy loop indices
	319	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
[12489]	320	REAL(wp) :: z_tmin, z_tmax ! local scalars
[6140]	321	LOGICAL :: ll_do_bclinic ! local logical
[9019]	322	REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv
[13458]	323	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ze3t
	324	LOGICAL , DIMENSION(:,:,:), ALLOCATABLE :: llmsk
[4292]	325	!!----------------------------------------------------------------------
[6140]	326	!
	327	IF( ln_linssh ) RETURN ! No calculation in linear free surface
	328	!
[9019]	329	IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt')
[6140]	330	!
	331	IF( kt == nit000 ) THEN
[4292]	332	IF(lwp) WRITE(numout,*)
	333	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
	334	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
	335	ENDIF
	336
[4338]	337	ll_do_bclinic = .TRUE.
	338	IF( PRESENT(kcall) ) THEN
[6140]	339	IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
[4338]	340	ENDIF
	341
[4292]	342	! ******************************* !
	343	! After acale factors at t-points !
	344	! ******************************* !
[4338]	345	! ! --------------------------------------------- !
[6140]	346	! ! z_star coordinate and barotropic z-tilde part !
[4338]	347	! ! --------------------------------------------- !
[6140]	348	!
[12377]	349	z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
[4338]	350	DO jk = 1, jpkm1
[12377]	351	! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
	352	e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
[4338]	353	END DO
[6140]	354	!
[11415]	355	IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
	356	! ! ------baroclinic part------ !
[4292]	357	! I - initialization
	358	! ==================
	359
	360	! 1 - barotropic divergence
	361	! -------------------------
[7753]	362	zhdiv(:,:) = 0._wp
	363	zht(:,:) = 0._wp
[4292]	364	DO jk = 1, jpkm1
[12377]	365	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
	366	zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
[592]	367	END DO
[7753]	368	zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
[2528]	369
[4292]	370	! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only)
	371	! --------------------------------------------------
	372	IF( ln_vvl_ztilde ) THEN
[6140]	373	IF( kt > nit000 ) THEN
[4292]	374	DO jk = 1, jpkm1
[12489]	375	hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) &
[12377]	376	& * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
[4292]	377	END DO
	378	ENDIF
[6140]	379	ENDIF
[3294]	380
[4292]	381	! II - after z_tilde increments of vertical scale factors
	382	! =======================================================
[7753]	383	tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms
[4292]	384
	385	! 1 - High frequency divergence term
	386	! ----------------------------------
	387	IF( ln_vvl_ztilde ) THEN ! z_tilde case
	388	DO jk = 1, jpkm1
[12377]	389	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
[4292]	390	END DO
	391	ELSE ! layer case
	392	DO jk = 1, jpkm1
[12377]	393	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
[4292]	394	END DO
[6140]	395	ENDIF
[4292]	396
	397	! 2 - Restoring term (z-tilde case only)
	398	! ------------------
	399	IF( ln_vvl_ztilde ) THEN
	400	DO jk = 1, jpk
[7753]	401	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
[4292]	402	END DO
[6140]	403	ENDIF
[4292]	404
	405	! 3 - Thickness diffusion term
	406	! ----------------------------
[7753]	407	zwu(:,:) = 0._wp
	408	zwv(:,:) = 0._wp
[13497]	409	DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! a - first derivative: diffusive fluxes
[12377]	410	un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) &
	411	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )
[14098]	412	vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) &
[12377]	413	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )
	414	zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
	415	zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
	416	END_3D
[13497]	417	DO_2D( 1, 1, 1, 1 ) ! b - correction for last oceanic u-v points
[12377]	418	un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
	419	vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
	420	END_2D
[13497]	421	DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! c - second derivative: divergence of diffusive fluxes
[12377]	422	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &
	423	& + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) &
	424	& ) * r1_e1e2t(ji,jj)
	425	END_3D
[13497]	426	! ! d - thickness diffusion transport: boundary conditions
[6140]	427	! (stored for tracer advction and continuity equation)
[10425]	428	CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
[4292]	429	! 4 - Time stepping of baroclinic scale factors
	430	! ---------------------------------------------
[10425]	431	CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
[12489]	432	tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
[4292]	433
	434	! Maximum deformation control
	435	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
[13458]	436	ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) )
	437	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	438	ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
	439	END_3D
	440	!
	441	llmsk( 1:Nis1,:,:) = .FALSE. ! exclude halos from the checked region
	442	llmsk(Nie1: jpi,:,:) = .FALSE.
	443	llmsk(:, 1:Njs1,:) = .FALSE.
	444	llmsk(:,Nje1: jpj,:) = .FALSE.
	445	!
	446	llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp ! define only the inner domain
	447	z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk ) ; CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
	448	z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk ) ; CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain
[4292]	449	! - ML - test: for the moment, stop simulation for too large e3_t variations
[6140]	450	IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
[13458]	451	CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max )
	452	CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min )
[4292]	453	IF (lwp) THEN
	454	WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
	455	WRITE(numout, *) 'at i, j, k=', ijk_max
	456	WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
[14098]	457	WRITE(numout, *) 'at i, j, k=', ijk_min
[10425]	458	CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
[4292]	459	ENDIF
	460	ENDIF
[13458]	461	DEALLOCATE( ze3t, llmsk )
[4292]	462	! - ML - end test
	463	! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
[7753]	464	tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )
	465	tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
[4292]	466
[4338]	467	!
	468	! "tilda" change in the after scale factor
[4292]	469	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[4338]	470	DO jk = 1, jpkm1
[7753]	471	dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
[4338]	472	END DO
[4292]	473	! III - Barotropic repartition of the sea surface height over the baroclinic profile
	474	! ==================================================================================
[4338]	475	! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
[4292]	476	! - ML - baroclinicity error should be better treated in the future
	477	! i.e. locally and not spread over the water column.
	478	! (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
[7753]	479	zht(:,:) = 0.
[4292]	480	DO jk = 1, jpkm1
[7753]	481	zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
[4292]	482	END DO
[12377]	483	z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
[4292]	484	DO jk = 1, jpkm1
[12377]	485	dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
[4292]	486	END DO
[7753]	487
[4292]	488	ENDIF
	489
[4338]	490	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
	491	! ! ---baroclinic part--------- !
	492	DO jk = 1, jpkm1
[12377]	493	e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
[4338]	494	END DO
	495	ENDIF
	496
	497	IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging
[4292]	498	!
	499	IF( lwp ) WRITE(numout, *) 'kt =', kt
	500	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
	501	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
[10425]	502	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[4292]	503	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
	504	END IF
	505	!
[7753]	506	zht(:,:) = 0.0_wp
[4292]	507	DO jk = 1, jpkm1
[12377]	508	zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
[4292]	509	END DO
[12377]	510	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
[10425]	511	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[12377]	512	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
[4292]	513	!
[7753]	514	zht(:,:) = 0.0_wp
[4292]	515	DO jk = 1, jpkm1
[12377]	516	zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
[4292]	517	END DO
[12377]	518	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
[10425]	519	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[12377]	520	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
[4292]	521	!
[7753]	522	zht(:,:) = 0.0_wp
[4292]	523	DO jk = 1, jpkm1
[12377]	524	zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
[4292]	525	END DO
[12377]	526	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
[10425]	527	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[12377]	528	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
[4292]	529	!
[12377]	530	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) )
[10425]	531	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[12377]	532	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
[4292]	533	!
[12377]	534	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) )
[10425]	535	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[12377]	536	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
[4292]	537	!
[12377]	538	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) )
[10425]	539	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
[12377]	540	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
[4292]	541	END IF
	542
	543	! *********************************** !
	544	! After scale factors at u- v- points !
	545	! *********************************** !
	546
[12377]	547	CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
	548	CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
[4292]	549
[4370]	550	! *********************************** !
	551	! After depths at u- v points !
	552	! *********************************** !
	553
[12377]	554	hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
	555	hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
[6140]	556	DO jk = 2, jpkm1
[12377]	557	hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
	558	hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
[4370]	559	END DO
	560	! ! Inverse of the local depth
[6140]	561	!!gm BUG ? don't understand the use of umask_i here .....
[12377]	562	r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
	563	r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
[6140]	564	!
[9019]	565	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt')
[6140]	566	!
[4292]	567	END SUBROUTINE dom_vvl_sf_nxt
	568
	569
[12377]	570	SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
[3294]	571	!!----------------------------------------------------------------------
[12377]	572	!! * ROUTINE dom_vvl_sf_update *
[14098]	573	!!
	574	!! ** Purpose : for z tilde case: compute time filter and swap of scale factors
[4292]	575	!! compute all depths and related variables for next time step
	576	!! write outputs and restart file
[3294]	577	!!
[12377]	578	!! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
[4292]	579	!! - reconstruct scale factor at other grid points (interpolate)
	580	!! - recompute depths and water height fields
	581	!!
[12377]	582	!! ** Action : - tilde_e3t_(b/n) ready for next time step
[4292]	583	!! - Recompute:
[14098]	584	!! e3(u/v)_b
	585	!! e3w(:,:,:,Kmm)
	586	!! e3(u/v)w_b
	587	!! e3(u/v)w_n
[12377]	588	!! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
[4292]	589	!! h(u/v) and h(u/v)r
	590	!!
	591	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	592	!! Leclair, M., and G. Madec, 2011, Ocean Modelling.
[3294]	593	!!----------------------------------------------------------------------
[12377]	594	INTEGER, INTENT( in ) :: kt ! time step
	595	INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
[6140]	596	!
	597	INTEGER :: ji, jj, jk ! dummy loop indices
	598	REAL(wp) :: zcoef ! local scalar
[3294]	599	!!----------------------------------------------------------------------
[6140]	600	!
	601	IF( ln_linssh ) RETURN ! No calculation in linear free surface
	602	!
[12377]	603	IF( ln_timing ) CALL timing_start('dom_vvl_sf_update')
[3294]	604	!
[4292]	605	IF( kt == nit000 ) THEN
	606	IF(lwp) WRITE(numout,*)
[12377]	607	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
	608	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
[3294]	609	ENDIF
[4292]	610	!
	611	! Time filter and swap of scale factors
	612	! =====================================
[6140]	613	! - ML - e3(t/u/v)_b are allready computed in dynnxt.
[4292]	614	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[12489]	615	IF( l_1st_euler ) THEN
[7753]	616	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
[4292]	617	ELSE
[14098]	618	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
[12489]	619	& + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
[4292]	620	ENDIF
[7753]	621	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
[4292]	622	ENDIF
[4488]	623
[4292]	624	! Compute all missing vertical scale factor and depths
	625	! ====================================================
	626	! Horizontal scale factor interpolations
	627	! --------------------------------------
[12377]	628	! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
	629	! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
[14098]	630
[12377]	631	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
[14098]	632
[4292]	633	! Vertical scale factor interpolations
[12377]	634	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )
	635	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
	636	CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
	637	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' )
	638	CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
	639	CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
[5120]	640
[6140]	641	! t- and w- points depth (set the isf depth as it is in the initial step)
[12377]	642	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
	643	gdepw(:,:,1,Kmm) = 0.0_wp
	644	gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
[13295]	645	DO_3D( 1, 1, 1, 1, 2, jpk )
[12377]	646	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	647	! 1 for jk = mikt
	648	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
	649	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
	650	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) &
[14098]	651	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) )
[12377]	652	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
	653	END_3D
[5120]	654
[6140]	655	! Local depth and Inverse of the local depth of the water
	656	! -------------------------------------------------------
[7753]	657	!
[12377]	658	ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
[6140]	659	DO jk = 2, jpkm1
[12377]	660	ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
[4370]	661	END DO
[7753]	662
[4292]	663	! write restart file
	664	! ==================
[12377]	665	IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
[4292]	666	!
[12377]	667	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update')
[6140]	668	!
[12377]	669	END SUBROUTINE dom_vvl_sf_update
[4292]	670
	671
	672	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
	673	!!---------------------------------------------------------------------
	674	!! * ROUTINE dom_vvl__interpol *
	675	!!
	676	!! ** Purpose : interpolate scale factors from one grid point to another
	677	!!
	678	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
	679	!! - horizontal interpolation: grid cell surface averaging
	680	!! - vertical interpolation: simple averaging
	681	!!----------------------------------------------------------------------
[5836]	682	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated
	683	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3
	684	CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors
	685	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
	686	!
[6152]	687	INTEGER :: ji, jj, jk ! dummy loop indices
[9023]	688	REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F
[4292]	689	!!----------------------------------------------------------------------
[5836]	690	!
[9023]	691	IF(ln_wd_il) THEN
[6152]	692	zlnwd = 1.0_wp
	693	ELSE
	694	zlnwd = 0.0_wp
	695	END IF
	696	!
[5836]	697	SELECT CASE ( pout ) !== type of interpolation ==!
[4292]	698	!
[5836]	699	CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean
[13295]	700	DO_3D( 1, 0, 1, 0, 1, jpk )
[12377]	701	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
	702	& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	703	& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
	704	END_3D
[10425]	705	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
[7753]	706	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
[5836]	707	!
	708	CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean
[13295]	709	DO_3D( 1, 0, 1, 0, 1, jpk )
[12377]	710	pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
	711	& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	712	& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
	713	END_3D
[10425]	714	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
[7753]	715	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
[5836]	716	!
	717	CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean
[13295]	718	DO_3D( 1, 0, 1, 0, 1, jpk )
[12377]	719	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	720	& * r1_e1e2f(ji,jj) &
	721	& * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
	722	& + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
	723	END_3D
[10425]	724	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
[7753]	725	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
[5836]	726	!
	727	CASE( 'W' ) !* from T- to W-point : vertical simple mean
	728	!
[7753]	729	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
[5836]	730	! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
	731	!!gm BUG? use here wmask in case of ISF ? to be checked
[4292]	732	DO jk = 2, jpk
[7753]	733	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
	734	& * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
	735	& + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	736	& * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
[4292]	737	END DO
[5836]	738	!
	739	CASE( 'UW' ) !* from U- to UW-point : vertical simple mean
	740	!
[7753]	741	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
[4292]	742	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
[5836]	743	!!gm BUG? use here wumask in case of ISF ? to be checked
[4292]	744	DO jk = 2, jpk
[7753]	745	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
	746	& * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
	747	& + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	748	& * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
[4292]	749	END DO
[5836]	750	!
	751	CASE( 'VW' ) !* from V- to VW-point : vertical simple mean
	752	!
[7753]	753	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
[4292]	754	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
[5836]	755	!!gm BUG? use here wvmask in case of ISF ? to be checked
[4292]	756	DO jk = 2, jpk
[7753]	757	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
	758	& * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
	759	& + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	760	& * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
[4292]	761	END DO
	762	END SELECT
	763	!
	764	END SUBROUTINE dom_vvl_interpol
	765
[5836]	766
[12377]	767	SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
[4292]	768	!!---------------------------------------------------------------------
	769	!! * ROUTINE dom_vvl_rst *
[14098]	770	!!
[4292]	771	!! ** Purpose : Read or write VVL file in restart file
	772	!!
[14066]	773	!! ** Method : * restart comes from a linear ssh simulation :
	774	!! an attempt to read e3t_n stops simulation
	775	!! * restart comes from a z-star, z-tilde, or layer :
	776	!! read e3t_n and e3t_b
	777	!! * restart comes from a z-star :
	778	!! set tilde_e3t_n, tilde_e3t_n, and hdiv_lf to 0
	779	!! * restart comes from layer :
	780	!! read tilde_e3t_n and tilde_e3t_b
	781	!! set hdiv_lf to 0
	782	!! * restart comes from a z-tilde:
	783	!! read tilde_e3t_n, tilde_e3t_b, and hdiv_lf
	784	!!
	785	!! NB: if l_1st_euler = T (ln_1st_euler or ssh_b not found)
	786	!! Kbb fields set to Kmm ones
[4292]	787	!!----------------------------------------------------------------------
[12377]	788	INTEGER , INTENT(in) :: kt ! ocean time-step
	789	INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
	790	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[5836]	791	!
[14066]	792	INTEGER :: ji, jj, jk ! dummy loop indices
	793	INTEGER :: id3, id4, id5 ! local integers
[4292]	794	!!----------------------------------------------------------------------
	795	!
[14066]	796	! !=====================!
	797	IF( TRIM(cdrw) == 'READ' ) THEN ! Read / initialise !
	798	! !=====================!
	799	!
	800	IF( ln_rstart ) THEN !== Read the restart file ==!
[4366]	801	!
[14098]	802	#if defined key_agrif
	803	IF ( (.NOT.Agrif_root()).AND.(ln_init_chfrpar) ) THEN
	804	! skip reading restart if initialized from parent:
	805	id3 = -1 ; id4 = -1 ; id5 = -1
	806	ELSE
	807	#endif
[14066]	808	CALL rst_read_open !* open the restart file if necessary
	809	! ! --------- !
	810	! ! all cases !
	811	! ! --------- !
	812	!
	813	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) !* check presence
[4292]	814	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
[14066]	815	id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. )
[12377]	816	!
[14098]	817	#if defined key_agrif
	818	ENDIF
	819	#endif
[14066]	820	! !* scale factors
	821	IF(lwp) WRITE(numout,*) ' Kmm scale factor read in the restart file'
	822	CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
[14098]	823	WHERE ( tmask(:,:,:) == 0.0_wp )
[14066]	824	e3t(:,:,:,Kmm) = e3t_0(:,:,:)
	825	END WHERE
	826	IF( l_1st_euler ) THEN ! euler
	827	IF(lwp) WRITE(numout,*) ' Euler first time step : e3t(Kbb) = e3t(Kmm)'
	828	e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
	829	ELSE ! leap frog
	830	IF(lwp) WRITE(numout,*) ' Kbb scale factor read in the restart file'
	831	CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) )
[14098]	832	WHERE ( tmask(:,:,:) == 0.0_wp )
[12377]	833	e3t(:,:,:,Kbb) = e3t_0(:,:,:)
[4990]	834	END WHERE
[4292]	835	ENDIF
[14066]	836	! ! ------------ !
	837	IF( ln_vvl_zstar ) THEN ! z_star case !
	838	! ! ------------ !
[4292]	839	IF( MIN( id3, id4 ) > 0 ) THEN
	840	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
	841	ENDIF
[14066]	842	! ! ------------------------ !
	843	ELSE ! z_tilde and layer cases !
	844	! ! ------------------------ !
	845	!
	846	IF( id4 > 0 ) THEN !* scale factor increments
	847	IF(lwp) WRITE(numout,*) ' Kmm scale factor increments read in the restart file'
	848	CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
	849	IF( l_1st_euler ) THEN ! euler
	850	IF(lwp) WRITE(numout,*) ' Euler first time step : tilde_e3t(Kbb) = tilde_e3t(Kmm)'
	851	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
	852	ELSE ! leap frog
	853	IF(lwp) WRITE(numout,*) ' Kbb scale factor increments read in the restart file'
	854	CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
	855	ENDIF
[14098]	856	ELSE
	857	tilde_e3t_b(:,:,:) = 0.0_wp
	858	tilde_e3t_n(:,:,:) = 0.0_wp
[4292]	859	ENDIF
[14066]	860	! ! ------------ !
	861	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	862	! ! ------------ !
[4292]	863	IF( id5 > 0 ) THEN ! required array exists
[14066]	864	CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:) )
[4292]	865	ELSE ! array is missing
[14098]	866	hdiv_lf(:,:,:) = 0.0_wp
[4292]	867	ENDIF
	868	ENDIF
	869	ENDIF
	870	!
[14066]	871	ELSE !== Initialize at "rest" with ssh ==!
[7646]	872	!
[14066]	873	DO jk = 1, jpk
	874	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) )
	875	END DO
	876	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
	877	!
[4292]	878	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
[7646]	879	tilde_e3t_b(:,:,:) = 0._wp
	880	tilde_e3t_n(:,:,:) = 0._wp
	881	IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
[14066]	882	ENDIF
[4292]	883	ENDIF
[14066]	884	! !=======================!
	885	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file !
	886	! !=======================!
[14098]	887	#if defined key_agrif
	888	IF ( .NOT.Agrif_root().AND.(ln_init_chfrpar) ) THEN
	889	! Interpolate initial ssh from parent:
	890	CALL Agrif_istate_ssh( Kbb, Kmm )
	891	!
	892	DO jk = 1, jpk
	893	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
	894	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
	895	& + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
	896	END DO
	897	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
	898	ENDIF
	899	#endif
[5836]	900	!
[4292]	901	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
	902	! ! --------- !
	903	! ! all cases !
	904	! ! --------- !
[14066]	905	CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb) )
	906	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm) )
[4292]	907	! ! ----------------------- !
	908	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
	909	! ! ----------------------- !
[14066]	910	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:))
	911	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:))
[4292]	912	END IF
[14098]	913	! ! -------------!
[4292]	914	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	915	! ! ------------ !
[14066]	916	CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:))
[4292]	917	ENDIF
[5836]	918	!
[4292]	919	ENDIF
[5836]	920	!
[4292]	921	END SUBROUTINE dom_vvl_rst
	922
	923
	924	SUBROUTINE dom_vvl_ctl
	925	!!---------------------------------------------------------------------
	926	!! * ROUTINE dom_vvl_ctl *
[14098]	927	!!
[4292]	928	!! ** Purpose : Control the consistency between namelist options
	929	!! for vertical coordinate
	930	!!----------------------------------------------------------------------
[5836]	931	INTEGER :: ioptio, ios
	932	!!
[4292]	933	NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
[5836]	934	& ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &
	935	& rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe
[14098]	936	!!----------------------------------------------------------------------
[5836]	937	!
[4294]	938	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
[11536]	939	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
[4294]	940	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
[11536]	941	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
[4624]	942	IF(lwm) WRITE ( numond, nam_vvl )
[5836]	943	!
[4292]	944	IF(lwp) THEN ! Namelist print
	945	WRITE(numout,*)
	946	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
	947	WRITE(numout,*) '~~~~~~~~~~~'
[9168]	948	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
	949	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
	950	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
	951	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
	952	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
[4292]	953	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
[9168]	954	WRITE(numout,*) ' !'
	955	WRITE(numout,*) ' thickness diffusion coefficient rn_ahe3 = ', rn_ahe3
	956	WRITE(numout,*) ' maximum e3t deformation fractional change rn_zdef_max = ', rn_zdef_max
[4292]	957	IF( ln_vvl_ztilde_as_zstar ) THEN
[9168]	958	WRITE(numout,*) ' ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
	959	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
	960	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
	961	WRITE(numout,*) ' rn_rst_e3t = 0.e0'
	962	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
[12489]	963	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rn_Dt'
[4292]	964	ELSE
[9168]	965	WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t
	966	WRITE(numout,*) ' z-tilde cutoff frequency of low-pass filter (days) rn_lf_cutoff = ', rn_lf_cutoff
[4292]	967	ENDIF
[9168]	968	WRITE(numout,*) ' debug prints flag ln_vvl_dbg = ', ln_vvl_dbg
[4292]	969	ENDIF
[5836]	970	!
[4292]	971	ioptio = 0 ! Parameter control
[5836]	972	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
	973	IF( ln_vvl_zstar ) ioptio = ioptio + 1
	974	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
	975	IF( ln_vvl_layer ) ioptio = ioptio + 1
	976	!
[4292]	977	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
[5836]	978	!
[4292]	979	IF(lwp) THEN ! Print the choice
	980	WRITE(numout,*)
[9168]	981	IF( ln_vvl_zstar ) WRITE(numout,*) ' ==>>> zstar vertical coordinate is used'
	982	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ==>>> ztilde vertical coordinate is used'
	983	IF( ln_vvl_layer ) WRITE(numout,*) ' ==>>> layer vertical coordinate is used'
	984	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' ==>>> to emulate a zstar coordinate'
[4292]	985	ENDIF
[5836]	986	!
[4486]	987	#if defined key_agrif
[9190]	988	IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
[4486]	989	#endif
[5836]	990	!
[4292]	991	END SUBROUTINE dom_vvl_ctl
	992
[13237]	993	#endif
	994
[592]	995	!!======================================================================
	996	END MODULE domvvl

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

Download in other formats: