Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

vgrid.f90

Last change on this file was 13369, checked in by jpaul, 4 years ago
update: cf changelog inside documentation
File size: 32.3 KB

Rev	Line
[4213]	1	!----------------------------------------------------------------------
	2	! NEMO system team, System and Interface for oceanic RElocable Nesting
	3	!----------------------------------------------------------------------
	4	!
	5	! DESCRIPTION:
[5037]	6	!> @brief This module manage vertical grid.
[4213]	7	!>
	8	!> @details
[5037]	9	!> to set the depth of model levels and the resulting vertical scale
	10	!> factors:<br/>
	11	!> @code
[13369]	12	!> CALL vgrid_zgr_z(dd_gdepw(:), dd_gdept(:), dd_e3w(:), dd_e3t(:),
	13	!> dd_ppkth, dd_ppkth2, dd_ppacr, dd_ppacr2,
	14	!> dd_ppdzmin, dd_pphmax, dd_pp_to_be_computed,
[5037]	15	!> dd_ppa0, dd_ppa1, dd_ppa2, dd_ppsur)
	16	!> @endcode
	17	!> - dd_gdepw is array of depth value on W point
	18	!> - dd_gdept is array of depth value on T point
	19	!> - dd_e3w is array of vertical mesh size on W point
	20	!> - dd_e3t is array of vertical mesh size on T point
	21	!> - dd_ppkth see NEMO documentation
	22	!> - dd_ppkth2 see NEMO documentation
	23	!> - dd_ppacr see NEMO documentation
	24	!> - dd_ppdzmin see NEMO documentation
	25	!> - dd_pphmax see NEMO documentation
[12080]	26	!> - dd_pp_to_be_computed see NEMO documentation
[5037]	27	!> - dd_ppa1 see NEMO documentation
	28	!> - dd_ppa2 see NEMO documentation
	29	!> - dd_ppa0 see NEMO documentation
	30	!> - dd_ppsur see NEMO documentation
[13369]	31	!>
	32	!>
[5037]	33	!> to set the depth and vertical scale factor in partial step z-coordinate
	34	!> case:<br/>
	35	!> @code
[13369]	36	!> CALL vgrid_zgr_zps(id_mbathy(:,:), dd_bathy(:,:), id_jpkmax, dd_gdepw(:),
[5037]	37	!> dd_e3t(:), dd_e3zps_min, dd_e3zps_rat)
	38	!> @endcode
	39	!> - id_mbathy is array of bathymetry level
	40	!> - dd_bathy is array of bathymetry
	41	!> - id_jpkmax is the maximum number of level to be used
	42	!> - dd_gdepw is array of vertical mesh size on W point
	43	!> - dd_e3t is array of vertical mesh size on T point
	44	!> - dd_e3zps_min see NEMO documentation
	45	!> - dd_e3zps_rat see NEMO documentation
	46	!>
	47	!> to check the bathymetry in levels:<br/>
	48	!> @code
	49	!> CALL vgrid_zgr_bat_ctl(id_mbathy, id_jpkmax, id_jpk)
	50	!> @endcode
	51	!> - id_mbathy is array of bathymetry level
	52	!> - id_jpkmax is the maximum number of level to be used
	53	!> - id_jpk is the number of level
[13369]	54	!>
[5037]	55	!> to compute bathy level in T,U,V,F point from Bathymetry file:<br/>
	56	!> @code
	57	!> tl_level(:)=vgrid_get_level(td_bathy, [cd_namelist,] [td_dom,] [id_nlevel])
	58	!> @endcode
	59	!> - td_bathy is Bathymetry file structure
	60	!> - cd_namelist is namelist [optional]
	61	!> - td_dom is domain structure [optional]
	62	!> - id_nlevel is number of lelvel to be used [optional]
[13369]	63	!>
[4213]	64	!> @author
	65	!> J.Paul
[12080]	66	!>
[5037]	67	!> @date November, 2013 - Initial Version
	68	!> @date Spetember, 2014
	69	!> - add header
[5609]	70	!> @date June, 2015 - update subroutine with NEMO 3.6
[5037]	71	!>
[12080]	72	!> @todo
	73	!> - fusionner vgrid et grid_zgr
	74	!>
	75	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[4213]	76	!----------------------------------------------------------------------
	77	MODULE vgrid
[12080]	78
[5037]	79	USE netcdf ! nf90 library
[4213]	80	USE kind ! F90 kind parameter
	81	USE fct ! basic usefull function
	82	USE global ! global parameter
	83	USE phycst ! physical constant
[5037]	84	USE logger ! log file manager
[4213]	85	USE file ! file manager
	86	USE var ! variable manager
	87	USE dim ! dimension manager
	88	USE dom ! domain manager
[5037]	89	USE grid ! grid manager
[4213]	90	USE iom ! I/O manager
	91	USE mpp ! MPP manager
	92	USE iom_mpp ! I/O MPP manager
	93	IMPLICIT NONE
	94	! NOTE_avoid_public_variables_if_possible
	95
	96	! type and variable
	97
	98	! function and subroutine
[13369]	99	PUBLIC :: vgrid_zgr_z
[4213]	100	PUBLIC :: vgrid_zgr_zps
	101	PUBLIC :: vgrid_zgr_bat_ctl
	102	PUBLIC :: vgrid_get_level
	103
	104	CONTAINS
