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.

domzgr.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago
Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
Property svn:keywords set to `Id`
File size: 88.4 KB

Rev	Line
[3]	1	MODULE domzgr
	2	!!==============================================================================
	3	!! * MODULE domzgr *
	4	!! Ocean initialization : domain initialization
	5	!!==============================================================================
[1566]	6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
	7	!! ! 1997-07 (G. Madec) lbc_lnk call
	8	!! ! 1997-04 (J.-O. Beismann)
[2528]	9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
	10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
[1566]	11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
	12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
	13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
	14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
	15	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
[2528]	16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
[1099]	17	!!----------------------------------------------------------------------
[3]	18
	19	!!----------------------------------------------------------------------
[1099]	20	!! dom_zgr : defined the ocean vertical coordinate system
[3]	21	!! zgr_bat : bathymetry fields (levels and meters)
	22	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
	23	!! zgr_bat_ctl : check the bathymetry files
[2528]	24	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
[3]	25	!! zgr_z : reference z-coordinate
[454]	26	!! zgr_zco : z-coordinate
[3]	27	!! zgr_zps : z-coordinate with partial steps
[454]	28	!! zgr_sco : s-coordinate
[1083]	29	!! fssig : sigma coordinate non-dimensional function
[1099]	30	!! dfssig : derivative of the sigma coordinate function !!gm (currently missing!)
[3]	31	!!---------------------------------------------------------------------
[2528]	32	USE oce ! ocean variables
	33	USE dom_oce ! ocean domain
	34	USE closea ! closed seas
	35	USE c1d ! 1D vertical configuration
	36	USE in_out_manager ! I/O manager
	37	USE iom ! I/O library
	38	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	39	USE lib_mpp ! distributed memory computing library
[3]	40
	41	IMPLICIT NONE
	42	PRIVATE
	43
[2715]	44	PUBLIC dom_zgr ! called by dom_init.F90
[3837]	45	PUBLIC zgr_z, zgr_bat, zgr_zco, zgr_zps ! called by nemogcm::recursive_partition
	46	PUBLIC fssig1 ! called by partition_mod::smooth_bathy
[3]	47
[2528]	48	! !!* Namelist namzgr_sco *
	49	REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m)
	50	REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
	51	REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20)
	52	REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1)
	53	REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1)
	54	LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
	55	REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter for song and haidvogel stretching
	56	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
	57	REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates
[3837]	58	PUBLIC rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, rn_rmax, &
	59	ln_s_sigma, rn_bb, rn_hc
	60	PUBLIC ln_zco, ln_zps, ln_sco
	61
[3211]	62	!! * Control permutation of array indices
	63	# include "oce_ftrans.h90"
	64	# include "dom_oce_ftrans.h90"
	65
[2715]	66	!! * Substitutions
[3]	67	# include "domzgr_substitute.h90"
	68	# include "vectopt_loop_substitute.h90"
[3837]	69
	70	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
	71	NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, &
	72	rn_rmax, ln_s_sigma, rn_bb, rn_hc
	73	PUBLIC namzgr, namzgr_sco
[3]	74	!!----------------------------------------------------------------------
[2715]	75	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
[1146]	76	!! $Id$
[2528]	77	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	78	!!----------------------------------------------------------------------
	79	CONTAINS
	80
	81	SUBROUTINE dom_zgr
	82	!!----------------------------------------------------------------------
	83	!! * ROUTINE dom_zgr *
	84	!!
	85	!! ** Purpose : set the depth of model levels and the resulting
	86	!! vertical scale factors.
	87	!!
[1099]	88	!! ** Method : - reference 1D vertical coordinate (gdep._0, e3._0)
	89	!! - read/set ocean depth and ocean levels (bathy, mbathy)
	90	!! - vertical coordinate (gdep., e3.) depending on the
	91	!! coordinate chosen :
[2528]	92	!! ln_zco=T z-coordinate
[1099]	93	!! ln_zps=T z-coordinate with partial steps
[3432]	94	!! ln_zco=T s-coordinate
[3]	95	!!
[1099]	96	!! ** Action : define gdep., e3., mbathy and bathy
	97	!!----------------------------------------------------------------------
	98	INTEGER :: ioptio = 0 ! temporary integer
[2528]	99	!
[3837]	100	!NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
[3]	101	!!----------------------------------------------------------------------
	102
[3837]	103	REWIND( numnam ) ! Read Namelist namzgr : vertical coordinate
[2528]	104	READ ( numnam, namzgr )
[454]	105
[1099]	106	IF(lwp) THEN ! Control print
[454]	107	WRITE(numout,*)
	108	WRITE(numout,*) 'dom_zgr : vertical coordinate'
	109	WRITE(numout,*) '~~~~~~~'
[1601]	110	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
[454]	111	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
	112	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
	113	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
	114	ENDIF
	115
[1099]	116	ioptio = 0 ! Check Vertical coordinate options
[454]	117	IF( ln_zco ) ioptio = ioptio + 1
	118	IF( ln_zps ) ioptio = ioptio + 1
	119	IF( ln_sco ) ioptio = ioptio + 1
[2528]	120	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
	121	!
[3]	122	! Build the vertical coordinate system
	123	! ------------------------------------
[2528]	124	CALL zgr_z ! Reference z-coordinate system (always called)
	125	CALL zgr_bat ! Bathymetry fields (levels and meters)
	126	IF( ln_zco ) CALL zgr_zco ! z-coordinate
	127	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
	128	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
[2465]	129	!
[2528]	130	! final adjustment of mbathy & check
	131	! -----------------------------------
	132	IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain
	133	IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isoated ocean points
	134	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
	135	!
	136	!
[2715]	137
[1348]	138	IF( nprint == 1 .AND. lwp ) THEN
	139	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	140	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	141	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
	142	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), &
	143	& ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), &
	144	& ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), &
	145	& ' w ', MINVAL( fse3w(:,:,:) )
	146
	147	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	148	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
	149	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), &
	150	& ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), &
	151	& ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), &
	152	& ' w ', MAXVAL( fse3w(:,:,:) )
	153	ENDIF
[2528]	154	!
[3]	155	END SUBROUTINE dom_zgr
	156
	157
	158	SUBROUTINE zgr_z
	159	!!----------------------------------------------------------------------
	160	!! * ROUTINE zgr_z *
	161	!!
	162	!! ** Purpose : set the depth of model levels and the resulting
	163	!! vertical scale factors.
	164	!!
	165	!! ** Method : z-coordinate system (use in all type of coordinate)
	166	!! The depth of model levels is defined from an analytical
	167	!! function the derivative of which gives the scale factors.
	168	!! both depth and scale factors only depend on k (1d arrays).
[454]	169	!! w-level: gdepw_0 = fsdep(k)
	170	!! e3w_0(k) = dk(fsdep)(k) = fse3(k)
	171	!! t-level: gdept_0 = fsdep(k+0.5)
	172	!! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
[3]	173	!!
[454]	174	!! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m)
[1099]	175	!! - e3t_0 , e3w_0 : scale factors at T- and W-levels (m)
[3]	176	!!
[1099]	177	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
[3]	178	!!----------------------------------------------------------------------
	179	INTEGER :: jk ! dummy loop indices
	180	REAL(wp) :: zt, zw ! temporary scalars
[1099]	181	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
	182	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
[1577]	183	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
[2528]	184	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
[3]	185	!!----------------------------------------------------------------------
	186
	187	! Set variables from parameters
	188	! ------------------------------
	189	zkth = ppkth ; zacr = ppacr
	190	zdzmin = ppdzmin ; zhmax = pphmax
[2528]	191	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
[3]	192
	193	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
	194	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
[1099]	195	IF( ppa1 == pp_to_be_computed .AND. &
[3]	196	& ppa0 == pp_to_be_computed .AND. &
	197	& ppsur == pp_to_be_computed ) THEN
[1099]	198	!
	199	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
	200	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
	201	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
	202	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
	203	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
	204	ELSE
[3]	205	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
[2528]	206	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
[1099]	207	ENDIF
[3]	208
[1099]	209	IF(lwp) THEN ! Parameter print
[3]	210	WRITE(numout,*)
	211	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
	212	WRITE(numout,*) ' ~~~~~~~'
[2528]	213	IF( ppkth == 0._wp ) THEN
[250]	214	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
	215	WRITE(numout,*) ' Total depth :', zhmax
	216	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
	217	ELSE
[2528]	218	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
[250]	219	WRITE(numout,*) ' zsur, za0, za1 computed from '
	220	WRITE(numout,*) ' zdzmin = ', zdzmin
	221	WRITE(numout,*) ' zhmax = ', zhmax
	222	ENDIF
	223	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
	224	WRITE(numout,*) ' zsur = ', zsur
	225	WRITE(numout,*) ' za0 = ', za0
	226	WRITE(numout,*) ' za1 = ', za1
	227	WRITE(numout,*) ' zkth = ', zkth
	228	WRITE(numout,*) ' zacr = ', zacr
