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

Last change on this file since 1707 was 1707, checked in by rblod, 14 years ago

Fix bug for ORCA2 at Danish Straits, see ticket #586

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 28.5 KB
Line 
1MODULE domhgr
2   !!==============================================================================
3   !!                       ***  MODULE domhgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :       !  88-03  (G. Madec)
7   !!                 !  91-11  (G. Madec)
8   !!                 !  92-06  (M. Imbard)
9   !!                 !  96-01  (G. Madec)  terrain following coordinates
10   !!                 !  97-02  (G. Madec)  print mesh informations
11   !!                 !  99-11  (M. Imbard) NetCDF format with IO-IPSL
12   !!                 !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
13   !!                 !  01-09  (M. Levy)  eel config: grid in km, beta-plane
14   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module, namelist
15   !!            9.0  !  04-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
16   !!                           use of parameters in par_CONFIG-Rxx.h90, not in namelist
17   !!                 !  04-05  (A. Koch-Larrouy) Add Gyre configuration
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dom_hgr        : initialize the horizontal mesh
22   !!   hgr_read       : read "coordinate" NetCDF file
23   !!----------------------------------------------------------------------
24   !! * Modules used
25   USE dom_oce         ! ocean space and time domain
26   USE phycst          ! physical constants
27   USE in_out_manager  ! I/O manager
28   USE lib_mpp
29
30   IMPLICIT NONE
31   PRIVATE
32
33   !! * Module variables
34   REAL(wp) ::   glam0, gphi0           ! variables corresponding to parameters
35      !                                 ! ppglam0 ppgphi0 set in par_oce
36
37   !! * Routine accessibility
38   PUBLIC dom_hgr        ! called by domain.F90
39   !!----------------------------------------------------------------------
40   !!   OPA 9.0 , LOCEAN-IPSL (2005)
41   !! $Id$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE dom_hgr
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dom_hgr  ***
50      !!
51      !! ** Purpose :   Compute the geographical position (in degre) of the
52      !!      model grid-points,  the horizontal scale factors (in meters) and
53      !!      the Coriolis factor (in s-1).
54      !!
55      !! ** Method  :   The geographical position of the model grid-points is
56      !!      defined from analytical functions, fslam and fsphi, the deriva-
57      !!      tives of which gives the horizontal scale factors e1,e2.
58      !!      Defining two function fslam and fsphi and their derivatives in
59      !!      the two horizontal directions (fse1 and fse2), the model grid-
60      !!      point position and scale factors are given by:
61      !!         t-point:
62      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
63      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
64      !!         u-point:
65      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
66      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
67      !!         v-point:
68      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
69      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
70      !!            f-point:
71      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
72      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
73      !!      Where fse1 and fse2 are defined by:
74      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
75      !!                                     +          di(fsphi) **2 )(i,j)
76      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
77      !!                                     +          dj(fsphi) **2 )(i,j)
78      !!
79      !!        The coriolis factor is given at z-point by:
80      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
81      !!
82      !!        This routine is given as an example, it must be modified
83      !!      following the user s desiderata. nevertheless, the output as
84      !!      well as the way to compute the model grid-point position and
85      !!      horizontal scale factors must be respected in order to insure
86      !!      second order accuracy schemes.
87      !!
88      !! N.B. If the domain is periodic, verify that scale factors are also
89      !!      periodic, and the coriolis term again.
90      !!
91      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
92      !!                u-, v- and f-points (in degre)
93      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
94      !!               u-, v-  and f-points (in degre)
95      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
96      !!      scale factors (in meters) at t-, u-, v-, and f-points.
97      !!        define ff: coriolis factor at f-point
98      !!
99      !! References :   Marti, Madec and Delecluse, 1992, JGR
100      !!                Madec, Imbard, 1996, Clim. Dyn.
101      !!----------------------------------------------------------------------
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               !
147               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
148               ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
149                                               e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3
150               IF(lwp) WRITE(numout,*)
151               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km'
152               IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km'
153            ENDIF
154
155            ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km)
156            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
157            IF(lwp) WRITE(numout,*)
158            IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km'
159            !
160         ENDIF
161
162         !                                                ! ======================
163         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
164            !                                             ! ======================
165            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
166            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
167            IF(lwp) WRITE(numout,*)
168            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
169            !
170            ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km)
171            ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
172            IF(lwp) WRITE(numout,*)
173            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait'
174            !
175            ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km)
176            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3
177            IF(lwp) WRITE(numout,*)
178            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait'
179            !
180            ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km)
181            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
182            IF(lwp) WRITE(numout,*)
183            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait'
184            !
185            ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km)
186            ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
187            IF(lwp) WRITE(numout,*)
188            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait'
189            !
190            ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km)
191            ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
192            IF(lwp) WRITE(numout,*)
193            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait'
194            !
195            !
196            ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km)
197            ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3
198            IF(lwp) WRITE(numout,*)
199            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb'
200            !
201         ENDIF
202
203
204         ! N.B. :  General case, lat and long function of both i and j indices:
205         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
206         !                                  + (                           fsdiph( zti, ztj ) )**2  )
207         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
208         !                                  + (                           fsdiph( zui, zuj ) )**2  )
209         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
210         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
211         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
212         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
213         !
214         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
215         !                                  + (                           fsdjph( zti, ztj ) )**2  )
216         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
217         !                                  + (                           fsdjph( zui, zuj ) )**2  )
218         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
219         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
220         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
221         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
222
223
224      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
225
226         IF(lwp) WRITE(numout,*)
227         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
228         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
229
230         DO jj = 1, jpj
231            DO ji = 1, jpi
232               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
233               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
234               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
235               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
236         ! Longitude
237               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
238               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
239               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
240               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
241         ! Latitude
242               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
243               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
244               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
245               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
246         ! e1
247               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
248               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
249               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
250               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
251         ! e2
252               e2t(ji,jj) = ra * rad * ppe2_deg
253               e2u(ji,jj) = ra * rad * ppe2_deg
254               e2v(ji,jj) = ra * rad * ppe2_deg
255               e2f(ji,jj) = ra * rad * ppe2_deg
256            END DO
257         END DO
258
259
260      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
261
262         IF(lwp) WRITE(numout,*)
263         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
264         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
265
266         ! Position coordinates (in kilometers)
267         !                          ==========
268         glam0 = 0.e0
269         gphi0 = - ppe2_m * 1.e-3
270         
271#if defined key_agrif && defined key_eel_r6
272         IF (.Not.Agrif_Root()) THEN
273           glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
274           gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
275           ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
276           ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
277         ENDIF
278#endif         
279         DO jj = 1, jpj
280            DO ji = 1, jpi
281               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
282               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
283               glamv(ji,jj) = glamt(ji,jj)
284               glamf(ji,jj) = glamu(ji,jj)
285   
286               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
287               gphiu(ji,jj) = gphit(ji,jj)
288               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
289               gphif(ji,jj) = gphiv(ji,jj)
290            END DO
291         END DO
292
293         ! Horizontal scale factors (in meters)
294         !                              ======
295         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
296         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
297         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
298         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
299
300      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
301
302         IF(lwp) WRITE(numout,*)
303         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
304         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
305         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
306
307         !  Find index corresponding to the equator, given the grid spacing e1_deg
308         !  and the (approximate) southern latitude ppgphi0.
309         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
310         !  The formula should work even if the equator is outside the domain.
311         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
312         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
313         IF(  ppgphi0 > 0 )  ijeq = -ijeq
314
315         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
316
317         DO jj = 1, jpj
318            DO ji = 1, jpi
319               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
320               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
321               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
322               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
323         ! Longitude
324               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
325               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
326               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
327               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
328         ! Latitude
329               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
330               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
331               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
332               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
333         ! e1
334               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
335               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
336               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
337               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
338         ! e2
339               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
340               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
341               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
342               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
343            END DO
344         END DO
345
346      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
347
348         IF(lwp) WRITE(numout,*)
349         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
350         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
351
352         ! Position coordinates (in kilometers)
353         !                          ==========
354
355         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
356         zlam1 = -85
357         zphi1 = 29
358         ! resolution in meters
359         ze1 = 106000. / FLOAT(jp_cfg)           
360         ! benchmark: forced the resolution to be about 100 km
361         IF( nbench /= 0 )   ze1 = 106000.e0     
362         zsin_alpha = - SQRT( 2. ) / 2.
363         zcos_alpha =   SQRT( 2. ) / 2.
364         ze1deg = ze1 / (ra * rad)
365         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
366         !                                                          ! at the right jp_cfg resolution
367         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
368         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
369
370         IF( nprint==1 .AND. lwp )   THEN
371            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
372            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
373         ENDIF
374
375         DO jj = 1, jpj
376           DO ji = 1, jpi
377             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
378             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
379
380             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
381             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
382
383             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
384             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
385
386             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
387             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
388
389             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
390             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
391           END DO
392          END DO
393
394         ! Horizontal scale factors (in meters)
395         !                              ======
396         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
397         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
398         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
399         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
400
401      CASE DEFAULT
402         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
403         CALL ctl_stop( ctmp1 )
404
405      END SELECT
406
407
408      ! Control printing : Grid informations (if not restart)
409      ! ----------------
410
411      IF( lwp .AND. .NOT.ln_rstart ) THEN
412         WRITE(numout,*)
413         WRITE(numout,*) '          longitude and e1 scale factors'
414         WRITE(numout,*) '          ------------------------------'
415         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
416            glamv(ji,1), glamf(ji,1),   &
417            e1t(ji,1), e1u(ji,1),   &
418            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4199300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
420            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
421         
422         WRITE(numout,*)
423         WRITE(numout,*) '          latitude and e2 scale factors'
424         WRITE(numout,*) '          -----------------------------'
425         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
426            &                     gphiv(1,jj), gphif(1,jj),   &
427            &                     e2t  (1,jj), e2u  (1,jj),   &
428            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
429      ENDIF
430
431     
432      IF( nprint == 1 .AND. lwp ) THEN
433         WRITE(numout,*) '          e1u e2u '
434         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
435         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
436         WRITE(numout,*) '          e1v e2v  '
437         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
438         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
439         WRITE(numout,*) '          e1f e2f  '
440         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
441         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
442      ENDIF
443
444
445      ! ================= !
446      !  Coriolis factor  !
447      ! ================= !
448
449      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
450
451      CASE ( 0, 1, 4 )               ! mesh on the sphere
452
453         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
454
455      CASE ( 2 )                     ! f-plane at ppgphi0
456
457         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
458
459         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
460
461      CASE ( 3 )                     ! beta-plane
462
463         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
464         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
465         
466#if defined key_agrif && defined key_eel_r6
467         IF (.Not.Agrif_Root()) THEN
468           zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
469         ENDIF
470#endif         
471         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
472
473         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
474         
475         IF(lwp) THEN
476            WRITE(numout,*) 
477            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
478            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
479         ENDIF
480         IF( lk_mpp ) THEN
481            zminff=ff(nldi,nldj)
482            zmaxff=ff(nldi,nlej)
483            CALL mpp_min( zminff )   ! min over the global domain
484            CALL mpp_max( zmaxff )   ! max over the global domain
485            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
486         END IF
487
488      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
489
490         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
491         zphi0 = 15.e0                                                      ! latitude of the first row F-points
492         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
493
494         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
495
496         IF(lwp) THEN
497            WRITE(numout,*) 
498            WRITE(numout,*) '          Beta-plane and rotated domain : '
499            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
500         ENDIF
501
502         IF( lk_mpp ) THEN
503            zminff=ff(nldi,nldj)
504            zmaxff=ff(nldi,nlej)
505            CALL mpp_min( zminff )   ! min over the global domain
506            CALL mpp_max( zmaxff )   ! max over the global domain
507            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
508         END IF
509
510      END SELECT
511
512
513      ! Control of domain for symetrical condition
514      ! ------------------------------------------
515      ! The equator line must be the latitude coordinate axe
516
517      IF( nperio == 2 ) THEN
518         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
519         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
520      ENDIF
521
522   END SUBROUTINE dom_hgr
523
524
525   SUBROUTINE hgr_read
526      !!---------------------------------------------------------------------
527      !!              ***  ROUTINE hgr_read  ***
528      !!
529      !! ** Purpose :   Read a coordinate file in NetCDF format
530      !!
531      !! ** Method  :   The mesh file has been defined trough a analytical
532      !!      or semi-analytical method. It is read in a NetCDF file.
533      !!     
534      !!----------------------------------------------------------------------
535      USE iom
536
537      INTEGER ::   inum   ! temporary logical unit
538      !!----------------------------------------------------------------------
539
540      IF(lwp) THEN
541         WRITE(numout,*)
542         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
543         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
544      ENDIF
545     
546      CALL iom_open( 'coordinates', inum )
547     
548      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
549      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
550      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
551      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
552     
553      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
554      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
555      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
556      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
557     
558      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
559      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
560      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
561      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
562     
563      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
564      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
565      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
566      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
567     
568      CALL iom_close( inum )
569     
570    END SUBROUTINE hgr_read
571   
572   !!======================================================================
573END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.