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

Last change on this file since 1057 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.3 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   !! $Header$
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               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            ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km)
168            ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
169            IF(lwp) WRITE(numout,*)
170            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait'
171            !
172            ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km)
173            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3
174            IF(lwp) WRITE(numout,*)
175            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait'
176            !
177            ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km)
178            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
179            IF(lwp) WRITE(numout,*)
180            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait'
181            !
182            ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km)
183            ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
184            IF(lwp) WRITE(numout,*)
185            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait'
186            !
187            ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km)
188            ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
189            IF(lwp) WRITE(numout,*)
190            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait'
191            !
192            !
193            ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km)
194            ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3
195            IF(lwp) WRITE(numout,*)
196            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb'
197            !
198         ENDIF
199
200
201         ! N.B. :  General case, lat and long function of both i and j indices:
202         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
203         !                                  + (                           fsdiph( zti, ztj ) )**2  )
204         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
205         !                                  + (                           fsdiph( zui, zuj ) )**2  )
206         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
207         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
208         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
209         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
210         !
211         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
212         !                                  + (                           fsdjph( zti, ztj ) )**2  )
213         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
214         !                                  + (                           fsdjph( zui, zuj ) )**2  )
215         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
216         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
217         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
218         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
219
220
221      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
222
223         IF(lwp) WRITE(numout,*)
224         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
225         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
226
227         DO jj = 1, jpj
228            DO ji = 1, jpi
229               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
230               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
231               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
232               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
233         ! Longitude
234               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
235               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
236               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
237               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
238         ! Latitude
239               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
240               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
241               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
242               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
243         ! e1
244               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
245               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
246               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
247               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
248         ! e2
249               e2t(ji,jj) = ra * rad * ppe2_deg
250               e2u(ji,jj) = ra * rad * ppe2_deg
251               e2v(ji,jj) = ra * rad * ppe2_deg
252               e2f(ji,jj) = ra * rad * ppe2_deg
253            END DO
254         END DO
255
256
257      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
258
259         IF(lwp) WRITE(numout,*)
260         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
261         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
262
263         ! Position coordinates (in kilometers)
264         !                          ==========
265         glam0 = 0.e0
266         gphi0 = - ppe2_m * 1.e-3
267         
268#if defined key_agrif && defined key_eel_r6
269         IF (.Not.Agrif_Root()) THEN
270           glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
271           gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
272           ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
273           ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
274         ENDIF
275#endif         
276         DO jj = 1, jpj
277            DO ji = 1, jpi
278               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
279               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
280               glamv(ji,jj) = glamt(ji,jj)
281               glamf(ji,jj) = glamu(ji,jj)
282   
283               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
284               gphiu(ji,jj) = gphit(ji,jj)
285               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
286               gphif(ji,jj) = gphiv(ji,jj)
287            END DO
288         END DO
289
290         ! Horizontal scale factors (in meters)
291         !                              ======
292         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
293         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
294         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
295         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
296
297      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
298
299         IF(lwp) WRITE(numout,*)
300         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
301         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
302         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
303
304         !  Find index corresponding to the equator, given the grid spacing e1_deg
305         !  and the (approximate) southern latitude ppgphi0.
306         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
307         !  The formula should work even if the equator is outside the domain.
308         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
309         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
310         IF(  ppgphi0 > 0 )  ijeq = -ijeq
311
312         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
313
314         DO jj = 1, jpj
315            DO ji = 1, jpi
316               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
317               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
318               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
319               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
320         ! Longitude
321               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
322               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
323               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
324               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
325         ! Latitude
326               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
327               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
328               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
329               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
330         ! e1
331               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
332               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
333               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
334               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
335         ! e2
336               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
337               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
338               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
339               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
340            END DO
341         END DO
342
343      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
344
345         IF(lwp) WRITE(numout,*)
346         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
347         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
348
349         ! Position coordinates (in kilometers)
350         !                          ==========
351
352         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
353         zlam1 = -85
354         zphi1 = 29
355         ! resolution in meters
356         ze1 = 106000. / FLOAT(jp_cfg)           
357         ! benchmark: forced the resolution to be about 100 km
358         IF( nbench /= 0 )   ze1 = 106000.e0     
359         zsin_alpha = - SQRT( 2. ) / 2.
360         zcos_alpha =   SQRT( 2. ) / 2.
361         ze1deg = ze1 / (ra * rad)
362         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
363         !                                                          ! at the right jp_cfg resolution
364         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
365         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
366
367         IF( nprint==1 .AND. lwp )   THEN
368            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
369            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
370         ENDIF
371
372         DO jj = 1, jpj
373           DO ji = 1, jpi
374             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
375             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
376
377             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
378             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
379
380             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
381             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
382
383             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
384             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
385
386             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
387             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
388           END DO
389          END DO
390
391         ! Horizontal scale factors (in meters)
392         !                              ======
393         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
394         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
395         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
396         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
397
398      CASE DEFAULT
399         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
400         CALL ctl_stop( ctmp1 )
401
402      END SELECT
403
404
405      ! Control printing : Grid informations (if not restart)
406      ! ----------------
407
408      IF( lwp .AND. .NOT.ln_rstart ) THEN
409         WRITE(numout,*)
410         WRITE(numout,*) '          longitude and e1 scale factors'
411         WRITE(numout,*) '          ------------------------------'
412         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
413            glamv(ji,1), glamf(ji,1),   &
414            e1t(ji,1), e1u(ji,1),   &
415            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4169300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
417            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
418         
419         WRITE(numout,*)
420         WRITE(numout,*) '          latitude and e2 scale factors'
421         WRITE(numout,*) '          -----------------------------'
422         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
423            &                     gphiv(1,jj), gphif(1,jj),   &
424            &                     e2t  (1,jj), e2u  (1,jj),   &
425            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
426      ENDIF
427
428     
429      IF( nprint == 1 .AND. lwp ) THEN
430         WRITE(numout,*) '          e1u e2u '
431         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
432         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
433         WRITE(numout,*) '          e1v e2v  '
434         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
435         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
436         WRITE(numout,*) '          e1f e2f  '
437         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
438         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
439      ENDIF
440
441
442      ! ================= !
443      !  Coriolis factor  !
444      ! ================= !
445
446      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
447
448      CASE ( 0, 1, 4 )               ! mesh on the sphere
449
450         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
451
452      CASE ( 2 )                     ! f-plane at ppgphi0
453
454         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
455
456         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
457
458      CASE ( 3 )                     ! beta-plane
459
460         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
461         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
462         
463#if defined key_agrif && defined key_eel_r6
464         IF (.Not.Agrif_Root()) THEN
465           zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
466         ENDIF
467#endif         
468         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
469
470         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
471         
472         IF(lwp) THEN
473            WRITE(numout,*) 
474            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
475            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
476         ENDIF
477         IF( lk_mpp ) THEN
478            zminff=ff(nldi,nldj)
479            zmaxff=ff(nldi,nlej)
480            CALL mpp_min( zminff )   ! min over the global domain
481            CALL mpp_max( zmaxff )   ! max over the global domain
482            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
483         END IF
484
485      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
486
487         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
488         zphi0 = 15.e0                                                      ! latitude of the first row F-points
489         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
490
491         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
492
493         IF(lwp) THEN
494            WRITE(numout,*) 
495            WRITE(numout,*) '          Beta-plane and rotated domain : '
496            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
497         ENDIF
498
499         IF( lk_mpp ) THEN
500            zminff=ff(nldi,nldj)
501            zmaxff=ff(nldi,nlej)
502            CALL mpp_min( zminff )   ! min over the global domain
503            CALL mpp_max( zmaxff )   ! max over the global domain
504            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
505         END IF
506
507      END SELECT
508
509
510      ! Control of domain for symetrical condition
511      ! ------------------------------------------
512      ! The equator line must be the latitude coordinate axe
513
514      IF( nperio == 2 ) THEN
515         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
516         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
517      ENDIF
518
519   END SUBROUTINE dom_hgr
520
521
522   SUBROUTINE hgr_read
523      !!---------------------------------------------------------------------
524      !!              ***  ROUTINE hgr_read  ***
525      !!
526      !! ** Purpose :   Read a coordinate file in NetCDF format
527      !!
528      !! ** Method  :   The mesh file has been defined trough a analytical
529      !!      or semi-analytical method. It is read in a NetCDF file.
530      !!     
531      !!----------------------------------------------------------------------
532      USE iom
533
534      INTEGER ::   inum   ! temporary logical unit
535      !!----------------------------------------------------------------------
536
537      IF(lwp) THEN
538         WRITE(numout,*)
539         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
540         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
541      ENDIF
542     
543      CALL iom_open( 'coordinates', inum )
544     
545      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
546      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
547      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
548      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
549     
550      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
551      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
552      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
553      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
554     
555      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
556      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
557      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
558      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
559     
560      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
561      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
562      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
563      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
564     
565      CALL iom_close( inum )
566     
567    END SUBROUTINE hgr_read
568   
569   !!======================================================================
570END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.