[2528]	229	IF( ldbletanh ) THEN
	230	WRITE(numout,*) ' (Double tanh za2 = ', za2
	231	WRITE(numout,*) ' parameters) zkth2= ', zkth2
	232	WRITE(numout,*) ' zacr2= ', zacr2
	233	ENDIF
[3]	234	ENDIF
	235	ENDIF
	236
	237
	238	! Reference z-coordinate (depth - scale factor at T- and W-points)
	239	! ======================
[2528]	240	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
[454]	241	za1 = zhmax / FLOAT(jpk-1)
[250]	242	DO jk = 1, jpk
	243	zw = FLOAT( jk )
[2528]	244	zt = FLOAT( jk ) + 0.5_wp
[454]	245	gdepw_0(jk) = ( zw - 1 ) * za1
	246	gdept_0(jk) = ( zt - 1 ) * za1
	247	e3w_0 (jk) = za1
	248	e3t_0 (jk) = za1
[250]	249	END DO
[1099]	250	ELSE ! Madec & Imbard 1996 function
[2528]	251	IF( .NOT. ldbletanh ) THEN
	252	DO jk = 1, jpk
	253	zw = REAL( jk , wp )
	254	zt = REAL( jk , wp ) + 0.5_wp
	255	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
	256	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
	257	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
	258	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
	259	END DO
	260	ELSE
	261	DO jk = 1, jpk
	262	zw = FLOAT( jk )
	263	zt = FLOAT( jk ) + 0.5_wp
	264	! Double tanh function
	265	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
	266	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
	267	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
	268	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
	269	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
	270	& + za2 * TANH( (zw-zkth2) / zacr2 )
	271	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
	272	& + za2 * TANH( (zt-zkth2) / zacr2 )
	273	END DO
	274	ENDIF
	275	gdepw_0(1) = 0._wp ! force first w-level to be exactly at zero
[250]	276	ENDIF
	277
[1601]	278	!!gm BUG in s-coordinate this does not work!
[2528]	279	! deepest/shallowest W level Above/Below ~10m
	280	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 ) ! ref. depth with tolerance (10% of minimum layer thickness)
	281	nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
	282	nla10 = nlb10 - 1 ! deepest W level Above ~10m
[1601]	283	!!gm end bug
[1577]	284
[1099]	285	IF(lwp) THEN ! control print
[3]	286	WRITE(numout,*)
	287	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
	288	WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" )
