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

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 6404

Last change on this file since 6404 was 6404, checked in by timgraham, 8 years ago

First attempt at upgrading branch to the head of the trunk. This should include all of the simplification branch from the merge in Dec 2015.

  • Property svn:keywords set to Id
File size: 25.5 KB
RevLine 
[3]1MODULE domhgr
2   !!==============================================================================
[93]3   !!                       ***  MODULE domhgr   ***
[3]4   !! Ocean initialization : domain initialization
5   !!==============================================================================
[2528]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
[6404]16   !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse
17   !!                                       add optional read of e1e2u & e1e2v
[473]18   !!----------------------------------------------------------------------
[3]19
20   !!----------------------------------------------------------------------
[2528]21   !!   dom_hgr       : initialize the horizontal mesh
22   !!   hgr_read      : read "coordinate" NetCDF file
[3]23   !!----------------------------------------------------------------------
[2528]24   USE dom_oce        ! ocean space and time domain
25   USE phycst         ! physical constants
[6404]26   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
27   !
[2528]28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
[3294]30   USE timing         ! Timing
[3]31
32   IMPLICIT NONE
33   PRIVATE
34
[2528]35   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
[3]36
[2528]37   PUBLIC   dom_hgr   ! called by domain.F90
38
[3]39   !!----------------------------------------------------------------------
[6404]40   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
[5819]41   !! $Id$
[2528]42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]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:
[81]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)
[3]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      !!
[473]98      !! References :   Marti, Madec and Delecluse, 1992, JGR
99      !!                Madec, Imbard, 1996, Clim. Dyn.
[3]100      !!----------------------------------------------------------------------
[2528]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
[5385]110      INTEGER  ::   isrow                ! index for ORCA1 starting row
[6404]111      INTEGER  ::   ie1e2u_v             ! fag for u- & v-surface read in coordinate file or not
[3]112      !!----------------------------------------------------------------------
[3294]113      !
114      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
115      !
[3]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
[6404]127      !
128      !
129      SELECT CASE( jphgr_msh )   !  type of horizontal mesh 
130      !
131      CASE ( 0 )                     !==  read in coordinate.nc file  ==!
132         !
[3]133         IF(lwp) WRITE(numout,*)
[81]134         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
[6404]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            e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
143            e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
[81]144         ENDIF
[3]145         !
[6404]146      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==!
147         !
[3]148         IF(lwp) WRITE(numout,*)
149         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
150         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
[6404]151         !
[3]152         DO jj = 1, jpj
153            DO ji = 1, jpi
[6404]154               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 )
155               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 )
156               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5
157               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5
[3]158         ! Longitude
159               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
160               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
161               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
162               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
163         ! Latitude
164               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
165               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
166               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
167               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
168         ! e1
169               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
170               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
171               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
172               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
173         ! e2
174               e2t(ji,jj) = ra * rad * ppe2_deg
175               e2u(ji,jj) = ra * rad * ppe2_deg
176               e2v(ji,jj) = ra * rad * ppe2_deg
177               e2f(ji,jj) = ra * rad * ppe2_deg
178            END DO
179         END DO
[6404]180         !
181      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==!
182         !
[3]183         IF(lwp) WRITE(numout,*)
184         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
185         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
[6404]186         !
[3]187         ! Position coordinates (in kilometers)
188         !                          ==========
[6404]189         glam0 = 0._wp
[3]190         gphi0 = - ppe2_m * 1.e-3
[6404]191         !
[4147]192#if defined key_agrif 
193         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
194            IF( .NOT. Agrif_Root() ) THEN
195              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
196              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
197              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
198              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
199            ENDIF
[389]200         ENDIF
201#endif         
[3]202         DO jj = 1, jpj
203            DO ji = 1, jpi
[6404]204               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
205               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
[3]206               glamv(ji,jj) = glamt(ji,jj)
207               glamf(ji,jj) = glamu(ji,jj)
[6404]208               !
209               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
[3]210               gphiu(ji,jj) = gphit(ji,jj)
[6404]211               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
[3]212               gphif(ji,jj) = gphiv(ji,jj)
213            END DO
214         END DO
[6404]215         !
[3]216         ! Horizontal scale factors (in meters)
217         !                              ======
218         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
219         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
220         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
221         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
[6404]222         !
223      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
224         !
[3]225         IF(lwp) WRITE(numout,*)
226         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
227         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
[473]228         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
[6404]229         !
[3]230         !  Find index corresponding to the equator, given the grid spacing e1_deg
231         !  and the (approximate) southern latitude ppgphi0.
232         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
233         !  The formula should work even if the equator is outside the domain.
234         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
[29]235         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
[224]236         IF(  ppgphi0 > 0 )  ijeq = -ijeq
[6404]237         !
[29]238         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
[6404]239         !
[3]240         DO jj = 1, jpj
241            DO ji = 1, jpi
[6404]242               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
243               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
244               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
245               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
[3]246         ! Longitude
247               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
248               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
249               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
250               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
251         ! Latitude
252               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
[93]253               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
254               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
255               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
[3]256         ! e1
257               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
258               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
259               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
260               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
261         ! e2
262               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
263               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
264               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
265               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
266            END DO
267         END DO
[6404]268         !
269      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
270         !
[93]271         IF(lwp) WRITE(numout,*)
272         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
273         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
[6404]274         !
[93]275         ! Position coordinates (in kilometers)
276         !                          ==========
[6404]277         !
[167]278         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
[6404]279         zlam1 = -85._wp
280         zphi1 =  29._wp
[167]281         ! resolution in meters
[6404]282         ze1 = 106000. / REAL( jp_cfg , wp )           
[167]283         ! benchmark: forced the resolution to be about 100 km
[6404]284         IF( nbench /= 0 )   ze1 = 106000._wp     
285         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
286         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
[93]287         ze1deg = ze1 / (ra * rad)
[6404]288         IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
289         !                                                           ! at the right jp_cfg resolution
290         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
291         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
292         !
[516]293         IF( nprint==1 .AND. lwp )   THEN
294            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
295            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
296         ENDIF
[6404]297         !
[93]298         DO jj = 1, jpj
[6404]299            DO ji = 1, jpi
300               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
301               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
302               !
303               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
304               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
305               !
306               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
307               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
308               !
309               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
310               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
311               !
312               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
313               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
314            END DO
315         END DO
316         !
[93]317         ! Horizontal scale factors (in meters)
318         !                              ======
319         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
320         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
321         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
322         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
[6404]323         !
[3]324      CASE DEFAULT
[473]325         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
326         CALL ctl_stop( ctmp1 )
[6404]327         !
[3]328      END SELECT
[2715]329     
[6404]330      ! associated horizontal metrics
331      ! -----------------------------
332      !
333      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
334      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
335      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
336      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
337      !
338      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
339      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
340      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
341         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
342         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
343      ENDIF
344      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
345      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
346      !   
347      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
348      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
[3]349
[6404]350      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
[3]351         WRITE(numout,*)
352         WRITE(numout,*) '          longitude and e1 scale factors'
353         WRITE(numout,*) '          ------------------------------'
354         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
355            glamv(ji,1), glamf(ji,1),   &
356            e1t(ji,1), e1u(ji,1),   &
357            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3589300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
359            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
[6404]360            !
[3]361         WRITE(numout,*)
362         WRITE(numout,*) '          latitude and e2 scale factors'
363         WRITE(numout,*) '          -----------------------------'
364         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
365            &                     gphiv(1,jj), gphif(1,jj),   &
366            &                     e2t  (1,jj), e2u  (1,jj),   &
367            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
368      ENDIF
369
370
371      ! ================= !
372      !  Coriolis factor  !
373      ! ================= !
374
375      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
[6404]376      !
[3]377      CASE ( 0, 1, 4 )               ! mesh on the sphere
[6404]378         !
[3]379         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
[6404]380         !
[3]381      CASE ( 2 )                     ! f-plane at ppgphi0
[6404]382         !
[3]383         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
[6404]384         !
[3]385         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
[6404]386         !
[3]387      CASE ( 3 )                     ! beta-plane
[6404]388         !
[187]389         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
[6404]390         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
391         !
[4147]392#if defined key_agrif
[6404]393         IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN       ! for EEL6 configuration only
394            IF( .NOT.Agrif_Root() ) THEN
395              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
[4147]396            ENDIF
[389]397         ENDIF
398#endif         
[187]399         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
[6404]400         !
[187]401         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
[6404]402         !
[516]403         IF(lwp) THEN
404            WRITE(numout,*) 
405            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
406            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
407         ENDIF
[434]408         IF( lk_mpp ) THEN
409            zminff=ff(nldi,nldj)
410            zmaxff=ff(nldi,nlej)
411            CALL mpp_min( zminff )   ! min over the global domain
412            CALL mpp_max( zmaxff )   ! max over the global domain
[516]413            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
[434]414         END IF
[6404]415         !
[434]416      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
[6404]417         !
[93]418         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
[6404]419         zphi0 = 15._wp                                                     ! latitude of the first row F-points
[93]420         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
[6404]421         !
[93]422         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
[6404]423         !
[516]424         IF(lwp) THEN
425            WRITE(numout,*) 
426            WRITE(numout,*) '          Beta-plane and rotated domain : '
427            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
428         ENDIF
[6404]429         !
[434]430         IF( lk_mpp ) THEN
431            zminff=ff(nldi,nldj)
432            zmaxff=ff(nldi,nlej)
433            CALL mpp_min( zminff )   ! min over the global domain
434            CALL mpp_max( zmaxff )   ! max over the global domain
[516]435            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
[434]436         END IF
[6404]437         !
[3]438      END SELECT
439
440
441      ! Control of domain for symetrical condition
442      ! ------------------------------------------
443      ! The equator line must be the latitude coordinate axe
444
445      IF( nperio == 2 ) THEN
[6404]446         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
[473]447         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
[3]448      ENDIF
[3294]449      !
450      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
451      !
[3]452   END SUBROUTINE dom_hgr
453
454
[6404]455   SUBROUTINE hgr_read( ke1e2u_v )
[3]456      !!---------------------------------------------------------------------
[81]457      !!              ***  ROUTINE hgr_read  ***
[3]458      !!
[6404]459      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
[3]460      !!
461      !!----------------------------------------------------------------------
[473]462      USE iom
[6404]463      !!
464      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0)
465      !
[473]466      INTEGER ::   inum   ! temporary logical unit
[3]467      !!----------------------------------------------------------------------
[6404]468      !
[3]469      IF(lwp) THEN
470         WRITE(numout,*)
[81]471         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
[473]472         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
[3]473      ENDIF
[6404]474      !
[473]475      CALL iom_open( 'coordinates', inum )
[6404]476      !
[5118]477      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
478      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
479      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
480      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
[6404]481      !
[5118]482      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
483      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
484      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
485      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
[6404]486      !
487      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
488      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
489      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
490      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
491      !
492      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
493      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
494      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
495      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
496      !
497      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
498         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
499         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
500         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
501         ke1e2u_v = 1
502      ELSE
503         ke1e2u_v = 0
504      ENDIF
505      !
[473]506      CALL iom_close( inum )
507     
508    END SUBROUTINE hgr_read
509   
[3]510   !!======================================================================
511END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.