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

Last change on this file since 239 was 239, checked in by opalod, 19 years ago

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