[454]	289	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
[3]	290	ENDIF
[1099]	291	DO jk = 1, jpk ! control positivity
[2528]	292	IF( e3w_0 (jk) <= 0._wp .OR. e3t_0 (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 ' )
	293	IF( gdepw_0(jk) < 0._wp .OR. gdept_0(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' )
[3]	294	END DO
[1099]	295	!
[3]	296	END SUBROUTINE zgr_z
	297
	298
[3837]	299	SUBROUTINE zgr_bat(global_domain)
[3]	300	!!----------------------------------------------------------------------
	301	!! * ROUTINE zgr_bat *
	302	!!
	303	!! ** Purpose : set bathymetry both in levels and meters
	304	!!
	305	!! ** Method : read or define mbathy and bathy arrays
	306	!! * level bathymetry:
	307	!! The ocean basin geometry is given by a two-dimensional array,
	308	!! mbathy, which is defined as follow :
	309	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
	310	!! at t-point (ji,jj).
	311	!! = 0 over the continental t-point.
	312	!! The array mbathy is checked to verified its consistency with
	313	!! model option. in particular:
	314	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	315	!! along closed boundary.
	316	!! mbathy must be cyclic IF jperio=1.
	317	!! mbathy must be lower or equal to jpk-1.
	318	!! isolated ocean grid points are suppressed from mbathy
	319	!! since they are only connected to remaining
	320	!! ocean through vertical diffusion.
	321	!! ntopo=-1 : rectangular channel or bassin with a bump
	322	!! ntopo= 0 : flat rectangular channel or basin
[128]	323	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
[3]	324	!! bathy is read in 'bathy_meter.nc' NetCDF file
	325	!!
	326	!! ** Action : - mbathy: level bathymetry (in level index)
	327	!! - bathy : meter bathymetry (in meters)
	328	!!----------------------------------------------------------------------
[3837]	329	LOGICAL, OPTIONAL, INTENT(in) :: global_domain ! Whether dealing with
	330	! whole domain (T) or a
	331	! sub-domain after domain
	332	! decomposition
	333	! Locals
[1099]	334	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	335	INTEGER :: inum ! temporary logical unit
[1348]	336	INTEGER :: ii_bump, ij_bump, ih ! bump center position
[2528]	337	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
[1099]	338	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
[2528]	339	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
[1099]	340	INTEGER , DIMENSION(jpidta,jpjdta) :: idta ! global domain integer data
	341	REAL(wp), DIMENSION(jpidta,jpjdta) :: zdta ! global domain scalar data
[3837]	342	LOGICAL :: is_global
[3]	343	!!----------------------------------------------------------------------
	344
	345	IF(lwp) WRITE(numout,*)
	346	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
	347	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
	348
[3837]	349	! Set local flag to signal whether we're dealing with the global domain
	350	! (pre decomposition) or a local part of it. Required by the
	351	! recursive k-section partitioning.
	352	is_global = .FALSE.
	353	IF( PRESENT(global_domain) )THEN
	354	IF( global_domain ) is_global = .TRUE.
	355	END IF
	356
[1099]	357	! ! ================== !
	358	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
	359	! ! ================== !
	360	! ! global domain level and meter bathymetry (idta,zdta)
	361	!
[3]	362	IF( ntopo == 0 ) THEN ! flat basin
	363	IF(lwp) WRITE(numout,*)
	364	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
[1099]	365	idta(:,:) = jpkm1 ! before last level
	366	zdta(:,:) = gdepw_0(jpk) ! last w-point depth
[592]	367	h_oce = gdepw_0(jpk)
[1099]	368	ELSE ! bump centered in the basin
[3]	369	IF(lwp) WRITE(numout,*)
	370	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
[1099]	371	ii_bump = jpidta / 2 ! i-index of the bump center
	372	ij_bump = jpjdta / 2 ! j-index of the bump center
[3837]	373	r_bump = 0.165*MIN(jpidta,jpjdta) ! bump radius (grid cells)
	374	h_bump = 3000._wp ! bump height (meters)
[1099]	375	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
[3]	376	IF(lwp) WRITE(numout,*) ' bump characteristics: '
[3837]	377	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ij_bump
[3]	378	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
[3837]	379	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' cells'
[3]	380	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
[1099]	381	!
	382	DO jj = 1, jpjdta ! zdta :
[3]	383	DO ji = 1, jpidta
[3837]	384	!zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
	385	!zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
	386	zi = FLOAT( ji - ii_bump ) / r_bump
	387	zj = FLOAT( jj - ij_bump ) / r_bump
[3]	388	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
	389	END DO
	390	END DO
[1099]	391	! ! idta :
	392	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
[454]	393	idta(:,:) = jpkm1
[1099]	394	ELSE ! z-coordinate (zco or zps): step-like topography
[454]	395	idta(:,:) = jpkm1
	396	DO jk = 1, jpkm1
[2528]	397	WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) ) idta(:,:) = jk
[3]	398	END DO
[454]	399	ENDIF
[3]	400	ENDIF
[1099]	401	! ! set GLOBAL boundary conditions
	402	! ! Caution : idta on the global domain: use of jperio, not nperio
[3]	403	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
[2528]	404	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
	405	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
[3]	406	ELSEIF( jperio == 2 ) THEN
[30]	407	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
[2528]	408	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
	409	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
	410	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp
[3]	411	ELSE
[2528]	412	ih = 0 ; zh = 0._wp
[525]	413	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
[454]	414	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
	415	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
	416	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
	417	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
[3]	418	ENDIF
	419
[1099]	420	! ! local domain level and meter bathymetries (mbathy,bathy)
	421	mbathy(:,:) = 0 ! set to zero extra halo points
[2528]	422	bathy (:,:) = 0._wp ! (require for mpp case)
[3432]	423	#if defined key_mpp_rkpart
	424	DO jj = nldj, nlcj ! interior values
	425	DO ji = nldi, nlci
	426	#else
[1099]	427	DO jj = 1, nlcj ! interior values
[473]	428	DO ji = 1, nlci
[3432]	429	#endif
[473]	430	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
	431	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
	432	END DO
	433	END DO
[1099]	434	!
	435	! ! ================ !
	436	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
	437	! ! ================ !
	438	!
	439	IF( ln_zco ) THEN ! zco : read level bathymetry
[3837]	440	CALL iom_open ( 'bathy_level.nc', inum )
	441	IF(is_global)THEN
	442	CALL iom_get ( inum, jpdom_unknown, 'Bathymetry' , bathy, &
	443	kstart=(/jpizoom,jpjzoom/), &
	444	kcount=(/jpiglo,jpjglo/) )
	445	ELSE
	446	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
	447	END IF
	448
[2528]	449	CALL iom_close( inum )
[473]	450	mbathy(:,:) = INT( bathy(:,:) )
[1273]	451	! ! =====================
	452	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	453	! ! =====================
[2528]	454	IF( nn_cla == 0 ) THEN
[1273]	455	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	456	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	457	DO ji = mi0(ii0), mi1(ii1)
	458	DO jj = mj0(ij0), mj1(ij1)
	459	mbathy(ji,jj) = 15
	460	END DO
	461	END DO
	462	IF(lwp) WRITE(numout,*)
[2528]	463	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
[1273]	464	!
	465	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	466	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	467	DO ji = mi0(ii0), mi1(ii1)
	468	DO jj = mj0(ij0), mj1(ij1)
	469	mbathy(ji,jj) = 12
	470	END DO
	471	END DO
	472	IF(lwp) WRITE(numout,*)
[2528]	473	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
[1273]	474	ENDIF
	475	!
	476	ENDIF
	477	!
[454]	478	ENDIF
[1099]	479	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
[2528]	480	CALL iom_open ( 'bathy_meter.nc', inum )
[3837]	481	IF(is_global)THEN
	482	CALL iom_get ( inum, jpdom_unknown, 'Bathymetry' , bathy, &
	483	kstart=(/jpizoom,jpjzoom/), &
	484	kcount=(/jpiglo,jpjglo/) )
	485	ELSE
	486	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
	487	END IF
[2528]	488	CALL iom_close( inum )
[1348]	489	! ! =====================
[2528]	490	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	491	ii0 = 142 ; ii1 = 142 ! =====================
	492	ij0 = 51 ; ij1 = 53
	493	DO ji = mi0(ii0), mi1(ii1) ! Close Halmera Strait
	494	DO jj = mj0(ij0), mj1(ij1)
	495	bathy(ji,jj) = 0._wp
	496	END DO
	497	END DO
	498	IF(lwp) WRITE(numout,*)
	499	IF(lwp) WRITE(numout,*) ' orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
	500	ENDIF
	501	! ! =====================
	502	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	503	! ! =====================
	504	IF( nn_cla == 0 ) THEN
	505	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	506	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
[1273]	507	DO ji = mi0(ii0), mi1(ii1)
	508	DO jj = mj0(ij0), mj1(ij1)
[2528]	509	bathy(ji,jj) = 284._wp
[1273]	510	END DO
	511	END DO
	512	IF(lwp) WRITE(numout,*)
[2528]	513	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
[1273]	514	!
[2528]	515	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	516	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
[1273]	517	DO ji = mi0(ii0), mi1(ii1)
	518	DO jj = mj0(ij0), mj1(ij1)
[2528]	519	bathy(ji,jj) = 137._wp
[1273]	520	END DO
	521	END DO
	522	IF(lwp) WRITE(numout,*)
	523	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	524	ENDIF
	525	!
	526	ENDIF
[1348]	527	!
	528	ENDIF
[3]	529	! ! =============== !
	530	ELSE ! error !
	531	! ! =============== !
[1099]	532	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
[473]	533	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
[3]	534	ENDIF
[1099]	535	!
	536	! ! =========================== !
	537	IF( nclosea == 0 ) THEN ! NO closed seas or lakes !
	538	DO jl = 1, jpncs ! =========================== !
[3]	539	DO jj = ncsj1(jl), ncsj2(jl)
	540	DO ji = ncsi1(jl), ncsi2(jl)
[1099]	541	mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry
[2528]	542	bathy (ji,jj) = 0._wp
[3]	543	END DO
	544	END DO
	545	END DO
	546	ENDIF
[2528]	547	!
	548	! ! =========================== !
	549	! ! set a minimum depth !
	550	! ! =========================== !
[2712]	551	IF ( .not. ln_sco ) THEN
	552	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
	553	ELSE ; ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	554	ENDIF
	555	zhmin = gdepw_0(ik+1) ! minimum depth = ik+1 w-levels
	556	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
[3837]	557	ELSEWHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
[2712]	558	END WHERE
	559	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
[2528]	560	ENDIF
	561	!
[3]	562	END SUBROUTINE zgr_bat
	563
	564
	565	SUBROUTINE zgr_bat_zoom
	566	!!----------------------------------------------------------------------
	567	!! * ROUTINE zgr_bat_zoom *
	568	!!
	569	!! ** Purpose : - Close zoom domain boundary if necessary
	570	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
	571	!!
	572	!! ** Method :
	573	!!
	574	!! ** Action : - update mbathy: level bathymetry (in level index)
	575	!!----------------------------------------------------------------------
	576	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	577	!!----------------------------------------------------------------------
[1099]	578	!
[3]	579	IF(lwp) WRITE(numout,*)
	580	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
	581	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
[1099]	582	!
[3]	583	! Zoom domain
	584	! ===========
[1099]	585	!
[3]	586	! Forced closed boundary if required
[1099]	587	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
	588	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
	589	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
	590	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
	591	!
[3]	592	! Configuration specific domain modifications
	593	! (here, ORCA arctic configuration: suppress Med Sea)
	594	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
	595	SELECT CASE ( jp_cfg )
	596	! ! =======================
	597	CASE ( 2 ) ! ORCA_R2 configuration
	598	! ! =======================
	599	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
	600	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
	601	ij0 = 98 ; ij1 = 110
	602	! ! =======================
	603	CASE ( 05 ) ! ORCA_R05 configuration
	604	! ! =======================
	605	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
	606	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
	607	ij0 = 314 ; ij1 = 370
	608	END SELECT
	609	!
	610	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
	611	!
	612	ENDIF
[1099]	613	!
[3]	614	END SUBROUTINE zgr_bat_zoom
	615
	616
	617	SUBROUTINE zgr_bat_ctl
	618	!!----------------------------------------------------------------------
	619	!! * ROUTINE zgr_bat_ctl *
	620	!!
	621	!! ** Purpose : check the bathymetry in levels
	622	!!
	623	!! ** Method : The array mbathy is checked to verified its consistency
	624	!! with the model options. in particular:
	625	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	626	!! along closed boundary.
	627	!! mbathy must be cyclic IF jperio=1.
	628	!! mbathy must be lower or equal to jpk-1.
	629	!! isolated ocean grid points are suppressed from mbathy
	630	!! since they are only connected to remaining
	631	!! ocean through vertical diffusion.
	632	!! C A U T I O N : mbathy will be modified during the initializa-
	633	!! tion phase to become the number of non-zero w-levels of a water
	634	!! column, with a minimum value of 1.
	635	!!
	636	!! ** Action : - update mbathy: level bathymetry (in level index)
	637	!! - update bathy : meter bathymetry (in meters)
	638	!!----------------------------------------------------------------------
[3432]	639	USE mapcomm_mod, ONLY: trimmed, nidx,eidx,sidx,widx
	640	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	641	USE wrk_nemo, ONLY: zbathy => wrk_2d_1
[2715]	642	!!
[1099]	643	INTEGER :: ji, jj, jl ! dummy loop indices
	644	INTEGER :: icompt, ibtest, ikmax ! temporary integers
[3]	645	!!----------------------------------------------------------------------
	646
[2715]	647	IF( wrk_in_use(2, 1) ) THEN
	648	CALL ctl_stop('zgr_bat_ctl: requested workspace array unavailable') ; RETURN
	649	ENDIF
	650
[3]	651	IF(lwp) WRITE(numout,*)
	652	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
	653	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	654
[1099]	655	! ! Suppress isolated ocean grid points
	656	IF(lwp) WRITE(numout,*)
	657	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
	658	IF(lwp) WRITE(numout,*)' -----------------------------------'
	659	icompt = 0
	660	DO jl = 1, 2
	661	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	662	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	663	mbathy(jpi,:) = mbathy( 2 ,:)
	664	ENDIF
	665	DO jj = 2, jpjm1
	666	DO ji = 2, jpim1
	667	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
	668	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
	669	IF( ibtest < mbathy(ji,jj) ) THEN
	670	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
	671	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
	672	mbathy(ji,jj) = ibtest
	673	icompt = icompt + 1
	674	ENDIF
	675	END DO
	676	END DO
	677	END DO
	678	IF( icompt == 0 ) THEN
	679	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
	680	ELSE
	681	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
	682	ENDIF
[3432]	683
[1099]	684	IF( lk_mpp ) THEN
	685	zbathy(:,:) = FLOAT( mbathy(:,:) )
[2528]	686	CALL lbc_lnk( zbathy, 'T', 1._wp )
[1099]	687	mbathy(:,:) = INT( zbathy(:,:) )
	688	ENDIF
[3]	689
[1099]	690	! ! East-west cyclic boundary conditions
	691	IF( nperio == 0 ) THEN
	692	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
	693	IF( lk_mpp ) THEN
[3432]	694	IF( (nbondi == -1 .OR. nbondi == 2) .AND. (.NOT. trimmed(widx,narea) ) ) THEN
[1099]	695	IF( jperio /= 1 ) mbathy(1,:) = 0
[411]	696	ENDIF
[3432]	697	IF( (nbondi == 1 .OR. nbondi == 2) .AND. (.NOT. trimmed(eidx,narea) ) ) THEN
[1099]	698	IF( jperio /= 1 ) mbathy(nlci,:) = 0
	699	ENDIF
[411]	700	ELSE
[1099]	701	IF( ln_zco .OR. ln_zps ) THEN
	702	mbathy( 1 ,:) = 0
	703	mbathy(jpi,:) = 0
	704	ELSE
	705	mbathy( 1 ,:) = jpkm1
	706	mbathy(jpi,:) = jpkm1
	707	ENDIF
[411]	708	ENDIF
[1099]	709	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	710	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
	711	mbathy( 1 ,:) = mbathy(jpim1,:)
	712	mbathy(jpi,:) = mbathy( 2 ,:)
	713	ELSEIF( nperio == 2 ) THEN
	714	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
	715	ELSE
	716	IF(lwp) WRITE(numout,*) ' e r r o r'
	717	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
	718	! STOP 'dom_mba'
	719	ENDIF
	720
[1528]	721	! Boundary condition on mbathy
	722	IF( .NOT.lk_mpp ) THEN
	723	!!gm !!bug ??? think about it !
	724	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
	725	zbathy(:,:) = FLOAT( mbathy(:,:) )
[2528]	726	CALL lbc_lnk( zbathy, 'T', 1._wp )
[1528]	727	mbathy(:,:) = INT( zbathy(:,:) )
[3]	728	ENDIF
	729
	730	! Number of ocean level inferior or equal to jpkm1
	731	ikmax = 0
	732	DO jj = 1, jpj
	733	DO ji = 1, jpi
	734	ikmax = MAX( ikmax, mbathy(ji,jj) )
	735	END DO
	736	END DO
[1099]	737	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
[3]	738	IF( ikmax > jpkm1 ) THEN
	739	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
	740	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
	741	ELSE IF( ikmax < jpkm1 ) THEN
	742	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
	743	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
	744	ENDIF
	745
[1566]	746	IF( lwp .AND. nprint == 1 ) THEN ! control print
[3]	747	WRITE(numout,*)
[1099]	748	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
[3]	749	WRITE(numout,*) ' ------------------'
[1099]	750	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
[3]	751	WRITE(numout,*)
	752	ENDIF
[1099]	753	!
[2715]	754	IF( wrk_not_released(2, 1) ) CALL ctl_stop('zgr_bat_ctl: failed to release workspace array')
	755	!
[3]	756	END SUBROUTINE zgr_bat_ctl
	757
	758
[2528]	759	SUBROUTINE zgr_bot_level
	760	!!----------------------------------------------------------------------
	761	!! * ROUTINE zgr_bot_level *
	762	!!
	763	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
	764	!!
	765	!! ** Method : computes from mbathy with a minimum value of 1 over land
	766	!!
	767	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
	768	!! ocean level at t-, u- & v-points
	769	!! (min value = 1 over land)
	770	!!----------------------------------------------------------------------
[3837]	771	!USE arpdebugging, ONLY: dump_array
[2715]	772	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	773	USE wrk_nemo, ONLY: zmbk => wrk_2d_1
[4405]	774	USE par_oce, ONLY: jpkf, jpkfm1
[2715]	775	!!
[2528]	776	INTEGER :: ji, jj ! dummy loop indices
[4405]	777	LOGICAL :: domtrim_z
	778	CHARACTER(LEN=8) :: lstr ! Local string for reading env. var.
	779	INTEGER :: lztrim, ierr ! Local int for " " "
[2528]	780	!!----------------------------------------------------------------------
	781	!
[2715]	782	IF( wrk_in_use(2, 1) ) THEN
	783	CALL ctl_stop('zgr_bot_level: requested 2D workspace unavailable') ; RETURN
	784	ENDIF
	785	!
[2528]	786	IF(lwp) WRITE(numout,*)
	787	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
	788	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
	789	!
[3837]	790	!CALL dump_array(0, 'mbathy', mbathy, withHalos=.TRUE.)
	791
[2528]	792	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
	793	! ! bottom k-index of W-level = mbkt+1
	794	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
	795	DO ji = 1, jpim1
	796	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
	797	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
	798	END DO
	799	END DO
	800	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
	801	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	802	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	803	!
[4405]	804
	805	domtrim_z = .TRUE.
	806	CALL GET_ENVIRONMENT_VARIABLE(NAME='NEMO_TRIM_DOMZ', VALUE=lstr, &
	807	STATUS=ierr)
	808	IF( ierr == 0)THEN
[4409]	809	READ(lstr,FMT="(I10)",IOSTAT=ierr) lztrim
[4405]	810	IF(ierr == 0)THEN
	811	IF (lztrim == 0) domtrim_z = .FALSE.
	812	ELSE
	813	CALL ctl_warn('zgr_bot_level: failed to parse value of NEMO_TRIM_DOMZ environment variable: '//TRIM(lstr))
	814	END IF
	815	END IF
	816
	817	! Compute and store the deepest bottom level of any grid-type at each grid
	818	! point. For use in removing data below ocean floor from compute loops
[4409]	819	jpkorig = jpk
[4405]	820	IF( domtrim_z ) THEN
	821	mbkmax(:,:) = MAX(mbkt(:,:)+1, mbku(:,:), mbkv(:,:))
	822	! jpkf = MIN(jpk, 1 + MAXVAL( mbkmax(:,:) ) )
	823	jpkf = MAXVAL( mbkmax(:,:) )
	824	WRITE(,) narea,': ARPDBG: shallowest pt and jpkf = ', &
	825	MINVAL(mbkmax(:,:)), jpkf
[4409]	826	! Play rather fast and loose and just change the value of jpk
	827	! here...
	828	jpk = jpkf
	829	jpkm1 = jpk - 1
[4405]	830	ELSE
	831	WRITE(,) narea,': ARPDBG: NOT trimming domain in z'
	832	jpkf = jpk
	833	END IF
[4399]	834	jpkfm1 = jpkf - 1
	835
[3837]	836	!
[2715]	837	IF( wrk_not_released(2, 1) ) CALL ctl_stop('zgr_bot_level: failed to release workspace array')
	838	!
[2528]	839	END SUBROUTINE zgr_bot_level
	840
	841
[454]	842	SUBROUTINE zgr_zco
	843	!!----------------------------------------------------------------------
	844	!! * ROUTINE zgr_zco *
	845	!!
	846	!! ** Purpose : define the z-coordinate system
	847	!!
[2528]	848	!! ** Method : set 3D coord. arrays to reference 1D array
[454]	849	!!----------------------------------------------------------------------
[3211]	850	#if defined key_z_first
	851	INTEGER :: ji, jj ! Dummy loop indices
	852	#else
[454]	853	INTEGER :: jk
[3211]	854	#endif
[454]	855	!!----------------------------------------------------------------------
[1099]	856	!
[3211]	857	#if defined key_z_first
	858	DO jj = 1, jpj
	859	DO ji = 1, jpi
	860	fsdept(ji,jj,:) = gdept_0(:)
	861	fsdepw(ji,jj,:) = gdepw_0(:)
	862	fsde3w(ji,jj,:) = gdepw_0(:)
	863	fse3t (ji,jj,:) = e3t_0(:)
	864	fse3u (ji,jj,:) = e3t_0(:)
	865	fse3v (ji,jj,:) = e3t_0(:)
	866	fse3f (ji,jj,:) = e3t_0(:)
	867	fse3w (ji,jj,:) = e3w_0(:)
	868	fse3uw(ji,jj,:) = e3w_0(:)
	869	fse3vw(ji,jj,:) = e3w_0(:)
	870	END DO
	871	END DO
	872	#else
[2528]	873	DO jk = 1, jpk
	874	fsdept(:,:,jk) = gdept_0(jk)
	875	fsdepw(:,:,jk) = gdepw_0(jk)
	876	fsde3w(:,:,jk) = gdepw_0(jk)
	877	fse3t (:,:,jk) = e3t_0(jk)
	878	fse3u (:,:,jk) = e3t_0(jk)
	879	fse3v (:,:,jk) = e3t_0(jk)
	880	fse3f (:,:,jk) = e3t_0(jk)
	881	fse3w (:,:,jk) = e3w_0(jk)
	882	fse3uw(:,:,jk) = e3w_0(jk)
	883	fse3vw(:,:,jk) = e3w_0(jk)
	884	END DO
[3211]	885	#endif
[1099]	886	!
[454]	887	END SUBROUTINE zgr_zco
	888
	889
[3837]	890	SUBROUTINE zgr_zps(pre_domain_decomp)
[1083]	891	!!----------------------------------------------------------------------
	892	!! * ROUTINE zgr_zps *
	893	!!
	894	!! ** Purpose : the depth and vertical scale factor in partial step
	895	!! z-coordinate case
	896	!!
	897	!! ** Method : Partial steps : computes the 3D vertical scale factors
	898	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	899	!! a partial step representation of bottom topography.
	900	!!
	901	!! The reference depth of model levels is defined from an analytical
	902	!! function the derivative of which gives the reference vertical
	903	!! scale factors.
	904	!! From depth and scale factors reference, we compute there new value
	905	!! with partial steps on 3d arrays ( i, j, k ).
	906	!!
	907	!! w-level: gdepw(i,j,k) = fsdep(k)
	908	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
	909	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
	910	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
	911	!!
	912	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
	913	!! we find the mbathy index of the depth at each grid point.
	914	!! This leads us to three cases:
	915	!!
	916	!! - bathy = 0 => mbathy = 0
	917	!! - 1 < mbathy < jpkm1
	918	!! - bathy > gdepw(jpk) => mbathy = jpkm1
	919	!!
	920	!! Then, for each case, we find the new depth at t- and w- levels
	921	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
	922	!! and f-points.
	923	!!
	924	!! This routine is given as an example, it must be modified
	925	!! following the user s desiderata. nevertheless, the output as
	926	!! well as the way to compute the model levels and scale factors
	927	!! must be respected in order to insure second order accuracy
	928	!! schemes.
	929	!!
	930	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
	931	!! - - - - - - - gdept, gdepw and e3. are positives
	932	!!
[1099]	933	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
[1083]	934	!!----------------------------------------------------------------------
[2715]	935	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	936	USE wrk_nemo, ONLY: zprt => wrk_3d_1
[3211]	937	!! DCSE_NEMO: wrk_3d_1 renamed, need additional directive
	938	!FTRANS zprt :I :I :z
[3837]	939	LOGICAL, INTENT(in), OPTIONAL :: pre_domain_decomp
[2715]	940	!!
[1099]	941	INTEGER :: ji, jj, jk ! dummy loop indices
	942	INTEGER :: ik, it ! temporary integers
	943	LOGICAL :: ll_print ! Allow control print for debugging
	944	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	945	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
[2528]	946	REAL(wp) :: zmax ! Maximum depth
[1099]	947	REAL(wp) :: zdiff ! temporary scalar
[2528]	948	REAL(wp) :: zrefdep ! temporary scalar
[1099]	949	!!---------------------------------------------------------------------
[2715]	950	!
	951	IF( wrk_in_use(3, 1) ) THEN
	952	CALL ctl_stop('zgr_zps: requested workspace unavailable.') ; RETURN
	953	ENDIF
[3]	954
[1099]	955	IF(lwp) WRITE(numout,*)
	956	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
	957	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
	958	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
[3]	959
[2528]	960	ll_print = .FALSE. ! Local variable for debugging
[1083]	961
[1099]	962	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
[1083]	963	WRITE(numout,*)
	964	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
	965	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
	966	ENDIF
	967
	968
	969	! bathymetry in level (from bathy_meter)
	970	! ===================
[2528]	971	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
	972	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
	973	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
[3837]	974	ELSEWHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
[2528]	975	END WHERE
[1083]	976
	977	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	978	! find the number of ocean levels such that the last level thickness
	979	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
	980	! e3t_0 is the reference level thickness
	981	DO jk = jpkm1, 1, -1
	982	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
[2528]	983	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
[1083]	984	END DO
	985
[3837]	986	! If we've been called before domain decomposition then we only want to compute
	987	! mbathy and the return.
	988	IF( PRESENT(pre_domain_decomp) )THEN
	989	IF( pre_domain_decomp )RETURN
	990	ENDIF
	991
[1099]	992	! Scale factors and depth at T- and W-points
[3211]	993	#if defined key_z_first
	994	DO jj = 1, jpj
	995	DO ji = 1, jpi ! intitialization to the reference z-coordinate
	996	gdept(ji,jj,:) = gdept_0(:)
	997	gdepw(ji,jj,:) = gdepw_0(:)
	998	e3t (ji,jj,:) = e3t_0 (:)
	999	e3w (ji,jj,:) = e3w_0 (:)
	1000	END DO
	1001	END DO
	1002	#else
[1099]	1003	DO jk = 1, jpk ! intitialization to the reference z-coordinate
	1004	gdept(:,:,jk) = gdept_0(jk)
	1005	gdepw(:,:,jk) = gdepw_0(jk)
	1006	e3t (:,:,jk) = e3t_0 (jk)
	1007	e3w (:,:,jk) = e3w_0 (jk)
[1083]	1008	END DO
[3211]	1009	#endif
[1099]	1010	!
	1011	DO jj = 1, jpj
	1012	DO ji = 1, jpi
	1013	ik = mbathy(ji,jj)
	1014	IF( ik > 0 ) THEN ! ocean point only
	1015	! max ocean level case
	1016	IF( ik == jpkm1 ) THEN
	1017	zdepwp = bathy(ji,jj)
	1018	ze3tp = bathy(ji,jj) - gdepw_0(ik)
[2528]	1019	ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) )
[1099]	1020	e3t(ji,jj,ik ) = ze3tp
	1021	e3t(ji,jj,ik+1) = ze3tp
	1022	e3w(ji,jj,ik ) = ze3wp
	1023	e3w(ji,jj,ik+1) = ze3tp
	1024	gdepw(ji,jj,ik+1) = zdepwp
	1025	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
	1026	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
	1027	!
	1028	ELSE ! standard case
	1029	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
	1030	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
	1031	ENDIF
	1032	!gm Bug? check the gdepw_0
	1033	! ... on ik
	1034	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	1035	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
	1036	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
	1037	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	1038	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
[2528]	1039	e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) &
	1040	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
