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 branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_hgr.f90 @ 7233

Last change on this file since 7233 was 7233, checked in by jpaul, 7 years ago
see ticket #1781
File size: 59.1 KB

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

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

Download in other formats: