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 utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/domhgr.F90 @ 14698

Last change on this file since 14698 was 14630, checked in by jchanut, 3 years ago

AGFdomcfg: 1) restore use of analytical grids (i.e. case where jphgr_msh>0). That may be useful to set up test cases with AGRIF. 2) Add random topography over a flat bottom if nn_bathy = -1 (in place of a Gaussian bump). This illustrates well where the interface matching and update are done. #2638
File size: 25.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   !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse
17   !!                                       add optional read of e1e2u & e1e2v
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dom_hgr       : initialize the horizontal mesh
22   !!   hgr_read      : read "coordinate" NetCDF file
23   !!----------------------------------------------------------------------
24   USE dom_oce        ! ocean space and time domain
25   USE phycst         ! physical constants
26   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
27   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30
31   IMPLICIT NONE
32   PRIVATE
33
34   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
35
36   PUBLIC   dom_hgr   ! called by domain.F90
37
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
40   !! $Id: domhgr.F90 6140 2015-12-21 11:35:23Z timgraham $
41   !! Software governed by the CeCILL licence     (./LICENSE)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE dom_hgr
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE dom_hgr  ***
48      !!
49      !! ** Purpose :   Compute the geographical position (in degre) of the
50      !!      model grid-points,  the horizontal scale factors (in meters) and
51      !!      the Coriolis factor (in s-1).
52      !!
53      !! ** Method  :   The geographical position of the model grid-points is
54      !!      defined from analytical functions, fslam and fsphi, the deriva-
55      !!      tives of which gives the horizontal scale factors e1,e2.
56      !!      Defining two function fslam and fsphi and their derivatives in
57      !!      the two horizontal directions (fse1 and fse2), the model grid-
58      !!      point position and scale factors are given by:
59      !!         t-point:
60      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
61      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
62      !!         u-point:
63      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
64      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
65      !!         v-point:
66      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
67      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
68      !!            f-point:
69      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
70      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
71      !!      Where fse1 and fse2 are defined by:
72      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
73      !!                                     +          di(fsphi) **2 )(i,j)
74      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
75      !!                                     +          dj(fsphi) **2 )(i,j)
76      !!
77      !!        The coriolis factor is given at z-point by:
78      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
79      !!
80      !!        This routine is given as an example, it must be modified
81      !!      following the user s desiderata. nevertheless, the output as
82      !!      well as the way to compute the model grid-point position and
83      !!      horizontal scale factors must be respected in order to insure
84      !!      second order accuracy schemes.
85      !!
86      !! N.B. If the domain is periodic, verify that scale factors are also
87      !!      periodic, and the coriolis term again.
88      !!
89      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
90      !!                u-, v- and f-points (in degre)
91      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
92      !!               u-, v-  and f-points (in degre)
93      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
94      !!      scale factors (in meters) at t-, u-, v-, and f-points.
95      !!        define ff: coriolis factor at f-point
96      !!
97      !! References :   Marti, Madec and Delecluse, 1992, JGR
98      !!                Madec, Imbard, 1996, Clim. Dyn.
99      !!----------------------------------------------------------------------
100      INTEGER  ::   ji, jj               ! dummy loop indices
101      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
102      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
103      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
104      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
105      REAL(wp) ::   zphi0, zbeta, znorme !
106      REAL(wp) ::   zarg, zf0, zminff, zmaxff
107      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
108      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
109      INTEGER  ::   isrow                ! index for ORCA1 starting row
110      INTEGER  ::   ie1e2u_v             ! fag for u- & v-surface read in coordinate file or not
111      !!----------------------------------------------------------------------
112      !
113      IF(lwp) THEN
114         WRITE(numout,*)
115         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
116         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
117         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
118         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
119         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
120         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
121         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
122         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
123      ENDIF
124      !
125      !
126#if defined key_agrif
127      IF (agrif_root()) THEN
128#endif
129      !
130      SELECT CASE( jphgr_msh )   !  type of horizontal mesh 
131      !
132      CASE ( 0 )                     !==  read in coordinate.nc file  ==!
133         !
134         IF(lwp) WRITE(numout,*)
135         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
136         !
137         ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v
138         !
139         CALL hgr_read( ie1e2u_v )     ! read the coordinate.nc file
140         !
141         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them
142            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction
143            e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
144            e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
145         ENDIF
146         !
147      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==!
148         !
149         IF(lwp) WRITE(numout,*)
150         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
151         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
152         !
153         DO jj = 1, jpj
154            DO ji = 1, jpi
155               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 )
156               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 )
157               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5
158               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5
159         ! Longitude
160               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
161               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
162               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
163               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
164         ! Latitude
165               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
166               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
167               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
168               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
169         ! e1
170               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
171               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
172               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
173               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
174         ! e2
175               e2t(ji,jj) = ra * rad * ppe2_deg
176               e2u(ji,jj) = ra * rad * ppe2_deg
177               e2v(ji,jj) = ra * rad * ppe2_deg
178               e2f(ji,jj) = ra * rad * ppe2_deg
179            END DO
180         END DO
181         !
182      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==!
183         !
184         IF(lwp) WRITE(numout,*)
185         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
186         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
187         !
188         ! Position coordinates (in kilometers)
189         !                          ==========
190         glam0 = 0._wp
191         gphi0 = - ppe2_m * 1.e-3
192         !
193         DO jj = 1, jpj
194            DO ji = 1, jpi
195               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
196               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
197               glamv(ji,jj) = glamt(ji,jj)
198               glamf(ji,jj) = glamu(ji,jj)
199               !
200               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
201               gphiu(ji,jj) = gphit(ji,jj)
202               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
203               gphif(ji,jj) = gphiv(ji,jj)
204            END DO
205         END DO
206         !
207         ! Horizontal scale factors (in meters)
208         !                              ======
209         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
210         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
211         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
212         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
213         !
214      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
215         !
216         IF(lwp) WRITE(numout,*)
217         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
218         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
219         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
220         !
221         !  Find index corresponding to the equator, given the grid spacing e1_deg
222         !  and the (approximate) southern latitude ppgphi0.
223         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
224         !  The formula should work even if the equator is outside the domain.
225         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
226         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
227         IF(  ppgphi0 > 0 )  ijeq = -ijeq
228         !
229         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
230         !
231         DO jj = 1, jpj
232            DO ji = 1, jpi
233               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
234               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
235               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
236               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
237         ! Longitude
238               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
239               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
240               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
241               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
242         ! Latitude
243               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
244               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
245               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
246               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
247         ! e1
248               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
249               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
250               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
251               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
252         ! e2
253               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
254               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
255               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
256               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
257            END DO
258         END DO
259         !
260      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
261         !
262         IF(lwp) WRITE(numout,*)
263         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
264         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
265         !
266         ! Position coordinates (in kilometers)
267         !                          ==========
268         !
269         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
270         zlam1 = -85._wp
271         zphi1 =  29._wp
272         ! resolution in meters
273         ze1 = 106000. / REAL( jp_cfg , wp )           
274         ! benchmark: forced the resolution to be about 100 km
275       !  IF( nbench /= 0 )   ze1 = 106000._wp     
276         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
277         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
278         ze1deg = ze1 / (ra * rad)
279       !  IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
280         !                                                           ! at the right jp_cfg resolution
281         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
282         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
283         !
284    !        WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
285    !        WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
286         !
287         DO jj = 1, jpj
288            DO ji = 1, jpi
289               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
290               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
291               !
292               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
293               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
294               !
295               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
296               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
297               !
298               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
299               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
300               !
301               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
302               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
303            END DO
304         END DO
305         !
306         ! Horizontal scale factors (in meters)
307         !                              ======
308         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
309         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
310         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
311         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
312         !
313      CASE DEFAULT
314         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
315         CALL ctl_stop( ctmp1 )
316         !
317      END SELECT
318     
319#if defined key_agrif
320      ELSE
321         CALL Agrif_InitValues_cont()
322      ENDIF
323#endif
324      ! associated horizontal metrics
325      ! -----------------------------
326      !
327      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
328      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
329      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
330      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
331      !
332      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
333      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
334      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
335         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
336         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
337      ENDIF
338      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
339      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
340      !   
341      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
342      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
343
344      IF( lwp ) THEN      ! Control print : Grid informations (if not restart)
345         WRITE(numout,*)
346         WRITE(numout,*) '          longitude and e1 scale factors'
347         WRITE(numout,*) '          ------------------------------'
348         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
349            glamv(ji,1), glamf(ji,1),   &
350            e1t(ji,1), e1u(ji,1),   &
351            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3529300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
353            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
354            !
355         WRITE(numout,*)
356         WRITE(numout,*) '          latitude and e2 scale factors'
357         WRITE(numout,*) '          -----------------------------'
358         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
359            &                     gphiv(1,jj), gphif(1,jj),   &
360            &                     e2t  (1,jj), e2u  (1,jj),   &
361            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
362      ENDIF
363
364
365      ! ================= !
366      !  Coriolis factor  !
367      ! ================= !
368
369      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
370      !
371      CASE ( 0, 1, 4 )               ! mesh on the sphere
372         !
373         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
374         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       - at t-point
375         !
376      CASE ( 2 )                     ! f-plane at ppgphi0
377         !
378         ff_f(:,:) = 2. * omega * SIN( rad * ppgphi0 )
379         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 )
380         !
381         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff_f(1,1)
382         !
383      CASE ( 3 )                     ! beta-plane
384         !
385         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
386         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
387         !
388         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
389         !
390         ff_f(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
391         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
392         !
393         IF(lwp) THEN
394            WRITE(numout,*) 
395            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff_f(Nis0,Njs0)
396            WRITE(numout,*) '          Coriolis parameter varies from ', ff_f(Nis0,Njs0),' to ', ff_f(Nis0,Nje0)
397         ENDIF
398         IF( lk_mpp ) THEN
399            zminff=ff_f(Nis0,Njs0)
400            zmaxff=ff_f(Nis0,Nje0)
401            CALL mpp_min( 'toto',zminff )   ! min over the global domain
402            CALL mpp_max( 'toto',zmaxff )   ! max over the global domain
403            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
404         END IF
405         !
406      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
407         !
408         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
409         zphi0 = 15._wp                                                     ! latitude of the first row F-points
410         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
411         !
412         ff_f(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
413         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
414         !
415         IF(lwp) THEN
416            WRITE(numout,*) 
417            WRITE(numout,*) '          Beta-plane and rotated domain : '
418            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff_f(Nis0,Njs0),' to ', ff_f(Nis0,Nje0)
419         ENDIF
420         !
421         IF( lk_mpp ) THEN
422            zminff=ff_f(Nis0,Njs0)
423            zmaxff=ff_f(Nis0,Nje0)
424            CALL mpp_min('toto', zminff )   ! min over the global domain
425            CALL mpp_max( 'toto',zmaxff )   ! max over the global domain
426            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
427         END IF
428         !
429      END SELECT
430
431
432      ! Control of domain for symetrical condition
433      ! ------------------------------------------
434      ! The equator line must be the latitude coordinate axe
435 !(PM) be carefull with nperio/jperio
436      IF( jperio == 2 ) THEN
437         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
438         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
439      ENDIF
440      !
441   END SUBROUTINE dom_hgr
442
443
444   SUBROUTINE hgr_read( ke1e2u_v )
445      !!---------------------------------------------------------------------
446      !!              ***  ROUTINE hgr_read  ***
447      !!
448      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
449      !!
450      !!----------------------------------------------------------------------
451      USE iom
452      !!
453      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0)
454      !
455      INTEGER ::   inum   ! temporary logical unit
456      CHARACTER(LEN=135) :: coordinate_filename
457      !!----------------------------------------------------------------------
458      !
459      IF(lwp) THEN
460         WRITE(numout,*)
461         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
462         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
463      ENDIF
464      !
465
466      IF (ln_read_cfg) THEN
467         coordinate_filename=TRIM(cn_domcfg)
468      ELSE
469         coordinate_filename=cn_fcoord
470      ENDIF
471      CALL iom_open( coordinate_filename, inum )
472      !
473      CALL iom_get( inum, jpdom_global, 'glamt', glamt, cd_type = 'T', psgn = 1._wp )
474      CALL iom_get( inum, jpdom_global, 'glamu', glamu, cd_type = 'U', psgn = 1._wp )
475      CALL iom_get( inum, jpdom_global, 'glamv', glamv, cd_type = 'V', psgn = 1._wp )
476      CALL iom_get( inum, jpdom_global, 'glamf', glamf, cd_type = 'F', psgn = 1._wp )
477      !
478      CALL iom_get( inum, jpdom_global, 'gphit', gphit, cd_type = 'T', psgn = 1._wp )
479      CALL iom_get( inum, jpdom_global, 'gphiu', gphiu, cd_type = 'U', psgn = 1._wp )
480      CALL iom_get( inum, jpdom_global, 'gphiv', gphiv, cd_type = 'V', psgn = 1._wp )
481      CALL iom_get( inum, jpdom_global, 'gphif', gphif, cd_type = 'F', psgn = 1._wp )
482      !
483      CALL iom_get( inum, jpdom_global, 'e1t'  , e1t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
484      CALL iom_get( inum, jpdom_global, 'e1u'  , e1u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
485      CALL iom_get( inum, jpdom_global, 'e1v'  , e1v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
486      CALL iom_get( inum, jpdom_global, 'e1f'  , e1f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
487      !
488      CALL iom_get( inum, jpdom_global, 'e2t'  , e2t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
489      CALL iom_get( inum, jpdom_global, 'e2u'  , e2u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
490      CALL iom_get( inum, jpdom_global, 'e2v'  , e2v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
491      CALL iom_get( inum, jpdom_global, 'e2f'  , e2f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
492      !
493      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
494         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
495         CALL iom_get( inum, jpdom_global, 'e1e2u', e1e2u, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
496         CALL iom_get( inum, jpdom_global, 'e1e2v', e1e2v, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
497         ke1e2u_v = 1
498      ELSE
499         ke1e2u_v = 0
500      ENDIF
501      !
502      CALL iom_close( inum )
503     
504    END SUBROUTINE hgr_read
505   
506   !!======================================================================
507END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.