[1099]	1041	! ... on ik+1
	1042	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
	1043	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
	1044	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
	1045	ENDIF
	1046	ENDIF
	1047	END DO
	1048	END DO
	1049	!
	1050	it = 0
	1051	DO jj = 1, jpj
	1052	DO ji = 1, jpi
	1053	ik = mbathy(ji,jj)
	1054	IF( ik > 0 ) THEN ! ocean point only
	1055	e3tp (ji,jj) = e3t(ji,jj,ik )
	1056	e3wp (ji,jj) = e3w(ji,jj,ik )
	1057	! test
	1058	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
[2528]	1059	IF( zdiff <= 0._wp .AND. lwp ) THEN
[1099]	1060	it = it + 1
	1061	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	1062	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	1063	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
	1064	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
	1065	ENDIF
	1066	ENDIF
	1067	END DO
	1068	END DO
[1083]	1069
	1070	! Scale factors and depth at U-, V-, UW and VW-points
[3211]	1071	#if defined key_z_first
	1072	DO jj = 1, jpj ! initialisation to z-scale factors
	1073	DO ji = 1, jpi
	1074	e3u (ji,jj,:) = e3t_0(:)
	1075	e3v (ji,jj,:) = e3t_0(:)
	1076	e3uw(ji,jj,:) = e3w_0(:)
	1077	e3vw(ji,jj,:) = e3w_0(:)
