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/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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