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 NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domhgr.F90 @ 11129

Last change on this file since 11129 was 11129, checked in by mathiot, 5 years ago

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