[3226]	1078	END DO
[3211]	1079	END DO
	1080	#else
[1099]	1081	DO jk = 1, jpk ! initialisation to z-scale factors
[2528]	1082	e3u (:,:,jk) = e3t_0(jk)
	1083	e3v (:,:,jk) = e3t_0(jk)
	1084	e3uw(:,:,jk) = e3w_0(jk)
	1085	e3vw(:,:,jk) = e3w_0(jk)
[1083]	1086	END DO
[3211]	1087	#endif
	1088	#if defined key_z_first
	1089	DO jj = 1, jpjm1
	1090	DO ji = 1, jpim1
	1091	DO jk = 1, jpk ! Computed as the minimum of neighbouring scale factors
	1092	#else
	1093	DO jk = 1,jpk ! Computed as the minimum of neighbouring scale factors
[1099]	1094	DO jj = 1, jpjm1
	1095	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	1096	#endif
[2528]	1097	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
	1098	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
[1099]	1099	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
	1100	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
	1101	END DO
	1102	END DO
	1103	END DO
[2528]	1104	CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions
	1105	CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp )
[1099]	1106	!
	1107	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
[2528]	1108	WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk)
	1109	WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk)
	1110	WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk)
	1111	WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk)
[1099]	1112	END DO
	1113
	1114	! Scale factor at F-point
