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

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 6434

Last change on this file since 6434 was 6434, checked in by flavoni, 8 years ago

update usr_def module: add explicit inputs and outputs arguments, see ticket 1692
  • Property svn:keywords set to Id
File size: 24.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   !!             -   ! 2016-04  (S. Flavoni) change configuration's interface:
19   !!                            read file or CALL usr_def module to compute
20   !!                            horizontal grid (example given for GYRE)
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dom_hgr       : initialize the horizontal mesh
25   !!   hgr_read      : read "coordinate" NetCDF file
26   !!----------------------------------------------------------------------
27   USE dom_oce        ! ocean space and time domain
28   USE par_oce        ! ocean space and time domain
29   USE phycst         ! physical constants
30   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
31   !
32   USE in_out_manager ! I/O manager
33   USE lib_mpp        ! MPP library
34   USE timing         ! Timing
35   USE usrdef        ! User defined routine
36
37   IMPLICIT NONE
38   PRIVATE
39
40   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in namelist !SF modify some routines????
41
42   PUBLIC   dom_hgr   ! called by domain.F90
43
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE dom_hgr
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE dom_hgr  ***
54      !!
55      !! ** Purpose :   Compute the geographical position (in degre) of the
56      !!      model grid-points,  the horizontal scale factors (in meters) and
57      !!      the Coriolis factor (in s-1).
58      !!
59      !! ** Method  :   The geographical position of the model grid-points is
60      !!      defined from analytical functions, fslam and fsphi, the deriva-
61      !!      tives of which gives the horizontal scale factors e1,e2.
62      !!      Defining two function fslam and fsphi and their derivatives in
63      !!      the two horizontal directions (fse1 and fse2), the model grid-
64      !!      point position and scale factors are given by:
65      !!         t-point:
66      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
67      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
68      !!         u-point:
69      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
70      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
71      !!         v-point:
72      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
73      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
74      !!            f-point:
75      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
76      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
77      !!      Where fse1 and fse2 are defined by:
78      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
79      !!                                     +          di(fsphi) **2 )(i,j)
80      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
81      !!                                     +          dj(fsphi) **2 )(i,j)
82      !!
83      !!        The coriolis factor is given at z-point by:
84      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
85      !!
86      !!        This routine is given as an example, it must be modified
87      !!      following the user s desiderata. nevertheless, the output as
88      !!      well as the way to compute the model grid-point position and
89      !!      horizontal scale factors must be respected in order to insure
90      !!      second order accuracy schemes.
91      !!
92      !! N.B. If the domain is periodic, verify that scale factors are also
93      !!      periodic, and the coriolis term again.
94      !!
95      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
96      !!                u-, v- and f-points (in degre)
97      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
98      !!               u-, v-  and f-points (in degre)
99      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
100      !!      scale factors (in meters) at t-, u-, v-, and f-points.
101      !!        define ff: coriolis factor at f-point
102      !!
103      !! References :   Marti, Madec and Delecluse, 1992, JGR
104      !!                Madec, Imbard, 1996, Clim. Dyn.
105      !!----------------------------------------------------------------------
106      INTEGER  ::   ji, jj                   ! dummy loop indices
107      INTEGER  ::   ii0, ii1, ij0, ij1, iff  ! temporary integers
108      INTEGER  ::   ijeq                     ! index of equator T point (used in case 4)
109      REAL(wp) ::   zti, zui, zvi, zfi       ! local scalars
110      REAL(wp) ::   ztj, zuj, zvj, zfj       !   -      -
111      REAL(wp) ::   zphi0, zbeta, znorme     !
112      REAL(wp) ::   zarg, zf0, zminff, zmaxff
113      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
114      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
115      INTEGER  ::   isrow                    ! index for ORCA1 starting row
116      INTEGER  ::   ie1e2u_v                 ! fag for u- & v-surface read in coordinate file or not
117      !!----------------------------------------------------------------------
118      !
119      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
120      !
121      IF(lwp) THEN
122         WRITE(numout,*)
123         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
124         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
125         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
126         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
127         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
128         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
129         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
130         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
131      ENDIF
132      !
133      !
134      SELECT CASE( jphgr_msh )   !  type of horizontal mesh 
135      !
136      CASE ( 0 )                     !==  read in coordinate.nc file  ==!
137         !
138         IF(lwp) WRITE(numout,*)
139         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
140         !
141         ie1e2u_v = 0             ! set to unread e1e2u and e1e2v
142         iff = 0                  ! set to unread iff
143         !
144         CALL hgr_read( ie1e2u_v, iff )     ! read the coordinate.nc file !SF to be changed in mesh_mask
145         !
146         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them
147            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction
148            e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
149            e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
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#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         DO jj = 1, jpj
209            DO ji = 1, jpi
210               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
211               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
212               glamv(ji,jj) = glamt(ji,jj)
213               glamf(ji,jj) = glamu(ji,jj)
214               !
215               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
216               gphiu(ji,jj) = gphit(ji,jj)
217               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
218               gphif(ji,jj) = gphiv(ji,jj)
219            END DO
220         END DO
221         !
222         ! Horizontal scale factors (in meters)
223         !                              ======
224         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
225         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
226         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
227         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
228         !
229      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
230         !
231         IF(lwp) WRITE(numout,*)
232         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
233         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
234         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
235         !
236         !  Find index corresponding to the equator, given the grid spacing e1_deg
237         !  and the (approximate) southern latitude ppgphi0.
238         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
239         !  The formula should work even if the equator is outside the domain.
240         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
241         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
242         IF(  ppgphi0 > 0 )  ijeq = -ijeq
243         !
244         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
245         !
246         DO jj = 1, jpj
247            DO ji = 1, jpi
248               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
249               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
250               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
251               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
252         ! Longitude
253               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
254               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
255               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
256               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
257         ! Latitude
258               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
259               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
260               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
261               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
262         ! e1
263               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
264               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
265               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
266               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
267         ! e2
268               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
269               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
270               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
271               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
272            END DO
273         END DO
274         !
275      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
276         !
277         IF(lwp) WRITE(numout,*)
278         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
279         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
280         !
281         IF(lwp) THEN
282            WRITE(numout,*) '               control print of value of ln_read_cfg     = ', ln_read_cfg
283         ENDIF
284         IF (ln_read_cfg) THEN
285            IF(lwp) WRITE(numout,*) '          ln_read_cfg read input files'
286            ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v
287            iff = 0                       ! set to unread ff
288            !
289            CALL hgr_read( ie1e2u_v, iff )     ! read the coordinate.nc file !SF to be changed in mesh_mask
290            !
291            IF( ie1e2u_v == 0 ) THEN   ! e1e2u and e1e2v have not been read: compute them
292               !                          ! e2u and e1v does not include a
293               !                          reduction in some strait: apply reduction
294               e1e2u (:,:) = e1u(:,:) * e2u(:,:)
295               e1e2v (:,:) = e1v(:,:) * e2v(:,:)
296            ENDIF
297         ELSE
298            IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module'
299            CALL usr_def_hgr( nbench , jp_cfg  ,                    &
300            &                 ff,                         &
301            &                 glamt  , glamu   , glamv   , glamf   ,         &
302            &                 gphit  , gphiu   , gphiv   , gphif   ,         &
303            &                 e1t    , e1u     , e1v     , e1f     ,         &
304            &                 e2t    , e2u     , e2v     , e2f     ,         & 
305            &                 e1e2u  , e1e2v   , e2_e1u  , e1_e2v     )
306         !
307         ENDIF
308         !
309      CASE DEFAULT
310         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
311         CALL ctl_stop( ctmp1 )
312         !
313      END SELECT
314     
315      ! associated horizontal metrics
316      ! -----------------------------
317      !
318      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
319      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
320      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
321      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
322      !
323      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
324      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
325      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
326         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
327         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
328      ENDIF
329      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
330      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
331      !   
332      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
333      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
334
335      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
336         WRITE(numout,*)
337         WRITE(numout,*) '          longitude and e1 scale factors'
338         WRITE(numout,*) '          ------------------------------'
339         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
340            glamv(ji,1), glamf(ji,1),   &
341            e1t(ji,1), e1u(ji,1),   &
342            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3439300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
344            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
345            !
346         WRITE(numout,*)
347         WRITE(numout,*) '          latitude and e2 scale factors'
348         WRITE(numout,*) '          -----------------------------'
349         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
350            &                     gphiv(1,jj), gphif(1,jj),   &
351            &                     e2t  (1,jj), e2u  (1,jj),   &
352            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
353      ENDIF
354
355
356      ! ================= !
357      !  Coriolis factor  !
358      ! ================= !
359
360      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
361      !
362      CASE ( 0, 1, 4 )               ! mesh on the sphere
363         !
364         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
365         !
366      CASE ( 2 )                     ! f-plane at ppgphi0
367         !
368         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
369         !
370         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
371         !
372      CASE ( 3 )                     ! beta-plane
373         !
374         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
375         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
376         !
377#if defined key_agrif
378         IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN       ! for EEL6 configuration only
379            IF( .NOT.Agrif_Root() ) THEN
380              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
381            ENDIF
382         ENDIF
383#endif         
384         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
385         !
386         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
387         !
388         IF(lwp) THEN
389            WRITE(numout,*) 
390            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
391            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
392         ENDIF
393         IF( lk_mpp ) THEN
394            zminff=ff(nldi,nldj)
395            zmaxff=ff(nldi,nlej)
396            CALL mpp_min( zminff )   ! min over the global domain
397            CALL mpp_max( zmaxff )   ! max over the global domain
398            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
399         END IF
400         !
401      END SELECT
402
403
404      ! Control of domain for symetrical condition
405      ! ------------------------------------------
406      ! The equator line must be the latitude coordinate axe
407
408      IF( nperio == 2 ) THEN
409         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
410         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
411      ENDIF
412      !
413      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
414      !
415   END SUBROUTINE dom_hgr
416
417
418   SUBROUTINE hgr_read( ke1e2u_v, kff )
419      !!---------------------------------------------------------------------
420      !!              ***  ROUTINE hgr_read  ***
421      !!
422      !! ** Purpose :   Read a mesh_mask file in NetCDF format using IOM
423      !!
424      !!----------------------------------------------------------------------
425      USE iom
426      !!
427      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in mesh_mask file (=1) or not (=0)
428      INTEGER, INTENT( inout ) ::   kff        ! fag: kff read in mesh_mask file (=1) or not (=0)
429      !
430      INTEGER ::   inum   ! temporary logical unit
431      !!----------------------------------------------------------------------
432      !
433      IF(lwp) THEN
434         WRITE(numout,*)
435         WRITE(numout,*) 'hgr_read : read the horizontal coordinates in mesh_mask'
436         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
437      ENDIF
438      !
439      !SF CALL iom_open( 'coordinates', inum )
440      CALL iom_open( 'mesh_mask', inum )
441      !
442      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
443      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
444      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
445      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
446      !
447      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
448      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
449      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
450      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
451      !
452      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
453      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
454      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
455      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
456      !
457      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
458      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
459      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
460      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
461      !SF add read coriolis in mesh_mask file
462      CALL iom_get( inum, jpdom_data, 'ff'   , ff   , lrowattr=ln_use_jattr )
463      !
464      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
465         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in mesh_mask file'
466         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
467         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
468         ke1e2u_v = 1
469      ELSE
470         ke1e2u_v = 0
471      ENDIF
472      !
473      IF( iom_varid( inum, 'kff', ldstop = .FALSE. ) > 0 ) THEN
474         IF(lwp) WRITE(numout,*) 'hgr_read : ff read in mesh_mask file'
475         CALL iom_get( inum, jpdom_data, 'ff'  , ff  , lrowattr=ln_use_jattr )
476         kff = 1
477      ELSE
478         kff = 0
479      ENDIF
480      !
481      CALL iom_close( inum )
482     
483    END SUBROUTINE hgr_read
484   
485   !!======================================================================
486END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.