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 utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src – NEMO

source: utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src/domhgr.F90 @ 10727

Last change on this file since 10727 was 10727, checked in by rblod, 5 years ago

new nesting tools (attempt) and brutal cleaning of DOMAINcfg, see ticket #2129
File size: 25.4 KB
Line 
1MODULE domhgr
2   !!==============================================================================
3   !!                       ***  MODULE domhgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1988-03  (G. Madec) Original code
7   !!            7.0  ! 1996-01  (G. Madec)  terrain following coordinates
8   !!            8.0  ! 1997-02  (G. Madec)  print mesh informations
9   !!            8.1  ! 1999-11  (M. Imbard) NetCDF format with IO-IPSL
10   !!            8.2  ! 2000-08  (D. Ludicone) Reduced section at Bab el Mandeb
11   !!             -   ! 2001-09  (M. Levy)  eel config: grid in km, beta-plane
12   !!  NEMO      1.0  ! 2002-08  (G. Madec)  F90: Free form and module, namelist
13   !!             -   ! 2004-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
14   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist
15   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration
16   !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse
17   !!                                       add optional read of e1e2u & e1e2v
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dom_hgr       : initialize the horizontal mesh
22   !!   hgr_read      : read "coordinate" NetCDF file
23   !!----------------------------------------------------------------------
24   USE dom_oce        ! ocean space and time domain
25   USE phycst         ! physical constants
26   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
27   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
36
37   PUBLIC   dom_hgr   ! called by domain.F90
38
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
41   !! $Id: domhgr.F90 6140 2015-12-21 11:35:23Z timgraham $
42   !! Software governed by the CeCILL licence     (./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE dom_hgr
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dom_hgr  ***
49      !!
50      !! ** Purpose :   Compute the geographical position (in degre) of the
51      !!      model grid-points,  the horizontal scale factors (in meters) and
52      !!      the Coriolis factor (in s-1).
53      !!
54      !! ** Method  :   The geographical position of the model grid-points is
55      !!      defined from analytical functions, fslam and fsphi, the deriva-
56      !!      tives of which gives the horizontal scale factors e1,e2.
57      !!      Defining two function fslam and fsphi and their derivatives in
58      !!      the two horizontal directions (fse1 and fse2), the model grid-
59      !!      point position and scale factors are given by:
60      !!         t-point:
61      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
62      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
63      !!         u-point:
64      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
65      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
66      !!         v-point:
67      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
68      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
69      !!            f-point:
70      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
71      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
72      !!      Where fse1 and fse2 are defined by:
73      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
74      !!                                     +          di(fsphi) **2 )(i,j)
75      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
76      !!                                     +          dj(fsphi) **2 )(i,j)
77      !!
78      !!        The coriolis factor is given at z-point by:
79      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
80      !!
81      !!        This routine is given as an example, it must be modified
82      !!      following the user s desiderata. nevertheless, the output as
83      !!      well as the way to compute the model grid-point position and
84      !!      horizontal scale factors must be respected in order to insure
85      !!      second order accuracy schemes.
86      !!
87      !! N.B. If the domain is periodic, verify that scale factors are also
88      !!      periodic, and the coriolis term again.
89      !!
90      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
91      !!                u-, v- and f-points (in degre)
92      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
93      !!               u-, v-  and f-points (in degre)
94      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
95      !!      scale factors (in meters) at t-, u-, v-, and f-points.
96      !!        define ff: coriolis factor at f-point
97      !!
98      !! References :   Marti, Madec and Delecluse, 1992, JGR
99      !!                Madec, Imbard, 1996, Clim. Dyn.
100      !!----------------------------------------------------------------------
101      INTEGER  ::   ji, jj               ! dummy loop indices
102      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
103      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
104      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
105      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
106      REAL(wp) ::   zphi0, zbeta, znorme !
107      REAL(wp) ::   zarg, zf0, zminff, zmaxff
108      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
109      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
110      INTEGER  ::   isrow                ! index for ORCA1 starting row
111      INTEGER  ::   ie1e2u_v             ! fag for u- & v-surface read in coordinate file or not
112      !!----------------------------------------------------------------------
113      !
114  !    IF( 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 defined key_agrif
134         IF (agrif_root()) THEN
135#endif
136         IF(lwp) WRITE(numout,*)
137         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
138         !
139         ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v
140         !
141         CALL hgr_read( ie1e2u_v )     ! read the coordinate.nc file
142         !
143         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them
144            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction
145            e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
146            e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
147         ENDIF
148#if defined key_agrif
149         ELSE
150           CALL Agrif_InitValues_cont()
151         ENDIF
152#endif
153         !
154      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==!
155         !
156         IF(lwp) WRITE(numout,*)
157         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
158         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
159         !
160         DO jj = 1, jpj
161            DO ji = 1, jpi
162               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 )
163               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 )
164               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5
165               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5
166         ! Longitude
167               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
168               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
169               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
170               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
171         ! Latitude
172               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
173               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
174               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
175               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
176         ! e1
177               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
178               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
179               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
180               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
181         ! e2
182               e2t(ji,jj) = ra * rad * ppe2_deg
183               e2u(ji,jj) = ra * rad * ppe2_deg
184               e2v(ji,jj) = ra * rad * ppe2_deg
185               e2f(ji,jj) = ra * rad * ppe2_deg
186            END DO
187         END DO
188         !
189      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==!
190         !
191         IF(lwp) WRITE(numout,*)
192         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
193         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
194         !
195         ! Position coordinates (in kilometers)
196         !                          ==========
197         glam0 = 0._wp
198         gphi0 = - ppe2_m * 1.e-3
199         !
200         DO jj = 1, jpj
201            DO ji = 1, jpi
202               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
203               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
204               glamv(ji,jj) = glamt(ji,jj)
205               glamf(ji,jj) = glamu(ji,jj)
206               !
207               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
208               gphiu(ji,jj) = gphit(ji,jj)
209               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
210               gphif(ji,jj) = gphiv(ji,jj)
211            END DO
212         END DO
213         !
214         ! Horizontal scale factors (in meters)
215         !                              ======
216         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
217         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
218         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
219         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
220         !
221      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
222         !
223         IF(lwp) WRITE(numout,*)
224         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
225         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
226         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
227         !
228         !  Find index corresponding to the equator, given the grid spacing e1_deg
229         !  and the (approximate) southern latitude ppgphi0.
230         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
231         !  The formula should work even if the equator is outside the domain.
232         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
233         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
234         IF(  ppgphi0 > 0 )  ijeq = -ijeq
235         !
236         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
237         !
238         DO jj = 1, jpj
239            DO ji = 1, jpi
240               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
241               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
242               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
243               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
244         ! Longitude
245               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
246               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
247               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
248               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
249         ! Latitude
250               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
251               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
252               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
253               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
254         ! e1
255               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
256               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
257               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
258               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
259         ! e2
260               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
261               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
262               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
263               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
264            END DO
265         END DO
266         !
267      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
268         !
269         IF(lwp) WRITE(numout,*)
270         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
271         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
272         !
273         ! Position coordinates (in kilometers)
274         !                          ==========
275         !
276         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
277         zlam1 = -85._wp
278         zphi1 =  29._wp
279         ! resolution in meters
280         ze1 = 106000. / REAL( jp_cfg , wp )           
281         ! benchmark: forced the resolution to be about 100 km
282       !  IF( nbench /= 0 )   ze1 = 106000._wp     
283         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
284         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
285         ze1deg = ze1 / (ra * rad)
286       !  IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
287         !                                                           ! at the right jp_cfg resolution
288         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
289         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
290         !
291         IF( nprint==1 .AND. lwp )   THEN
292            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
293            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
294         ENDIF
295         !
296         DO jj = 1, jpj
297            DO ji = 1, jpi
298               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
299               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
300               !
301               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
302               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
303               !
304               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
305               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
306               !
307               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
308               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
309               !
310               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
311               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
312            END DO
313         END DO
314         !
315         ! Horizontal scale factors (in meters)
316         !                              ======
317         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
318         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
319         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
320         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
321         !
322      CASE DEFAULT
323         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
324         CALL ctl_stop( ctmp1 )
325         !
326      END SELECT
327     
328      ! associated horizontal metrics
329      ! -----------------------------
330      !
331      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
332      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
333      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
334      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
335      !
336      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
337      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
338      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
339         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
340         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
341      ENDIF
342      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
343      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
344      !   
345      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
346      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
347
348      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
349         WRITE(numout,*)
350         WRITE(numout,*) '          longitude and e1 scale factors'
351         WRITE(numout,*) '          ------------------------------'
352         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
353            glamv(ji,1), glamf(ji,1),   &
354            e1t(ji,1), e1u(ji,1),   &
355            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3569300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
357            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
358            !
359         WRITE(numout,*)
360         WRITE(numout,*) '          latitude and e2 scale factors'
361         WRITE(numout,*) '          -----------------------------'
362         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
363            &                     gphiv(1,jj), gphif(1,jj),   &
364            &                     e2t  (1,jj), e2u  (1,jj),   &
365            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
366      ENDIF
367
368
369      ! ================= !
370      !  Coriolis factor  !
371      ! ================= !
372
373      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
374      !
375      CASE ( 0, 1, 4 )               ! mesh on the sphere
376         !
377         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
378         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       - at t-point
379         !
380      CASE ( 2 )                     ! f-plane at ppgphi0
381         !
382         ff_f(:,:) = 2. * omega * SIN( rad * ppgphi0 )
383         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 )
384         !
385         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff_f(1,1)
386         !
387      CASE ( 3 )                     ! beta-plane
388         !
389         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
390         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
391         !
392         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
393         !
394         ff_f(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
395         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
396         !
397         IF(lwp) THEN
398            WRITE(numout,*) 
399            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff_f(nldi,nldj)
400            WRITE(numout,*) '          Coriolis parameter varies from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej)
401         ENDIF
402         IF( lk_mpp ) THEN
403            zminff=ff_f(nldi,nldj)
404            zmaxff=ff_f(nldi,nlej)
405            CALL mpp_min( 'toto',zminff )   ! min over the global domain
406            CALL mpp_max( 'toto',zmaxff )   ! max over the global domain
407            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
408         END IF
409         !
410      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
411         !
412         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
413         zphi0 = 15._wp                                                     ! latitude of the first row F-points
414         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
415         !
416         ff_f(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
417         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
418         !
419         IF(lwp) THEN
420            WRITE(numout,*) 
421            WRITE(numout,*) '          Beta-plane and rotated domain : '
422            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej)
423         ENDIF
424         !
425         IF( lk_mpp ) THEN
426            zminff=ff_f(nldi,nldj)
427            zmaxff=ff_f(nldi,nlej)
428            CALL mpp_min('toto', zminff )   ! min over the global domain
429            CALL mpp_max( 'toto',zmaxff )   ! max over the global domain
430            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
431         END IF
432         !
433      END SELECT
434
435
436      ! Control of domain for symetrical condition
437      ! ------------------------------------------
438      ! The equator line must be the latitude coordinate axe
439
440!      IF( nperio == 2 ) THEN
441!         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
442!         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
443!      ENDIF
444      !
445    !  IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
446      !
447   END SUBROUTINE dom_hgr
448
449
450   SUBROUTINE hgr_read( ke1e2u_v )
451      !!---------------------------------------------------------------------
452      !!              ***  ROUTINE hgr_read  ***
453      !!
454      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
455      !!
456      !!----------------------------------------------------------------------
457      USE iom
458      !!
459      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0)
460      !
461      INTEGER ::   inum   ! temporary logical unit
462      !!----------------------------------------------------------------------
463      !
464      IF(lwp) THEN
465         WRITE(numout,*)
466         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
467         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
468      ENDIF
469      !
470      CALL iom_open( 'coordinates', inum )
471      !
472      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
473      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
474      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
475      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
476      !
477      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
478      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
479      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
480      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
481      !
482      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
483      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
484      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
485      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
486      !
487      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
488      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
489      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
490      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
491      !
492      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
493         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
494         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
495         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
496         ke1e2u_v = 1
497      ELSE
498         ke1e2u_v = 0
499      ENDIF
500      !
501      CALL iom_close( inum )
502     
503    END SUBROUTINE hgr_read
504   
505   !!======================================================================
506END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.