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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago

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