[3211]	1115	#if defined key_z_first
	1116	DO jj = 1, jpj
	1117	DO ji = 1, jpi ! initialisation to z-scale factors
	1118	e3f(ji,jj,:) = e3t_0(:)
	1119	END DO
	1120	END DO
	1121	DO jj = 1, jpjm1
	1122	DO ji = 1, jpim1 ! NO vector opt.
	1123	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	1124	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
	1125	END DO
	1126	END DO
	1127	END DO
	1128	#else
[1099]	1129	DO jk = 1, jpk ! initialisation to z-scale factors
[2528]	1130	e3f(:,:,jk) = e3t_0(jk)
[1099]	1131	END DO
	1132	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	1133	DO jj = 1, jpjm1
	1134	DO ji = 1, fs_jpim1 ! vector opt.
	1135	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
	1136	END DO
	1137	END DO
	1138	END DO
[3211]	1139	#endif
[2528]	1140	CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions
[1099]	1141	!
	1142	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
[2528]	1143	WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk)
[1099]	1144	END DO
	1145	!!gm bug ? : must be a do loop with mj0,mj1
	1146	!
	1147	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
	1148	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
	1149	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
	1150	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
	1151	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
[1083]	1152
[1099]	1153	! Control of the sign
[2528]	1154	IF( MINVAL( e3t (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
	1155	IF( MINVAL( e3w (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
	1156	IF( MINVAL( gdept(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
	1157	IF( MINVAL( gdepw(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
[1083]	1158
[1099]	1159	! Compute gdep3w (vertical sum of e3w)
[2528]	1160	gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
[1099]	1161	DO jk = 2, jpk
	1162	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
	1163	END DO
[2528]	1164
[1099]	1165	! ! ================= !
	1166	IF(lwp .AND. ll_print) THEN ! Control print !
	1167	! ! ================= !
[1083]	1168	DO jj = 1,jpj
	1169	DO ji = 1, jpi
[1099]	1170	ik = MAX( mbathy(ji,jj), 1 )
	1171	zprt(ji,jj,1) = e3t (ji,jj,ik)
	1172	zprt(ji,jj,2) = e3w (ji,jj,ik)
	1173	zprt(ji,jj,3) = e3u (ji,jj,ik)
	1174	zprt(ji,jj,4) = e3v (ji,jj,ik)
	1175	zprt(ji,jj,5) = e3f (ji,jj,ik)
	1176	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
[1083]	1177	END DO
	1178	END DO
	1179	WRITE(numout,*)
[1099]	1180	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1181	WRITE(numout,*)
[1099]	1182	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1183	WRITE(numout,*)
[1099]	1184	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1185	WRITE(numout,*)
[1099]	1186	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1187	WRITE(numout,*)
[1099]	1188	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1189	WRITE(numout,*)
[1099]	1190	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1191	ENDIF
[2528]	1192	!
[2715]	1193	IF( wrk_not_released(3, 1) ) CALL ctl_stop('zgr_zps: failed to release workspace')
	1194	!
[1083]	1195	END SUBROUTINE zgr_zps
	1196
	1197
	1198	FUNCTION fssig( pk ) RESULT( pf )
	1199	!!----------------------------------------------------------------------
	1200	!! * ROUTINE eos_init *
	1201	!!
	1202	!! ** Purpose : provide the analytical function in s-coordinate
	1203	!!
	1204	!! ** Method : the function provide the non-dimensional position of
	1205	!! T and W (i.e. between 0 and 1)
	1206	!! T-points at integer values (between 1 and jpk)
	1207	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1208	!!----------------------------------------------------------------------
[2528]	1209	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
	1210	REAL(wp) :: pf ! sigma value
[1083]	1211	!!----------------------------------------------------------------------
	1212	!
[2528]	1213	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
[1601]	1214	& - TANH( rn_thetb * rn_theta ) ) &
[2528]	1215	& * ( COSH( rn_theta ) &
	1216	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
	1217	& / ( 2._wp * SINH( rn_theta ) )
[1083]	1218	!
	1219	END FUNCTION fssig
	1220
	1221
[1601]	1222	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
[1348]	1223	!!----------------------------------------------------------------------
	1224	!! * ROUTINE eos_init *
	1225	!!
	1226	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
	1227	!!
	1228	!! ** Method : the function provides the non-dimensional position of
	1229	!! T and W (i.e. between 0 and 1)
	1230	!! T-points at integer values (between 1 and jpk)
	1231	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1232	!!----------------------------------------------------------------------
[2528]	1233	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
	1234	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
	1235	REAL(wp) :: pf1 ! sigma value
[1348]	1236	!!----------------------------------------------------------------------
	1237	!
[1601]	1238	IF ( rn_theta == 0 ) then ! uniform sigma
[2528]	1239	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
[1566]	1240	ELSE ! stretched sigma
[2528]	1241	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
	1242	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
	1243	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
[1348]	1244	ENDIF
[1566]	1245	!
[1348]	1246	END FUNCTION fssig1
	1247
	1248
[454]	1249	SUBROUTINE zgr_sco
	1250	!!----------------------------------------------------------------------
	1251	!! * ROUTINE zgr_sco *
	1252	!!
	1253	!! ** Purpose : define the s-coordinate system
	1254	!!
	1255	!! ** Method : s-coordinate
	1256	!! The depth of model levels is defined as the product of an
	1257	!! analytical function by the local bathymetry, while the vertical
	1258	!! scale factors are defined as the product of the first derivative
	1259	!! of the analytical function by the bathymetry.
	1260	!! (this solution save memory as depth and scale factors are not
	1261	!! 3d fields)
	1262	!! - Read bathymetry (in meters) at t-point and compute the
	1263	!! bathymetry at u-, v-, and f-points.
	1264	!! hbatu = mi( hbatt )
	1265	!! hbatv = mj( hbatt )
	1266	!! hbatf = mi( mj( hbatt ) )
	1267	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
[1083]	1268	!! function and its derivative given as function.
	1269	!! gsigt(k) = fssig (k )
	1270	!! gsigw(k) = fssig (k-0.5)
	1271	!! esigt(k) = fsdsig(k )
	1272	!! esigw(k) = fsdsig(k-0.5)
[454]	1273	!! This routine is given as an example, it must be modified
	1274	!! following the user s desiderata. nevertheless, the output as
	1275	!! well as the way to compute the model levels and scale factors
	1276	!! must be respected in order to insure second order a!!uracy
	1277	!! schemes.
	1278	!!
[1099]	1279	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
	1280	!!----------------------------------------------------------------------
[2715]	1281	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	1282	USE wrk_nemo, ONLY: zenv => wrk_2d_1 , ztmp => wrk_2d_2 , zmsk => wrk_2d_3
[3432]	1283	USE wrk_nemo, ONLY: ztmp2 => wrk_2d_4 , zhbat => wrk_2d_5
[2715]	1284	USE wrk_nemo, ONLY: gsigw3 => wrk_3d_1
	1285	USE wrk_nemo, ONLY: gsigt3 => wrk_3d_2
	1286	USE wrk_nemo, ONLY: gsi3w3 => wrk_3d_3
	1287	USE wrk_nemo, ONLY: esigt3 => wrk_3d_4
	1288	USE wrk_nemo, ONLY: esigw3 => wrk_3d_5
	1289	USE wrk_nemo, ONLY: esigtu3 => wrk_3d_6
	1290	USE wrk_nemo, ONLY: esigtv3 => wrk_3d_7
	1291	USE wrk_nemo, ONLY: esigtf3 => wrk_3d_8
	1292	USE wrk_nemo, ONLY: esigwu3 => wrk_3d_9
	1293	USE wrk_nemo, ONLY: esigwv3 => wrk_3d_10
[3432]	1294	USE mapcomm_mod, ONLY: trimmed, cyclic_bc
	1295	USE mapcomm_mod, ONLY: nidx, eidx, sidx, widx
[3837]	1296
[3211]	1297	!! DCSE_NEMO: wrk_nemo module variables renamed, need additional directives
	1298	!FTRANS gsigw3 :I :I :z
	1299	!FTRANS gsigt3 :I :I :z
	1300	!FTRANS gsi3w3 :I :I :z
	1301	!FTRANS esigt3 :I :I :z
	1302	!FTRANS esigw3 :I :I :z
	1303	!FTRANS esigtu3 :I :I :z
	1304	!FTRANS esigtv3 :I :I :z
	1305	!FTRANS esigtf3 :I :I :z
	1306	!FTRANS esigwu3 :I :I :z
	1307	!FTRANS esigwv3 :I :I :z
[2715]	1308	!
[1099]	1309	INTEGER :: ji, jj, jk, jl ! dummy loop argument
	1310	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
[3432]	1311	REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper, zri, zrj ! temporary scalars
	1312	REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than this assumed zero
[2715]	1313	!
	1314
[3837]	1315	! NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
[454]	1316	!!----------------------------------------------------------------------
	1317
[3432]	1318	IF( wrk_in_use(2, 1,2,3,4,5) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN
[2715]	1319	CALL ctl_stop('zgr_sco: ERROR - requested workspace arrays unavailable') ; RETURN
	1320	ENDIF
	1321
	1322	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
[1601]	1323	READ ( numnam, namzgr_sco )
[454]	1324
[2715]	1325	IF(lwp) THEN ! control print
[454]	1326	WRITE(numout,*)
	1327	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
	1328	WRITE(numout,*) '~~~~~~~~~~~'
[1601]	1329	WRITE(numout,*) ' Namelist namzgr_sco'
	1330	WRITE(numout,*) ' sigma-stretching coeffs '
	1331	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max
	1332	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min
	1333	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta
	1334	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb
	1335	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax
	1336	WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma
	1337	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb
	1338	WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc
[454]	1339	ENDIF
	1340
[2528]	1341	gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp
	1342	esigt3 = 0._wp ; esigw3 = 0._wp
	1343	esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp
	1344	esigwu3 = 0._wp ; esigwv3 = 0._wp
[1629]	1345
[1601]	1346	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
	1347	hifu(:,:) = rn_sbot_min
	1348	hifv(:,:) = rn_sbot_min
	1349	hiff(:,:) = rn_sbot_min
[1348]	1350
	1351	! ! set maximum ocean depth
[1601]	1352	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
[454]	1353
[1461]	1354	DO jj = 1, jpj
	1355	DO ji = 1, jpi
[2715]	1356	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
[1461]	1357	END DO
	1358	END DO
[1099]	1359	! ! =============================
	1360	! ! Define the envelop bathymetry (hbatt)
	1361	! ! =============================
[454]	1362	! use r-value to create hybrid coordinates
	1363	DO jj = 1, jpj
	1364	DO ji = 1, jpi
[1601]	1365	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
[454]	1366	END DO
	1367	END DO
[3432]	1368
	1369	CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. )
[1639]	1370	!
	1371	! Smooth the bathymetry (if required)
[2528]	1372	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
[1639]	1373	scobot(:,:) = bathy(:,:) ! ocean bottom depth
	1374	!
[454]	1375	jl = 0
[2528]	1376	zrmax = 1._wp
[454]	1377	! ! ================ !
[2528]	1378	DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
[454]	1379	! ! ================ !
	1380	jl = jl + 1
[2528]	1381	zrmax = 0._wp
	1382	zmsk(:,:) = 0._wp
[3432]	1383
[454]	1384	DO jj = 1, nlcj
	1385	DO ji = 1, nlci
	1386	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
[3432]	1387	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last row (jj=nclj+1 to jpj)
	1388	zri = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
	1389	zrj = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
	1390	zrmax = MAX( zrmax, zri, zrj )
	1391	IF( zri > rn_rmax ) zmsk(ji ,jj ) = 1._wp
	1392	IF( zri > rn_rmax ) zmsk(iip1,jj ) = 1._wp
	1393	IF( zrj > rn_rmax ) zmsk(ji ,jj ) = 1._wp
	1394	IF( zrj > rn_rmax ) zmsk(ji ,ijp1) = 1._wp
[454]	1395	END DO
	1396	END DO
[3432]	1397
	1398	! lateral boundary condition on zmsk: retain any 1's along closed
	1399	! boundary (use of lzero flag to lbc_lnk)
	1400	CALL lbc_lnk( zmsk, 'T', 1._wp, lzero=.FALSE. )
	1401
[454]	1402	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
[1348]	1403	!
[3432]	1404	IF(lwp)WRITE(numout,"('zgr_sco : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") &
	1405	jl, zrmax, INT( SUM(zmsk(:,:) ) )
[1099]	1406	!
[3432]	1407
	1408	! Copy current surface before next smoothing iteration
	1409	ztmp(:,:) = zenv(:,:)
	1410
	1411	DO jj = nldj, nlcj
	1412	DO ji = nldi, nlci
[454]	1413	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
	1414	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
	1415	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
	1416	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
[2528]	1417	IF( zmsk(ji,jj) == 1._wp ) THEN
[454]	1418	ztmp(ji,jj) = ( &
	1419	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
[2528]	1420	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) &
[454]	1421	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
	1422	& ) / ( &
	1423	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
[2528]	1424	& + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) &
[454]	1425	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
	1426	& )
	1427	ENDIF
	1428	END DO
	1429	END DO
[1348]	1430	!
[3432]	1431
	1432	! Need to update halos of ztmp here but do not zero halos on closed
	1433	! boundaries
	1434	CALL lbc_lnk( ztmp, 'T', 1._wp, lzero=.FALSE.)
	1435
	1436	DO jj = 1,nlcj
	1437	DO ji = 1,nlci
	1438	IF( zmsk(ji,jj) >= 1._wp-TOL_ZERO ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
[454]	1439	END DO
	1440	END DO
[1099]	1441	!
[3432]	1442	! Apply lateral boundary condition but do not zero on closed boundaries
	1443	CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. )
	1444
[454]	1445	! ! ================ !
	1446	END DO ! End loop !
	1447	! ! ================ !
[1099]	1448	!
	1449	! ! envelop bathymetry saved in hbatt
[454]	1450	hbatt(:,:) = zenv(:,:)
[2528]	1451	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
[1099]	1452	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
	1453	DO jj = 1, jpj
	1454	DO ji = 1, jpi
[2528]	1455	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
	1456	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
[1099]	1457	END DO
	1458	END DO
[516]	1459	ENDIF
[1099]	1460	!
	1461	IF(lwp) THEN ! Control print
[454]	1462	WRITE(numout,*)
	1463	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
	1464	WRITE(numout,*)
[2528]	1465	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
[1099]	1466	IF( nprint == 1 ) THEN
	1467	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
	1468	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
	1469	ENDIF
[454]	1470	ENDIF
	1471
[3432]	1472	! CALL dump_array(0, 'hbatt', hbatt, withHalos=.FALSE.)
	1473
[1099]	1474	! ! ==============================
	1475	! ! hbatu, hbatv, hbatf fields
	1476	! ! ==============================
[454]	1477	IF(lwp) THEN
	1478	WRITE(numout,*)
[1601]	1479	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
[454]	1480	ENDIF
[1601]	1481	hbatu(:,:) = rn_sbot_min
	1482	hbatv(:,:) = rn_sbot_min
	1483	hbatf(:,:) = rn_sbot_min
[454]	1484	DO jj = 1, jpjm1
[1694]	1485	DO ji = 1, jpim1 ! NO vector opt.
[2528]	1486	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
	1487	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
	1488	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
	1489	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
[454]	1490	END DO
	1491	END DO
[1099]	1492	!
[454]	1493	! Apply lateral boundary condition
[1099]	1494	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
[2528]	1495	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
[454]	1496	DO jj = 1, jpj
	1497	DO ji = 1, jpi
[2528]	1498	IF( hbatu(ji,jj) == 0._wp ) THEN
	1499	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
	1500	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
[454]	1501	ENDIF
	1502	END DO
	1503	END DO
[2528]	1504	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
[454]	1505	DO jj = 1, jpj
	1506	DO ji = 1, jpi
[2528]	1507	IF( hbatv(ji,jj) == 0._wp ) THEN
	1508	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
	1509	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
[454]	1510	ENDIF
	1511	END DO
	1512	END DO
[2528]	1513	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
[454]	1514	DO jj = 1, jpj
	1515	DO ji = 1, jpi
[2528]	1516	IF( hbatf(ji,jj) == 0._wp ) THEN
	1517	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
	1518	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
[454]	1519	ENDIF
	1520	END DO
	1521	END DO
	1522
	1523	!!bug: key_helsinki a verifer
	1524	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
	1525	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
	1526	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
	1527	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
	1528
[516]	1529	IF( nprint == 1 .AND. lwp ) THEN
[1099]	1530	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
	1531	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
	1532	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
	1533	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
[516]	1534	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
	1535	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
	1536	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
	1537	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
	1538	ENDIF
[454]	1539	!! helsinki
	1540
[1099]	1541	! ! =======================
	1542	! ! s-ccordinate fields (gdep., e3.)
	1543	! ! =======================
	1544	!
	1545	! non-dimensional "sigma" for model level depth at w- and t-levels
[1348]	1546
[1601]	1547	IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths
	1548	! ! below rn_hc, with uniform sigma in shallower waters
	1549	DO ji = 1, jpi
	1550	DO jj = 1, jpj
[1348]	1551
[2528]	1552	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
	1553	DO jk = 1, jpk
	1554	gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
	1555	gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
	1556	END DO
	1557	ELSE ! shallow water, uniform sigma
	1558	DO jk = 1, jpk
	1559	gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
	1560	gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
	1561	END DO
	1562	ENDIF
	1563	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
	1564	!
	1565	DO jk = 1, jpkm1
	1566	esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
	1567	esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
	1568	END DO
	1569	esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) )
	1570	esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
	1571	!
	1572	! Coefficients for vertical depth as the sum of e3w scale factors
	1573	gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
	1574	DO jk = 2, jpk
	1575	gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
	1576	END DO
	1577	!
[1348]	1578	DO jk = 1, jpk
[2528]	1579	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	1580	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
	1581	gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)gsigt3(ji,jj,jk)+rn_hczcoeft )
	1582	gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)gsigw3(ji,jj,jk)+rn_hczcoefw )
	1583	gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)gsi3w3(ji,jj,jk)+rn_hczcoeft )
[1348]	1584	END DO
[2528]	1585	!
	1586	END DO ! for all jj's
	1587	END DO ! for all ji's
	1588
[2715]	1589	DO ji = 1, jpim1
	1590	DO jj = 1, jpjm1
[1348]	1591	DO jk = 1, jpk
[2528]	1592	esigtu3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) ) &
	1593	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1594	esigtv3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji,jj+1)esigt3(ji,jj+1,jk) ) &
	1595	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1596	esigtf3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) &
	1597	& + hbatt(ji,jj+1)esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)esigt3(ji+1,jj+1,jk) ) &
	1598	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
	1599	esigwu3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji+1,jj)esigw3(ji+1,jj,jk) ) &
	1600	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1601	esigwv3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji,jj+1)esigw3(ji,jj+1,jk) ) &
	1602	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1603	!
	1604	e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
	1605	e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
	1606	e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
	1607	e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
	1608	!
	1609	e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
	1610	e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
	1611	e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
