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

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

Last change on this file since 7153 was 7153, checked in by jpaul, 7 years ago

see ticket #1781

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