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.
domhgr.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 7158

Last change on this file since 7158 was 7158, checked in by clem, 7 years ago

debug branch
  • Property svn:keywords set to Id
File size: 37.9 KB
Line 
1MODULE domhgr
2   !!==============================================================================
3   !!                       ***  MODULE domhgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1988-03  (G. Madec) Original code
7   !!            7.0  ! 1996-01  (G. Madec)  terrain following coordinates
8   !!            8.0  ! 1997-02  (G. Madec)  print mesh informations
9   !!            8.1  ! 1999-11  (M. Imbard) NetCDF format with IO-IPSL
10   !!            8.2  ! 2000-08  (D. Ludicone) Reduced section at Bab el Mandeb
11   !!             -   ! 2001-09  (M. Levy)  eel config: grid in km, beta-plane
12   !!  NEMO      1.0  ! 2002-08  (G. Madec)  F90: Free form and module, namelist
13   !!             -   ! 2004-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
14   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist
15   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration
16   !!            4.0  ! 2011-02  (G. Madec) add cell surface (e1e2t)
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   dom_hgr       : initialize the horizontal mesh
21   !!   hgr_read      : read "coordinate" NetCDF file
22   !!----------------------------------------------------------------------
23   USE dom_oce        ! ocean space and time domain
24   USE phycst         ! physical constants
25   USE in_out_manager ! I/O manager
26   USE lib_mpp        ! MPP library
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
33
34   PUBLIC   dom_hgr   ! called by domain.F90
35
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dom_hgr
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE dom_hgr  ***
46      !!
47      !! ** Purpose :   Compute the geographical position (in degre) of the
48      !!      model grid-points,  the horizontal scale factors (in meters) and
49      !!      the Coriolis factor (in s-1).
50      !!
51      !! ** Method  :   The geographical position of the model grid-points is
52      !!      defined from analytical functions, fslam and fsphi, the deriva-
53      !!      tives of which gives the horizontal scale factors e1,e2.
54      !!      Defining two function fslam and fsphi and their derivatives in
55      !!      the two horizontal directions (fse1 and fse2), the model grid-
56      !!      point position and scale factors are given by:
57      !!         t-point:
58      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
59      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
60      !!         u-point:
61      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
62      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
63      !!         v-point:
64      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
65      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
66      !!            f-point:
67      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
68      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
69      !!      Where fse1 and fse2 are defined by:
70      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
71      !!                                     +          di(fsphi) **2 )(i,j)
72      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
73      !!                                     +          dj(fsphi) **2 )(i,j)
74      !!
75      !!        The coriolis factor is given at z-point by:
76      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
77      !!
78      !!        This routine is given as an example, it must be modified
79      !!      following the user s desiderata. nevertheless, the output as
80      !!      well as the way to compute the model grid-point position and
81      !!      horizontal scale factors must be respected in order to insure
82      !!      second order accuracy schemes.
83      !!
84      !! N.B. If the domain is periodic, verify that scale factors are also
85      !!      periodic, and the coriolis term again.
86      !!
87      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
88      !!                u-, v- and f-points (in degre)
89      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
90      !!               u-, v-  and f-points (in degre)
91      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
92      !!      scale factors (in meters) at t-, u-, v-, and f-points.
93      !!        define ff: coriolis factor at f-point
94      !!
95      !! References :   Marti, Madec and Delecluse, 1992, JGR
96      !!                Madec, Imbard, 1996, Clim. Dyn.
97      !!----------------------------------------------------------------------
98      INTEGER  ::   ji, jj               ! dummy loop indices
99      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
100      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
101      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
102      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
103      REAL(wp) ::   zphi0, zbeta, znorme !
104      REAL(wp) ::   zarg, zf0, zminff, zmaxff
105      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
106      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
107      INTEGER  ::   isrow                ! index for ORCA1 starting row
108
109      !!----------------------------------------------------------------------
110      !
111      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
112      !
113      IF(lwp) THEN
114         WRITE(numout,*)
115         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
116         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
117         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
118         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
119         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
120         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
121         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
122         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
123      ENDIF
124
125
126      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
127
128      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
129
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
132
133         CALL hgr_read           ! Defaultl option  :   NetCDF file
134
135         !                                                ! =====================
136         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
137            !                                             ! =====================
138            IF( nn_cla == 0 ) THEN
139               !
140               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
141               ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
142               IF(lwp) WRITE(numout,*)
143               IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km'
144               !
145               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
146               ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
147                                               e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3
148               IF(lwp) WRITE(numout,*)
149               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km'
150               IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km'
151            ENDIF
152
153            ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km)
154            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
155            IF(lwp) WRITE(numout,*)
156            IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km'
157            !
158         ENDIF
159
160            !                                             ! =====================
161         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
162            !                                             ! =====================
163            ! This dirty section will be suppressed by simplification process: all this will come back in input files
164            ! Currently these hard-wired indices relate to configuration with
165            ! extend grid (jpjglo=332)
166            ! which had a grid-size of 362x292.
167            !
168            isrow = 332 - jpjglo
169            !
170            ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km)
171            ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
172            IF(lwp) WRITE(numout,*)
173            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km'
174
175            ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km)
176            ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
177            IF(lwp) WRITE(numout,*)
178            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km'
179
180            ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km)
181            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
182            IF(lwp) WRITE(numout,*)
183            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km'
184
185            ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
186            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
187            IF(lwp) WRITE(numout,*)
188            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km'
189
190            ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km)
191            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
192            IF(lwp) WRITE(numout,*)
193            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km'
194
195            ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km)
196            ij0 = 164 - isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
197            IF(lwp) WRITE(numout,*)
198            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km'
199
200            ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km)
201            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
202            IF(lwp) WRITE(numout,*)
203            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km'
204
205            ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km)
206            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
207            IF(lwp) WRITE(numout,*)
208            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km'
209            !
210            !
211         ENDIF
212
213         !                                                ! ======================
214         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
215            !                                             ! ======================
216            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
217            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
218            IF(lwp) WRITE(numout,*)
219            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
220            !
221            ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km)
222            ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
223            IF(lwp) WRITE(numout,*)
224            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait'
225            !
226            ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km)
227            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3
228            IF(lwp) WRITE(numout,*)
229            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait'
230            !
231            ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km)
232            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
233            IF(lwp) WRITE(numout,*)
234            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait'
235            !
236            ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km)
237            ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
238            IF(lwp) WRITE(numout,*)
239            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait'
240            !
241            ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km)
242            ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
243            IF(lwp) WRITE(numout,*)
244            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait'
245            !
246            !
247            ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km)
248            ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3
249            IF(lwp) WRITE(numout,*)
250            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb'
251            !
252         ENDIF
253
254
255         ! N.B. :  General case, lat and long function of both i and j indices:
256         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
257         !                                  + (                           fsdiph( zti, ztj ) )**2  )
258         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
259         !                                  + (                           fsdiph( zui, zuj ) )**2  )
260         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
261         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
262         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
263         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
264         !
265         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
266         !                                  + (                           fsdjph( zti, ztj ) )**2  )
267         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
268         !                                  + (                           fsdjph( zui, zuj ) )**2  )
269         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
270         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
271         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
272         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
273
274
275      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
276
277         IF(lwp) WRITE(numout,*)
278         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
279         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
280
281         DO jj = 1, jpj
282            DO ji = 1, jpi
283               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
284               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
285               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
286               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
287         ! Longitude
288               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
289               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
290               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
291               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
292         ! Latitude
293               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
294               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
295               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
296               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
297         ! e1
298               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
299               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
300               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
301               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
302         ! e2
303               e2t(ji,jj) = ra * rad * ppe2_deg
304               e2u(ji,jj) = ra * rad * ppe2_deg
305               e2v(ji,jj) = ra * rad * ppe2_deg
306               e2f(ji,jj) = ra * rad * ppe2_deg
307            END DO
308         END DO
309
310
311      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
312
313         IF(lwp) WRITE(numout,*)
314         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
315         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
316
317         ! Position coordinates (in kilometers)
318         !                          ==========
319         glam0 = 0.e0
320         gphi0 = - ppe2_m * 1.e-3
321         
322#if defined key_agrif 
323         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
324            IF( .NOT. Agrif_Root() ) THEN
325              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
326              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
327              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
328              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
329            ENDIF
330         ENDIF
331#endif         
332         DO jj = 1, jpj
333            DO ji = 1, jpi
334               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
335               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
336               glamv(ji,jj) = glamt(ji,jj)
337               glamf(ji,jj) = glamu(ji,jj)
338   
339               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
340               gphiu(ji,jj) = gphit(ji,jj)
341               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
342               gphif(ji,jj) = gphiv(ji,jj)
343            END DO
344         END DO
345
346         ! Horizontal scale factors (in meters)
347         !                              ======
348         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
349         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
350         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
351         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
352
353      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
354
355         IF(lwp) WRITE(numout,*)
356         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
357         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
358         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
359
360         !  Find index corresponding to the equator, given the grid spacing e1_deg
361         !  and the (approximate) southern latitude ppgphi0.
362         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
363         !  The formula should work even if the equator is outside the domain.
364         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
365         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
366         IF(  ppgphi0 > 0 )  ijeq = -ijeq
367
368         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
369
370         DO jj = 1, jpj
371            DO ji = 1, jpi
372               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
373               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
374               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
375               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
376         ! Longitude
377               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
378               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
379               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
380               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
381         ! Latitude
382               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
383               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
384               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
385               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
386         ! e1
387               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
388               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
389               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
390               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
391         ! e2
392               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
393               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
394               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
395               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
396            END DO
397         END DO
398
399      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
400
401         IF(lwp) WRITE(numout,*)
402         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
403         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
404
405         ! Position coordinates (in kilometers)
406         !                          ==========
407
408         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
409         zlam1 = -85
410         zphi1 = 29
411         ! resolution in meters
412         ze1 = 106000. / FLOAT(jp_cfg)           
413         ! benchmark: forced the resolution to be about 100 km
414         IF( nbench /= 0 )   ze1 = 106000.e0     
415         zsin_alpha = - SQRT( 2. ) / 2.
416         zcos_alpha =   SQRT( 2. ) / 2.
417         ze1deg = ze1 / (ra * rad)
418         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
419         !                                                          ! at the right jp_cfg resolution
420         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
421         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
422
423         IF( nprint==1 .AND. lwp )   THEN
424            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
425            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
426         ENDIF
427
428         DO jj = 1, jpj
429           DO ji = 1, jpi
430             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
431             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
432
433             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
434             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
435
436             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
437             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
438
439             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
440             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
441
442             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
443             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
444           END DO
445          END DO
446
447         ! Horizontal scale factors (in meters)
448         !                              ======
449         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
450         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
451         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
452         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
453
454      CASE ( 6 )                   ! clem: f-plane with irregular grid-spacing
455
456         IF(lwp) WRITE(numout,*)
457         IF(lwp) WRITE(numout,*) '          f-plane with irregular grid-spacing (+- 10%)'
458         IF(lwp) WRITE(numout,*) '          the max is given by ppe1_m and ppe2_m' 
459
460         ! Position coordinates (in kilometers)
461         !                          ==========
462         glam0 = 0._wp
463         gphi0 = 0._wp
464         
465#if defined key_agrif 
466         IF( .NOT. Agrif_Root() ) THEN
467            glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-5
468            gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-5
469            ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
470            ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
471         ENDIF
472#endif         
473
474         DO jj = 1, jpj
475            DO ji = 1, jpi
476               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
477               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
478               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
479               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
480
481               glamt(ji,jj) = glam0 + ppe1_m * 1.e-5 * zti
482               glamu(ji,jj) = glam0 + ppe1_m * 1.e-5 * zui
483               glamv(ji,jj) = glam0 + ppe1_m * 1.e-5 * zvi
484               glamf(ji,jj) = glam0 + ppe1_m * 1.e-5 * zfi
485   
486               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-5 * ztj
487               gphiu(ji,jj) = gphi0 + ppe2_m * 1.e-5 * zuj
488               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-5 * zvj
489               gphif(ji,jj) = gphi0 + ppe2_m * 1.e-5 * zfj
490            END DO
491         END DO
492         
493         ! Horizontal scale factors (in meters)
494         !                              ======
495!! ==> EITHER 1) variable scale factors
496         DO jj = 1, jpj
497            DO ji = 1, jpi
498               !!e1t(ji,jj) = ppe1_m * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji)-REAL(jpiglo+1)*0.5)**2 )  ! gaussian shape
499               !!e2t(ji,jj) = ppe2_m * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj)-REAL(jpjglo+1)*0.5)**2 )  ! gaussian shape
500               e1t(ji,jj) = ppe1_m * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
501               e2t(ji,jj) = ppe2_m * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
502            END DO
503         END DO
504#if defined key_agrif 
505         IF( .NOT. Agrif_Root() ) THEN ! only works if the zoom is positioned at the center of the parent grid
506            DO jj = 1, jpj
507               DO ji = 1, jpi
508                  e1t(ji,jj) = ppe1_m * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5)  &
509                     &                            * REAL(jpiglo) / REAL(Agrif_Parent(jpiglo) * Agrif_Rhox()) )       ! factor to match parent grid
510                  e2t(ji,jj) = ppe2_m * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5)  &
511                     &                            * REAL(jpjglo) / REAL(Agrif_Parent(jpjglo) * Agrif_Rhoy()) )       ! factor to match parent grid
512               END DO
513            END DO
514         ENDIF
515#endif
516!! ==> OR 2) constant scale factors
517         e1t(:,:) = ppe1_m
518         e2t(:,:) = ppe2_m
519         
520         e1u(:,:) = e1t(:,:)      ;      e2u(:,:) = e2t(:,:)
521         e1v(:,:) = e1t(:,:)      ;      e2v(:,:) = e2t(:,:)
522         e1f(:,:) = e1t(:,:)      ;      e2f(:,:) = e2t(:,:)
523
524      CASE DEFAULT
525         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
526         CALL ctl_stop( ctmp1 )
527
528      END SELECT
529     
530      ! T-cell surface
531      ! --------------
532      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
533   
534      ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning)
535      ! ---------------------------------------------------------------------------
536      e12t    (:,:) = e1t(:,:) * e2t(:,:)
537      e12u    (:,:) = e1u(:,:) * e2u(:,:)
538      e12v    (:,:) = e1v(:,:) * e2v(:,:)
539      e12f    (:,:) = e1f(:,:) * e2f(:,:)
540      r1_e12t (:,:) = 1._wp    / e12t(:,:)
541      r1_e12u (:,:) = 1._wp    / e12u(:,:)
542      r1_e12v (:,:) = 1._wp    / e12v(:,:)
543      r1_e12f (:,:) = 1._wp    / e12f(:,:)
544      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
545      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
546      r1_e1t  (:,:) = 1._wp    / e1t(:,:)
547      r1_e1u  (:,:) = 1._wp    / e1u(:,:)
548      r1_e1v  (:,:) = 1._wp    / e1v(:,:)
549      r1_e1f  (:,:) = 1._wp    / e1f(:,:)
550      r1_e2t  (:,:) = 1._wp    / e2t(:,:)
551      r1_e2u  (:,:) = 1._wp    / e2u(:,:)
552      r1_e2v  (:,:) = 1._wp    / e2v(:,:)
553      r1_e2f  (:,:) = 1._wp    / e2f(:,:)
554
555      ! Control printing : Grid informations (if not restart)
556      ! ----------------
557
558      IF( lwp .AND. .NOT.ln_rstart ) THEN
559         WRITE(numout,*)
560         WRITE(numout,*) '          longitude and e1 scale factors'
561         WRITE(numout,*) '          ------------------------------'
562         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
563            glamv(ji,1), glamf(ji,1),   &
564            e1t(ji,1), e1u(ji,1),   &
565            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
5669300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
567            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
568         
569         WRITE(numout,*)
570         WRITE(numout,*) '          latitude and e2 scale factors'
571         WRITE(numout,*) '          -----------------------------'
572         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
573            &                     gphiv(1,jj), gphif(1,jj),   &
574            &                     e2t  (1,jj), e2u  (1,jj),   &
575            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
576      ENDIF
577
578     
579      IF( nprint == 1 .AND. lwp ) THEN
580         WRITE(numout,*) '          e1u e2u '
581         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
582         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
583         WRITE(numout,*) '          e1v e2v  '
584         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
585         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
586         WRITE(numout,*) '          e1f e2f  '
587         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
588         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
589      ENDIF
590
591
592      ! ================= !
593      !  Coriolis factor  !
594      ! ================= !
595
596      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
597
598      CASE ( 0, 1, 4 )               ! mesh on the sphere
599
600         ff  (:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
601         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )  ! clem: coriolis at T-point
602
603      CASE ( 2 )                     ! f-plane at ppgphi0
604
605         ff  (:,:) = 2. * omega * SIN( rad * ppgphi0 )
606         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 )  ! clem: coriolis at T-point
607
608         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
609
610      CASE ( 3 )                     ! beta-plane
611
612         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
613         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
614         
615#if defined key_agrif
616         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
617            IF( .NOT. Agrif_Root() ) THEN
618              zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m)   & 
619                    &           / (ra * rad)
620            ENDIF
621         ENDIF
622#endif         
623         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
624
625         ff  (:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
626         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! clem: coriolis at T-point
627         
628         IF(lwp) THEN
629            WRITE(numout,*) 
630            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
631            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
632         ENDIF
633         IF( lk_mpp ) THEN
634            zminff=ff(nldi,nldj)
635            zmaxff=ff(nldi,nlej)
636            CALL mpp_min( zminff )   ! min over the global domain
637            CALL mpp_max( zmaxff )   ! max over the global domain
638            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
639         END IF
640
641      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
642
643         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
644         zphi0 = 15.e0                                                      ! latitude of the first row F-points
645         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
646
647         ff  (:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
648         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! clem: coriolis at T-point
649
650         IF(lwp) THEN
651            WRITE(numout,*) 
652            WRITE(numout,*) '          Beta-plane and rotated domain : '
653            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
654         ENDIF
655
656         IF( lk_mpp ) THEN
657            zminff=ff(nldi,nldj)
658            zmaxff=ff(nldi,nlej)
659            CALL mpp_min( zminff )   ! min over the global domain
660            CALL mpp_max( zmaxff )   ! max over the global domain
661            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
662         END IF
663
664      END SELECT
665
666
667      ! Control of domain for symetrical condition
668      ! ------------------------------------------
669      ! The equator line must be the latitude coordinate axe
670
671      IF( nperio == 2 ) THEN
672         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
673         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
674      ENDIF
675      !
676      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
677      !
678   END SUBROUTINE dom_hgr
679
680
681   SUBROUTINE hgr_read
682      !!---------------------------------------------------------------------
683      !!              ***  ROUTINE hgr_read  ***
684      !!
685      !! ** Purpose :   Read a coordinate file in NetCDF format
686      !!
687      !! ** Method  :   The mesh file has been defined trough a analytical
688      !!      or semi-analytical method. It is read in a NetCDF file.
689      !!     
690      !!----------------------------------------------------------------------
691      USE iom
692
693      INTEGER ::   inum   ! temporary logical unit
694      !!----------------------------------------------------------------------
695
696      IF(lwp) THEN
697         WRITE(numout,*)
698         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
699         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
700      ENDIF
701     
702      CALL iom_open( 'coordinates', inum )
703     
704      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
705      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
706      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
707      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
708     
709      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
710      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
711      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
712      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
713     
714      CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr )
715      CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr )
716      CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr )
717      CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr )
718     
719      CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr )
720      CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr )
721      CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr )
722      CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr )
723     
724      CALL iom_close( inum )
725     
726! need to be define for the extended grid south of -80S
727! some point are undefined but you need to have e1 and e2 .NE. 0
728      WHERE (e1t==0.0_wp)
729         e1t=1.0e2
730      END WHERE
731      WHERE (e1v==0.0_wp)
732         e1v=1.0e2
733      END WHERE
734      WHERE (e1u==0.0_wp)
735         e1u=1.0e2
736      END WHERE
737      WHERE (e1f==0.0_wp)
738         e1f=1.0e2
739      END WHERE
740      WHERE (e2t==0.0_wp)
741         e2t=1.0e2
742      END WHERE
743      WHERE (e2v==0.0_wp)
744         e2v=1.0e2
745      END WHERE
746      WHERE (e2u==0.0_wp)
747         e2u=1.0e2
748      END WHERE
749      WHERE (e2f==0.0_wp)
750         e2f=1.0e2
751      END WHERE
752       
753    END SUBROUTINE hgr_read
754   
755   !!======================================================================
756END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.