[1348]	1612	END DO
[2528]	1613	END DO
	1614	END DO
[2715]	1615
	1616	CALL lbc_lnk( e3t , 'T', 1._wp )
	1617	CALL lbc_lnk( e3u , 'U', 1._wp )
	1618	CALL lbc_lnk( e3v , 'V', 1._wp )
	1619	CALL lbc_lnk( e3f , 'F', 1._wp )
	1620	CALL lbc_lnk( e3w , 'W', 1._wp )
	1621	CALL lbc_lnk( e3uw, 'U', 1._wp )
	1622	CALL lbc_lnk( e3vw, 'V', 1._wp )
	1623
[2528]	1624	!
[1348]	1625	ELSE ! not ln_s_sigma
[2528]	1626	!
[1348]	1627	DO jk = 1, jpk
	1628	gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
[2528]	1629	gsigt(jk) = -fssig( REAL(jk,wp) )
[1348]	1630	END DO
[2528]	1631	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
[1348]	1632	!
[2528]	1633	! Coefficients for vertical scale factors at w-, t- levels
[1099]	1634	!!gm bug : define it from analytical function, not like juste bellow....
	1635	!!gm or betteroffer the 2 possibilities....
[1348]	1636	DO jk = 1, jpkm1
	1637	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
	1638	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
	1639	END DO