[12080]	105	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	106	SUBROUTINE vgrid_zgr_z(dd_gdepw, dd_gdept, dd_e3w, dd_e3t, &
	107	& dd_e3w_1d, dd_e3t_1d, &
	108	& dd_ppkth, dd_ppkth2, dd_ppacr, dd_ppacr2, &
	109	& dd_ppdzmin, dd_pphmax, dd_pp_to_be_computed, &
	110	& dd_ppa0, dd_ppa1, dd_ppa2, dd_ppsur )
[4213]	111	!-------------------------------------------------------------------
[13369]	112	!> @brief This subroutine set the depth of model levels and the resulting
[4213]	113	!> vertical scale factors.
[12080]	114	!>
[4213]	115	!> @details
	116	!> ** Method : z-coordinate system (use in all type of coordinate)
	117	!> The depth of model levels is defined from an analytical
	118	!> function the derivative of which gives the scale factors.
[12080]	119	!> both depth and scale factors only depend on k (1d arrays). <br/>
	120	!> w-level: gdepw = fsdep(k) <br/>
	121	!> e3w(k) = dk(fsdep)(k) = fse3(k) <br/>
	122	!> t-level: gdept = fsdep(k+0.5) <br/>
	123	!> e3t(k) = dk(fsdep)(k+0.5) = fse3(k+0.5) <br/>
[4213]	124	!>
[12080]	125	!> ** Action : - gdept, gdepw : depth of T- and W-point (m) <br/>
	126	!> - e3t, e3w : scale factors at T- and W-levels (m) <br/>
[4213]	127	!>
	128	!> @author G. Madec
[5617]	129	!> @date Marsh,2008 - F90: Free form and module
[12080]	130	!>
[4213]	131	!> @note Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
	132	!>
	133	!> @param[inout] dd_gdepw
	134	!> @param[inout] dd_gedpt
	135	!> @param[inout] dd_e3w
	136	!> @param[inout] dd_e2t
[13369]	137	!> @param[in] dd_ppkth
[4213]	138	!> @param[in] dd_ppkth2
	139	!> @param[in] dd_ppacr
	140	!> @param[in] dd_ppacr2
	141	!> @param[in] dd_ppdzmin
[13369]	142	!> @param[in] dd_pphmax
[12080]	143	!> @param[in] dd_pp_to_be_computed
[4213]	144	!> @param[in] dd_ppa1
[13369]	145	!> @param[in] dd_ppa2
[4213]	146	!> @param[in] dd_ppa0
	147	!> @param[in] dd_ppsur
	148	!-------------------------------------------------------------------
[12080]	149
[4213]	150	IMPLICIT NONE
[12080]	151
[13369]	152	! Argument
[4213]	153	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_gdepw
	154	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_gdept
	155	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3w
	156	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3t
[5609]	157	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3w_1d
	158	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3t_1d
[4213]	159
	160	REAL(dp) , INTENT(IN ) :: dd_ppkth
	161	REAL(dp) , INTENT(IN ) :: dd_ppkth2
	162	REAL(dp) , INTENT(IN ) :: dd_ppacr
	163	REAL(dp) , INTENT(IN ) :: dd_ppacr2
	164
	165	REAL(dp) , INTENT(IN ) :: dd_ppdzmin
	166	REAL(dp) , INTENT(IN ) :: dd_pphmax
[12080]	167	REAL(dp) , INTENT(IN ) :: dd_pp_to_be_computed
[4213]	168
	169	REAL(dp) , INTENT(IN ) :: dd_ppa0
	170	REAL(dp) , INTENT(IN ) :: dd_ppa1
	171	REAL(dp) , INTENT(IN ) :: dd_ppa2
	172	REAL(dp) , INTENT(IN ) :: dd_ppsur
	173
	174	! local variable
	175	REAL(dp) :: dl_zkth
	176	REAL(dp) :: dl_zkth2
	177	REAL(dp) :: dl_zdzmin
	178	REAL(dp) :: dl_zhmax
	179	REAL(dp) :: dl_zacr
	180	REAL(dp) :: dl_zacr2
	181
	182	REAL(dp) :: dl_ppacr
	183	REAL(dp) :: dl_ppacr2
	184
	185	REAL(dp) :: dl_za0
	186	REAL(dp) :: dl_za1
	187	REAL(dp) :: dl_za2
	188	REAL(dp) :: dl_zsur
	189	REAL(dp) :: dl_zw
	190	REAL(dp) :: dl_zt
	191
	192	INTEGER(i4) :: il_jpk
	193
	194	! loop indices
	195	INTEGER(i4) :: jk
	196	!----------------------------------------------------------------
	197
	198	dl_ppacr = dd_ppacr
	199	IF( dd_ppa1 == 0._dp ) dl_ppacr =1.0
	200	dl_ppacr2= dd_ppacr2
	201	IF( dd_ppa2 == 0._dp ) dl_ppacr2=1.0
	202
	203	! Set variables from parameters
	204	! ------------------------------
	205	dl_zkth = dd_ppkth ; dl_zacr = dl_ppacr
	206	dl_zdzmin = dd_ppdzmin ; dl_zhmax = dd_pphmax
	207	dl_zkth2 = dd_ppkth2 ; dl_zacr2 = dl_ppacr2
	208
	209	il_jpk = SIZE(dd_gdepw(:))
	210
	211	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
	212	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
	213	!
[12080]	214	IF( dd_ppa1 == dd_pp_to_be_computed .AND. &
	215	& dd_ppa0 == dd_pp_to_be_computed .AND. &
	216	& dd_ppsur == dd_pp_to_be_computed ) THEN
[4213]	217	dl_za1 = ( dl_zdzmin - dl_zhmax / REAL((il_jpk-1),dp) ) &
	218	& / ( TANH((1-dl_zkth)/dl_zacr) - dl_zacr/REAL((il_jpk-1),dp) &
	219	& * ( LOG( COSH( (REAL(il_jpk,dp) - dl_zkth) / dl_zacr) ) &
	220	& - LOG( COSH(( 1.0 - dl_zkth) / dl_zacr) ) ) )
	221
	222	dl_za0 = dl_zdzmin - dl_za1 * TANH( (1.0-dl_zkth) / dl_zacr )
	223	dl_zsur = - dl_za0 - dl_za1 * dl_zacr * LOG( COSH( (1-dl_zkth) / dl_zacr ) )
	224
	225	ELSE
	226	dl_za1 = dd_ppa1 ; dl_za0 = dd_ppa0 ; dl_zsur = dd_ppsur
	227	dl_za2 = dd_ppa2
	228	ENDIF
	229
	230	! Reference z-coordinate (depth - scale factor at T- and W-points)
	231	! ======================
[13369]	232	IF( dd_ppkth == 0. )THEN ! uniform vertical grid
[4213]	233
[13369]	234	dl_za1 = dl_zhmax/REAL((il_jpk-1),dp)
[4213]	235	DO jk = 1, il_jpk
	236	dl_zw = REAL(jk,dp)
[5609]	237	dl_zt = REAL(jk,dp) + 0.5_dp
[4213]	238	dd_gdepw(jk) = ( dl_zw - 1.0 ) * dl_za1
	239	dd_gdept(jk) = ( dl_zt - 1.0 ) * dl_za1
	240	dd_e3w (jk) = dl_za1
	241	dd_e3t (jk) = dl_za1
	242	END DO
	243
	244	ELSE
	245
	246	DO jk = 1, il_jpk
	247	dl_zw = REAL( jk,dp)
[5609]	248	dl_zt = REAL( jk,dp) + 0.5_dp
[4213]	249	dd_gdepw(jk) = ( dl_zsur + dl_za0 * dl_zw + &
	250	& dl_za1 * dl_zacr * LOG( COSH( (dl_zw-dl_zkth)/dl_zacr ) ) + &
	251	& dl_za2 * dl_zacr2* LOG( COSH( (dl_zw-dl_zkth2)/dl_zacr2 ) ) )
	252	dd_gdept(jk) = ( dl_zsur + dl_za0 * dl_zt + &
	253	& dl_za1 * dl_zacr * LOG( COSH( (dl_zt-dl_zkth)/dl_zacr ) ) + &
	254	& dl_za2 * dl_zacr2* LOG( COSH( (dl_zt-dl_zkth2)/dl_zacr2 ) ) )
	255	dd_e3w (jk) = dl_za0 + &
	256	& dl_za1 * TANH( (dl_zw-dl_zkth)/dl_zacr ) + &
	257	& dl_za2 * TANH( (dl_zw-dl_zkth2)/dl_zacr2 )
	258	dd_e3t (jk) = dl_za0 + &
	259	& dl_za1 * TANH( (dl_zt-dl_zkth)/dl_zacr ) + &
	260	& dl_za2 * TANH( (dl_zt-dl_zkth2)/dl_zacr2 )
	261	END DO
	262	dd_gdepw(1) = 0.e0 ! force first w-level to be exactly at zero
	263
	264	ENDIF
	265
[5609]	266	! need to be like this to compute the pressure gradient with ISF.
	267	! If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
	268	! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
	269	DO jk = 1, il_jpk-1
[13369]	270	dd_e3t_1d(jk) = dd_gdepw(jk+1)-dd_gdepw(jk)
[5609]	271	END DO
	272	dd_e3t_1d(il_jpk) = dd_e3t_1d(il_jpk-1) ! we don't care because this level is masked in NEMO
	273
	274	DO jk = 2, il_jpk
[13369]	275	dd_e3w_1d(jk) = dd_gdept(jk) - dd_gdept(jk-1)
[5609]	276	END DO
	277	dd_e3w_1d(1 ) = 2._dp * (dd_gdept(1) - dd_gdepw(1))
	278
[4213]	279	! Control and print
	280	! ==================
	281
	282	DO jk = 1, il_jpk
	283	IF( dd_e3w(jk) <= 0. .OR. dd_e3t(jk) <= 0. )then
[5609]	284	CALL logger_debug("VGRID ZGR Z: e3w or e3t <= 0 ")
[13369]	285	ENDIF
[4213]	286
[5609]	287	IF( dd_e3w_1d(jk) <= 0. .OR. dd_e3t_1d(jk) <= 0. )then
	288	CALL logger_debug("VGRID ZGR Z: e3w_1d or e3t_1d <= 0 ")
	289	ENDIF
	290
[4213]	291	IF( dd_gdepw(jk) < 0. .OR. dd_gdept(jk) < 0. )then
	292	CALL logger_debug("VGRID ZGR Z: gdepw or gdept < 0 ")
	293	ENDIF
	294	END DO
	295
	296	END SUBROUTINE vgrid_zgr_z
[12080]	297	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	298	SUBROUTINE vgrid_zgr_bat(dd_bathy, dd_gdepw, dd_hmin, dd_fill)
[4213]	299	!-------------------------------------------------------------------
[6393]	300	!> @brief This subroutine
	301	!>
	302	!> @todo add subroutine description
	303	!>
	304	!> @param[inout] dd_bathy
	305	!> @param[in] dd_gdepw
	306	!> @param[in] dd_hmin
	307	!> @param[in] dd_fill
