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 @ 13204

Last change on this file since 13204 was 13204, checked in by smasson, 4 years ago

tools: update with tools_dev_r12970_AGRIF_CMEMS
File size: 25.3 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      SELECT CASE( jphgr_msh )   !  type of horizontal mesh 
127      !
128      CASE ( 0 )                     !==  read in coordinate.nc file  ==!
129         !
130#if defined key_agrif
131         IF (agrif_root()) THEN
132#endif
133         IF(lwp) WRITE(numout,*)
134         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
135         !
136         ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v
137         !
138         CALL hgr_read( ie1e2u_v )     ! read the coordinate.nc file
139         !
140         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them
141            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction
142            e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
143            e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
144         ENDIF
145#if defined key_agrif
146         ELSE
147           CALL Agrif_InitValues_cont()
148         ENDIF
149#endif
150         !
151      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==!
152         !
153         IF(lwp) WRITE(numout,*)
154         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
155         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
156         !
157         DO jj = 1, jpj
158            DO ji = 1, jpi
159               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 )
160               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 )
161               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5
162               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5
163         ! Longitude
164               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
165               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
166               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
167               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
168         ! Latitude
169               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
170               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
171               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
172               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
173         ! e1
174               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
175               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
176               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
177               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
178         ! e2
179               e2t(ji,jj) = ra * rad * ppe2_deg
180               e2u(ji,jj) = ra * rad * ppe2_deg
181               e2v(ji,jj) = ra * rad * ppe2_deg
182               e2f(ji,jj) = ra * rad * ppe2_deg
183            END DO
184         END DO
185         !
186      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==!
187         !
188         IF(lwp) WRITE(numout,*)
189         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
190         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
191         !
192         ! Position coordinates (in kilometers)
193         !                          ==========
194         glam0 = 0._wp
195         gphi0 = - ppe2_m * 1.e-3
196         !
197         DO jj = 1, jpj
198            DO ji = 1, jpi
199               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
200               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
201               glamv(ji,jj) = glamt(ji,jj)
202               glamf(ji,jj) = glamu(ji,jj)
203               !
204               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
205               gphiu(ji,jj) = gphit(ji,jj)
206               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
207               gphif(ji,jj) = gphiv(ji,jj)
208            END DO
209         END DO
210         !
211         ! Horizontal scale factors (in meters)
212         !                              ======
213         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
214         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
215         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
216         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
217         !
218      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
219         !
220         IF(lwp) WRITE(numout,*)
221         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
222         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
223         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
224         !
225         !  Find index corresponding to the equator, given the grid spacing e1_deg
226         !  and the (approximate) southern latitude ppgphi0.
227         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
228         !  The formula should work even if the equator is outside the domain.
229         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
230         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
231         IF(  ppgphi0 > 0 )  ijeq = -ijeq
232         !
233         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
234         !
235         DO jj = 1, jpj
236            DO ji = 1, jpi
237               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
238               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
239               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
240               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
241         ! Longitude
242               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
243               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
244               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
245               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
246         ! Latitude
247               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
248               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
249               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
250               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
251         ! e1
252               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
253               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
254               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
255               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
256         ! e2
257               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
258               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
259               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
260               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
261            END DO
262         END DO
263         !
264      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
265         !
266         IF(lwp) WRITE(numout,*)
267         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
268         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
269         !
270         ! Position coordinates (in kilometers)
271         !                          ==========
272         !
273         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
274         zlam1 = -85._wp
275         zphi1 =  29._wp
276         ! resolution in meters
277         ze1 = 106000. / REAL( jp_cfg , wp )           
278         ! benchmark: forced the resolution to be about 100 km
279       !  IF( nbench /= 0 )   ze1 = 106000._wp     
280         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
281         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
282         ze1deg = ze1 / (ra * rad)
283       !  IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
284         !                                                           ! at the right jp_cfg resolution
285         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
286         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
287         !
288    !        WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
289    !        WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
290         !
291         DO jj = 1, jpj
292            DO ji = 1, jpi
293               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
294               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
295               !
296               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
297               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
298               !
299               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
300               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
301               !
302               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
303               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
304               !
305               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
306               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
307            END DO
308         END DO
309         !
310         ! Horizontal scale factors (in meters)
311         !                              ======
312         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
313         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
314         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
315         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
316         !
317      CASE DEFAULT
318         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
319         CALL ctl_stop( ctmp1 )
320         !
321      END SELECT
322     
323      ! associated horizontal metrics
324      ! -----------------------------
325      !
326      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
327      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
328      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
329      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
330      !
331      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
332      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
333      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
334         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
335         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
336      ENDIF
337      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
338      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
339      !   
340      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
341      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
342
343      IF( lwp ) THEN      ! Control print : Grid informations (if not restart)
344         WRITE(numout,*)
345         WRITE(numout,*) '          longitude and e1 scale factors'
346         WRITE(numout,*) '          ------------------------------'
347         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
348            glamv(ji,1), glamf(ji,1),   &
349            e1t(ji,1), e1u(ji,1),   &
350            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3519300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
352            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
353            !
354         WRITE(numout,*)
355         WRITE(numout,*) '          latitude and e2 scale factors'
356         WRITE(numout,*) '          -----------------------------'
357         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
358            &                     gphiv(1,jj), gphif(1,jj),   &
359            &                     e2t  (1,jj), e2u  (1,jj),   &
360            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
361      ENDIF
362
363
364      ! ================= !
365      !  Coriolis factor  !
366      ! ================= !
367
368      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
369      !
370      CASE ( 0, 1, 4 )               ! mesh on the sphere
371         !
372         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
373         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       - at t-point
374         !
375      CASE ( 2 )                     ! f-plane at ppgphi0
376         !
377         ff_f(:,:) = 2. * omega * SIN( rad * ppgphi0 )
378         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 )
379         !
380         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff_f(1,1)
381         !
382      CASE ( 3 )                     ! beta-plane
383         !
384         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
385         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
386         !
387         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
388         !
389         ff_f(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
390         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
391         !
392         IF(lwp) THEN
393            WRITE(numout,*) 
394            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff_f(nldi,nldj)
395            WRITE(numout,*) '          Coriolis parameter varies from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej)
396         ENDIF
397         IF( lk_mpp ) THEN
398            zminff=ff_f(nldi,nldj)
399            zmaxff=ff_f(nldi,nlej)
400            CALL mpp_min( 'toto',zminff )   ! min over the global domain
401            CALL mpp_max( 'toto',zmaxff )   ! max over the global domain
402            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
403         END IF
404         !
405      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
406         !
407         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
408         zphi0 = 15._wp                                                     ! latitude of the first row F-points
409         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
410         !
411         ff_f(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
412         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
413         !
414         IF(lwp) THEN
415            WRITE(numout,*) 
416            WRITE(numout,*) '          Beta-plane and rotated domain : '
417            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej)
418         ENDIF
419         !
420         IF( lk_mpp ) THEN
421            zminff=ff_f(nldi,nldj)
422            zmaxff=ff_f(nldi,nlej)
423            CALL mpp_min('toto', zminff )   ! min over the global domain
424            CALL mpp_max( 'toto',zmaxff )   ! max over the global domain
425            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
426         END IF
427         !
428      END SELECT
429
430
431      ! Control of domain for symetrical condition
432      ! ------------------------------------------
433      ! The equator line must be the latitude coordinate axe
434 !(PM) be carefull with nperio/jperio
435      IF( jperio == 2 ) THEN
436         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
437         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
438      ENDIF
439      !
440   END SUBROUTINE dom_hgr
441
442
443   SUBROUTINE hgr_read( ke1e2u_v )
444      !!---------------------------------------------------------------------
445      !!              ***  ROUTINE hgr_read  ***
446      !!
447      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
448      !!
449      !!----------------------------------------------------------------------
450      USE iom
451      !!
452      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0)
453      !
454      INTEGER ::   inum   ! temporary logical unit
455      CHARACTER(LEN=135) :: coordinate_filename
456      !!----------------------------------------------------------------------
457      !
458      IF(lwp) THEN
459         WRITE(numout,*)
460         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
461         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
462      ENDIF
463      !
464      IF (ln_read_cfg) THEN
465         coordinate_filename=TRIM(cn_domcfg)
466      ELSE
467         coordinate_filename='coordinates'
468      ENDIF
469      CALL iom_open( coordinate_filename, inum )
470      !
471      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
472      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
473      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
474      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
475      !
476      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
477      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
478      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
479      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
480      !
481      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
482      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
483      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
484      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
485      !
486      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
487      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
488      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
489      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
490      !
491      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
492         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
493         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
494         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
495         ke1e2u_v = 1
496      ELSE
497         ke1e2u_v = 0
498      ENDIF
499      !
500      CALL iom_close( inum )
501     
502    END SUBROUTINE hgr_read
503   
504   !!======================================================================
505END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.