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/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 5755

Last change on this file since 5755 was 5755, checked in by flavoni, 9 years ago

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