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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

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