source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 4018

Last change on this file since 4018 was 4018, checked in by clevy, 8 years ago

Configuration setting, bugfixes for AGRIF, see ticket:#1074

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