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

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

RB:nemo_v1_update_038: first integration of Agrif :

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