source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 2380

Last change on this file since 2380 was 2380, checked in by acc, 10 years ago

nemo_v3_3beta. ORCA_R1 settings (step 2, see ticket #758). Introduces key_orca_r1 (46 level default, 75 level if key_orca_r1=75)

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