source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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