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 @ 7025

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

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

Download in other formats: