Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grid_hgr.f90 @ 11290

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

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

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

Context Navigation

source: branches/UKMO/r8395_mix-lyr_diag/NEMOGCM/TOOLS/SIREN/src/grid_hgr.f90 @ 11290

Download in other formats: