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 branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5974

Last change on this file since 5974 was 5974, checked in by timgraham, 8 years ago
Upgrade to head of trunk (r5936)
Property svn:keywords set to `Id`
File size: 51.9 KB

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

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

Download in other formats: