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.

grid_hgr.f90 in utils/tools/SIREN/src – NEMO

Context Navigation

source: utils/tools/SIREN/src/grid_hgr.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 4 years ago
update nemo trunk
File size: 59.7 KB

Rev	Line
[7314]	1	!----------------------------------------------------------------------
	2	! NEMO system team, System and Interface for oceanic RElocable Nesting
	3	!----------------------------------------------------------------------
	4	!
	5	! DESCRIPTION:
	6	!> @brief This module manage Horizontal grid.
	7	!>
	8	!> @details
	9	!> ** Purpose : Compute the geographical position (in degre) of the
	10	!> model grid-points, the horizontal scale factors (in meters) and
	11	!> the Coriolis factor (in s-1).
	12	!>
	13	!> ** Method : The geographical position of the model grid-points is
	14	!> defined from analytical functions, fslam and fsphi, the derivatives of which gives the horizontal scale factors e1,e2.
	15	!> Defining two function fslam and fsphi and their derivatives in the two horizontal directions (fse1 and fse2),
	16	!> the model grid-point position and scale factors are given by:
	17	!> - t-point:
	18	!> - glamt(i,j) = fslam(i ,j ) e1t(i,j) = fse1(i ,j )
	19	!> - gphit(i,j) = fsphi(i ,j ) e2t(i,j) = fse2(i ,j )
	20	!> - u-point:
	21	!> - glamu(i,j) = fslam(i+1/2,j ) e1u(i,j) = fse1(i+1/2,j )
	22	!> - gphiu(i,j) = fsphi(i+1/2,j ) e2u(i,j) = fse2(i+1/2,j )
	23	!> - v-point:
	24	!> - glamv(i,j) = fslam(i ,j+1/2) e1v(i,j) = fse1(i ,j+1/2)
	25	!> - gphiv(i,j) = fsphi(i ,j+1/2) e2v(i,j) = fse2(i ,j+1/2)
	26	!> - f-point:
	27	!> - glamf(i,j) = fslam(i+1/2,j+1/2) e1f(i,j) = fse1(i+1/2,j+1/2)
	28	!> - gphif(i,j) = fsphi(i+1/2,j+1/2) e2f(i,j) = fse2(i+1/2,j+1/2)
	29	!>
	30	!> Where fse1 and fse2 are defined by:
	31	!> - fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
	32	!> + di(fsphi) **2 )(i,j)
	33	!> - fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
	34	!> + dj(fsphi) **2 )(i,j)
	35	!>
	36	!> The coriolis factor is given at z-point by:
	37	!> - ff = 2.omegasin(gphif) (in s-1)<br/>
	38	!>
	39	!> This routine is given as an example, it must be modified
	40	!> following the user s desiderata. nevertheless, the output as
	41	!> well as the way to compute the model grid-point position and
	42	!> horizontal scale factors must be respected in order to insure
	43	!> second order accuracy schemes.
	44	!>
	45	!> @note If the domain is periodic, verify that scale factors are also
	46	!> periodic, and the coriolis term again.
	47	!>
	48	!> ** Action :
	49	!> - define glamt, glamu, glamv, glamf: longitude of t-, u-, v- and f-points (in degre)
	50	!> - define gphit, gphiu, gphiv, gphit: latitude of t-, u-, v- and f-points (in degre)
	51	!> - define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
	52	!> - scale factors (in meters) at t-, u-, v-, and f-points.
	53	!> - define ff: coriolis factor at f-point
	54	!>
	55	!> References : Marti, Madec and Delecluse, 1992, JGR
	56	!> Madec, Imbard, 1996, Clim. Dyn.
	57	!>
	58	!> @author
	59	!> G, Madec
[12080]	60	!>
[7314]	61	!> @date March, 1988 - Original code
	62	!> @date January, 1996
	63	!> - terrain following coordinates
	64	!> @date February, 1997
	65	!> - print mesh informations
	66	!> @date November, 1999
	67	!> - M. Imbard : NetCDF format with IO-IPSL
	68	!> @date Augustr, 2000
	69	!> - D. Ludicone : Reduced section at Bab el Mandeb
	70	!> @date September, 2001
	71	!> - M. Levy : eel config: grid in km, beta-plane
	72	!> @date August, 2002
	73	!> - G. Madec : F90: Free form and module, namelist
	74	!> @date January, 2004
	75	!> - A.M. Treguier, J.M. Molines : Case 4 (Mercator mesh)
	76	!> use of parameters in par_CONFIG-Rxx.h90, not in namelist
	77	!> @date May, 2004
	78	!> - A. Koch-Larrouy : Add Gyre configuration
	79	!> @date February, 2011
	80	!> - G. Madec : add cell surface (e1e2t)
	81	!> @date September, 2015
	82	!> - J, Paul : rewrite to SIREN format from $Id: domhgr.F90 5506 2015-06-29 15:19:38Z clevy $
	83	!> @date October, 2016
	84	!> - J, Paul : update from trunk (revision 6961): add wetting and drying, ice sheet coupling..
	85	!> - J, Paul : compute coriolis factor at f-point and at t-point
	86	!> - J, Paul : do not use anymore special case for ORCA grid
	87	!>
[12080]	88	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[7314]	89	!----------------------------------------------------------------------
	90	MODULE grid_hgr
[12080]	91
[7314]	92	USE netcdf ! nf90 library
	93	USE kind ! F90 kind parameter
	94	USE fct ! basic usefull function
	95	USE global ! global parameter
	96	USE phycst ! physical constant
	97	USE logger ! log file manager
	98	USE file ! file manager
	99	USE var ! variable manager
	100	USE dim ! dimension manager
	101	USE dom ! domain manager
	102	USE grid ! grid manager
	103	USE iom ! I/O manager
	104	USE mpp ! MPP manager
	105	USE iom_mpp ! I/O MPP manager
	106	USE lbc ! lateral boundary conditions
[12080]	107
[7314]	108	IMPLICIT NONE
	109	! NOTE_avoid_public_variables_if_possible
	110
	111	! type and variable
	112	PUBLIC :: TNAMH
	113
	114	PUBLIC :: tg_tmask
	115	PUBLIC :: tg_umask
	116	PUBLIC :: tg_vmask
	117	PUBLIC :: tg_fmask
	118
	119	! PUBLIC :: tg_wmask
	120	! PUBLIC :: tg_wumask
	121	! PUBLIC :: tg_wvmask
	122
	123	PUBLIC :: tg_ssmask
	124	! PUBLIC :: tg_ssumask
	125	! PUBLIC :: tg_ssvmask
	126	! PUBLIC :: tg_ssfmask
	127
	128	PUBLIC :: tg_glamt
	129	PUBLIC :: tg_glamu
	130	PUBLIC :: tg_glamv
	131	PUBLIC :: tg_glamf
	132
	133	PUBLIC :: tg_gphit
	134	PUBLIC :: tg_gphiu
	135	PUBLIC :: tg_gphiv
	136	PUBLIC :: tg_gphif
	137
	138	PUBLIC :: tg_e1t
	139	PUBLIC :: tg_e1u
	140	PUBLIC :: tg_e1v
	141	PUBLIC :: tg_e1f
	142
	143	PUBLIC :: tg_e2t
	144	PUBLIC :: tg_e2u
	145	PUBLIC :: tg_e2v
	146	PUBLIC :: tg_e2f
	147
	148	PUBLIC :: tg_ff_t
	149	PUBLIC :: tg_ff_f
	150
	151	PUBLIC :: tg_gcost
	152	PUBLIC :: tg_gcosu
	153	PUBLIC :: tg_gcosv
	154	PUBLIC :: tg_gcosf
	155
	156	PUBLIC :: tg_gsint
	157	PUBLIC :: tg_gsinu
	158	PUBLIC :: tg_gsinv
	159	PUBLIC :: tg_gsinf
	160
	161	! function and subroutine
	162	PUBLIC :: grid_hgr_init
	163	PUBLIC :: grid_hgr_fill
	164	PUBLIC :: grid_hgr_clean
	165	PUBLIC :: grid_hgr_nam
	166
	167	PRIVATE :: grid_hgr__fill_curv
	168	PRIVATE :: grid_hgr__fill_reg
	169	PRIVATE :: grid_hgr__fill_plan
	170	PRIVATE :: grid_hgr__fill_merc
	171	PRIVATE :: grid_hgr__fill_gyre
	172	PRIVATE :: grid_hgr__fill_coriolis
	173	PRIVATE :: grid_hgr__angle
	174
	175	TYPE TNAMH
	176
	177	CHARACTER(LEN=lc) :: c_coord
	178	INTEGER(i4) :: i_perio
	179
	180	INTEGER(i4) :: i_mshhgr
	181	REAL(dp) :: d_ppglam0
	182	REAL(dp) :: d_ppgphi0
	183
	184	REAL(dp) :: d_ppe1_deg
	185	REAL(dp) :: d_ppe2_deg
	186	! REAL(dp) :: d_ppe1_m
	187	! REAL(dp) :: d_ppe2_m
	188
	189	! INTEGER(i4) :: i_cla
	190
	191	! CHARACTER(LEN=lc) :: c_cfg
	192	INTEGER(i4) :: i_cfg
	193	LOGICAL :: l_bench
	194
	195	END TYPE
	196
	197	TYPE(TVAR), SAVE :: tg_tmask
	198	TYPE(TVAR), SAVE :: tg_umask
	199	TYPE(TVAR), SAVE :: tg_vmask
	200	TYPE(TVAR), SAVE :: tg_fmask
	201	! TYPE(TVAR), SAVE :: tg_wmask
	202	! TYPE(TVAR), SAVE :: tg_wumask
	203	! TYPE(TVAR), SAVE :: tg_wvmask
	204
	205	TYPE(TVAR), SAVE :: tg_ssmask
	206	! TYPE(TVAR), SAVE :: tg_ssumask
	207	! TYPE(TVAR), SAVE :: tg_ssvmask
	208	! TYPE(TVAR), SAVE :: tg_ssfmask
	209
	210	TYPE(TVAR), SAVE :: tg_glamt
	211	TYPE(TVAR), SAVE :: tg_glamu
	212	TYPE(TVAR), SAVE :: tg_glamv
	213	TYPE(TVAR), SAVE :: tg_glamf
	214
	215	TYPE(TVAR), SAVE :: tg_gphit
	216	TYPE(TVAR), SAVE :: tg_gphiu
	217	TYPE(TVAR), SAVE :: tg_gphiv
	218	TYPE(TVAR), SAVE :: tg_gphif
	219
	220	TYPE(TVAR), SAVE :: tg_e1t
	221	TYPE(TVAR), SAVE :: tg_e1u
	222	TYPE(TVAR), SAVE :: tg_e1v
	223	TYPE(TVAR), SAVE :: tg_e1f
	224
	225	TYPE(TVAR), SAVE :: tg_e2t
	226	TYPE(TVAR), SAVE :: tg_e2u
	227	TYPE(TVAR), SAVE :: tg_e2v
	228	TYPE(TVAR), SAVE :: tg_e2f
	229
	230	TYPE(TVAR), SAVE :: tg_ff_t
	231	TYPE(TVAR), SAVE :: tg_ff_f
	232
	233	TYPE(TVAR), SAVE :: tg_gcost
	234	TYPE(TVAR), SAVE :: tg_gcosu
	235	TYPE(TVAR), SAVE :: tg_gcosv
	236	TYPE(TVAR), SAVE :: tg_gcosf
	237
	238	TYPE(TVAR), SAVE :: tg_gsint
	239	TYPE(TVAR), SAVE :: tg_gsinu
	240	TYPE(TVAR), SAVE :: tg_gsinv
	241	TYPE(TVAR), SAVE :: tg_gsinf
	242
	243	CONTAINS
[12080]	244	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	245	SUBROUTINE grid_hgr_init(jpi, jpj, jpk, ld_domcfg)
[7314]	246	!-------------------------------------------------------------------
[12080]	247	!> @brief This subroutine initialise hgr structure
[7314]	248	!>
	249	!> @author J.Paul
	250	!> @date September, 2015 - Initial version
	251	!>
	252	!> @param[in] jpi
	253	!> @param[in] jpj
	254	!-------------------------------------------------------------------
[12080]	255
[7314]	256	IMPLICIT NONE
[12080]	257
[7314]	258	! Argument
	259	INTEGER(i4), INTENT(IN) :: jpi
	260	INTEGER(i4), INTENT(IN) :: jpj
	261	INTEGER(i4), INTENT(IN) :: jpk
	262	LOGICAL , INTENT(IN) :: ld_domcfg
	263
	264	REAL(dp), DIMENSION(jpi,jpj) :: dl_tmp2D
	265	REAL(dp), DIMENSION(jpi,jpj,jpk) :: dl_tmp3D
	266	! loop indices
	267	!----------------------------------------------------------------
	268	! variable 2D
	269	dl_tmp2D(:,:) =dp_fill_i1
	270
	271	tg_ssmask = var_init('ssmask' ,dl_tmp2D(:,:) , dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	272	! tg_ssumask = var_init('ssumask',dl_tmp2D(:,:) , dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	273	! tg_ssvmask = var_init('ssvmask',dl_tmp2D(:,:) , dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	274	! tg_ssfmask = var_init('ssfmask',dl_tmp2D(:,:) , dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	275
	276	dl_tmp2D(:,:)=dp_fill
	277
	278	tg_glamt = var_init('glamt',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	279	tg_glamu = var_init('glamu',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	280	tg_glamv = var_init('glamv',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	281	tg_glamf = var_init('glamf',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	282
	283	tg_gphit = var_init('gphit',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	284	tg_gphiu = var_init('gphiu',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	285	tg_gphiv = var_init('gphiv',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	286	tg_gphif = var_init('gphif',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	287
	288	tg_e1t = var_init('e1t' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	289	tg_e1u = var_init('e1u' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	290	tg_e1v = var_init('e1v' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	291	tg_e1f = var_init('e1f' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	292
	293	tg_e2t = var_init('e2t' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	294	tg_e2u = var_init('e2u' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	295	tg_e2v = var_init('e2v' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	296	tg_e2f = var_init('e2f' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	297
	298	tg_ff_t = var_init('ff_t' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	299	tg_ff_f = var_init('ff_f' ,dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	300
	301	IF( .NOT. ld_domcfg )THEN
	302	tg_gcost =var_init('gcost',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	303	tg_gcosu =var_init('gcosu',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	304	tg_gcosv =var_init('gcosv',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	305	tg_gcosf =var_init('gcosf',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	306
	307	tg_gsint =var_init('gsint',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	308	tg_gsinu =var_init('gsinu',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	309	tg_gsinv =var_init('gsinv',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	310	tg_gsinf =var_init('gsinf',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
	311	ENDIF
	312
	313	! variable 3D
	314	dl_tmp3D(:,:,:)=dp_fill_i1
	315
	316	tg_tmask = var_init('tmask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	317	tg_umask = var_init('umask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	318	tg_vmask = var_init('vmask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	319	IF( .NOT. ld_domcfg )THEN
	320	tg_fmask = var_init('fmask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	321	ENDIF
	322
	323	! tg_wmask = var_init('wmask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	324	! tg_wumask = var_init('wumask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	325	! tg_wvmask = var_init('wvmask' ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
	326
	327	END SUBROUTINE grid_hgr_init
[12080]	328	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	329	SUBROUTINE grid_hgr_clean(ld_domcfg)
[7314]	330	!-------------------------------------------------------------------
[12080]	331	!> @brief This subroutine clean hgr structure
[7314]	332	!>
	333	!> @author J.Paul
	334	!> @date September, 2015 - Initial version
	335	!>
	336	!-------------------------------------------------------------------
[12080]	337
[7314]	338	IMPLICIT NONE
[12080]	339
[7314]	340	! Argument
	341	LOGICAL , INTENT(IN) :: ld_domcfg
	342
	343	! local variable
	344	! loop indices
	345	!----------------------------------------------------------------
	346	CALL var_clean(tg_ssmask )
	347
	348	CALL var_clean(tg_glamt)
	349	CALL var_clean(tg_glamu)
	350	CALL var_clean(tg_glamv)
	351	CALL var_clean(tg_glamf)
	352
	353	CALL var_clean(tg_gphit)
	354	CALL var_clean(tg_gphiu)
	355	CALL var_clean(tg_gphiv)
	356	CALL var_clean(tg_gphif)
	357
	358	CALL var_clean(tg_e1t )
	359	CALL var_clean(tg_e1u )
	360	CALL var_clean(tg_e1v )
	361	CALL var_clean(tg_e1f )
	362
	363	CALL var_clean(tg_e2t )
	364	CALL var_clean(tg_e2u )
	365	CALL var_clean(tg_e2v )
	366	CALL var_clean(tg_e2f )
	367
	368	CALL var_clean(tg_ff_t )
	369	CALL var_clean(tg_ff_f )
	370
	371	IF( .NOT. ld_domcfg )THEN
	372	CALL var_clean(tg_gcost )
	373	CALL var_clean(tg_gcosu )
	374	CALL var_clean(tg_gcosv )
	375	CALL var_clean(tg_gcosf )
	376
	377	CALL var_clean(tg_gsint )
	378	CALL var_clean(tg_gsinu )
	379	CALL var_clean(tg_gsinv )
	380	CALL var_clean(tg_gsinf )
	381	ENDIF
	382
	383	CALL var_clean(tg_tmask )
	384	CALL var_clean(tg_umask )
	385	CALL var_clean(tg_vmask )
	386	IF( .NOT. ld_domcfg )THEN
	387	CALL var_clean(tg_fmask )
	388	ENDIF
	389	END SUBROUTINE grid_hgr_clean
[12080]	390	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	391	FUNCTION grid_hgr_nam(cd_coord, id_perio, cd_namelist) &
	392	& RESULT (tf_namh)
[7314]	393	!-------------------------------------------------------------------
	394	!> @brief This function initialise hgr namelist structure
	395	!>
	396	!> @author J.Paul
	397	!> @date September, 2015 - Initial version
	398	!>
	399	!> @param[in] cd_coord
	400	!> @param[in] id_perio
	401	!> @param[in] cd_namelist
	402	!> @return hgr namelist structure
	403	!-------------------------------------------------------------------
[12080]	404
[7314]	405	IMPLICIT NONE
[12080]	406
[7314]	407	! Argument
	408	CHARACTER(LEN=*), INTENT(IN) :: cd_coord
	409	INTEGER(i4) , INTENT(IN) :: id_perio
	410	CHARACTER(LEN=*), INTENT(IN) :: cd_namelist
	411
	412	! function
[12080]	413	TYPE(TNAMH) :: tf_namh
[7314]	414
	415	! local variable
	416	INTEGER(i4) :: il_status
	417	INTEGER(i4) :: il_fileid
	418
	419	LOGICAL :: ll_exist
	420
	421	! loop indices
	422	! namelist
	423
	424	! namhgr
	425	INTEGER(i4) :: in_mshhgr = 0
	426	REAL(dp) :: dn_ppglam0 = NF90_FILL_DOUBLE
	427	REAL(dp) :: dn_ppgphi0 = NF90_FILL_DOUBLE
	428	REAL(dp) :: dn_ppe1_deg = NF90_FILL_DOUBLE
	429	REAL(dp) :: dn_ppe2_deg = NF90_FILL_DOUBLE
	430	! REAL(dp) :: dn_ppe1_m = NF90_FILL_DOUBLE
	431	! REAL(dp) :: dn_ppe2_m = NF90_FILL_DOUBLE
	432
	433	! ! namcla
	434	! INTEGER(i4) :: in_cla = 0
	435
	436	! namgrd
	437	! CHARACTER(LEN=lc) :: cn_cfg = ''
	438	INTEGER(i4) :: in_cfg = 0
	439	LOGICAL :: ln_bench = .FALSE.
	440
	441	!----------------------------------------------------------------
	442	NAMELIST /namhgr/ &
	443	& in_mshhgr, & !< type of horizontal mesh
	444	!< 0: curvilinear coordinate on the sphere read in coordinate.nc
	445	!< 1: geographical mesh on the sphere with regular grid-spacing
	446	!< 2: f-plane with regular grid-spacing
	447	!< 3: beta-plane with regular grid-spacing
	448	!< 4: Mercator grid with T/U point at the equator
	449	!< 5: beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
	450	& dn_ppglam0, & !< longitude of first raw and column T-point (in_mshhgr = 1 or 4)
	451	& dn_ppgphi0, & !< latitude of first raw and column T-point (in_mshhgr = 1 or 4)
	452	& dn_ppe1_deg, & !< zonal grid-spacing (degrees) (in_mshhgr = 1,2,3 or 4)
	453	& dn_ppe2_deg !< meridional grid-spacing (degrees) (in_mshhgr = 1,2,3 or 4)
	454	! & dn_ppe1_m, & !< zonal grid-spacing (degrees)
	455	! & dn_ppe2_m !< meridional grid-spacing (degrees)
	456
	457	! NAMELIST /namcla/ &
	458	! & in_cla !< =1 cross land advection for exchanges through some straits (ORCA2)
	459
	460	NAMELIST/namgrd/ & !< orca grid namelist
	461	! & cn_cfg, & !< name of the configuration (orca)
	462	& in_cfg, & !< resolution of the configuration (2,1,025..)
	463	& ln_bench !< benchmark parameter (in_mshhgr = 5 ).
	464
	465	!----------------------------------------------------------------
	466	! read namelist
	467	INQUIRE(FILE=TRIM(cd_namelist), EXIST=ll_exist)
	468	IF( ll_exist )THEN
	469
	470	il_fileid=fct_getunit()
	471
	472	OPEN( il_fileid, FILE=TRIM(cd_namelist), &
[12080]	473	& FORM='FORMATTED', &
	474	& ACCESS='SEQUENTIAL', &
	475	& STATUS='OLD', &
	476	& ACTION='READ', &
	477	& IOSTAT=il_status)
[7314]	478	CALL fct_err(il_status)
	479	IF( il_status /= 0 )THEN
	480	CALL logger_fatal("GRID HGR NAM: error opening "//&
	481	& TRIM(cd_namelist))
	482	ENDIF
	483
	484	READ( il_fileid, NML = namhgr )
	485	! READ( il_fileid, NML = namcla )
	486	! READ( il_fileid, NML = namgrd )
	487
	488	CLOSE( il_fileid, IOSTAT=il_status )
	489	CALL fct_err(il_status)
	490	IF( il_status /= 0 )THEN
	491	CALL logger_error("GRID HGR NAM: closing "//TRIM(cd_namelist))
	492	ENDIF
	493
[12080]	494	tf_namh%c_coord = TRIM(cd_coord)
	495	tf_namh%i_perio = id_perio
[7314]	496
[12080]	497	tf_namh%i_mshhgr = in_mshhgr
	498	tf_namh%d_ppglam0 = dn_ppglam0
	499	tf_namh%d_ppgphi0 = dn_ppgphi0
[7314]	500
[12080]	501	tf_namh%d_ppe1_deg= dn_ppe1_deg
	502	tf_namh%d_ppe2_deg= dn_ppe2_deg
	503	! tf_namh%d_ppe1_m = dn_ppe1_m
	504	! tf_namh%d_ppe2_m = dn_ppe2_m
[7314]	505
[12080]	506	! tf_namh%i_cla = in_cla
[7314]	507
[12080]	508	! tf_namh%c_cfg = TRIM(cn_cfg)
	509	tf_namh%i_cfg = in_cfg
	510	tf_namh%l_bench = ln_bench
[7314]	511
	512	ELSE
	513
	514	CALL logger_fatal(" GRID HGR NAM: can't find "//TRIM(cd_namelist))
	515
	516	ENDIF
	517
	518	END FUNCTION grid_hgr_nam
[12080]	519	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	520	SUBROUTINE grid_hgr_fill(td_nam, jpi, jpj, ld_domcfg)
[7314]	521	!-------------------------------------------------------------------
	522	!> @brief This subroutine fill horizontal mesh (hgr structure)
	523	!>
	524	!> @author J.Paul
	525	!> @date September, 2015 - Initial version
	526	!>
	527	!> @param[in] td_nam
	528	!> @param[in] jpi
	529	!> @param[in] jpj
	530	!-------------------------------------------------------------------
[12080]	531
[7314]	532	IMPLICIT NONE
[12080]	533
[7314]	534	! Argument
	535	TYPE(TNAMH), INTENT(IN) :: td_nam
	536	INTEGER(i4), INTENT(IN) :: jpi
	537	INTEGER(i4), INTENT(IN) :: jpj
	538	LOGICAL , INTENT(IN) :: ld_domcfg
	539
	540	! local variable
	541	REAL(dp) :: znorme
	542	! loop indices
	543	!----------------------------------------------------------------
	544	CALL logger_info('GRID HGR FILL : define the horizontal mesh from the'//&
	545	& ' type of horizontal mesh mshhgr = '//TRIM(fct_str(td_nam%i_mshhgr)))
	546	IF( td_nam%i_mshhgr == 1 )THEN
	547	CALL logger_info(' position of the first row and ppglam0 = '//&
	548	& TRIM(fct_str(td_nam%d_ppglam0 )) )
	549	CALL logger_info(' column grid-point (degrees) ppgphi0 = '//&
	550	& TRIM(fct_str(td_nam%d_ppgphi0 )) )
	551	ELSEIF( td_nam%i_mshhgr == 2 .OR. td_nam%i_mshhgr == 3 )THEN
	552	CALL logger_info(' zonal grid-spacing (degrees) ppe1_deg = '//&
	553	& TRIM(fct_str(td_nam%d_ppe1_deg )) )
	554	CALL logger_info(' meridional grid-spacing (degrees) ppe2_deg = '//&
	555	& TRIM(fct_str(td_nam%d_ppe2_deg )) )
	556	! CALL logger_info(' zonal grid-spacing (meters) ppe1_m = '//&
	557	! & TRIM(fct_str(td_nam%d_ppe1_m )) )
	558	! CALL logger_info(' meridional grid-spacing (meters) ppe2_m = '//&
	559	! & TRIM(fct_str(td_nam%d_ppe2_m )) )
	560	ENDIF
	561
	562	SELECT CASE( td_nam%i_mshhgr ) ! type of horizontal mesh
	563
	564	CASE(0) ! curvilinear coordinate on the sphere read in coordinate.nc file
	565
	566	CALL grid_hgr__fill_curv(td_nam)!,jpi,jpj)
	567
	568	CASE(1) ! geographical mesh on the sphere with regular grid-spacing
	569
	570	CALL grid_hgr__fill_reg(td_nam,jpi,jpj)
	571
	572	CASE(2:3) ! f- or beta-plane with regular grid-spacing
	573
	574	CALL grid_hgr__fill_plan(td_nam,jpi,jpj)
	575
	576	CASE(4) ! geographical mesh on the sphere, isotropic MERCATOR type
	577
	578	CALL grid_hgr__fill_merc(td_nam,jpi,jpj)
	579
	580	CASE(5) ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
	581
	582	CALL grid_hgr__fill_gyre(td_nam,jpi,jpj)
	583
	584	CASE DEFAULT
	585
	586	CALL logger_fatal('GRID HGR FILL : bad flag value for mshhgr = '//&
	587	& TRIM(fct_str(td_nam%i_mshhgr)))
	588
	589	END SELECT
	590
	591	! No Useful associated horizontal metrics
	592	! ---------------------------------------
	593
	594	! create coriolis factor
	595	CALL grid_hgr__fill_coriolis(td_nam,jpi)!,jpj)
	596
	597	! Control of domain for symetrical condition
	598	! ------------------------------------------
	599	! The equator line must be the latitude coordinate axe
	600
	601	IF( td_nam%i_perio == 2 ) THEN
	602	znorme = SQRT( SUM(tg_gphiu%d_value(:,2,1,1)*tg_gphiu%d_value(:,2,1,1)) ) / FLOAT( jpi )
	603	IF( znorme > 1.e-13 )THEN
	604	CALL logger_fatal( ' ===>>>> : symmetrical condition: rerun with good equator line' )
	605	ENDIF
	606	ENDIF
	607
	608	! compute angles between model grid lines and the North direction
	609	! ---------------------------------------------------------------
	610	IF( .NOT. ld_domcfg )THEN
	611	CALL grid_hgr__angle(td_nam,jpi,jpj)
	612	ENDIF
	613
	614	END SUBROUTINE grid_hgr_fill
[12080]	615	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	616	SUBROUTINE grid_hgr__fill_curv(td_nam)!,jpi,jpj)
[7314]	617	!-------------------------------------------------------------------
	618	!> @brief This subroutine fill horizontal mesh (hgr structure)
	619	!> for case of curvilinear coordinate on the sphere read in coordinate.nc file
	620	!>
	621	!> @author J.Paul
	622	!> @date September, 2015 - Initial version
	623	!> @date October, 2016
	624	!> - do not use anymore special case for ORCA grid
	625	!>
	626	!> @param[in] td_nam
	627	! @param[in] jpi
	628	! @param[in] jpj
	629	!-------------------------------------------------------------------
[12080]	630
[7314]	631	IMPLICIT NONE
[12080]	632
[7314]	633	! Argument
	634	TYPE(TNAMH), INTENT(IN) :: td_nam
	635	! INTEGER(i4), INTENT(IN) :: jpi
	636	! INTEGER(i4), INTENT(IN) :: jpj
	637
	638	! local variable
	639	! INTEGER(i4) :: ii0, ii1, ij0, ij1 ! temporary integers
	640	! INTEGER(i4) :: isrow ! index for ORCA1 starting row
	641
	642	TYPE(TMPP) :: tl_coord
	643
	644	! loop indices
	645	!----------------------------------------------------------------
	646
	647	! read coordinates
	648	! open file
	649	IF( td_nam%c_coord /= '' )THEN
	650	tl_coord=mpp_init( file_init(TRIM(td_nam%c_coord)), id_perio=td_nam%i_perio)
	651	CALL grid_get_info(tl_coord)
	652	ELSE
	653	CALL logger_fatal("GRID HGR FILL: no input coordinates file found. "//&
	654	& "check namelist")
	655	ENDIF
	656
	657	CALL iom_mpp_open( tl_coord )
	658
	659	! read variable in coordinates
	660	tg_glamt=iom_mpp_read_var(tl_coord, 'glamt')
	661	tg_glamu=iom_mpp_read_var(tl_coord, 'glamu')
	662	tg_glamv=iom_mpp_read_var(tl_coord, 'glamv')
	663	tg_glamf=iom_mpp_read_var(tl_coord, 'glamf')
	664
	665	tg_gphit=iom_mpp_read_var(tl_coord, 'gphit')
	666	tg_gphiu=iom_mpp_read_var(tl_coord, 'gphiu')
	667	tg_gphiv=iom_mpp_read_var(tl_coord, 'gphiv')
	668	tg_gphif=iom_mpp_read_var(tl_coord, 'gphif')
	669
	670	! force output type
	671	tg_glamt%i_type=NF90_DOUBLE
	672	tg_glamu%i_type=NF90_DOUBLE
	673	tg_glamv%i_type=NF90_DOUBLE
	674	tg_glamf%i_type=NF90_DOUBLE
	675
	676	tg_gphit%i_type=NF90_DOUBLE
	677	tg_gphiu%i_type=NF90_DOUBLE
	678	tg_gphiv%i_type=NF90_DOUBLE
	679	tg_gphif%i_type=NF90_DOUBLE
	680
	681	tg_e1t =iom_mpp_read_var(tl_coord, 'e1t')
	682	tg_e1u =iom_mpp_read_var(tl_coord, 'e1u')
	683	tg_e1v =iom_mpp_read_var(tl_coord, 'e1v')
	684	tg_e1f =iom_mpp_read_var(tl_coord, 'e1f')
	685
	686	tg_e2t =iom_mpp_read_var(tl_coord, 'e2t')
	687	tg_e2u =iom_mpp_read_var(tl_coord, 'e2u')
	688	tg_e2v =iom_mpp_read_var(tl_coord, 'e2v')
	689	tg_e2f =iom_mpp_read_var(tl_coord, 'e2f')
	690
	691	CALL iom_mpp_close( tl_coord )
	692	! clean
	693	CALL mpp_clean(tl_coord)
	694
	695	!! WARNING extended grid have to be correctly fill
	696
	697	! !! special case for ORCA grid
	698	! ! ORCA R2 configuration
	699	! IF( TRIM(td_nam%c_cfg) == "orca" .AND. td_nam%i_cfg == 2 ) THEN
	700	! IF( td_nam%i_cla == 0 ) THEN
	701	! !
	702	! ! Gibraltar Strait (e2u = 20 km)
	703	! ii0 = 139 ; ii1 = 140
	704	! ij0 = 102 ; ij1 = 102
	705	! ! e2u = 20 km
	706	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 20.e3
	707	! CALL logger_info('orca_r2: Gibraltar : e2u reduced to 20 km')
	708	! !
	709	! ! Bab el Mandeb (e2u = 18 km)
	710	! ii0 = 160 ; ii1 = 160
	711	! ij0 = 88 ; ij1 = 88
	712	! ! e1v = 18 km
	713	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 18.e3
	714	! ! e2u = 30 km
	715	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 30.e3
	716	!
	717	! CALL logger_info('orca_r2: Bab el Mandeb: e2u reduced to 30 km')
	718	! CALL logger_info('e1v reduced to 18 km')
	719	! ENDIF
	720	! ! Danish Straits
	721	! ii0 = 145 ; ii1 = 146
	722	! ij0 = 116 ; ij1 = 116
	723	! ! e2u = 10 km
	724	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 10.e3
	725	! CALL logger_info('orca_r2: Danish Straits : e2u reduced to 10 km')
	726	! ENDIF
	727	!
	728	! ! ORCA R1 configuration
	729	! IF( TRIM(td_nam%c_cfg) == "orca" .AND. td_nam%i_cfg == 1 ) THEN
	730	! ! This dirty section will be suppressed by simplification process: all this will come back in input files
	731	! ! Currently these hard-wired indices relate to configuration with
	732	! ! extend grid (jpjglo=332)
	733	! ! which had a grid-size of 362x292.
	734	!
	735	! isrow = 332 - jpj
	736	!
	737	! ! Gibraltar Strait (e2u = 20 km)
	738	! ii0 = 282 ; ii1 = 283
	739	! ij0 = 201 + isrow ; ij1 = 241 - isrow
	740	! ! e2u = 20 km
	741	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 20.e3
	742	! CALL logger_info('orca_r1: Gibraltar : e2u reduced to 20 km')
	743	!
	744	! ! Bhosporus Strait (e2u = 10 km)
	745	! ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u = 10 km)
	746	! ij0 = 208 + isrow ; ij1 = 248 - isrow
	747	! ! Bhosporus Strait (e2u = 10 km)
	748	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 10.e3
	749	! CALL logger_info('orca_r1: Bhosporus : e2u reduced to 10 km')
	750	!
	751	! ! Lombok Strait (e1v = 13 km)
	752	! ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v = 13 km)
	753	! ij0 = 124 + isrow ; ij1 = 165 - isrow
	754	! ! Lombok Strait (e1v = 13 km)
	755	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 13.e3
	756	! CALL logger_info('orca_r1: Lombok : e1v reduced to 10 km')
	757	!
	758	! ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
	759	! ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
	760	! ij0 = 124 + isrow ; ij1 = 165 - isrow
	761	! ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
	762	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 8.e3
	763	! CALL logger_info('orca_r1: Sumba : e1v reduced to 8 km')
	764	!
	765	! ! Ombai Strait (e1v = 13 km)
	766	! ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v = 13 km)
	767	! ij0 = 124 + isrow ; ij1 = 165 - isrow
	768	! ! Ombai Strait (e1v = 13 km)
	769	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 13.e3
	770	! CALL logger_info('orca_r1: Ombai : e1v reduced to 13 km')
	771	!
	772	! ! Timor Passage (e1v = 20 km)
	773	! ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v = 20 km)
	774	! ij0 = 124 + isrow ; ij1 = 145 - isrow
	775	! ! Timor Passage (e1v = 20 km)
	776	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 20.e3
	777	! CALL logger_info('orca_r1: Timor Passage : e1v reduced to 20 km')
	778	!
	779	! ! West Halmahera Strait (e1v = 30 km)
	780	! ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v = 30 km)
	781	! ij0 = 141 + isrow ; ij1 = 182 - isrow
	782	! ! West Halmahera Strait (e1v = 30 km)
	783	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 30.e3
	784	! CALL logger_info('orca_r1: W Halmahera : e1v reduced to 30 km')
	785	!
	786	! ! East Halmahera Strait (e1v = 50 km)
	787	! ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v = 50 km)
	788	! ij0 = 141 + isrow ; ij1 = 182 - isrow
	789	! ! East Halmahera Strait (e1v = 50 km)
	790	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 50.e3
	791	! CALL logger_info('orca_r1: E Halmahera : e1v reduced to 50 km')
	792	!
	793	! ENDIF
	794	!
	795	! ! ORCA R05 configuration
	796	! IF( TRIM(td_nam%c_cfg) == "orca" .AND. td_nam%i_cfg == 05 ) THEN
	797	!
	798	! ! Gibraltar Strait (e2u = 20 km)
	799	! ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u = 20 km)
	800	! ij0 = 327 ; ij1 = 327
	801	! ! Gibraltar Strait (e2u = 20 km)
	802	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 20.e3
	803	! CALL logger_info('orca_r05: Reduced e2u at the Gibraltar Strait')
	804	! !
	805	! ! Bosphore Strait (e2u = 10 km)
	806	! ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u = 10 km)
	807	! ij0 = 343 ; ij1 = 343
	808	! ! Bosphore Strait (e2u = 10 km)
	809	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 10.e3
	810	! CALL logger_info('orca_r05: Reduced e2u at the Bosphore Strait')
	811	! !
	812	! ! Sumba Strait (e2u = 40 km)
	813	! ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u = 40 km)
	814	! ij0 = 232 ; ij1 = 232
	815	! ! Sumba Strait (e2u = 40 km)
	816	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 40.e3
	817	! CALL logger_info('orca_r05: Reduced e2u at the Sumba Strait')
	818	! !
	819	! ! Ombai Strait (e2u = 15 km)
	820	! ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u = 15 km)
	821	! ij0 = 232 ; ij1 = 232
	822	! ! Ombai Strait (e2u = 15 km)
	823	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 15.e3
	824	! CALL logger_info('orca_r05: Reduced e2u at the Ombai Strait')
	825	! !
	826	! ! Palk Strait (e2u = 10 km)
	827	! ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u = 10 km)
	828	! ij0 = 270 ; ij1 = 270
	829	! ! Palk Strait (e2u = 10 km)
	830	! tg_e2u%d_value(ii0:ii1,ij0:ij1,1,1) = 10.e3
	831	! CALL logger_info('orca_r05: Reduced e2u at the Palk Strait')
	832	! !
	833	! ! Lombok Strait (e1v = 10 km)
	834	! ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v = 10 km)
	835	! ij0 = 232 ; ij1 = 233
	836	! ! Lombok Strait (e1v = 10 km)
	837	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 10.e3
	838	! CALL logger_info('orca_r05: Reduced e1v at the Lombok Strait')
	839	! !
	840	! !
	841	! ! Bab el Mandeb (e1v = 25 km)
	842	! ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v = 25 km)
	843	! ij0 = 276 ; ij1 = 276
	844	! ! Bab el Mandeb (e1v = 25 km)
	845	! tg_e1v%d_value(ii0:ii1,ij0:ij1,1,1) = 25.e3
	846	! CALL logger_info('orca_r05: Reduced e1v at the Bab el Mandeb')
	847	!
	848	! ENDIF
	849
	850	END SUBROUTINE grid_hgr__fill_curv
[12080]	851	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	852	SUBROUTINE grid_hgr__fill_reg(td_nam, jpi, jpj)
[7314]	853	!-------------------------------------------------------------------
	854	!> @brief This subroutine fill horizontal mesh (hgr structure)
	855	!> for case of geographical mesh on the sphere with regular grid-spacing
	856	!>
	857	!> @author J.Paul
	858	!> @date September, 2015 - Initial version
	859	!>
	860	!> @param[in] td_nam
	861	!> @param[in] jpi
	862	!> @param[in] jpj
	863	!-------------------------------------------------------------------
[12080]	864
[7314]	865	IMPLICIT NONE
[12080]	866
[7314]	867	! Argument
	868	TYPE(TNAMH), INTENT(IN) :: td_nam
	869	INTEGER(i4), INTENT(IN) :: jpi
	870	INTEGER(i4), INTENT(IN) :: jpj
	871
	872	! local variable
	873	REAL(dp) :: zti, zui, zvi, zfi ! local scalars
	874	REAL(dp) :: ztj, zuj, zvj, zfj !
	875
	876	! loop indices
	877	INTEGER(i4) :: ji
	878	INTEGER(i4) :: jj
	879	!----------------------------------------------------------------
	880
	881	CALL logger_info('GRID HGR FILL : geographical mesh on the sphere with'//&
	882	& ' regular grid-spacing given by ppe1_deg and ppe2_deg')
	883
	884	DO jj = 1, jpj
	885	DO ji = 1, jpi
	886	zti = FLOAT( ji - 1 ) ; ztj = FLOAT( jj - 1 )
	887	zui = FLOAT( ji - 1 ) + 0.5 ; zuj = FLOAT( jj - 1 )
	888	zvi = FLOAT( ji - 1 ) ; zvj = FLOAT( jj - 1 ) + 0.5
	889	zfi = FLOAT( ji - 1 ) + 0.5 ; zfj = FLOAT( jj - 1 ) + 0.5
	890	! Longitude
	891	tg_glamt%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zti
	892	tg_glamu%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zui
	893	tg_glamv%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zvi
	894	tg_glamf%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zfi
	895	! Latitude
	896	tg_gphit%d_value(ji,jj,1,1) = td_nam%d_ppgphi0 + td_nam%d_ppe2_deg * ztj
	897	tg_gphiu%d_value(ji,jj,1,1) = td_nam%d_ppgphi0 + td_nam%d_ppe2_deg * zuj
	898	tg_gphiv%d_value(ji,jj,1,1) = td_nam%d_ppgphi0 + td_nam%d_ppe2_deg * zvj
	899	tg_gphif%d_value(ji,jj,1,1) = td_nam%d_ppgphi0 + td_nam%d_ppe2_deg * zfj
	900	! e1
	901	tg_e1t%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphit%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	902	tg_e1u%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphiu%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	903	tg_e1v%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphiv%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	904	tg_e1f%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphif%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	905	! e2
	906	tg_e2t%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * td_nam%d_ppe2_deg
	907	tg_e2u%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * td_nam%d_ppe2_deg
	908	tg_e2v%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * td_nam%d_ppe2_deg
	909	tg_e2f%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * td_nam%d_ppe2_deg
	910	END DO
	911	END DO
	912
	913	END SUBROUTINE grid_hgr__fill_reg
[12080]	914	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	915	SUBROUTINE grid_hgr__fill_plan(td_nam, jpi, jpj)
[7314]	916	!-------------------------------------------------------------------
	917	!> @brief This subroutine fill horizontal mesh (hgr structure)
	918	!> for case of f- or beta-plane with regular grid-spacing
	919	!>
	920	!> @author J.Paul
	921	!> @date September, 2015 - Initial version
	922	!>
	923	!> @param[in] td_nam
	924	!> @param[in] jpi
	925	!> @param[in] jpj
	926	!-------------------------------------------------------------------
[12080]	927
[7314]	928	IMPLICIT NONE
[12080]	929
[7314]	930	! Argument
	931	TYPE(TNAMH), INTENT(IN) :: td_nam
	932	INTEGER(i4), INTENT(IN) :: jpi
	933	INTEGER(i4), INTENT(IN) :: jpj
	934
	935	! local variable
	936	REAL(dp) :: dl_glam0
	937	REAL(dp) :: dl_gphi0
	938
	939	! loop indices
	940	INTEGER(i4) :: ji
	941	INTEGER(i4) :: jj
	942	!----------------------------------------------------------------
	943
	944	CALL logger_info('GRID HGR FILL : f- or beta-plane with regular'//&
	945	& ' grid-spacing given by ppe1_deg and ppe2_deg')
	946	! & ' grid-spacing given by ppe1_m and ppe2_m')
	947
	948	! Position coordinates (in kilometers)
	949	! ==========
	950	dl_glam0 = 0.e0
	951	dl_gphi0 = - td_nam%d_ppe2_deg * 1.e-3
	952	! dl_gphi0 = - td_nam%d_ppe2_m * 1.e-3
	953
	954	!
	955	DO jj = 1, jpj
	956	DO ji = 1, jpi
	957	! tg_glamt%d_value(ji,jj,1,1) = dl_glam0 + td_nam%d_ppe1_m * 1.e-3 * ( FLOAT( ji - 1 ) )
	958	! tg_glamu%d_value(ji,jj,1,1) = dl_glam0 + td_nam%d_ppe1_m * 1.e-3 * ( FLOAT( ji - 1 ) + 0.5 )
	959	tg_glamt%d_value(ji,jj,1,1) = dl_glam0 + td_nam%d_ppe1_deg * 1.e-3 * ( FLOAT( ji - 1 ) )
	960	tg_glamu%d_value(ji,jj,1,1) = dl_glam0 + td_nam%d_ppe1_deg * 1.e-3 * ( FLOAT( ji - 1 ) + 0.5 )
	961	tg_glamv%d_value(ji,jj,1,1) = tg_glamt%d_value(ji,jj,1,1)
	962	tg_glamf%d_value(ji,jj,1,1) = tg_glamu%d_value(ji,jj,1,1)
	963
	964	!tg_gphit%d_value(ji,jj,1,1) = dl_gphi0 + td_nam%d_ppe2_m * 1.e-3 * ( FLOAT( jj - 1 ) )
	965	tg_gphit%d_value(ji,jj,1,1) = dl_gphi0 + td_nam%d_ppe2_deg * 1.e-3 * ( FLOAT( jj - 1 ) )
	966	tg_gphiu%d_value(ji,jj,1,1) = tg_gphit%d_value(ji,jj,1,1)
	967	!tg_gphiv%d_value(ji,jj,1,1) = dl_gphi0 + td_nam%d_ppe2_m * 1.e-3 * ( FLOAT( jj - 1 ) + 0.5 )
	968	tg_gphiv%d_value(ji,jj,1,1) = dl_gphi0 + td_nam%d_ppe2_deg * 1.e-3 * ( FLOAT( jj - 1 ) + 0.5 )
	969	tg_gphif%d_value(ji,jj,1,1) = tg_gphiv%d_value(ji,jj,1,1)
	970	END DO
	971	END DO
	972
	973	! Horizontal scale factors (in meters)
	974	! ======
	975	! tg_e1t%d_value(:,:,1,1) = td_nam%d_ppe1_m
	976	! tg_e1u%d_value(:,:,1,1) = td_nam%d_ppe1_m
	977	! tg_e1v%d_value(:,:,1,1) = td_nam%d_ppe1_m
	978	! tg_e1f%d_value(:,:,1,1) = td_nam%d_ppe1_m
	979	tg_e1t%d_value(:,:,1,1) = td_nam%d_ppe1_deg
	980	tg_e1u%d_value(:,:,1,1) = td_nam%d_ppe1_deg
	981	tg_e1v%d_value(:,:,1,1) = td_nam%d_ppe1_deg
	982	tg_e1f%d_value(:,:,1,1) = td_nam%d_ppe1_deg
	983
	984	! tg_e2t%d_value(:,:,1,1) = td_nam%d_ppe2_m
	985	! tg_e2u%d_value(:,:,1,1) = td_nam%d_ppe2_m
	986	! tg_e2v%d_value(:,:,1,1) = td_nam%d_ppe2_m
	987	! tg_e2f%d_value(:,:,1,1) = td_nam%d_ppe2_m
	988	tg_e2t%d_value(:,:,1,1) = td_nam%d_ppe2_deg
	989	tg_e2u%d_value(:,:,1,1) = td_nam%d_ppe2_deg
	990	tg_e2v%d_value(:,:,1,1) = td_nam%d_ppe2_deg
	991	tg_e2f%d_value(:,:,1,1) = td_nam%d_ppe2_deg
	992
	993	END SUBROUTINE grid_hgr__fill_plan
[12080]	994	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	995	SUBROUTINE grid_hgr__fill_merc(td_nam, jpi, jpj)
[7314]	996	!-------------------------------------------------------------------
	997	!> @brief This subroutine fill horizontal mesh (hgr structure)
	998	!> for case of geographical mesh on the sphere, isotropic MERCATOR type
	999	!>
	1000	!> @author J.Paul
	1001	!> @date September, 2015 - Initial version
	1002	!>
	1003	!> @param[in] td_nam
	1004	!> @param[in] jpi
	1005	!> @param[in] jpj
	1006	!-------------------------------------------------------------------
[12080]	1007
[7314]	1008	IMPLICIT NONE
[12080]	1009
[7314]	1010	! Argument
	1011	TYPE(TNAMH), INTENT(IN) :: td_nam
	1012	INTEGER(i4), INTENT(IN) :: jpi
	1013	INTEGER(i4), INTENT(IN) :: jpj
	1014
	1015	! local variable
	1016	INTEGER :: ijeq ! index of equator T point (used in case 4)
	1017
	1018	REAL(dp) :: zti, zui, zvi, zfi ! local scalars
	1019	REAL(dp) :: ztj, zuj, zvj, zfj !
	1020	REAL(dp) :: zarg
	1021
	1022	! loop indices
	1023	INTEGER(i4) :: ji
	1024	INTEGER(i4) :: jj
	1025	!----------------------------------------------------------------
	1026
	1027	CALL logger_info('GRID HGR FILL : geographical mesh on the sphere, '//&
	1028	& 'MERCATOR type longitudinal/latitudinal spacing given by ppe1_deg')
	1029
	1030	IF( td_nam%d_ppgphi0 == -90 )THEN
	1031	CALL logger_fatal(' Mercator grid cannot start at south pole !!!! ' )
	1032	ENDIF
	1033
	1034	! Find index corresponding to the equator, given the grid spacing e1_deg
	1035	! and the (approximate) southern latitude ppgphi0.
	1036	! This way we ensure that the equator is at a "T / U" point, when in the domain.
	1037	! The formula should work even if the equator is outside the domain.
	1038	zarg = dp_pi / 4. - dp_pi / 180. * td_nam%d_ppgphi0 / 2.
	1039	ijeq = ABS( 180./dp_pi * LOG( COS( zarg ) / SIN( zarg ) ) / td_nam%d_ppe1_deg )
	1040	IF( td_nam%d_ppgphi0 > 0 ) ijeq = -ijeq
	1041
	1042	CALL logger_info('Index of the equator on the MERCATOR grid: '//TRIM(fct_str(ijeq)))
	1043
	1044	DO jj = 1, jpj
	1045	DO ji = 1, jpi
	1046	zti = FLOAT( ji - 1 ) ; ztj = FLOAT( jj - ijeq )
	1047	zui = FLOAT( ji - 1 ) + 0.5 ; zuj = FLOAT( jj - ijeq )
	1048	zvi = FLOAT( ji - 1 ) ; zvj = FLOAT( jj - ijeq ) + 0.5
	1049	zfi = FLOAT( ji - 1 ) + 0.5 ; zfj = FLOAT( jj - ijeq ) + 0.5
	1050	! Longitude
	1051	tg_glamt%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zti
	1052	tg_glamu%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zui
	1053	tg_glamv%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zvi
	1054	tg_glamf%d_value(ji,jj,1,1) = td_nam%d_ppglam0 + td_nam%d_ppe1_deg * zfi
	1055	! Latitude
	1056	tg_gphit%d_value(ji,jj,1,1) = 1./dp_deg2rad * ASIN ( TANH( td_nam%d_ppe1_deg dp_deg2rad ztj ) )
	1057	tg_gphiu%d_value(ji,jj,1,1) = 1./dp_deg2rad * ASIN ( TANH( td_nam%d_ppe1_deg dp_deg2rad zuj ) )
	1058	tg_gphiv%d_value(ji,jj,1,1) = 1./dp_deg2rad * ASIN ( TANH( td_nam%d_ppe1_deg dp_deg2rad zvj ) )
	1059	tg_gphif%d_value(ji,jj,1,1) = 1./dp_deg2rad * ASIN ( TANH( td_nam%d_ppe1_deg dp_deg2rad zfj ) )
	1060	! e1
	1061	tg_e1t%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphit%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1062	tg_e1u%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphiu%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1063	tg_e1v%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphiv%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1064	tg_e1f%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphif%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1065	! e2
	1066	tg_e2t%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphit%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1067	tg_e2u%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphiu%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1068	tg_e2v%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphiv%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1069	tg_e2f%d_value(ji,jj,1,1) = dp_rearth * dp_deg2rad * COS( dp_deg2rad * tg_gphif%d_value(ji,jj,1,1) ) * td_nam%d_ppe1_deg
	1070	END DO
	1071	END DO
	1072
	1073	END SUBROUTINE grid_hgr__fill_merc
[12080]	1074	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1075	SUBROUTINE grid_hgr__fill_gyre(td_nam, jpi, jpj)
[7314]	1076	!-------------------------------------------------------------------
	1077	!> @brief This subroutine fill horizontal mesh (hgr structure)
	1078	!> for case of beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
	1079	!>
	1080	!> @author J.Paul
	1081	!> @date September, 2015 - Initial version
	1082	!>
	1083	!> @param[in] td_nam
	1084	!> @param[in] jpi
	1085	!> @param[in] jpj
	1086	!-------------------------------------------------------------------
[12080]	1087
[7314]	1088	IMPLICIT NONE
[12080]	1089
[7314]	1090	! Argument
	1091	TYPE(TNAMH), INTENT(IN) :: td_nam
	1092	INTEGER(i4), INTENT(IN) :: jpi
	1093	INTEGER(i4), INTENT(IN) :: jpj
	1094
	1095	! local variable
	1096	REAL(dp) :: zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
	1097	REAL(dp) :: zphi1, zsin_alpha, zim05, zjm05
	1098
	1099	REAL(dp) :: dl_glam0
	1100	REAL(dp) :: dl_gphi0
	1101
	1102	! loop indices
	1103	INTEGER(i4) :: ji
	1104	INTEGER(i4) :: jj
	1105	!----------------------------------------------------------------
	1106
	1107	CALL logger_info('GRID HGR FILL : beta-plane with regular grid-spacing '//&
	1108	& 'and rotated domain (GYRE configuration)')
	1109
	1110	! Position coordinates (in kilometers)
	1111	!
	1112	! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
	1113	zlam1 = -85
	1114	zphi1 = 29
	1115	! resolution in meters
	1116	ze1 = 106000. / FLOAT(td_nam%i_cfg)
	1117	! benchmark: forced the resolution to be about 100 km
	1118	IF( td_nam%l_bench ) ze1 = 106000.e0
	1119	zsin_alpha = - SQRT( 2. ) / 2.
	1120	zcos_alpha = SQRT( 2. ) / 2.
	1121	ze1deg = ze1 / (dp_rearth * dp_deg2rad)
	1122	dl_glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpj-2 )
	1123	dl_gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpj-2 )
	1124
	1125	DO jj = 1, jpj
	1126	DO ji = 1, jpi
	1127	zim1 = FLOAT( ji - 1 ) ; zim05 = FLOAT( ji ) - 1.5
	1128	zjm1 = FLOAT( jj - 1 ) ; zjm05 = FLOAT( jj ) - 1.5
	1129
	1130	tg_glamf%d_value(ji,jj,1,1) = dl_glam0 &
	1131	& + zim1 * ze1deg * zcos_alpha &
	1132	& + zjm1 * ze1deg * zsin_alpha
	1133	tg_gphif%d_value(ji,jj,1,1) = dl_gphi0 &
	1134	& - zim1 * ze1deg * zsin_alpha &
	1135	& + zjm1 * ze1deg * zcos_alpha
	1136
	1137	tg_glamt%d_value(ji,jj,1,1) = dl_glam0 &
	1138	& + zim05 * ze1deg * zcos_alpha &
	1139	& + zjm05 * ze1deg * zsin_alpha
	1140	tg_gphit%d_value(ji,jj,1,1) = dl_gphi0 &
	1141	& - zim05 * ze1deg * zsin_alpha &
	1142	& + zjm05 * ze1deg * zcos_alpha
	1143
	1144	tg_glamu%d_value(ji,jj,1,1) = dl_glam0 &
	1145	& + zim1 * ze1deg * zcos_alpha &
	1146	& + zjm05 * ze1deg * zsin_alpha
	1147	tg_gphiu%d_value(ji,jj,1,1) = dl_gphi0 &
	1148	& - zim1 * ze1deg * zsin_alpha &
	1149	& + zjm05 * ze1deg * zcos_alpha
	1150
	1151	tg_glamv%d_value(ji,jj,1,1) = dl_glam0 &
	1152	& + zim05 * ze1deg * zcos_alpha &
	1153	& + zjm1 * ze1deg * zsin_alpha
	1154	tg_gphiv%d_value(ji,jj,1,1) = dl_gphi0 &
	1155	& - zim05 * ze1deg * zsin_alpha &
	1156	& + zjm1 * ze1deg * zcos_alpha
	1157
	1158	END DO
	1159	END DO
	1160
	1161	! Horizontal scale factors (in meters)
	1162	! ======
	1163	tg_e1t%d_value(:,:,1,1) = ze1
	1164	tg_e1u%d_value(:,:,1,1) = ze1
	1165	tg_e1v%d_value(:,:,1,1) = ze1
	1166	tg_e1f%d_value(:,:,1,1) = ze1
	1167
	1168	tg_e2t%d_value(:,:,1,1) = ze1
	1169	tg_e2u%d_value(:,:,1,1) = ze1
	1170	tg_e2v%d_value(:,:,1,1) = ze1
	1171	tg_e2f%d_value(:,:,1,1) = ze1
	1172
	1173	END SUBROUTINE grid_hgr__fill_gyre
[12080]	1174	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1175	SUBROUTINE grid_hgr__fill_coriolis(td_nam, jpi)!,jpj)
[7314]	1176	!-------------------------------------------------------------------
	1177	!> @brief This subroutine fill coriolis factor
	1178	!>
	1179	!> @author J.Paul
	1180	!> @date September, 2015 - Initial version
	1181	!> @date October, 2016
	1182	!> - compute coriolis factor at f-point and at t-point
	1183	!>
	1184	!> @param[in] td_nam
	1185	!> @param[in] jpi
	1186	! @param[in] jpj
	1187	!-------------------------------------------------------------------
[12080]	1188
[7314]	1189	IMPLICIT NONE
[12080]	1190
[7314]	1191	! Argument
	1192	TYPE(TNAMH), INTENT(IN) :: td_nam
	1193	INTEGER(i4), INTENT(IN) :: jpi
	1194	! INTEGER(i4), INTENT(IN) :: jpj
	1195
	1196	! local variable
	1197	REAL(dp) :: zbeta
	1198	REAL(dp) :: zphi0
	1199	REAL(dp) :: zf0
	1200
	1201	! loop indices
	1202	!----------------------------------------------------------------
	1203
	1204	! Coriolis factor
	1205	SELECT CASE( td_nam%i_mshhgr ) ! type of horizontal mesh
	1206
	1207	CASE ( 0, 1, 4 ) ! mesh on the sphere
	1208
	1209	tg_ff_f%d_value(:,:,1,1) = 2. * dp_omega * SIN(dp_deg2rad * tg_gphif%d_value(:,:,1,1))
	1210	tg_ff_t%d_value(:,:,1,1) = 2. * dp_omega * SIN(dp_deg2rad * tg_gphit%d_value(:,:,1,1)) ! at t-point
	1211
	1212	CASE ( 2 ) ! f-plane at ppgphi0
	1213
	1214	tg_ff_f%d_value(:,:,1,1) = 2. * dp_omega * SIN( dp_deg2rad * td_nam%d_ppgphi0 )
	1215	tg_ff_t%d_value(:,:,1,1) = 2. * dp_omega * SIN( dp_deg2rad * td_nam%d_ppgphi0 )
	1216	CALL logger_info('f-plane: Coriolis parameter = constant = '//&
	1217	& TRIM(fct_str(tg_ff_f%d_value(1,1,1,1))) )
	1218
	1219	CASE ( 3 ) ! beta-plane
	1220
	1221	! beta at latitude ppgphi0
	1222	zbeta = 2. * dp_omega * COS( dp_deg2rad * td_nam%d_ppgphi0 ) / dp_rearth
	1223	! latitude of the first row F-points
	1224	! zphi0 = td_nam%d_ppgphi0 - FLOAT( jpi/2 ) * td_nam%d_ppe2_m / ( dp_rearth * dp_deg2rad )
	1225	zphi0 = td_nam%d_ppgphi0 - FLOAT( jpi/2 ) * td_nam%d_ppe2_deg / ( dp_rearth * dp_deg2rad )
	1226
	1227	! compute f0 1st point south
	1228	zf0 = 2. * dp_omega * SIN( dp_deg2rad * zphi0 )
	1229	! f = f0 +beta* y ( y=0 at south)
	1230	tg_ff_f%d_value(:,:,1,1) = zf0 + zbeta * tg_gphif%d_value(:,:,1,1) * 1.e3
	1231	tg_ff_t%d_value(:,:,1,1) = zf0 + zbeta * tg_gphit%d_value(:,:,1,1) * 1.e3
	1232
	1233	CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration)
	1234
	1235	! beta at latitude ppgphi0
	1236	zbeta = 2. * dp_omega * COS( dp_deg2rad * td_nam%d_ppgphi0 ) / dp_rearth
	1237	! latitude of the first row F-points
	1238	zphi0 = 15.e0
	1239	! compute f0 1st point south
	1240	zf0 = 2. * dp_omega * SIN( dp_deg2rad * zphi0 )
	1241
	1242	! f = f0 +beta* y ( y=0 at south)
	1243	tg_ff_f%d_value(:,:,1,1) = ( zf0 + zbeta * ABS( tg_gphif%d_value(:,:,1,1) - zphi0 ) &
	1244	& * dp_deg2rad * dp_rearth )
	1245	tg_ff_t%d_value(:,:,1,1) = ( zf0 + zbeta * ABS( tg_gphit%d_value(:,:,1,1) - zphi0 ) &
	1246	& * dp_deg2rad * dp_rearth )
	1247
	1248	END SELECT
	1249
	1250	END SUBROUTINE grid_hgr__fill_coriolis
[12080]	1251	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1252	SUBROUTINE grid_hgr__angle(td_nam, jpi, jpj)
[7314]	1253	!!----------------------------------------------------------------------
[12080]	1254	!> @brief This subroutine compute angles between model grid lines and the North direction
[7314]	1255	!>
	1256	!> @details
	1257	!> ** Method :
	1258	!>
	1259	!> ** Action : Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays:
	1260	!> sinus and cosinus of the angle between the north-south axe and the
	1261	!> j-direction at t, u, v and f-points
	1262	!>
	1263	!> History :
	1264	!> 7.0 ! 96-07 (O. Marti ) Original code
	1265	!> 8.0 ! 98-06 (G. Madec )
	1266	!> 8.5 ! 98-06 (G. Madec ) Free form, F90 + opt.
	1267	!> 9.2 ! 07-04 (S. Masson) Add T, F points and bugfix in cos lateral boundary
	1268	!>
	1269	!> @author J.Paul
	1270	!> @date September, 2015 - rewrite from geo2ocean
	1271	!>
	1272	!> @param[in] td_nam
	1273	!> @param[in] jpi
	1274	!> @param[in] jpj
	1275	!!----------------------------------------------------------------------
[12080]	1276
[7314]	1277	IMPLICIT NONE
[12080]	1278
[7314]	1279	! Argument
	1280	TYPE(TNAMH), INTENT(IN) :: td_nam
	1281	INTEGER(i4), INTENT(IN) :: jpi
	1282	INTEGER(i4), INTENT(IN) :: jpj
	1283
	1284	! local variable
	1285	REAL(dp) :: zlam, zphi
	1286	REAL(dp) :: zlan, zphh
	1287	REAL(dp) :: zxnpt, zynpt, znnpt ! x,y components and norm of the vector: T point to North Pole
	1288	REAL(dp) :: zxnpu, zynpu, znnpu ! x,y components and norm of the vector: U point to North Pole
	1289	REAL(dp) :: zxnpv, zynpv, znnpv ! x,y components and norm of the vector: V point to North Pole
	1290	REAL(dp) :: zxnpf, zynpf, znnpf ! x,y components and norm of the vector: F point to North Pole
	1291	REAL(dp) :: zxvvt, zyvvt, znvvt ! x,y components and norm of the vector: between V points below and above a T point
	1292	REAL(dp) :: zxffu, zyffu, znffu ! x,y components and norm of the vector: between F points below and above a U point
	1293	REAL(dp) :: zxffv, zyffv, znffv ! x,y components and norm of the vector: between F points left and right a V point
	1294	REAL(dp) :: zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point
	1295
	1296	! loop indices
	1297	INTEGER(i4) :: ji
	1298	INTEGER(i4) :: jj
	1299	!!----------------------------------------------------------------------
	1300
	1301	! ============================= !
	1302	! Compute the cosinus and sinus !
	1303	! ============================= !
	1304	! (computation done on the north stereographic polar plane)
	1305
	1306	DO jj = 2, jpj-1
	1307	!CDIR NOVERRCHK
	1308	DO ji = 2, jpi ! vector opt.
	1309
	1310	! north pole direction & modulous (at t-point)
	1311	zlam = tg_glamt%d_value(ji,jj,1,1)
	1312	zphi = tg_gphit%d_value(ji,jj,1,1)
	1313	zxnpt = 0._dp - 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1314	zynpt = 0._dp - 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1315	znnpt = zxnptzxnpt + zynptzynpt
	1316
	1317	! north pole direction & modulous (at u-point)
	1318	zlam = tg_glamu%d_value(ji,jj,1,1)
	1319	zphi = tg_gphiu%d_value(ji,jj,1,1)
	1320	zxnpu = 0._dp - 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1321	zynpu = 0._dp - 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1322	znnpu = zxnpuzxnpu + zynpuzynpu
	1323
	1324	! north pole direction & modulous (at v-point)
	1325	zlam = tg_glamv%d_value(ji,jj,1,1)
	1326	zphi = tg_gphiv%d_value(ji,jj,1,1)
	1327	zxnpv = 0._dp - 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1328	zynpv = 0._dp - 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1329	znnpv = zxnpvzxnpv + zynpvzynpv
	1330
	1331	! north pole direction & modulous (at f-point)
	1332	zlam = tg_glamf%d_value(ji,jj,1,1)
	1333	zphi = tg_gphif%d_value(ji,jj,1,1)
	1334	zxnpf = 0._dp - 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1335	zynpf = 0._dp - 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp )
	1336	znnpf = zxnpfzxnpf + zynpfzynpf
	1337
	1338	! j-direction: v-point segment direction (around t-point)
	1339	zlam = tg_glamv%d_value(ji,jj ,1,1)
	1340	zphi = tg_gphiv%d_value(ji,jj ,1,1)
	1341	zlan = tg_glamv%d_value(ji,jj-1,1,1)
	1342	zphh = tg_gphiv%d_value(ji,jj-1,1,1)
	1343	zxvvt = 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1344	& - 2._dp * COS( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1345	zyvvt = 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1346	& - 2._dp * SIN( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1347	znvvt = SQRT( znnpt * ( zxvvtzxvvt + zyvvtzyvvt ) )
	1348	znvvt = MAX( znvvt, dp_eps )
	1349
	1350	! j-direction: f-point segment direction (around u-point)
	1351	zlam = tg_glamf%d_value(ji,jj ,1,1)
	1352	zphi = tg_gphif%d_value(ji,jj ,1,1)
	1353	zlan = tg_glamf%d_value(ji,jj-1,1,1)
	1354	zphh = tg_gphif%d_value(ji,jj-1,1,1)
	1355	zxffu = 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1356	& - 2._dp * COS( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1357	zyffu = 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1358	& - 2._dp * SIN( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1359	znffu = SQRT( znnpu * ( zxffuzxffu + zyffuzyffu ) )
	1360	znffu = MAX( znffu, dp_eps )
	1361
	1362	! i-direction: f-point segment direction (around v-point)
	1363	zlam = tg_glamf%d_value(ji ,jj,1,1)
	1364	zphi = tg_gphif%d_value(ji ,jj,1,1)
	1365	zlan = tg_glamf%d_value(ji-1,jj,1,1)
	1366	zphh = tg_gphif%d_value(ji-1,jj,1,1)
	1367	zxffv = 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1368	& - 2._dp * COS( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1369	zyffv = 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1370	& - 2._dp * SIN( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1371	znffv = SQRT( znnpv * ( zxffvzxffv + zyffvzyffv ) )
	1372	znffv = MAX( znffv, dp_eps )
	1373
	1374	! j-direction: u-point segment direction (around f-point)
	1375	zlam = tg_glamu%d_value(ji,jj+1,1,1)
	1376	zphi = tg_gphiu%d_value(ji,jj+1,1,1)
	1377	zlan = tg_glamu%d_value(ji,jj ,1,1)
	1378	zphh = tg_gphiu%d_value(ji,jj ,1,1)
	1379	zxuuf = 2._dp * COS( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1380	& - 2._dp * COS( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1381	zyuuf = 2._dp * SIN( dp_deg2radzlam ) TAN( dp_pi/4._dp - dp_deg2rad*zphi/2._dp ) &
	1382	& - 2._dp * SIN( dp_deg2radzlan ) TAN( dp_pi/4._dp - dp_deg2rad*zphh/2._dp )
	1383	znuuf = SQRT( znnpf * ( zxuufzxuuf + zyuufzyuuf ) )
	1384	znuuf = MAX( znuuf, dp_eps )
	1385
	1386	! cosinus and sinus using scalar and vectorial products
	1387	tg_gsint%d_value(ji,jj,1,1) = ( zxnptzyvvt - zynptzxvvt ) / znvvt
	1388	tg_gcost%d_value(ji,jj,1,1) = ( zxnptzxvvt + zynptzyvvt ) / znvvt
	1389
	1390	tg_gsinu%d_value(ji,jj,1,1) = ( zxnpuzyffu - zynpuzxffu ) / znffu
	1391	tg_gcosu%d_value(ji,jj,1,1) = ( zxnpuzxffu + zynpuzyffu ) / znffu
	1392
	1393	tg_gsinf%d_value(ji,jj,1,1) = ( zxnpfzyuuf - zynpfzxuuf ) / znuuf
	1394	tg_gcosf%d_value(ji,jj,1,1) = ( zxnpfzxuuf + zynpfzyuuf ) / znuuf
	1395
	1396	! (caution, rotation of 90 degres)
	1397	tg_gsinv%d_value(ji,jj,1,1) = ( zxnpvzxffv + zynpvzyffv ) / znffv
	1398	tg_gcosv%d_value(ji,jj,1,1) =-( zxnpvzyffv - zynpvzxffv ) / znffv
	1399
	1400	END DO
	1401	END DO
	1402
	1403	! =============== !
	1404	! Geographic mesh !
	1405	! =============== !
	1406
	1407	DO jj = 2, jpj-1
	1408	DO ji = 2, jpi ! vector opt.
	1409	IF( MOD( ABS( tg_glamv%d_value(ji,jj,1,1) - tg_glamv%d_value(ji,jj-1,1,1) ), 360._dp ) < 1.e-8 ) THEN
	1410	tg_gsint%d_value(ji,jj,1,1) = 0._dp
	1411	tg_gcost%d_value(ji,jj,1,1) = 1._dp
	1412	ENDIF
	1413	IF( MOD( ABS( tg_glamf%d_value(ji,jj,1,1) - tg_glamf%d_value(ji,jj-1,1,1) ), 360._dp ) < 1.e-8 ) THEN
	1414	tg_gsinu%d_value(ji,jj,1,1) = 0._dp
	1415	tg_gcosu%d_value(ji,jj,1,1) = 1._dp
	1416	ENDIF
	1417	IF( ABS( tg_gphif%d_value(ji,jj,1,1) - tg_gphif%d_value(ji-1,jj,1,1) ) < 1.e-8 ) THEN
	1418	tg_gsinv%d_value(ji,jj,1,1) = 0._dp
	1419	tg_gcosv%d_value(ji,jj,1,1) = 1._dp
	1420	ENDIF
	1421	IF( MOD( ABS( tg_glamu%d_value(ji,jj,1,1) - tg_glamu%d_value(ji,jj+1,1,1) ), 360._dp ) < 1.e-8 ) THEN
	1422	tg_gsinf%d_value(ji,jj,1,1) = 0._dp
	1423	tg_gcosf%d_value(ji,jj,1,1) = 1._dp
	1424	ENDIF
	1425	END DO
	1426	END DO
	1427
	1428	! =========================== !
	1429	! Lateral boundary conditions !
	1430	! =========================== !
	1431
	1432	! lateral boundary cond.: T-, U-, V-, F-pts, sgn
	1433	CALL lbc_lnk( tg_gcost%d_value(:,:,1,1), 'T', td_nam%i_perio, -1._dp )
	1434	CALL lbc_lnk( tg_gcosu%d_value(:,:,1,1), 'U', td_nam%i_perio, -1._dp )
	1435	CALL lbc_lnk( tg_gcosv%d_value(:,:,1,1), 'V', td_nam%i_perio, -1._dp )
	1436	CALL lbc_lnk( tg_gcosf%d_value(:,:,1,1), 'F', td_nam%i_perio, -1._dp )
	1437
	1438	CALL lbc_lnk( tg_gsint%d_value(:,:,1,1), 'T', td_nam%i_perio, -1._dp )
	1439	CALL lbc_lnk( tg_gsinu%d_value(:,:,1,1), 'U', td_nam%i_perio, -1._dp )
	1440	CALL lbc_lnk( tg_gsinv%d_value(:,:,1,1), 'V', td_nam%i_perio, -1._dp )
	1441	CALL lbc_lnk( tg_gsinf%d_value(:,:,1,1), 'F', td_nam%i_perio, -1._dp )
	1442
	1443	END SUBROUTINE grid_hgr__angle
[12080]	1444	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[7314]	1445	END MODULE grid_hgr

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

Download in other formats: