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

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months ago

update nemo trunk

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