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

Last change on this file since 7512 was 7508, checked in by mocavero, 7 years ago

changes on code duplication and workshare construct
  • Property svn:keywords set to Id
File size: 27.0 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 DO schedule(static) private(jj, ji)
143            DO jj = 1, jpj
144               DO ji = 1, jpi
145                  e1e2u (ji,jj) = e1u(ji,jj) * e2u(ji,jj)   
146                  e1e2v (ji,jj) = e1v(ji,jj) * e2v(ji,jj) 
147               END DO
148            END DO
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!$OMP PARALLEL DO schedule(static) private(jj, ji, zti, zui, zvi, zfi, ztj, zuj, zvj, zfj)
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#if defined key_agrif 
199         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
200            IF( .NOT. Agrif_Root() ) THEN
201              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
202              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
203              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
204              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
205            ENDIF
206         ENDIF
207#endif         
208!$OMP PARALLEL
209!$OMP DO schedule(static) private(jj, ji)
210         DO jj = 1, jpj
211            DO ji = 1, jpi
212               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
213               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
214               glamv(ji,jj) = glamt(ji,jj)
215               glamf(ji,jj) = glamu(ji,jj)
216               !
217               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
218               gphiu(ji,jj) = gphit(ji,jj)
219               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
220               gphif(ji,jj) = gphiv(ji,jj)
221            END DO
222         END DO
223!$OMP END DO NOWAIT
224         !
225         ! Horizontal scale factors (in meters)
226         !                              ======
227!$OMP DO schedule(static) private(jj, ji)
228         DO jj = 1, jpj
229            DO ji = 1, jpi
230               e1t(ji,jj) = ppe1_m      ;      e2t(ji,jj) = ppe2_m
231               e1u(ji,jj) = ppe1_m      ;      e2u(ji,jj) = ppe2_m
232               e1v(ji,jj) = ppe1_m      ;      e2v(ji,jj) = ppe2_m
233               e1f(ji,jj) = ppe1_m      ;      e2f(ji,jj) = ppe2_m
234            END DO
235         END DO
236!$OMP END DO NOWAIT
237!$OMP END PARALLEL
238         !
239      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
240         !
241         IF(lwp) WRITE(numout,*)
242         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
243         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
244         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
245         !
246         !  Find index corresponding to the equator, given the grid spacing e1_deg
247         !  and the (approximate) southern latitude ppgphi0.
248         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
249         !  The formula should work even if the equator is outside the domain.
250         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
251         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
252         IF(  ppgphi0 > 0 )  ijeq = -ijeq
253         !
254         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
255         !
256!$OMP PARALLEL DO schedule(static) private(jj, ji, zti, zui, zvi, zfi, ztj, zuj, zvj, zfj)
257         DO jj = 1, jpj
258            DO ji = 1, jpi
259               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
260               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
261               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
262               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
263         ! Longitude
264               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
265               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
266               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
267               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
268         ! Latitude
269               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
270               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
271               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
272               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
273         ! e1
274               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
275               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
276               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
277               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
278         ! e2
279               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
280               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
281               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
282               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
283            END DO
284         END DO
285         !
286      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
287         !
288         IF(lwp) WRITE(numout,*)
289         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
290         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
291         !
292         ! Position coordinates (in kilometers)
293         !                          ==========
294         !
295         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
296         zlam1 = -85._wp
297         zphi1 =  29._wp
298         ! resolution in meters
299         ze1 = 106000. / REAL( jp_cfg , wp )           
300         ! benchmark: forced the resolution to be about 100 km
301         IF( nbench /= 0 )   ze1 = 106000._wp     
302         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
303         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
304         ze1deg = ze1 / (ra * rad)
305         IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
306         !                                                           ! at the right jp_cfg resolution
307         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
308         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
309         !
310         IF( nprint==1 .AND. lwp )   THEN
311            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
312            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
313         ENDIF
314         !
315!$OMP PARALLEL
316!$OMP DO schedule(static) private(jj, ji, zim1, zjm1, zim05, zjm05)
317         DO jj = 1, jpj
318            DO ji = 1, jpi
319               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
320               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
321               !
322               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
323               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
324               !
325               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
326               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
327               !
328               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
329               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
330               !
331               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
332               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
333            END DO
334         END DO
335!$OMP END DO NOWAIT
336         !
337         ! Horizontal scale factors (in meters)
338         !                              ======
339!$OMP DO schedule(static) private(jj, ji)
340         DO jj = 1, jpj
341            DO ji = 1, jpi
342               e1t(ji,jj) =  ze1     ;      e2t(ji,jj) = ze1
343               e1u(ji,jj) =  ze1     ;      e2u(ji,jj) = ze1
344               e1v(ji,jj) =  ze1     ;      e2v(ji,jj) = ze1
345               e1f(ji,jj) =  ze1     ;      e2f(ji,jj) = ze1
346            END DO
347         END DO
348!$OMP END PARALLEL
349         !
350      CASE DEFAULT
351         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
352         CALL ctl_stop( ctmp1 )
353         !
354      END SELECT
355     
356      ! associated horizontal metrics
357      ! -----------------------------
358      !
359!$OMP PARALLEL DO schedule(static) private(jj, ji)
360      DO jj = 1, jpj
361         DO ji = 1, jpi
362            r1_e1t(ji,jj) = 1._wp / e1t(ji,jj)   ;   r1_e2t (ji,jj) = 1._wp / e2t(ji,jj)
363            r1_e1u(ji,jj) = 1._wp / e1u(ji,jj)   ;   r1_e2u (ji,jj) = 1._wp / e2u(ji,jj)
364            r1_e1v(ji,jj) = 1._wp / e1v(ji,jj)   ;   r1_e2v (ji,jj) = 1._wp / e2v(ji,jj)
365            r1_e1f(ji,jj) = 1._wp / e1f(ji,jj)   ;   r1_e2f (ji,jj) = 1._wp / e2f(ji,jj)
366      !
367            e1e2t (ji,jj) = e1t(ji,jj) * e2t(ji,jj)   ;   r1_e1e2t(ji,jj) = 1._wp / e1e2t(ji,jj)
368            e1e2f (ji,jj) = e1f(ji,jj) * e2f(ji,jj)   ;   r1_e1e2f(ji,jj) = 1._wp / e1e2f(ji,jj)
369         END DO
370      END DO
371
372      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
373!$OMP PARALLEL DO schedule(static) private(jj, ji)
374         DO jj = 1, jpj
375            DO ji = 1, jpi
376               e1e2u (ji,jj) = e1u(ji,jj) * e2u(ji,jj)   
377               e1e2v (ji,jj) = e1v(ji,jj) * e2v(ji,jj) 
378            END DO
379         END DO
380      ENDIF
381!$OMP PARALLEL DO schedule(static) private(jj, ji)
382         DO jj = 1, jpj
383            DO ji = 1, jpi
384               r1_e1e2u(ji,jj) = 1._wp / e1e2u(ji,jj)     ! compute their invert in both cases
385               r1_e1e2v(ji,jj) = 1._wp / e1e2v(ji,jj)
386               !   
387               e2_e1u(ji,jj) = e2u(ji,jj) / e1u(ji,jj)
388               e1_e2v(ji,jj) = e1v(ji,jj) / e2v(ji,jj)
389            END DO
390         END DO
391
392      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
393         WRITE(numout,*)
394         WRITE(numout,*) '          longitude and e1 scale factors'
395         WRITE(numout,*) '          ------------------------------'
396         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
397            glamv(ji,1), glamf(ji,1),   &
398            e1t(ji,1), e1u(ji,1),   &
399            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4009300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
401            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
402            !
403         WRITE(numout,*)
404         WRITE(numout,*) '          latitude and e2 scale factors'
405         WRITE(numout,*) '          -----------------------------'
406         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
407            &                     gphiv(1,jj), gphif(1,jj),   &
408            &                     e2t  (1,jj), e2u  (1,jj),   &
409            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
410      ENDIF
411
412
413      ! ================= !
414      !  Coriolis factor  !
415      ! ================= !
416
417      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
418      !
419      CASE ( 0, 1, 4 )               ! mesh on the sphere
420         !
421         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
422         !
423      CASE ( 2 )                     ! f-plane at ppgphi0
424         !
425         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
426         !
427         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
428         !
429      CASE ( 3 )                     ! beta-plane
430         !
431         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
432         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
433         !
434#if defined key_agrif
435         IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN       ! for EEL6 configuration only
436            IF( .NOT.Agrif_Root() ) THEN
437              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
438            ENDIF
439         ENDIF
440#endif         
441         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
442         !
443         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
444         !
445         IF(lwp) THEN
446            WRITE(numout,*) 
447            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
448            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
449         ENDIF
450         IF( lk_mpp ) THEN
451            zminff=ff(nldi,nldj)
452            zmaxff=ff(nldi,nlej)
453            CALL mpp_min( zminff )   ! min over the global domain
454            CALL mpp_max( zmaxff )   ! max over the global domain
455            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
456         END IF
457         !
458      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
459         !
460         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
461         zphi0 = 15._wp                                                     ! latitude of the first row F-points
462         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
463         !
464         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
465         !
466         IF(lwp) THEN
467            WRITE(numout,*) 
468            WRITE(numout,*) '          Beta-plane and rotated domain : '
469            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
470         ENDIF
471         !
472         IF( lk_mpp ) THEN
473            zminff=ff(nldi,nldj)
474            zmaxff=ff(nldi,nlej)
475            CALL mpp_min( zminff )   ! min over the global domain
476            CALL mpp_max( zmaxff )   ! max over the global domain
477            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
478         END IF
479         !
480      END SELECT
481
482
483      ! Control of domain for symetrical condition
484      ! ------------------------------------------
485      ! The equator line must be the latitude coordinate axe
486
487      IF( nperio == 2 ) THEN
488         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
489         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
490      ENDIF
491      !
492      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
493      !
494   END SUBROUTINE dom_hgr
495
496
497   SUBROUTINE hgr_read( ke1e2u_v )
498      !!---------------------------------------------------------------------
499      !!              ***  ROUTINE hgr_read  ***
500      !!
501      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
502      !!
503      !!----------------------------------------------------------------------
504      USE iom
505      !!
506      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0)
507      !
508      INTEGER ::   inum   ! temporary logical unit
509      !!----------------------------------------------------------------------
510      !
511      IF(lwp) THEN
512         WRITE(numout,*)
513         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
514         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
515      ENDIF
516      !
517      CALL iom_open( 'coordinates', inum )
518      !
519      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
520      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
521      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
522      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
523      !
524      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
525      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
526      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
527      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
528      !
529      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
530      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
531      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
532      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
533      !
534      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
535      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
536      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
537      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
538      !
539      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
540         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
541         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
542         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
543         ke1e2u_v = 1
544      ELSE
545         ke1e2u_v = 0
546      ENDIF
547      !
548      CALL iom_close( inum )
549     
550    END SUBROUTINE hgr_read
551   
552   !!======================================================================
553END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.