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

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

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