[2528]	1640	esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) )
	1641	esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
[454]	1642
[1099]	1643	!!gm original form
	1644	!!org DO jk = 1, jpk
	1645	!!org esigt(jk)=fsdsig( FLOAT(jk) )
	1646	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
	1647	!!org END DO
	1648	!!gm
[1348]	1649	!
	1650	! Coefficients for vertical depth as the sum of e3w scale factors
[2528]	1651	gsi3w(1) = 0.5_wp * esigw(1)
[1348]	1652	DO jk = 2, jpk
	1653	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
	1654	END DO
[1099]	1655	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
[1348]	1656	DO jk = 1, jpk
[2528]	1657	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	1658	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
	1659	gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))gsigt(jk) + hift(:,:)zcoeft )
	1660	gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))gsigw(jk) + hift(:,:)zcoefw )
	1661	gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))gsi3w(jk) + hift(:,:)zcoeft )
[1348]	1662	END DO
[1099]	1663	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
[1348]	1664	DO jj = 1, jpj
	1665	DO ji = 1, jpi
	1666	DO jk = 1, jpk
[2528]	1667	e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	1668	e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	1669	e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
	1670	e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
[1348]	1671	!
[2528]	1672	e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	1673	e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	1674	e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
[1348]	1675	END DO
[454]	1676	END DO
	1677	END DO
[2528]	1678	!
[1348]	1679	ENDIF ! ln_s_sigma
	1680
	1681
[1099]	1682	!
[1461]	1683	!! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
[3432]	1684	where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0
	1685	where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.0
	1686	where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.0
	1687	where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.0
	1688	where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.0
	1689	where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.0
	1690	where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.0
[1461]	1691
[3432]	1692
[1461]	1693	fsdept(:,:,:) = gdept (:,:,:)
	1694	fsdepw(:,:,:) = gdepw (:,:,:)
	1695	fsde3w(:,:,:) = gdep3w(:,:,:)
	1696	fse3t (:,:,:) = e3t (:,:,:)
	1697	fse3u (:,:,:) = e3u (:,:,:)
	1698	fse3v (:,:,:) = e3v (:,:,:)
	1699	fse3f (:,:,:) = e3f (:,:,:)
	1700	fse3w (:,:,:) = e3w (:,:,:)
	1701	fse3uw(:,:,:) = e3uw (:,:,:)
	1702	fse3vw(:,:,:) = e3vw (:,:,:)
	1703	!!
[1099]	1704	! HYBRID :
[454]	1705	DO jj = 1, jpj
	1706	DO ji = 1, jpi
	1707	DO jk = 1, jpkm1
	1708	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
[2528]	1709	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
[454]	1710	END DO
	1711	END DO
	1712	END DO
[3432]	1713
[1099]	1714	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
	1715	& ' MAX ', MAXVAL( mbathy(:,:) )
[454]	1716
[1632]	1717	! ! =============
	1718	IF(lwp) THEN ! Control print
	1719	! ! =============
[454]	1720	WRITE(numout,*)
	1721	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
[1099]	1722	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
	1723	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
[454]	1724	ENDIF
[1099]	1725	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
	1726	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	1727	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	1728	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
	1729	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
	1730	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
	1731	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
	1732	& ' w ', MINVAL( fse3w (:,:,:) )
[454]	1733
[1099]	1734	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	1735	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
	1736	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
	1737	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
	1738	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
	1739	& ' w ', MAXVAL( fse3w (:,:,:) )
	1740	ENDIF
	1741	!
	1742	IF(lwp) THEN ! selected vertical profiles
[454]	1743	WRITE(numout,*)
	1744	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
	1745	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1746	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1747	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
	1748	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
[473]	1749	DO jj = mj0(20), mj1(20)
	1750	DO ji = mi0(20), mi1(20)
	1751	WRITE(numout,*)
	1752	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1753	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1754	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1755	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1756	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1757	END DO
	1758	END DO
[3432]	1759	!!$ ARPDBG - out of bounds if jpj < 74 or jpi < 100, e.g. default GYRE
	1760	!!$ DO jj = mj0(74), mj1(74)
	1761	!!$ DO ji = mi0(100), mi1(100)
	1762	!!$ WRITE(numout,*)
	1763	!!$ WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1764	!!$ WRITE(numout,*) ' ~~~~~~ --------------------'
	1765	!!$ WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1766	!!$ WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1767	!!$ & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
	1768	!!$ END DO
	1769	!!$ END DO
[454]	1770	ENDIF
	1771
[1099]	1772	!!gm bug? no more necessary? if ! defined key_helsinki
[3211]	1773	#if defined key_z_first
	1774	DO jj = 1, jpj
	1775	DO ji = 1, jpi
	1776	DO jk = 1, jpk
	1777	#else
[454]	1778	DO jk = 1, jpk
	1779	DO jj = 1, jpj
	1780	DO ji = 1, jpi
[3211]	1781	#endif
[2528]	1782	IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
[1099]	1783	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	1784	CALL ctl_stop( ctmp1 )
[454]	1785	ENDIF
[3837]	1786	#if defined key_vvl
[3432]	1787	IF( gdepw_1(ji,jj,jk) < 0._wp .OR. gdept_1(ji,jj,jk) < 0._wp ) THEN
[1099]	1788	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	1789	CALL ctl_stop( ctmp1 )
[454]	1790	ENDIF
[3837]	1791	#endif
[454]	1792	END DO
	1793	END DO
	1794	END DO
[1099]	1795	!!gm bug #endif
	1796	!
[3432]	1797	IF( wrk_not_released(2, 1,2,3,4,5) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) ) &
[2715]	1798	& CALL ctl_stop('dom:zgr_sco: failed to release workspace arrays')
	1799	!
[454]	1800	END SUBROUTINE zgr_sco
	1801
[3]	1802	!!======================================================================
	1803	END MODULE domzgr

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

Download in other formats: