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

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

nemo_v1_update_060: SM: IOM + 301 levels + CORE + begining of ctl_stop
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.1 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(lwp) WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
368         IF(lwp) WRITE(numout,*) 'ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
369
370         DO jj = 1, jpj
371           DO ji = 1, jpi
372             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
373             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
374
375             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
376             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
377
378             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
379             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
380
381             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
382             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
383
384             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
385             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
386           END DO
387          END DO
388
389         ! Horizontal scale factors (in meters)
390         !                              ======
391         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
392         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
393         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
394         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
395
396      CASE DEFAULT
397         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
398         CALL ctl_stop( ctmp1 )
399
400      END SELECT
401
402
403      ! Control printing : Grid informations (if not restart)
404      ! ----------------
405
406      IF(lwp .AND. .NOT.ln_rstart ) THEN
407         WRITE(numout,*)
408         WRITE(numout,*) '          longitude and e1 scale factors'
409         WRITE(numout,*) '          ------------------------------'
410         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
411            glamv(ji,1), glamf(ji,1),   &
412            e1t(ji,1), e1u(ji,1),   &
413            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4149300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
415            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
416         
417         WRITE(numout,*)
418         WRITE(numout,*) '          latitude and e2 scale factors'
419         WRITE(numout,*) '          -----------------------------'
420         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
421            &                     gphiv(1,jj), gphif(1,jj),   &
422            &                     e2t  (1,jj), e2u  (1,jj),   &
423            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
424      ENDIF
425
426     
427      IF( nprint == 1 .AND. lwp ) THEN
428         WRITE(numout,*) '          e1u e2u '
429         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
430         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
431         WRITE(numout,*) '          e1v e2v  '
432         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
433         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
434         WRITE(numout,*) '          e1f e2f  '
435         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
436         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
437      ENDIF
438
439
440      ! ================= !
441      !  Coriolis factor  !
442      ! ================= !
443
444      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
445
446      CASE ( 0, 1, 4 )               ! mesh on the sphere
447
448         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
449
450      CASE ( 2 )                     ! f-plane at ppgphi0
451
452         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
453
454         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
455
456      CASE ( 3 )                     ! beta-plane
457
458         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
459         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
460         
461#if defined key_agrif && defined key_eel_r6
462         IF (.Not.Agrif_Root()) THEN
463           zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
464         ENDIF
465#endif         
466         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
467
468         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
469         
470         IF(lwp) WRITE(numout,*) 
471         IF(lwp) WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
472         IF(lwp) WRITE(numout,*) ' Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
473         IF( lk_mpp ) THEN
474            zminff=ff(nldi,nldj)
475            zmaxff=ff(nldi,nlej)
476            CALL mpp_min( zminff )   ! min over the global domain
477            CALL mpp_max( zmaxff )   ! max over the global domain
478            IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff
479         END IF
480
481      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
482
483         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
484         zphi0 = 15.e0                                                      ! latitude of the first row F-points
485         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
486
487         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
488
489         IF(lwp) WRITE(numout,*) ' Beta-plane and rotated domain : '
490         IF(lwp) WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
491         IF( lk_mpp ) THEN
492            zminff=ff(nldi,nldj)
493            zmaxff=ff(nldi,nlej)
494            CALL mpp_min( zminff )   ! min over the global domain
495            CALL mpp_max( zmaxff )   ! max over the global domain
496            IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff
497         END IF
498
499      END SELECT
500
501
502      ! Control of domain for symetrical condition
503      ! ------------------------------------------
504      ! The equator line must be the latitude coordinate axe
505
506      IF( nperio == 2 ) THEN
507         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
508         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
509      ENDIF
510
511   END SUBROUTINE dom_hgr
512
513
514   SUBROUTINE hgr_read
515      !!---------------------------------------------------------------------
516      !!              ***  ROUTINE hgr_read  ***
517      !!
518      !! ** Purpose :   Read a coordinate file in NetCDF format
519      !!
520      !! ** Method  :   The mesh file has been defined trough a analytical
521      !!      or semi-analytical method. It is read in a NetCDF file.
522      !!     
523      !!----------------------------------------------------------------------
524      USE iom
525
526      INTEGER ::   inum   ! temporary logical unit
527      !!----------------------------------------------------------------------
528
529      IF(lwp) THEN
530         WRITE(numout,*)
531         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
532         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
533      ENDIF
534     
535      CALL iom_open( 'coordinates', inum )
536     
537      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
538      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
539      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
540      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
541     
542      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
543      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
544      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
545      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
546     
547      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
548      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
549      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
550      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
551     
552      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
553      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
554      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
555      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
556     
557      CALL iom_close( inum )
558     
559    END SUBROUTINE hgr_read
560   
561   !!======================================================================
562END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.