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 trunk/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMO/OPA_SRC/DOM/domhgr.F90 @ 418

Last change on this file since 418 was 418, checked in by opalod, 18 years ago

nemo_v1_compil_012 : CT : remove missing declaration or warning compilation
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.8 KB
Line 
1MODULE domhgr
2   !!==============================================================================
3   !!                       ***  MODULE domhgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_hgr        : initialize the horizontal mesh
9   !!   hgr_read       : read "coordinate" NetCDF file
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE dom_oce         ! ocean space and time domain
13   USE phycst          ! physical constants
14   USE in_out_manager  ! I/O manager
15
16   IMPLICIT NONE
17   PRIVATE
18
19   !! * Module variables
20   REAL(wp) ::   glam0, gphi0           ! variables corresponding to parameters
21      !                                 ! ppglam0 ppgphi0 set in par_oce
22
23   !! * Routine accessibility
24   PUBLIC dom_hgr        ! called by domain.F90
25   !!----------------------------------------------------------------------
26   !!   OPA 9.0 , LOCEAN-IPSL (2005)
27   !! $Header$
28   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33   SUBROUTINE dom_hgr
34      !!----------------------------------------------------------------------
35      !!                  ***  ROUTINE dom_hgr  ***
36      !!
37      !! ** Purpose :   Compute the geographical position (in degre) of the
38      !!      model grid-points,  the horizontal scale factors (in meters) and
39      !!      the Coriolis factor (in s-1).
40      !!
41      !! ** Method  :   The geographical position of the model grid-points is
42      !!      defined from analytical functions, fslam and fsphi, the deriva-
43      !!      tives of which gives the horizontal scale factors e1,e2.
44      !!      Defining two function fslam and fsphi and their derivatives in
45      !!      the two horizontal directions (fse1 and fse2), the model grid-
46      !!      point position and scale factors are given by:
47      !!         t-point:
48      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
49      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
50      !!         u-point:
51      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
52      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
53      !!         v-point:
54      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
55      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
56      !!            f-point:
57      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
58      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
59      !!      Where fse1 and fse2 are defined by:
60      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
61      !!                                     +          di(fsphi) **2 )(i,j)
62      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
63      !!                                     +          dj(fsphi) **2 )(i,j)
64      !!
65      !!        The coriolis factor is given at z-point by:
66      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
67      !!
68      !!        This routine is given as an example, it must be modified
69      !!      following the user s desiderata. nevertheless, the output as
70      !!      well as the way to compute the model grid-point position and
71      !!      horizontal scale factors must be respected in order to insure
72      !!      second order accuracy schemes.
73      !!
74      !! N.B. If the domain is periodic, verify that scale factors are also
75      !!      periodic, and the coriolis term again.
76      !!
77      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
78      !!                u-, v- and f-points (in degre)
79      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
80      !!               u-, v-  and f-points (in degre)
81      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
82      !!      scale factors (in meters) at t-, u-, v-, and f-points.
83      !!        define ff: coriolis factor at f-point
84      !!
85      !! References :
86      !!      Marti, Madec and Delecluse, 1992, j. geophys. res., in press.
87      !!
88      !! History :
89      !!        !  88-03  (G. Madec)
90      !!        !  91-11  (G. Madec)
91      !!        !  92-06  (M. Imbard)
92      !!        !  96-01  (G. Madec)  terrain following coordinates
93      !!        !  97-02  (G. Madec)  print mesh informations
94      !!        !  01-09  (M. Levy)  eel config: grid in km, beta-plane
95      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, namelist
96      !!   9.0  !  04-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
97      !!                  use of parameters in par_CONFIG-Rxx.h90, not in namelist
98      !!        !  04-05  (A. Koch-Larrouy) Add Gyre configuration
99      !!----------------------------------------------------------------------
100      !! * local declarations
101      INTEGER  ::   ji, jj              ! dummy loop indices
102      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
103      INTEGER  ::   ijeq                ! index of equator T point (used in case 4)
104      REAL(wp) ::   &
105         zti, zui, zvi, zfi,         &  ! temporary scalars
106         ztj, zuj, zvj, zfj,         &  !
107         zphi0, zbeta, znorme,       &  !
108         zarg, zf0
109      REAL(wp) ::   &
110         zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg,   &
111         zphi1, zsin_alpha, zim05, zjm05
112      !!----------------------------------------------------------------------
113
114      IF(lwp) THEN
115         WRITE(numout,*)
116         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
117         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
118         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
119         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
120         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
121         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
122         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
123         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
124      ENDIF
125
126
127      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
128
129      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
130
131         IF(lwp) WRITE(numout,*)
132         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
133
134         CALL hgr_read           ! Defaultl option  :   NetCDF file
135
136         !                                                ! =====================
137         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
138            !                                             ! =====================
139            IF( n_cla == 0 ) THEN
140
141               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
142               ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
143               IF(lwp) WRITE(numout,*)
144               IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar : e2u reduced to 20 km'
145               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
146               ij0 =  88   ;   ij1 =  88   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
147               IF(lwp) WRITE(numout,*)
148               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 18 km'
149            ENDIF
150
151            ii0 = 145   ;   ii1 = 146        ! Sound Strait (e2u = 15 km)
152            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
153            IF(lwp) WRITE(numout,*)
154            IF(lwp) WRITE(numout,*) '             orca_r2: Reduced e2u at the Sound Strait'
155            !
156         ENDIF
157
158         !                                                ! ======================
159         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
160            !                                             ! ======================
161            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
162            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
163            IF(lwp) WRITE(numout,*)
164            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
165            !
166         ENDIF
167
168
169         ! N.B. :  General case, lat and long function of both i and j indices:
170         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
171         !                                  + (                           fsdiph( zti, ztj ) )**2  )
172         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
173         !                                  + (                           fsdiph( zui, zuj ) )**2  )
174         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
175         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
176         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
177         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
178         !
179         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
180         !                                  + (                           fsdjph( zti, ztj ) )**2  )
181         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
182         !                                  + (                           fsdjph( zui, zuj ) )**2  )
183         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
184         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
185         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
186         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
187
188
189      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
190
191         IF(lwp) WRITE(numout,*)
192         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
193         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
194
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
198               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
199               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
200               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
201         ! Longitude
202               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
203               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
204               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
205               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
206         ! Latitude
207               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
208               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
209               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
210               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
211         ! e1
212               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
213               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
214               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
215               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
216         ! e2
217               e2t(ji,jj) = ra * rad * ppe2_deg
218               e2u(ji,jj) = ra * rad * ppe2_deg
219               e2v(ji,jj) = ra * rad * ppe2_deg
220               e2f(ji,jj) = ra * rad * ppe2_deg
221            END DO
222         END DO
223
224
225      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
226
227         IF(lwp) WRITE(numout,*)
228         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
229         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
230
231         ! Position coordinates (in kilometers)
232         !                          ==========
233         glam0 = 0.e0
234         gphi0 = - ppe2_m * 1.e-3
235         
236#if defined key_agrif && defined key_eel_r6
237         IF (.Not.Agrif_Root()) THEN
238           glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
239           gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
240           ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
241           ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
242         ENDIF
243#endif         
244         DO jj = 1, jpj
245            DO ji = 1, jpi
246               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
247               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
248               glamv(ji,jj) = glamt(ji,jj)
249               glamf(ji,jj) = glamu(ji,jj)
250   
251               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
252               gphiu(ji,jj) = gphit(ji,jj)
253               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
254               gphif(ji,jj) = gphiv(ji,jj)
255            END DO
256         END DO
257
258         ! Horizontal scale factors (in meters)
259         !                              ======
260         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
261         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
262         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
263         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
264
265      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
266
267         IF(lwp) WRITE(numout,*)
268         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
269         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
270         IF ( ppgphi0 == -90 ) THEN
271                IF(lwp) WRITE(numout,*) ' Mercator grid cannot start at south pole !!!! '
272                IF(lwp) WRITE(numout,*) ' We stop '
273                STOP
274         ENDIF
275
276         !  Find index corresponding to the equator, given the grid spacing e1_deg
277         !  and the (approximate) southern latitude ppgphi0.
278         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
279         !  The formula should work even if the equator is outside the domain.
280         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
281         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
282         IF(  ppgphi0 > 0 )  ijeq = -ijeq
283
284         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
285
286         DO jj = 1, jpj
287            DO ji = 1, jpi
288               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
289               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
290               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
291               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
292         ! Longitude
293               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
294               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
295               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
296               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
297         ! Latitude
298               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
299               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
300               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
301               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
302         ! e1
303               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
304               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
305               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
306               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
307         ! e2
308               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
309               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
310               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
311               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
312            END DO
313         END DO
314
315      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
316
317         IF(lwp) WRITE(numout,*)
318         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
319         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
320
321         ! Position coordinates (in kilometers)
322         !                          ==========
323
324         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
325         zlam1 = -85
326         zphi1 = 29
327         ! resolution in meters
328         ze1 = 106000. / FLOAT(jp_cfg)           
329         ! benchmark: forced the resolution to be about 100 km
330         IF( nbench /= 0 )   ze1 = 106000.e0     
331         zsin_alpha = - SQRT( 2. ) / 2.
332         zcos_alpha =   SQRT( 2. ) / 2.
333         ze1deg = ze1 / (ra * rad)
334         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
335         !                                                          ! at the right jp_cfg resolution
336         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
337         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
338
339         IF(lwp) WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
340         IF(lwp) WRITE(numout,*) 'ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
341
342         DO jj = 1, jpj
343           DO ji = 1, jpi
344             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
345             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
346
347             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
348             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
349
350             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
351             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
352
353             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
354             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
355
356             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
357             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
358           END DO
359          END DO
360
361         ! Horizontal scale factors (in meters)
362         !                              ======
363         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
364         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
365         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
366         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
367
368      CASE DEFAULT
369         IF(lwp) WRITE(numout,cform_err)
370         IF(lwp) WRITE(numout,*) '          bad flag value for jphgr_msh = ', jphgr_msh
371         nstop = nstop + 1
372
373      END SELECT
374
375
376      ! Control printing : Grid informations (if not restart)
377      ! ----------------
378
379      IF(lwp .AND. .NOT.ln_rstart ) THEN
380         WRITE(numout,*)
381         WRITE(numout,*) '          longitude and e1 scale factors'
382         WRITE(numout,*) '          ------------------------------'
383         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
384            glamv(ji,1), glamf(ji,1),   &
385            e1t(ji,1), e1u(ji,1),   &
386            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3879300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
388            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
389         
390         WRITE(numout,*)
391         WRITE(numout,*) '          latitude and e2 scale factors'
392         WRITE(numout,*) '          -----------------------------'
393         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
394            &                     gphiv(1,jj), gphif(1,jj),   &
395            &                     e2t  (1,jj), e2u  (1,jj),   &
396            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
397      ENDIF
398
399     
400      IF( nprint == 1 .AND. lwp ) THEN
401         WRITE(numout,*) '          e1u e2u '
402         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
403         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
404         WRITE(numout,*) '          e1v e2v  '
405         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
406         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
407         WRITE(numout,*) '          e1f e2f  '
408         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
409         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
410      ENDIF
411
412
413      ! ================= !
414      !  Coriolis factor  !
415      ! ================= !
416
417      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
418
419      CASE ( 0, 1, 4 )               ! mesh on the sphere
420
421         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
422
423      CASE ( 2 )                     ! f-plane at ppgphi0
424
425         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
426
427         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
428
429      CASE ( 3 )                     ! beta-plane
430
431         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
432         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
433         
434#if defined key_agrif && defined key_eel_r6
435         IF (.Not.Agrif_Root()) THEN
436           zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
437         ENDIF
438#endif         
439         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
440
441         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
442         
443         IF(lwp) WRITE(numout,*) 
444         IF(lwp) WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff(1,1)
445         IF(lwp) WRITE(numout,*) ' Coriolis parameter varies from ', ff(1,1),' to ', ff(1,jpj)
446
447      CASE ( 5 )                     ! beta-plane and rotated domain
448
449         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
450         zphi0 = 15.e0                                                      ! latitude of the first row F-points
451         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
452
453         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
454
455         IF(lwp) WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(1,1)
456         IF(lwp) WRITE(numout,*) '                      Coriolis parameter varies from ', ff(1,1),' to ', ff(1,jpj)
457
458      END SELECT
459
460
461      ! Control of domain for symetrical condition
462      ! ------------------------------------------
463      ! The equator line must be the latitude coordinate axe
464
465      IF( nperio == 2 ) THEN
466         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
467         IF( znorme > 1.e-13 ) THEN
468            IF(lwp) WRITE(numout,cform_err)
469            IF(lwp) WRITE(numout,*) ' ===>>>> : symmetrical condition: rerun with good equator line'
470            nstop = nstop + 1
471         ENDIF
472      ENDIF
473
474   END SUBROUTINE dom_hgr
475
476
477   SUBROUTINE hgr_read
478      !!---------------------------------------------------------------------
479      !!              ***  ROUTINE hgr_read  ***
480      !!
481      !! ** Purpose :   Read a coordinate file in NetCDF format
482      !!
483      !! ** Method  :   The mesh file has been defined trough a analytical
484      !!      or semi-analytical method. It is read in a NetCDF file.
485      !!     
486      !! References :
487      !!      Marti, Madec and Delecluse, 1992, JGR, 97, 12,763-12,766.
488      !!      Madec, Imbard, 1996, Clim. Dyn., 12, 381-388.
489      !!
490      !! History :
491      !!        !         (O. Marti)  Original code
492      !!        !  91-03  (G. Madec)
493      !!        !  92-07  (M. Imbard)
494      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
495      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
496      !!   8.5  !  02-06  (G. Madec)  F90: Free form
497      !!----------------------------------------------------------------------
498      !! * Modules used
499      USE ioipsl
500
501      !! * Local declarations
502      LOGICAL ::   llog = .FALSE.
503      CHARACTER(len=21) ::   clname
504      INTEGER  ::   ji, jj              ! dummy loop indices
505      INTEGER  ::   inum                ! temporary logical unit
506      INTEGER  ::   ilev, itime         ! temporary integers
507      REAL(wp) ::   zdt, zdate0         ! temporary scalars
508      REAL(wp) ::   zdept(1)            ! temporary workspace
509      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
510         zlamt, zphit, zdta             ! temporary workspace (NetCDF read)
511      !!----------------------------------------------------------------------
512      clname = 'coordinates'
513#if defined key_agrif
514      if ( .NOT. Agrif_Root() ) then
515         clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
516      endif
517#endif         
518
519
520      ! 1. Read of the grid coordinates and scale factors
521      ! -------------------------------------------------
522
523      IF(lwp) THEN
524         WRITE(numout,*)
525         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
526         WRITE(numout,*) '~~~~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
527      ENDIF
528
529      ! read the file
530      itime = 0
531      ilev = 1   
532      zlamt(:,:) = 0.e0
533      zphit(:,:) = 0.e0
534      CALL restini( clname, jpidta, jpjdta, zlamt , zphit,   &
535         &                  ilev  , zdept , 'NONE',          &
536         &                  itime , zdate0, zdt   , inum, domain_id=nidom )
537
538      CALL restget( inum, 'glamt', jpidta, jpjdta, 1, itime, llog, zdta )
539      DO jj = 1, nlcj
540         DO ji = 1, nlci
541            glamt(ji,jj) = zdta(mig(ji),mjg(jj))
542         END DO
543      END DO
544      CALL restget( inum, 'glamu', jpidta, jpjdta, 1, itime, llog, zdta )
545      DO jj = 1, nlcj
546         DO ji = 1, nlci
547            glamu(ji,jj) = zdta(mig(ji),mjg(jj))                   
548         END DO
549      END DO
550      CALL restget( inum, 'glamv', jpidta, jpjdta, 1, itime, llog, zdta )
551      DO jj = 1, nlcj
552         DO ji = 1, nlci
553            glamv(ji,jj) = zdta(mig(ji),mjg(jj))                   
554         END DO
555      END DO
556      CALL restget( inum, 'glamf', jpidta, jpjdta, 1, itime, llog, zdta )
557      DO jj = 1, nlcj
558         DO ji = 1, nlci
559            glamf(ji,jj) = zdta(mig(ji),mjg(jj))                   
560         END DO
561      END DO
562      CALL restget( inum, 'gphit', jpidta, jpjdta, 1, itime, llog, zdta )
563      DO jj = 1, nlcj
564         DO ji = 1, nlci
565            gphit(ji,jj) = zdta(mig(ji),mjg(jj))                   
566         END DO
567      END DO
568      CALL restget( inum, 'gphiu', jpidta, jpjdta, 1, itime, llog, zdta )
569      DO jj = 1, nlcj
570         DO ji = 1, nlci
571            gphiu(ji,jj) = zdta(mig(ji),mjg(jj))                   
572         END DO
573      END DO
574      CALL restget( inum, 'gphiv', jpidta, jpjdta, 1, itime, llog, zdta )
575      DO jj = 1, nlcj
576         DO ji = 1, nlci
577            gphiv(ji,jj) = zdta(mig(ji),mjg(jj))                   
578         END DO
579      END DO
580      CALL restget( inum, 'gphif', jpidta, jpjdta, 1, itime, llog, zdta )
581      DO jj = 1, nlcj
582         DO ji = 1, nlci
583            gphif(ji,jj) = zdta(mig(ji),mjg(jj))                   
584         END DO
585      END DO
586      CALL restget( inum, 'e1t', jpidta, jpjdta, 1, itime, llog, zdta )
587      DO jj = 1, nlcj
588         DO ji = 1, nlci
589            e1t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
590         END DO
591      END DO
592      CALL restget( inum, 'e1u', jpidta, jpjdta, 1, itime, llog, zdta )
593      DO jj = 1, nlcj
594         DO ji = 1, nlci
595            e1u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
596         END DO
597      END DO
598      CALL restget( inum, 'e1v', jpidta, jpjdta, 1, itime, llog, zdta )
599      DO jj = 1, nlcj
600         DO ji = 1, nlci
601            e1v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
602         END DO
603      END DO
604      CALL restget( inum, 'e1f', jpidta, jpjdta, 1, itime, llog, zdta )
605      DO jj = 1, nlcj
606         DO ji = 1, nlci
607            e1f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
608         END DO
609      END DO
610      CALL restget( inum, 'e2t', jpidta, jpjdta, 1, itime, llog, zdta )
611      DO jj = 1, nlcj
612         DO ji = 1, nlci
613            e2t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
614         END DO
615      END DO
616      CALL restget( inum, 'e2u', jpidta, jpjdta, 1, itime, llog, zdta )
617      DO jj = 1, nlcj
618         DO ji = 1, nlci
619            e2u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
620         END DO
621      END DO
622      CALL restget( inum, 'e2v', jpidta, jpjdta, 1, itime, llog, zdta )
623      DO jj = 1, nlcj
624         DO ji = 1, nlci
625            e2v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
626         END DO
627      END DO
628      CALL restget( inum, 'e2f', jpidta, jpjdta, 1, itime, llog, zdta )
629      DO jj = 1, nlcj
630         DO ji = 1, nlci
631            e2f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
632         END DO
633      END DO
634
635      CALL restclo( inum )
636
637      ! set extra rows add in mpp to none zero values
638      DO jj = nlcj+1, jpj
639         DO ji = 1, nlci
640            glamt(ji,jj) = glamt(ji,1)   ;   gphit(ji,jj) = gphit(ji,1)
641            glamu(ji,jj) = glamu(ji,1)   ;   gphiu(ji,jj) = gphiu(ji,1)
642            glamv(ji,jj) = glamv(ji,1)   ;   gphiv(ji,jj) = gphiv(ji,1)
643            glamf(ji,jj) = glamf(ji,1)   ;   gphif(ji,jj) = gphif(ji,1)
644            e1t  (ji,jj) = e1t  (ji,1)   ;   e2t  (ji,jj) = e2t  (ji,1)
645            e1u  (ji,jj) = e1u  (ji,1)   ;   e2u  (ji,jj) = e2u  (ji,1)
646            e1v  (ji,jj) = e1v  (ji,1)   ;   e2v  (ji,jj) = e2v  (ji,1)
647            e1f  (ji,jj) = e1f  (ji,1)   ;   e2f  (ji,jj) = e2f  (ji,1)
648         END DO
649      END DO
650
651      ! set extra columns add in mpp to none zero values
652      DO ji = nlci+1, jpi
653         glamt(ji,:) = glamt(1,:)   ;   gphit(ji,:) = gphit(1,:)
654         glamu(ji,:) = glamu(1,:)   ;   gphiu(ji,:) = gphiu(1,:)
655         glamv(ji,:) = glamv(1,:)   ;   gphiv(ji,:) = gphiv(1,:)
656         glamf(ji,:) = glamf(1,:)   ;   gphif(ji,:) = gphif(1,:)
657         e1t  (ji,:) = e1t  (1,:)   ;   e2t  (ji,:) = e2t  (1,:)
658         e1u  (ji,:) = e1u  (1,:)   ;   e2u  (ji,:) = e2u  (1,:)
659         e1v  (ji,:) = e1v  (1,:)   ;   e2v  (ji,:) = e2v  (1,:)
660         e1f  (ji,:) = e1f  (1,:)   ;   e2f  (ji,:) = e2f  (1,:)
661      END DO
662
663   END SUBROUTINE hgr_read
664
665   !!======================================================================
666END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.