[5609]	308	!-------------------------------------------------------------------
[12080]	309
[5609]	310	IMPLICIT NONE
[12080]	311
[5609]	312	! Argument
[13369]	313	REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: dd_bathy
	314	REAL(dp), DIMENSION(:) , INTENT(IN ) :: dd_gdepw
[5609]	315	REAL(dp) , INTENT(IN ) :: dd_hmin
	316	REAL(dp) , INTENT(IN ), OPTIONAL :: dd_fill
	317
	318	! local
	319	INTEGER(i4) :: il_jpk
[13369]	320
[5609]	321	REAL(dp) :: dl_hmin
	322	REAL(dp) :: dl_fill
	323
	324	! loop indices
	325	INTEGER(i4) :: jk
	326	!----------------------------------------------------------------
	327	il_jpk = SIZE(dd_gdepw(:))
	328
	329	dl_fill=0._dp
	330	IF( PRESENT(dd_fill) ) dl_fill=dd_fill
	331
	332	IF( dd_hmin < 0._dp ) THEN
	333	jk = - INT( dd_hmin ) ! from a nb of level
	334	ELSE
	335	jk = MINLOC( dd_gdepw, mask = dd_gdepw > dd_hmin, dim = 1 ) ! from a depth
	336	ENDIF
[13369]	337
	338	dl_hmin = dd_gdepw(jk+1) ! minimum depth = ik+1 w-levels
[5609]	339	WHERE( dd_bathy(:,:) <= 0._wp .OR. dd_bathy(:,:) == dl_fill )
	340	dd_bathy(:,:) = dl_fill ! min=0 over the lands
	341	ELSE WHERE
	342	dd_bathy(:,:) = MAX( dl_hmin , dd_bathy(:,:) ) ! min=dl_hmin over the oceans
	343	END WHERE
[13369]	344	WRITE(,) 'Minimum ocean depth: ', dl_hmin, ' minimum number of ocean levels : ', jk
[5609]	345
	346	END SUBROUTINE vgrid_zgr_bat
[12080]	347	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	348	SUBROUTINE vgrid_zgr_zps(id_mbathy, dd_bathy, id_jpkmax, &
	349	& dd_gdepw, dd_e3t, &
	350	& dd_e3zps_min, dd_e3zps_rat, &
	351	& dd_fill )
[5609]	352	!-------------------------------------------------------------------
[4213]	353	!> @brief This subroutine set the depth and vertical scale factor in partial step
[13369]	354	!> z-coordinate case
[12080]	355	!>
[4213]	356	!> @details
	357	!> ** Method : Partial steps : computes the 3D vertical scale factors
	358	!> of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	359	!> a partial step representation of bottom topography.
	360	!>
	361	!> The reference depth of model levels is defined from an analytical
	362	!> function the derivative of which gives the reference vertical
	363	!> scale factors.
[5037]	364	!> From depth and scale factors reference, we compute there new value
[4213]	365	!> with partial steps on 3d arrays ( i, j, k ).
	366	!>
