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

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

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.0 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         DO jj = 1, jpj
236            DO ji = 1, jpi
237               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
238               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
239               glamv(ji,jj) = glamt(ji,jj)
240               glamf(ji,jj) = glamu(ji,jj)
241   
242               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
243               gphiu(ji,jj) = gphit(ji,jj)
244               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
245               gphif(ji,jj) = gphiv(ji,jj)
246            END DO
247         END DO
248
249         ! Horizontal scale factors (in meters)
250         !                              ======
251         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
252         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
253         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
254         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
255
256      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
257
258         IF(lwp) WRITE(numout,*)
259         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
260         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
261         IF ( ppgphi0 == -90 ) THEN
262                IF(lwp) WRITE(numout,*) ' Mercator grid cannot start at south pole !!!! '
263                IF(lwp) WRITE(numout,*) ' We stop '
264                STOP
265         ENDIF
266
267         !  Find index corresponding to the equator, given the grid spacing e1_deg
268         !  and the (approximate) southern latitude ppgphi0.
269         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
270         !  The formula should work even if the equator is outside the domain.
271         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
272         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
273         IF(  ppgphi0 > 0 )  ijeq = -ijeq
274
275         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
276
277         DO jj = 1, jpj
278            DO ji = 1, jpi
279               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
280               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
281               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
282               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
283         ! Longitude
284               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
285               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
286               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
287               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
288         ! Latitude
289               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
290               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
291               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
292               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
293         ! e1
294               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
295               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
296               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
297               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
298         ! e2
299               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
300               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
301               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
302               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
303            END DO
304         END DO
305
306      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
307
308         IF(lwp) WRITE(numout,*)
309         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
310         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
311
312         ! Position coordinates (in kilometers)
313         !                          ==========
314
315         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
316         zlam1 = -85
317         zphi1 = 29
318         ! resolution in meters
319         ze1 = 106000. / FLOAT(jp_cfg)           
320         ! benchmark: forced the resolution to be about 100 km
321         IF( nbench /= 0 )   ze1 = 106000.e0     
322         zsin_alpha = - SQRT( 2. ) / 2.
323         zcos_alpha =   SQRT( 2. ) / 2.
324         ze1deg = ze1 / (ra * rad)
325         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
326         !                                                          ! at the right jp_cfg resolution
327         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
328         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
329
330         IF(lwp) WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
331         IF(lwp) WRITE(numout,*) 'ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
332
333         DO jj = 1, jpj
334           DO ji = 1, jpi
335             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
336             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
337
338             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
339             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
340
341             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
342             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
343
344             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
345             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
346
347             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
348             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
349           END DO
350          END DO
351
352         ! Horizontal scale factors (in meters)
353         !                              ======
354         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
355         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
356         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
357         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
358
359      CASE DEFAULT
360         IF(lwp) WRITE(numout,cform_err)
361         IF(lwp) WRITE(numout,*) '          bad flag value for jphgr_msh = ', jphgr_msh
362         nstop = nstop + 1
363
364      END SELECT
365
366
367      ! Control printing : Grid informations (if not restart)
368      ! ----------------
369
370      IF(lwp .AND. .NOT.ln_rstart ) THEN
371         WRITE(numout,*)
372         WRITE(numout,*) '          longitude and e1 scale factors'
373         WRITE(numout,*) '          ------------------------------'
374         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
375            glamv(ji,1), glamf(ji,1),   &
376            e1t(ji,1), e1u(ji,1),   &
377            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3789300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
379            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
380         
381         WRITE(numout,*)
382         WRITE(numout,*) '          latitude and e2 scale factors'
383         WRITE(numout,*) '          -----------------------------'
384         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
385            &                     gphiv(1,jj), gphif(1,jj),   &
386            &                     e2t  (1,jj), e2u  (1,jj),   &
387            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
388      ENDIF
389
390     
391      IF( nprint == 1 .AND. lwp ) THEN
392         WRITE(numout,*) '          e1u e2u '
393         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
394         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
395         WRITE(numout,*) '          e1v e2v  '
396         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
397         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
398         WRITE(numout,*) '          e1f e2f  '
399         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
400         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
401      ENDIF
402
403
404      ! ================= !
405      !  Coriolis factor  !
406      ! ================= !
407
408      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
409
410      CASE ( 0, 1, 4 )               ! mesh on the sphere
411
412         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
413
414      CASE ( 2 )                     ! f-plane at ppgphi0
415
416         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
417
418         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
419
420      CASE ( 3 )                     ! beta-plane
421
422         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
423         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
424         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
425
426         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
427       
428         IF(lwp) WRITE(numout,*) 
429         IF(lwp) WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff(1,1)
430         IF(lwp) WRITE(numout,*) ' Coriolis parameter varies from ', ff(1,1),' to ', ff(1,jpj)
431
432      CASE ( 5 )                     ! beta-plane and rotated domain
433
434         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
435         zphi0 = 15.e0                                                      ! latitude of the first row F-points
436         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
437
438         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
439
440         IF(lwp) WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(1,1)
441         IF(lwp) WRITE(numout,*) '                      Coriolis parameter varies from ', ff(1,1),' to ', ff(1,jpj)
442
443      END SELECT
444
445
446      ! Control of domain for symetrical condition
447      ! ------------------------------------------
448      ! The equator line must be the latitude coordinate axe
449
450      IF( nperio == 2 ) THEN
451         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
452         IF( znorme > 1.e-13 ) THEN
453            IF(lwp) WRITE(numout,cform_err)
454            IF(lwp) WRITE(numout,*) ' ===>>>> : symmetrical condition: rerun with good equator line'
455            nstop = nstop + 1
456         ENDIF
457      ENDIF
458
459   END SUBROUTINE dom_hgr
460
461
462   SUBROUTINE hgr_read
463      !!---------------------------------------------------------------------
464      !!              ***  ROUTINE hgr_read  ***
465      !!
466      !! ** Purpose :   Read a coordinate file in NetCDF format
467      !!
468      !! ** Method  :   The mesh file has been defined trough a analytical
469      !!      or semi-analytical method. It is read in a NetCDF file.
470      !!     
471      !! References :
472      !!      Marti, Madec and Delecluse, 1992, JGR, 97, 12,763-12,766.
473      !!      Madec, Imbard, 1996, Clim. Dyn., 12, 381-388.
474      !!
475      !! History :
476      !!        !         (O. Marti)  Original code
477      !!        !  91-03  (G. Madec)
478      !!        !  92-07  (M. Imbard)
479      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
480      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
481      !!   8.5  !  02-06  (G. Madec)  F90: Free form
482      !!----------------------------------------------------------------------
483      !! * Modules used
484      USE ioipsl
485
486      !! * Local declarations
487      LOGICAL ::   llog = .FALSE.
488      CHARACTER(len=21) ::   clname = 'coordinates'
489      INTEGER  ::   ji, jj              ! dummy loop indices
490      INTEGER  ::   inum                ! temporary logical unit
491      INTEGER  ::   ilev, itime         ! temporary integers
492      REAL(wp) ::   zdt, zdate0         ! temporary scalars
493      REAL(wp) ::   zdept(1)            ! temporary workspace
494      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
495         zlamt, zphit, zdta             ! temporary workspace (NetCDF read)
496      !!----------------------------------------------------------------------
497
498
499      ! 1. Read of the grid coordinates and scale factors
500      ! -------------------------------------------------
501
502      IF(lwp) THEN
503         WRITE(numout,*)
504         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
505         WRITE(numout,*) '~~~~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
506      ENDIF
507
508      ! read the file
509      itime = 0
510      ilev = 1   
511      zlamt(:,:) = 0.e0
512      zphit(:,:) = 0.e0
513      CALL restini( clname, jpidta, jpjdta, zlamt , zphit,   &
514         &                  ilev  , zdept , clname,          &
515         &                  itime , zdate0, zdt   , inum )
516
517      CALL restget( inum, 'glamt', jpidta, jpjdta, 1, 0, llog, zdta )
518      DO jj = 1, nlcj
519         DO ji = 1, nlci
520            glamt(ji,jj) = zdta(mig(ji),mjg(jj))
521         END DO
522      END DO
523      CALL restget( inum, 'glamu', jpidta, jpjdta, 1, 0, llog, zdta )
524      DO jj = 1, nlcj
525         DO ji = 1, nlci
526            glamu(ji,jj) = zdta(mig(ji),mjg(jj))                   
527         END DO
528      END DO
529      CALL restget( inum, 'glamv', jpidta, jpjdta, 1, 0, llog, zdta )
530      DO jj = 1, nlcj
531         DO ji = 1, nlci
532            glamv(ji,jj) = zdta(mig(ji),mjg(jj))                   
533         END DO
534      END DO
535      CALL restget( inum, 'glamf', jpidta, jpjdta, 1, 0, llog, zdta )
536      DO jj = 1, nlcj
537         DO ji = 1, nlci
538            glamf(ji,jj) = zdta(mig(ji),mjg(jj))                   
539         END DO
540      END DO
541      CALL restget( inum, 'gphit', jpidta, jpjdta, 1, 0, llog, zdta )
542      DO jj = 1, nlcj
543         DO ji = 1, nlci
544            gphit(ji,jj) = zdta(mig(ji),mjg(jj))                   
545         END DO
546      END DO
547      CALL restget( inum, 'gphiu', jpidta, jpjdta, 1, 0, llog, zdta )
548      DO jj = 1, nlcj
549         DO ji = 1, nlci
550            gphiu(ji,jj) = zdta(mig(ji),mjg(jj))                   
551         END DO
552      END DO
553      CALL restget( inum, 'gphiv', jpidta, jpjdta, 1, 0, llog, zdta )
554      DO jj = 1, nlcj
555         DO ji = 1, nlci
556            gphiv(ji,jj) = zdta(mig(ji),mjg(jj))                   
557         END DO
558      END DO
559      CALL restget( inum, 'gphif', jpidta, jpjdta, 1, 0, llog, zdta )
560      DO jj = 1, nlcj
561         DO ji = 1, nlci
562            gphif(ji,jj) = zdta(mig(ji),mjg(jj))                   
563         END DO
564      END DO
565      CALL restget( inum, 'e1t', jpidta, jpjdta, 1, 0, llog, zdta )
566      DO jj = 1, nlcj
567         DO ji = 1, nlci
568            e1t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
569         END DO
570      END DO
571      CALL restget( inum, 'e1u', jpidta, jpjdta, 1, 0, llog, zdta )
572      DO jj = 1, nlcj
573         DO ji = 1, nlci
574            e1u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
575         END DO
576      END DO
577      CALL restget( inum, 'e1v', jpidta, jpjdta, 1, 0, llog, zdta )
578      DO jj = 1, nlcj
579         DO ji = 1, nlci
580            e1v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
581         END DO
582      END DO
583      CALL restget( inum, 'e1f', jpidta, jpjdta, 1, 0, llog, zdta )
584      DO jj = 1, nlcj
585         DO ji = 1, nlci
586            e1f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
587         END DO
588      END DO
589      CALL restget( inum, 'e2t', jpidta, jpjdta, 1, 0, llog, zdta )
590      DO jj = 1, nlcj
591         DO ji = 1, nlci
592            e2t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
593         END DO
594      END DO
595      CALL restget( inum, 'e2u', jpidta, jpjdta, 1, 0, llog, zdta )
596      DO jj = 1, nlcj
597         DO ji = 1, nlci
598            e2u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
599         END DO
600      END DO
601      CALL restget( inum, 'e2v', jpidta, jpjdta, 1, 0, llog, zdta )
602      DO jj = 1, nlcj
603         DO ji = 1, nlci
604            e2v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
605         END DO
606      END DO
607      CALL restget( inum, 'e2f', jpidta, jpjdta, 1, 0, llog, zdta )
608      DO jj = 1, nlcj
609         DO ji = 1, nlci
610            e2f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
611         END DO
612      END DO
613
614      CALL restclo( inum )
615
616      ! set extra rows add in mpp to none zero values
617      DO jj = nlcj+1, jpj
618         DO ji = 1, nlci
619            glamt(ji,jj) = glamt(ji,1)   ;   gphit(ji,jj) = gphit(ji,1)
620            glamu(ji,jj) = glamu(ji,1)   ;   gphiu(ji,jj) = gphiu(ji,1)
621            glamv(ji,jj) = glamv(ji,1)   ;   gphiv(ji,jj) = gphiv(ji,1)
622            glamf(ji,jj) = glamf(ji,1)   ;   gphif(ji,jj) = gphif(ji,1)
623            e1t  (ji,jj) = e1t  (ji,1)   ;   e2t  (ji,jj) = e2t  (ji,1)
624            e1u  (ji,jj) = e1u  (ji,1)   ;   e2u  (ji,jj) = e2u  (ji,1)
625            e1v  (ji,jj) = e1v  (ji,1)   ;   e2v  (ji,jj) = e2v  (ji,1)
626            e1f  (ji,jj) = e1f  (ji,1)   ;   e2f  (ji,jj) = e2f  (ji,1)
627         END DO
628      END DO
629
630      ! set extra columns add in mpp to none zero values
631      DO ji = nlci+1, jpi
632         glamt(ji,:) = glamt(1,:)   ;   gphit(ji,:) = gphit(1,:)
633         glamu(ji,:) = glamu(1,:)   ;   gphiu(ji,:) = gphiu(1,:)
634         glamv(ji,:) = glamv(1,:)   ;   gphiv(ji,:) = gphiv(1,:)
635         glamf(ji,:) = glamf(1,:)   ;   gphif(ji,:) = gphif(1,:)
636         e1t  (ji,:) = e1t  (1,:)   ;   e2t  (ji,:) = e2t  (1,:)
637         e1u  (ji,:) = e1u  (1,:)   ;   e2u  (ji,:) = e2u  (1,:)
638         e1v  (ji,:) = e1v  (1,:)   ;   e2v  (ji,:) = e2v  (1,:)
639         e1f  (ji,:) = e1f  (1,:)   ;   e2f  (ji,:) = e2f  (1,:)
640      END DO
641
642   END SUBROUTINE hgr_read
643
644   !!======================================================================
645END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.