[13369]	367	!> w-level:
[5037]	368	!> - gdepw_ps(i,j,k) = fsdep(k)
	369	!> - e3w_ps(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
[13369]	370	!> t-level:
[5037]	371	!> - gdept_ps(i,j,k) = fsdep(k+0.5)
	372	!> - e3t_ps(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
[4213]	373	!>
[5037]	374	!> With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
[4213]	375	!> we find the mbathy index of the depth at each grid point.
	376	!> This leads us to three cases:
[5037]	377	!> - bathy = 0 => mbathy = 0
[13369]	378	!> - 1 < mbathy < jpkm1
	379	!> - bathy > gdepw(jpk) => mbathy = jpkm1
[4213]	380	!>
[5037]	381	!> Then, for each case, we find the new depth at t- and w- levels
[13369]	382	!> and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
[4213]	383	!> and f-points.
[13369]	384	!>
[4213]	385	!> This routine is given as an example, it must be modified
	386	!> following the user s desiderata. nevertheless, the output as
	387	!> well as the way to compute the model levels and scale factors
	388	!> must be respected in order to insure second order accuracy
	389	!> schemes.
	390	!>
[13369]	391	!> @warning
[5037]	392	!> - gdept, gdepw and e3 are positives
	393	!> - gdept_ps, gdepw_ps and e3_ps are positives
[6393]	394	!>
[4213]	395	!> @author A. Bozec, G. Madec
[5617]	396	!> @date February, 2009 - F90: Free form and module
[13369]	397	!> @date February, 2009
[5617]	398	!> - A. de Miranda : rigid-lid + islands
[4213]	399	!>
	400	!> @note Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
	401	!>
	402	!> @param[inout] id_mbathy
	403	!> @param[inout] dd_bathy
[13369]	404	!> @param[inout] id_jpkmax
[4213]	405	!> @param[in] dd_gdepw
[13369]	406	!> @param[in] dd_e3t
[4213]	407	!> @param[in] dd_e3zps_min
	408	!> @param[in] dd_e3zps_rat
[6393]	409	!> @param[in] dd_fill
[4213]	410	!-------------------------------------------------------------------
[12080]	411
[4213]	412	IMPLICIT NONE
[12080]	413
[13369]	414	! Argument
[4213]	415	INTEGER(i4), DIMENSION(:,:), INTENT( OUT) :: id_mbathy
	416	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_bathy
	417	INTEGER(i4) , INTENT(INOUT) :: id_jpkmax
	418	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_gdepw
	419	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_e3t
[5609]	420	REAL(dp) , INTENT(IN ) :: dd_e3zps_min
	421	REAL(dp) , INTENT(IN ) :: dd_e3zps_rat
	422	REAL(dp) , INTENT(IN ), OPTIONAL :: dd_fill
[4213]	423
	424	! local variable
	425	REAL(dp) :: dl_zmax ! Maximum depth
[5609]	426	!REAL(dp) :: dl_zmin ! Minimum depth
[13369]	427	REAL(dp) :: dl_zdepth ! Ajusted ocean depth to avoid too small e3t
	428	REAL(dp) :: dl_fill
[4213]	429
	430	INTEGER(i4) :: il_jpk
	431	INTEGER(i4) :: il_jpkm1
	432	INTEGER(i4) :: il_jpiglo
	433	INTEGER(i4) :: il_jpjglo
	434
	435	! loop indices
	436	INTEGER(i4) :: ji
	437	INTEGER(i4) :: jj
	438	INTEGER(i4) :: jk
	439	!----------------------------------------------------------------
	440
	441	il_jpk=SIZE(dd_gdepw(:))
	442	il_jpiglo=SIZE(id_mbathy(:,:),DIM=1)
	443	il_jpjglo=SIZE(id_mbathy(:,:),DIM=2)
	444
[5609]	445	dl_fill=0._dp
	446	IF( PRESENT(dd_fill) ) dl_fill=dd_fill
	447
[4213]	448	! Initialization of constant
[5609]	449	dl_zmax = dd_gdepw(il_jpk) + dd_e3t(il_jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
[4213]	450
[5609]	451	! bounded value of bathy (min already set at the end of zgr_bat)
	452	WHERE( dd_bathy(:,:) /= dl_fill )
	453	dd_bathy(:,:) = MIN( dl_zmax , dd_bathy(:,:) )
	454	END WHERE
	455
[4213]	456	! bathymetry in level (from bathy_meter)
	457	! ===================
	458	il_jpkm1=il_jpk-1
	459	! initialize mbathy to the maximum ocean level available
	460	id_mbathy(:,:) = il_jpkm1
	461
	462	! storage of land and island's number (zera and negative values) in mbathy
	463	DO jj = 1, il_jpjglo
	464	DO ji= 1, il_jpiglo
[5609]	465	IF( dd_bathy(ji,jj) <= 0._dp )THEN
	466	id_mbathy(ji,jj) = INT(dd_bathy(ji,jj),i4)
	467	ELSEIF( dd_bathy(ji,jj) == dl_fill )THEN
	468	id_mbathy(ji,jj) = 0_i4
[4213]	469	ENDIF
	470	END DO
	471	END DO
	472
	473	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	474	! find the number of ocean levels such that the last level thickness
	475	! is larger than the minimum of e3zps_min and e3zps_rat * e3t (where
	476	! e3t is the reference level thickness
	477
	478	DO jk = il_jpkm1, 1, -1
	479	dl_zdepth = dd_gdepw(jk) + MIN( dd_e3zps_min, dd_e3t(jk)*dd_e3zps_rat )
	480
	481	DO jj = 1, il_jpjglo
	482	DO ji = 1, il_jpiglo
[5609]	483	IF( dd_bathy(ji,jj) /= dl_fill )THEN
	484	IF( 0. < dd_bathy(ji,jj) .AND. &
	485	& dd_bathy(ji,jj) <= dl_zdepth ) id_mbathy(ji,jj) = jk-1
	486	ENDIF
[4213]	487	END DO
	488	END DO
	489	END DO
	490
	491	! ================
	492	! Bathymetry check
	493	! ================
	494
	495	CALL vgrid_zgr_bat_ctl( id_mbathy, id_jpkmax, il_jpk)
	496
	497	END SUBROUTINE vgrid_zgr_zps
[12080]	498	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	499	SUBROUTINE vgrid_zgr_bat_ctl(id_mbathy, id_jpkmax, id_jpk)
[4213]	500	!-------------------------------------------------------------------
[13369]	501	!> @brief This subroutine check the bathymetry in levels
[12080]	502	!>
[4213]	503	!> @details
	504	!> ** Method : The array mbathy is checked to verified its consistency
	505	!> with the model options. in particular:
	506	!> mbathy must have at least 1 land grid-points (mbathy<=0)
	507	!> along closed boundary.
	508	!> mbathy must be cyclic IF jperio=1.
	509	!> mbathy must be lower or equal to jpk-1.
	510	!> isolated ocean grid points are suppressed from mbathy
	511	!> since they are only connected to remaining
	512	!> ocean through vertical diffusion.
	513	!> C A U T I O N : mbathy will be modified during the initializa-
	514	!> tion phase to become the number of non-zero w-levels of a water
	515	!> column, with a minimum value of 1.
	516	!>
	517	!> ** Action : - update mbathy: level bathymetry (in level index)
	518	!> - update bathy : meter bathymetry (in meters)
[6393]	519	!>
[4213]	520	!> @author G.Madec
[5617]	521	!> @date Marsh, 2008 - Original code
[12080]	522	!>
[13369]	523	!> @param[in] id_mbathy
[5037]	524	!> @param[in] id_jpkmax
	525	!> @param[in] id_jpk
[4213]	526	!-------------------------------------------------------------------
[12080]	527
[4213]	528	IMPLICIT NONE
[12080]	529
[13369]	530	! Argument
[4213]	531	INTEGER(i4), DIMENSION(:,:), INTENT(INOUT) :: id_mbathy
	532	INTEGER(i4) , INTENT(INOUT) :: id_jpkmax
	533	INTEGER(i4) , INTENT(INOUT) :: id_jpk
	534
	535	! local variable
	536	INTEGER(i4) :: il_jpiglo
	537	INTEGER(i4) :: il_jpjglo
	538
	539	INTEGER(i4) :: il_icompt
	540	INTEGER(i4) :: il_ibtest
	541	INTEGER(i4) :: il_ikmax
	542	INTEGER(i4) :: il_jpkm1
	543
	544	INTEGER(i4) :: il_jim
	545	INTEGER(i4) :: il_jip
	546	INTEGER(i4) :: il_jjm
	547	INTEGER(i4) :: il_jjp
	548
	549	! loop indices
	550	INTEGER(i4) :: jl
	551	INTEGER(i4) :: ji
	552	INTEGER(i4) :: jj
	553	!----------------------------------------------------------------
	554
	555	il_jpiglo=SIZE(id_mbathy(:,:),DIM=1)
	556	il_jpjglo=SIZE(id_mbathy(:,:),DIM=2)
	557
	558	! ================
	559	! Bathymetry check
	560	! ================
	561
	562	! suppress isolated ocean grid points'
	563	il_icompt = 0
	564
	565	DO jl = 1, 2
	566	DO jj = 1, il_jpjglo
	567	DO ji = 1, il_jpiglo
	568	il_jim=max(ji-1,1) ; il_jip=min(ji+1,il_jpiglo)
	569	il_jjm=max(jj-1,1) ; il_jjp=min(jj+1,il_jpjglo)
	570
	571	if(il_jim==ji) il_jim=il_jip ; if(il_jip==ji) il_jip=il_jim
	572	if(il_jjm==jj) il_jjm=il_jjp ; if(il_jjp==jj) il_jjp=il_jjm
	573
	574	il_ibtest = MAX( id_mbathy(il_jim,jj), id_mbathy(il_jip,jj), &
	575	id_mbathy(ji,il_jjm),id_mbathy(ji,il_jjp) )
	576
	577	IF( il_ibtest < id_mbathy(ji,jj) ) THEN
	578	id_mbathy(ji,jj) = il_ibtest
	579	il_icompt = il_icompt + 1
	580	ENDIF
	581	END DO
	582	END DO
	583
	584	END DO
	585	IF( il_icompt == 0 ) THEN
	586	CALL logger_info("VGRID ZGR BAT CTL: no isolated ocean grid points")
	587	ELSE
	588	CALL logger_info("VGRID ZGR BAT CTL:"//TRIM(fct_str(il_icompt))//&
	589	& " ocean grid points suppressed")
	590	ENDIF
	591
	592	id_mbathy(:,:) = MAX( 0, id_mbathy(:,:))
	593
	594	! Number of ocean level inferior or equal to jpkm1
	595
	596	il_ikmax = 0
	597	DO jj = 1, il_jpjglo
	598	DO ji = 1, il_jpiglo
	599	il_ikmax = MAX( il_ikmax, id_mbathy(ji,jj) )
	600	END DO
	601	END DO
	602
	603	id_jpkmax=id_jpk
	604
	605	il_jpkm1=id_jpk-1
	606	IF( il_ikmax > il_jpkm1 ) THEN
	607	CALL logger_error("VGRID ZGR BAT CTL: maximum number of ocean level = "//&
	608	& TRIM(fct_str(il_ikmax))//" > jpk-1."//&
	609	& " Change jpk to "//TRIM(fct_str(il_ikmax+1))//&
	610	& " to use the exact ead bathymetry" )
	611	ELSE IF( il_ikmax < il_jpkm1 ) THEN
	612	id_jpkmax=il_ikmax+1
	613	ENDIF
	614
	615	END SUBROUTINE vgrid_zgr_bat_ctl
[12080]	616	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	617	FUNCTION vgrid_get_level(td_bathy, cd_namelist, td_dom, id_nlevel) &
	618	& RESULT (tf_var)
[4213]	619	!-------------------------------------------------------------------
[13369]	620	!> @brief This function compute bathy level in T,U,V,F point, and return
[5037]	621	!> them as array of variable structure
[12080]	622	!>
[4213]	623	!> @details
[13369]	624	!> Bathymetry is read on Bathymetry file, then bathy level is computed
[5037]	625	!> on T point, and finally fit to U,V,F point.
	626	!>
	627	!> you could specify :<br/>
	628	!> - namelist where find parameter to set the depth of model levels
	629	!> (default use GLORYS 75 levels parameters)
[13369]	630	!> - domain structure to specify one area to work on
[5037]	631	!> - number of level to be used
	632	!>
[4213]	633	!> @author J.Paul
[5617]	634	!> @date November, 2013 - Initial Version
[12080]	635	!>
[13369]	636	!> @param[in] td_bathy Bathymetry file structure
	637	!> @param[in] cd_namelist namelist
[5037]	638	!> @param[in] td_dom domain structure
[13369]	639	!> @param[in] id_nlevel number of lelvel to be used
[5037]	640	!> @return array of level on T,U,V,F point (variable structure)
[4213]	641	!-------------------------------------------------------------------
[12080]	642
[4213]	643	IMPLICIT NONE
[12080]	644
[4213]	645	! Argument
[5037]	646	TYPE(TMPP) , INTENT(IN) :: td_bathy
	647	CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: cd_namelist
[4213]	648	TYPE(TDOM) , INTENT(IN), OPTIONAL :: td_dom
	649	INTEGER(i4) , INTENT(IN), OPTIONAL :: id_nlevel
	650
	651	! function
[12080]	652	TYPE(TVAR), DIMENSION(ip_npoint) :: tf_var
[4213]	653
	654	! local variable
[13369]	655	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_gdepw
	656	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_gdept
	657	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_e3w
[4213]	658	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_e3t
[13369]	659	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_e3w_1d
[5609]	660	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_e3t_1d
[4213]	661
	662	INTEGER(i4) :: il_status
	663	INTEGER(i4) :: il_fileid
	664	INTEGER(i4) :: il_jpkmax
[5037]	665	INTEGER(i4), DIMENSION(2,2) :: il_xghost
[4213]	666	INTEGER(i4), DIMENSION(:,:) , ALLOCATABLE :: il_mbathy
	667	INTEGER(i4), DIMENSION(:,:,:,:), ALLOCATABLE :: il_level
[13369]	668
[4213]	669	LOGICAL :: ll_exist
	670
[5037]	671	TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim
	672
	673	TYPE(TDOM) :: tl_dom
	674
	675	TYPE(TVAR) :: tl_var
	676
	677	TYPE(TMPP) :: tl_bathy
	678
[4213]	679	! loop indices
	680	INTEGER(i4) :: ji
	681	INTEGER(i4) :: jj
	682
	683	INTEGER(i4) :: jip
	684	INTEGER(i4) :: jjp
	685
	686	!namelist (intialise with GLORYS 75 levels parameters)
[12080]	687	REAL(dp) :: dn_pp_to_be_computed = 0._dp
[4213]	688	REAL(dp) :: dn_ppsur = -3958.951371276829_dp
	689	REAL(dp) :: dn_ppa0 = 103.9530096000000_dp
	690	REAL(dp) :: dn_ppa1 = 2.4159512690000_dp
	691	REAL(dp) :: dn_ppa2 = 100.7609285000000_dp
	692	REAL(dp) :: dn_ppkth = 15.3510137000000_dp
	693	REAL(dp) :: dn_ppkth2 = 48.0298937200000_dp
	694	REAL(dp) :: dn_ppacr = 7.0000000000000_dp
	695	REAL(dp) :: dn_ppacr2 = 13.000000000000_dp
	696	REAL(dp) :: dn_ppdzmin = 6._dp
	697	REAL(dp) :: dn_pphmax = 5750._dp
	698	INTEGER(i4) :: in_nlevel = 75
	699
	700	REAL(dp) :: dn_e3zps_min = 25._dp
	701	REAL(dp) :: dn_e3zps_rat = 0.2_dp
	702	!----------------------------------------------------------------
	703	NAMELIST /namzgr/ &
[12080]	704	& dn_pp_to_be_computed, &
[4213]	705	& dn_ppsur, &
	706	& dn_ppa0, &
	707	& dn_ppa1, &
	708	& dn_ppa2, &
	709	& dn_ppkth, &
	710	& dn_ppkth2, &
	711	& dn_ppacr, &
	712	& dn_ppacr2, &
	713	& dn_ppdzmin, &
	714	& dn_pphmax, &
	715	& in_nlevel
	716
	717	NAMELIST /namzps/ &
	718	& dn_e3zps_min, &
	719	& dn_e3zps_rat
	720	!----------------------------------------------------------------
	721
[5037]	722	IF( PRESENT(cd_namelist) )THEN
	723	!1- read namelist
	724	INQUIRE(FILE=TRIM(cd_namelist), EXIST=ll_exist)
	725	IF( ll_exist )THEN
[13369]	726
[5037]	727	il_fileid=fct_getunit()
[4213]	728
[5037]	729	OPEN( il_fileid, FILE=TRIM(cd_namelist), &
	730	& FORM='FORMATTED', &
	731	& ACCESS='SEQUENTIAL', &
	732	& STATUS='OLD', &
	733	& ACTION='READ', &
	734	& IOSTAT=il_status)
	735	CALL fct_err(il_status)
	736	IF( il_status /= 0 )THEN
	737	CALL logger_fatal("VGRID GET LEVEL: ERROR opening "//&
	738	& TRIM(cd_namelist))
	739	ENDIF
[4213]	740
[5037]	741	READ( il_fileid, NML = namzgr )
	742	READ( il_fileid, NML = namzps )
[4213]	743
[5037]	744	CLOSE( il_fileid, IOSTAT=il_status )
	745	CALL fct_err(il_status)
	746	IF( il_status /= 0 )THEN
	747	CALL logger_error("VGRID GET LEVELL: ERROR closing "//&
	748	& TRIM(cd_namelist))
	749	ENDIF
[4213]	750
[5037]	751	ELSE
[4213]	752
[5037]	753	CALL logger_fatal("VGRID GET LEVEL: ERROR. can not find "//&
	754	& TRIM(cd_namelist))
[4213]	755
[5037]	756	ENDIF
[4213]	757	ENDIF
	758
[5037]	759	! copy structure
	760	tl_bathy=mpp_copy(td_bathy)
	761
	762	! get domain
[4213]	763	IF( PRESENT(td_dom) )THEN
[5037]	764	tl_dom=dom_copy(td_dom)
[4213]	765	ELSE
	766	CALL logger_debug("VGRID GET LEVEL: get dom from "//&
	767	& TRIM(tl_bathy%c_name))
	768	tl_dom=dom_init(tl_bathy)
	769	ENDIF
	770
[12080]	771	! get ghost cell
[5037]	772	il_xghost(:,:)=grid_get_ghost(tl_bathy)
[4213]	773
[5037]	774	! open mpp files
	775	CALL iom_dom_open(tl_bathy, tl_dom)
[4213]	776
[5037]	777	! check namelist
[4213]	778	IF( PRESENT(id_nlevel) ) in_nlevel=id_nlevel
	779	IF( in_nlevel == 0 )THEN
	780	CALL logger_fatal("VGRID GET LEVEL: number of level to be used "//&
	781	& "is not specify. check namelist.")
	782	ENDIF
	783
[5037]	784	! read bathymetry
	785	tl_var=iom_dom_read_var(tl_bathy,'bathymetry',tl_dom)
	786	! clean
	787	CALL dom_clean(tl_dom)
[4213]	788
[5037]	789	! remove ghost cell
	790	CALL grid_del_ghost(tl_var, il_xghost(:,:))
	791
	792	! force _FillValue (land) to be 0
[4213]	793	WHERE( tl_var%d_value(:,:,1,1) == tl_var%d_fill )
	794	tl_var%d_value(:,:,1,1)=0
	795	END WHERE
	796
[5037]	797	! clean
	798	CALL iom_dom_close(tl_bathy)
	799	CALL mpp_clean(tl_bathy)
[4213]	800
[5037]	801	! compute vertical grid
[13369]	802	ALLOCATE( dl_gdepw(in_nlevel), dl_gdept(in_nlevel) )
	803	ALLOCATE( dl_e3w(in_nlevel), dl_e3t(in_nlevel) )
	804	ALLOCATE( dl_e3w_1d(in_nlevel), dl_e3t_1d(in_nlevel) )
[4213]	805	CALL vgrid_zgr_z( dl_gdepw(:), dl_gdept(:), dl_e3w(:), dl_e3t(:), &
[5609]	806	& dl_e3w_1d, dl_e3t_1d, &
[4213]	807	& dn_ppkth, dn_ppkth2, dn_ppacr, dn_ppacr2, &
[12080]	808	& dn_ppdzmin, dn_pphmax, dn_pp_to_be_computed, &
[4213]	809	& dn_ppa0, dn_ppa1, dn_ppa2, dn_ppsur )
	810
[5037]	811	! compute bathy level on T point
[4213]	812	ALLOCATE( il_mbathy(tl_var%t_dim(1)%i_len, &
	813	& tl_var%t_dim(2)%i_len ) )
	814	CALL vgrid_zgr_zps( il_mbathy(:,:), tl_var%d_value(:,:,1,1), il_jpkmax, &
	815	& dl_gdepw(:), dl_e3t(:), &
	816	& dn_e3zps_min, dn_e3zps_rat )
	817
[13369]	818	DEALLOCATE( dl_gdepw, dl_gdept )
[4213]	819	DEALLOCATE( dl_e3w, dl_e3t )
	820
[5037]	821	! compute bathy level in T,U,V,F point
[4213]	822	ALLOCATE( il_level(tl_var%t_dim(1)%i_len, &
	823	& tl_var%t_dim(2)%i_len, &
[5037]	824	& ip_npoint,1) )
[4213]	825
	826	DO jj=1,tl_var%t_dim(2)%i_len
	827	DO ji= 1,tl_var%t_dim(1)%i_len
	828
	829	jip=MIN(ji+1,tl_var%t_dim(1)%i_len)
	830	jjp=MIN(jj+1,tl_var%t_dim(2)%i_len)
	831
	832	! T point
	833	il_level(ji,jj,jp_T,1)=il_mbathy(ji,jj)
	834	! U point
	835	il_level(ji,jj,jp_U,1)=MIN( il_mbathy(ji, jj ), il_mbathy(jip, jj ))
	836	! V point
	837	il_level(ji,jj,jp_V,1)=MIN( il_mbathy(ji, jj ), il_mbathy(ji , jjp))
	838	! F point
	839	il_level(ji,jj,jp_F,1)=MIN( il_mbathy(ji, jj ), il_mbathy(jip, jj ), &
	840	& il_mbathy(ji, jjp), il_mbathy(jip, jjp))
	841
	842	ENDDO
	843	ENDDO
	844
	845	DEALLOCATE( il_mbathy )
	846
[5037]	847	tl_dim(:)=dim_copy(tl_var%t_dim(:))
	848	! clean
[4213]	849	CALL var_clean(tl_var)
	850
	851	! only 2 first dimension to be used
	852	tl_dim(3:4)%l_use=.FALSE.
	853
[12080]	854	tf_var(jp_T)=var_init('tlevel', il_level(:,:,jp_T:jp_T,:), td_dim=tl_dim(:))
	855	tf_var(jp_U)=var_init('ulevel', il_level(:,:,jp_U:jp_U,:), td_dim=tl_dim(:))
	856	tf_var(jp_V)=var_init('vlevel', il_level(:,:,jp_V:jp_V,:), td_dim=tl_dim(:))
	857	tf_var(jp_F)=var_init('flevel', il_level(:,:,jp_F:jp_F,:), td_dim=tl_dim(:))
[4213]	858
	859	DEALLOCATE( il_level )
	860
[12080]	861	CALL grid_add_ghost( tf_var(jp_T), il_xghost(:,:) )
	862	CALL grid_add_ghost( tf_var(jp_U), il_xghost(:,:) )
	863	CALL grid_add_ghost( tf_var(jp_V), il_xghost(:,:) )
	864	CALL grid_add_ghost( tf_var(jp_F), il_xghost(:,:) )
[4213]	865
[5037]	866	! clean
	867	CALL dim_clean(tl_dim(:))
[4213]	868
	869	END FUNCTION vgrid_get_level
[12080]	870	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[4213]	871	END MODULE vgrid
	872

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

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

Context Navigation

source: utils/tools/SIREN/src/vgrid.f90

Download in other formats: