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 tags/nemo_dev_x2/NEMO/OPA_SRC/DOM – NEMO

source: tags/nemo_dev_x2/NEMO/OPA_SRC/DOM/domhgr.F90 @ 4310

Last change on this file since 4310 was 81, checked in by opalod, 20 years ago

CT : UPDATE055 : Change the subroutine name dom_hgr_coo to hgr_read, and include the old domhgr_coo_fdir.h90 file in domhgr.F90

as a subroutine called hgr_read_fdir; remove the domhgr_coo_fdir.h90
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.3 KB
Line 
1MODULE domhgr
2   !!==============================================================================
3   !!                       ***  MODULE domain   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_hgr        : initialize the horizontal mesh
9   !!   hgr_read       : read "coordinate" NetCDF file
10   !!   hgr_read_fdir  : read "coordinate" direct access file
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE dom_oce         ! ocean space and time domain
14   USE phycst          ! physical constants
15   USE in_out_manager  ! I/O manager
16
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Module variables
21   REAL (wp)  :: glam0, gphi0        ! variables corresponding to parameters
22      !                              ! ppglam0 ppgphi0 set in par_oce
23
24   !! * Routine accessibility
25   PUBLIC dom_hgr        ! called by domain.F90
26   !!----------------------------------------------------------------------
27   !!   OPA 9.0 , LODYC-IPSL   (2003)
28   !!----------------------------------------------------------------------
29
30CONTAINS
31
32   SUBROUTINE dom_hgr
33      !!----------------------------------------------------------------------
34      !!                  ***  ROUTINE dom_hgr  ***
35      !!
36      !! ** Purpose :   Compute the geographical position (in degre) of the
37      !!      model grid-points,  the horizontal scale factors (in meters) and
38      !!      the Coriolis factor (in s-1).
39      !!
40      !! ** Method  :   The geographical position of the model grid-points is
41      !!      defined from analytical functions, fslam and fsphi, the deriva-
42      !!      tives of which gives the horizontal scale factors e1,e2.
43      !!      Defining two function fslam and fsphi and their derivatives in
44      !!      the two horizontal directions (fse1 and fse2), the model grid-
45      !!      point position and scale factors are given by:
46      !!         t-point:
47      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
48      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
49      !!         u-point:
50      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
51      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
52      !!         v-point:
53      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
54      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
55      !!            f-point:
56      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
57      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
58      !!      Where fse1 and fse2 are defined by:
59      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
60      !!                                     +          di(fsphi) **2 )(i,j)
61      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
62      !!                                     +          dj(fsphi) **2 )(i,j)
63      !!
64      !!        The coriolis factor is given at z-point by:
65      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
66      !!
67      !!        This routine is given as an example, it must be modified
68      !!      following the user s desiderata. nevertheless, the output as
69      !!      well as the way to compute the model grid-point position and
70      !!      horizontal scale factors must be respected in order to insure
71      !!      second order accuracy schemes.
72      !!
73      !! N.B. If the domain is periodic, verify that scale factors are also
74      !!      periodic, and the coriolis term again.
75      !!
76      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
77      !!                u-, v- and f-points (in degre)
78      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
79      !!               u-, v-  and f-points (in degre)
80      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
81      !!      scale factors (in meters) at t-, u-, v-, and f-points.
82      !!        define ff: coriolis factor at f-point
83      !!
84      !! References :
85      !!      Marti, Madec and Delecluse, 1992, j. geophys. res., in press.
86      !!
87      !! History :
88      !!        !  88-03  (G. Madec)
89      !!        !  91-11  (G. Madec)
90      !!        !  92-06  (M. Imbard)
91      !!        !  96-01  (G. Madec)  terrain following coordinates
92      !!        !  97-02  (G. Madec)  print mesh informations
93      !!        !  01-09  (M. Levy)  eel config: grid in km, beta-plane
94      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, namelist
95      !!   9.0  !  04-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
96      !!                  use of parameters in par_CONFIG-Rxx.h90, not in namelist
97      !!----------------------------------------------------------------------
98      !! * local declarations
99      INTEGER  ::   ji, jj              ! dummy loop indices
100      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
101      INTEGER  ::   ijeq                ! index of equator T point (used in case 4)
102      REAL(wp) ::   &
103         zti, zui, zvi, zfi,         &  ! temporary scalars
104         ztj, zuj, zvj, zfj,         &  !
105         zphi0, zbeta, znorme,       &  !
106         zarg, zf0
107      !!----------------------------------------------------------------------
108
109      IF(lwp) THEN
110         WRITE(numout,*)
111         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
112         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
113         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
114         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
115         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
116         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
117         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
118         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
119      ENDIF
120
121
122      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
123
124      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
125
126         IF(lwp) WRITE(numout,*)
127         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
128#if defined key_fdir
129         CALL hgr_read_fdir      ! 'key_fdir'       :   direct access file
130#else
131         CALL hgr_read           ! Defaultl option  :   NetCDF file
132#endif
133
134         !                                                ! =====================
135         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
136            !                                             ! =====================
137            IF( n_cla == 0 ) THEN
138               ii0 = 160   ;   ii1 = 161        ! Bab el Mandeb (e2u = 18 km)
139               ij0 =  88   ;   ij1 =  88   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
140               IF(lwp) WRITE(numout,*)
141               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 18 km'
142            ENDIF
143
144            ii0 = 145   ;   ii1 = 146        ! Sound Strait (e2u = 15 km)
145            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
146            IF(lwp) WRITE(numout,*)
147            IF(lwp) WRITE(numout,*) '             orca_r2: Reduced e2u at the Sound Strait'
148            !
149         ENDIF
150
151         !                                                ! ======================
152         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
153            !                                             ! ======================
154            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
155            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
156            IF(lwp) WRITE(numout,*)
157            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
158            !
159         ENDIF
160
161
162         ! N.B. :  General case, lat and long function of both i and j indices:
163         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
164         !                                  + (                           fsdiph( zti, ztj ) )**2  )
165         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
166         !                                  + (                           fsdiph( zui, zuj ) )**2  )
167         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
168         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
169         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
170         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
171         !
172         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
173         !                                  + (                           fsdjph( zti, ztj ) )**2  )
174         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
175         !                                  + (                           fsdjph( zui, zuj ) )**2  )
176         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
177         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
178         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
179         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
180
181
182      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
183
184         IF(lwp) WRITE(numout,*)
185         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
186         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
187
188         DO jj = 1, jpj
189            DO ji = 1, jpi
190               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
191               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
192               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
193               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
194         ! Longitude
195               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
196               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
197               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
198               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
199         ! Latitude
200               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
201               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
202               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
203               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
204         ! e1
205               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
206               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
207               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
208               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
209         ! e2
210               e2t(ji,jj) = ra * rad * ppe2_deg
211               e2u(ji,jj) = ra * rad * ppe2_deg
212               e2v(ji,jj) = ra * rad * ppe2_deg
213               e2f(ji,jj) = ra * rad * ppe2_deg
214            END DO
215         END DO
216
217
218      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
219
220         IF(lwp) WRITE(numout,*)
221         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
222         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
223
224         ! Position coordinates (in kilometers)
225         !                          ==========
226         glam0 = 0.e0
227         gphi0 = - ppe2_m * 1.e-3
228         DO jj = 1, jpj
229            DO ji = 1, jpi
230               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
231               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
232               glamv(ji,jj) = glamt(ji,jj)
233               glamf(ji,jj) = glamu(ji,jj)
234   
235               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
236               gphiu(ji,jj) = gphit(ji,jj)
237               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
238               gphif(ji,jj) = gphiv(ji,jj)
239            END DO
240         END DO
241
242         ! Horizontal scale factors (in meters)
243         !                              ======
244         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
245         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
246         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
247         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
248
249      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
250
251         IF(lwp) WRITE(numout,*)
252         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
253         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
254         IF ( ppgphi0 == -90 ) THEN
255                IF(lwp) WRITE(numout,*) ' Mercator grid cannot start at south pole !!!! '
256                IF(lwp) WRITE(numout,*) ' We stop '
257                STOP
258         ENDIF
259
260         !  Find index corresponding to the equator, given the grid spacing e1_deg
261         !  and the (approximate) southern latitude ppgphi0.
262         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
263         !  The formula should work even if the equator is outside the domain.
264         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
265         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
266
267         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
268
269         DO jj = 1, jpj
270            DO ji = 1, jpi
271               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
272               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
273               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
274               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
275         ! Longitude
276               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
277               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
278               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
279               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
280         ! Latitude
281               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
282               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
283               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
284               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
285         ! e1
286               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
287               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
288               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
289               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
290         ! e2
291               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
292               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
293               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
294               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
295            END DO
296         END DO
297
298      CASE DEFAULT
299         IF(lwp) WRITE(numout,cform_err)
300         IF(lwp) WRITE(numout,*) '          bad flag value for jphgr_msh = ', jphgr_msh
301         nstop = nstop + 1
302
303      END SELECT
304
305
306      ! Control printing : Grid informations (if not restart)
307      ! ----------------
308
309      IF(lwp .AND. .NOT.ln_rstart ) THEN
310         WRITE(numout,*)
311         WRITE(numout,*) '          longitude and e1 scale factors'
312         WRITE(numout,*) '          ------------------------------'
313         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
314            glamv(ji,1), glamf(ji,1),   &
315            e1t(ji,1), e1u(ji,1),   &
316            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
3179300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
318            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
319         
320         WRITE(numout,*)
321         WRITE(numout,*) '          latitude and e2 scale factors'
322         WRITE(numout,*) '          -----------------------------'
323         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
324            &                     gphiv(1,jj), gphif(1,jj),   &
325            &                     e2t  (1,jj), e2u  (1,jj),   &
326            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
327      ENDIF
328
329     
330      IF( nprint == 1 .AND. lwp ) THEN
331         WRITE(numout,*) '          e1u e2u '
332         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
333         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
334         WRITE(numout,*) '          e1v e2v  '
335         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
336         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
337         WRITE(numout,*) '          e1f e2f  '
338         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
339         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
340      ENDIF
341
342
343      ! ================= !
344      !  Coriolis factor  !
345      ! ================= !
346
347      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
348
349      CASE ( 0, 1, 4 )               ! mesh on the sphere
350
351         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
352
353      CASE ( 2 )                     ! f-plane at ppgphi0
354
355         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
356
357         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
358
359      CASE ( 3 )                     ! beta-plane
360
361         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
362         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m *1.e-3  / ( ra * rad ) ! latitude of the first row F-points
363         zf0     = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
364
365         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                      ! f = f0 +beta* y ( y=0 at south)
366
367         IF(lwp) WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(1,1)
368         IF(lwp) WRITE(numout,*) '                      Coriolis parameter varies from ', ff(1,1),' to ', ff(1,jpj)
369
370      END SELECT
371
372
373      ! Control of domain for symetrical condition
374      ! ------------------------------------------
375      ! The equator line must be the latitude coordinate axe
376
377      IF( nperio == 2 ) THEN
378         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
379         IF( znorme > 1.e-13 ) THEN
380            IF(lwp) WRITE(numout,cform_err)
381            IF(lwp) WRITE(numout,*) ' ===>>>> : symmetrical condition: rerun with good equator line'
382            nstop = nstop + 1
383         ENDIF
384      ENDIF
385
386   END SUBROUTINE dom_hgr
387
388
389   SUBROUTINE hgr_read
390      !!---------------------------------------------------------------------
391      !!              ***  ROUTINE hgr_read  ***
392      !!
393      !! ** Purpose :   Read a coordinate file in NetCDF format
394      !!
395      !! ** Method  :   The mesh file has been defined trough a analytical
396      !!      or semi-analytical method. It is read in a NetCDF file.
397      !!     
398      !! References :
399      !!      Marti, Madec and Delecluse, 1992, JGR, 97, 12,763-12,766.
400      !!      Madec, Imbard, 1996, Clim. Dyn., 12, 381-388.
401      !!
402      !! History :
403      !!        !         (O. Marti)  Original code
404      !!        !  91-03  (G. Madec)
405      !!        !  92-07  (M. Imbard)
406      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
407      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
408      !!   8.5  !  02-06  (G. Madec)  F90: Free form
409      !!----------------------------------------------------------------------
410      !! * Modules used
411      USE ioipsl
412
413      !! * Local declarations
414      LOGICAL ::   llog = .FALSE.
415      CHARACTER(len=21) ::   clname = 'coordinates'
416      INTEGER  ::   ji, jj              ! dummy loop indices
417      INTEGER  ::   inum                ! temporary logical unit
418      INTEGER  ::   ilev, itime         ! temporary integers
419      REAL(wp) ::   zdt, zdate0         ! temporary scalars
420      REAL(wp) ::   zdept(1)            ! temporary workspace
421      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
422         zlamt, zphit, zdta             ! temporary workspace (NetCDF read)
423      !!----------------------------------------------------------------------
424
425
426      ! 1. Read of the grid coordinates and scale factors
427      ! -------------------------------------------------
428
429      IF(lwp) THEN
430         WRITE(numout,*)
431         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
432         WRITE(numout,*) '~~~~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
433      ENDIF
434
435      ! read the file
436      itime = 0
437      ilev = 1   
438      zlamt(:,:) = 0.e0
439      zphit(:,:) = 0.e0
440      CALL restini( clname, jpidta, jpjdta, zlamt , zphit,   &
441         &                  ilev  , zdept , clname,          &
442         &                  itime , zdate0, zdt   , inum )
443
444      CALL restget( inum, 'glamt', jpidta, jpjdta, 1, 0, llog, zdta )
445      DO jj = 1, nlcj
446         DO ji = 1, nlci
447            glamt(ji,jj) = zdta(mig(ji),mjg(jj))
448         END DO
449      END DO
450      CALL restget( inum, 'glamu', jpidta, jpjdta, 1, 0, llog, zdta )
451      DO jj = 1, nlcj
452         DO ji = 1, nlci
453            glamu(ji,jj) = zdta(mig(ji),mjg(jj))                   
454         END DO
455      END DO
456      CALL restget( inum, 'glamv', jpidta, jpjdta, 1, 0, llog, zdta )
457      DO jj = 1, nlcj
458         DO ji = 1, nlci
459            glamv(ji,jj) = zdta(mig(ji),mjg(jj))                   
460         END DO
461      END DO
462      CALL restget( inum, 'glamf', jpidta, jpjdta, 1, 0, llog, zdta )
463      DO jj = 1, nlcj
464         DO ji = 1, nlci
465            glamf(ji,jj) = zdta(mig(ji),mjg(jj))                   
466         END DO
467      END DO
468      CALL restget( inum, 'gphit', jpidta, jpjdta, 1, 0, llog, zdta )
469      DO jj = 1, nlcj
470         DO ji = 1, nlci
471            gphit(ji,jj) = zdta(mig(ji),mjg(jj))                   
472         END DO
473      END DO
474      CALL restget( inum, 'gphiu', jpidta, jpjdta, 1, 0, llog, zdta )
475      DO jj = 1, nlcj
476         DO ji = 1, nlci
477            gphiu(ji,jj) = zdta(mig(ji),mjg(jj))                   
478         END DO
479      END DO
480      CALL restget( inum, 'gphiv', jpidta, jpjdta, 1, 0, llog, zdta )
481      DO jj = 1, nlcj
482         DO ji = 1, nlci
483            gphiv(ji,jj) = zdta(mig(ji),mjg(jj))                   
484         END DO
485      END DO
486      CALL restget( inum, 'gphif', jpidta, jpjdta, 1, 0, llog, zdta )
487      DO jj = 1, nlcj
488         DO ji = 1, nlci
489            gphif(ji,jj) = zdta(mig(ji),mjg(jj))                   
490         END DO
491      END DO
492      CALL restget( inum, 'e1t', jpidta, jpjdta, 1, 0, llog, zdta )
493      DO jj = 1, nlcj
494         DO ji = 1, nlci
495            e1t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
496         END DO
497      END DO
498      CALL restget( inum, 'e1u', jpidta, jpjdta, 1, 0, llog, zdta )
499      DO jj = 1, nlcj
500         DO ji = 1, nlci
501            e1u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
502         END DO
503      END DO
504      CALL restget( inum, 'e1v', jpidta, jpjdta, 1, 0, llog, zdta )
505      DO jj = 1, nlcj
506         DO ji = 1, nlci
507            e1v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
508         END DO
509      END DO
510      CALL restget( inum, 'e1f', jpidta, jpjdta, 1, 0, llog, zdta )
511      DO jj = 1, nlcj
512         DO ji = 1, nlci
513            e1f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
514         END DO
515      END DO
516      CALL restget( inum, 'e2t', jpidta, jpjdta, 1, 0, llog, zdta )
517      DO jj = 1, nlcj
518         DO ji = 1, nlci
519            e2t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
520         END DO
521      END DO
522      CALL restget( inum, 'e2u', jpidta, jpjdta, 1, 0, llog, zdta )
523      DO jj = 1, nlcj
524         DO ji = 1, nlci
525            e2u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
526         END DO
527      END DO
528      CALL restget( inum, 'e2v', jpidta, jpjdta, 1, 0, llog, zdta )
529      DO jj = 1, nlcj
530         DO ji = 1, nlci
531            e2v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
532         END DO
533      END DO
534      CALL restget( inum, 'e2f', jpidta, jpjdta, 1, 0, llog, zdta )
535      DO jj = 1, nlcj
536         DO ji = 1, nlci
537            e2f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
538         END DO
539      END DO
540
541      CALL restclo( inum )
542
543      ! set extra rows add in mpp to none zero values
544      DO jj = nlcj+1, jpj
545         DO ji = 1, nlci
546            glamt(ji,jj) = glamt(ji,1)   ;   gphit(ji,jj) = gphit(ji,1)
547            glamu(ji,jj) = glamu(ji,1)   ;   gphiu(ji,jj) = gphiu(ji,1)
548            glamv(ji,jj) = glamv(ji,1)   ;   gphiv(ji,jj) = gphiv(ji,1)
549            glamf(ji,jj) = glamf(ji,1)   ;   gphif(ji,jj) = gphif(ji,1)
550            e1t  (ji,jj) = e1t  (ji,1)   ;   e2t  (ji,jj) = e2t  (ji,1)
551            e1u  (ji,jj) = e1u  (ji,1)   ;   e2u  (ji,jj) = e2u  (ji,1)
552            e1v  (ji,jj) = e1v  (ji,1)   ;   e2v  (ji,jj) = e2v  (ji,1)
553            e1f  (ji,jj) = e1f  (ji,1)   ;   e2f  (ji,jj) = e2f  (ji,1)
554         END DO
555      END DO
556
557      ! set extra columns add in mpp to none zero values
558      DO ji = nlci+1, jpi
559         glamt(ji,:) = glamt(1,:)   ;   gphit(ji,:) = gphit(1,:)
560         glamu(ji,:) = glamu(1,:)   ;   gphiu(ji,:) = gphiu(1,:)
561         glamv(ji,:) = glamv(1,:)   ;   gphiv(ji,:) = gphiv(1,:)
562         glamf(ji,:) = glamf(1,:)   ;   gphif(ji,:) = gphif(1,:)
563         e1t  (ji,:) = e1t  (1,:)   ;   e2t  (ji,:) = e2t  (1,:)
564         e1u  (ji,:) = e1u  (1,:)   ;   e2u  (ji,:) = e2u  (1,:)
565         e1v  (ji,:) = e1v  (1,:)   ;   e2v  (ji,:) = e2v  (1,:)
566         e1f  (ji,:) = e1f  (1,:)   ;   e2f  (ji,:) = e2f  (1,:)
567      END DO
568
569   END SUBROUTINE hgr_read
570
571
572   SUBROUTINE hgr_read_fdir
573      !!----------------------------------------------------------------------
574      !!                 ***  ROUTINE hgr_read_fdir  ***
575      !!
576      !!----------------------------------------------------------------------
577      !! * Local declarations
578      CHARACTER (len=5) ::   clfield
579      CHARACTER(len=21) ::   clname = 'coordinates'
580      INTEGER ::   ji, jj         ! dummy loop indices
581      INTEGER ::   inumcoo = 11   ! logical unit for coordinate file
582      INTEGER ::   ijpi, ijpj     ! temporary integers
583      REAL(wp), DIMENSION(jpi,jpj) ::   zdta   ! temporary workspace
584      !!----------------------------------------------------------------------
585
586
587      ! 1. Read of the grid coordinates and scale factors
588      ! -------------------------------------------------
589
590      IF(lwp) THEN
591         WRITE(numout,*)
592         WRITE(numout,*) 'hgrcoo : read the horizontal coordinates'
593         WRITE(numout,*) '~~~~~~'
594         WRITE(numout,*) '         jpiglo jpjglo jpk : ', jpiglo, jpjglo, jpk
595      ENDIF
596
597      ! open the file
598          CALL ctlopn( inumcoo, clname, 'OLD', 'UNFORMATTED', 'SEQUENTIAL',   &
599                       1      , numout       , lwp  , 1                            )
600
601      ! read the file
602      READ(inumcoo) ijpi,ijpj
603      IF( (ijpi /= jpidta) .OR. (ijpj /= jpjdta) ) THEN
604         IF(lwp) THEN
605            WRITE(numout,*)
606            WRITE(numout,*) '         inconsitency in reading coordinate file, unit=',inumcoo
607            WRITE(numout,*) '            jpidta = ',jpidta  ,' jpi  read = ',ijpi
608            WRITE(numout,*) '            jpjdta = ',jpjdta  ,' jpj  read = ',ijpj
609            WRITE(numout,*)
610         ENDIF
611         nstop = nstop + 1
612      ENDIF
613
614      READ(inumcoo) clfield, zdta
615      IF( clfield /= 'GLAMT' ) THEN
616         IF(lwp) THEN
617            WRITE(numout,cform_err)
618            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GLAMT'
619         ENDIF
620         nstop = nstop + 1
621      ENDIF
622      DO jj = 1, nlcj
623         DO ji = 1, nlci
624            glamt(ji,jj) = zdta(mig(ji),mjg(jj))
625         END DO
626      END DO
627      READ(inumcoo) clfield, zdta
628      IF(clfield /= 'GLAMU') THEN
629         IF(lwp) THEN
630            WRITE(numout,cform_err)
631            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GLAMU'
632         ENDIF
633         nstop = nstop + 1
634      ENDIF
635      DO jj = 1, nlcj
636         DO ji = 1, nlci
637            glamu(ji,jj) = zdta(mig(ji),mjg(jj))                   
638         END DO
639      END DO
640      READ(inumcoo) clfield, zdta
641      IF(clfield /= 'GLAMV') THEN
642         IF(lwp) THEN
643            WRITE(numout,cform_err)
644            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GLAMV'
645         ENDIF
646         nstop = nstop + 1
647      ENDIF
648      DO jj = 1, nlcj
649         DO ji = 1, nlci
650            glamv(ji,jj) = zdta(mig(ji),mjg(jj))                   
651         END DO
652      END DO
653      READ(inumcoo) clfield, zdta
654      IF(clfield /= 'GLAMF') THEN
655         IF(lwp) THEN
656            WRITE(numout,cform_err)
657            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GLAMF'
658         ENDIF
659         nstop = nstop + 1
660      ENDIF
661      DO jj = 1, nlcj
662         DO ji = 1, nlci
663            glamf(ji,jj) = zdta(mig(ji),mjg(jj))                   
664         END DO
665      END DO
666      READ(inumcoo) clfield, zdta
667      IF(clfield /= 'GPHIT') THEN
668         IF(lwp) THEN
669            WRITE(numout,cform_err)
670            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GPHIT'
671         ENDIF
672         nstop = nstop + 1
673      ENDIF
674      DO jj = 1, nlcj
675         DO ji = 1, nlci
676            gphit(ji,jj) = zdta(mig(ji),mjg(jj))                   
677         END DO
678      END DO
679      READ(inumcoo) clfield, zdta
680      IF(clfield /= 'GPHIU') THEN
681         IF(lwp) THEN
682            WRITE(numout,cform_err)
683            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GPHIU'
684         ENDIF
685         nstop = nstop + 1
686      ENDIF
687      DO jj = 1, nlcj
688         DO ji = 1, nlci
689            gphiu(ji,jj) = zdta(mig(ji),mjg(jj))                   
690         END DO
691      END DO
692      READ(inumcoo) clfield, zdta
693      IF(clfield /= 'GPHIV') THEN
694         IF(lwp) THEN
695            WRITE(numout,cform_err)
696            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GPHIV'
697         ENDIF
698         nstop = nstop + 1
699      ENDIF
700      DO jj = 1, nlcj
701         DO ji = 1, nlci
702            gphiv(ji,jj) = zdta(mig(ji),mjg(jj))                   
703         END DO
704      END DO
705      READ(inumcoo) clfield, zdta
706      IF(clfield /= 'GPHIF') THEN
707         IF(lwp) THEN
708            WRITE(numout,cform_err)
709            WRITE(numout,*) 'hgrcoo: bad read',clfield,' GPHIF'
710         ENDIF
711         nstop = nstop + 1
712      ENDIF
713      DO jj = 1, nlcj
714         DO ji = 1, nlci
715            gphif(ji,jj) = zdta(mig(ji),mjg(jj))                   
716         END DO
717      END DO
718      READ(inumcoo) clfield, zdta
719      IF(clfield /= 'E1T  ') THEN
720         IF(lwp) THEN
721            WRITE(numout,cform_err)
722            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E1T  '
723         ENDIF
724         nstop = nstop + 1
725      ENDIF
726      DO jj = 1, nlcj
727         DO ji = 1, nlci
728            e1t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
729         END DO
730      END DO
731      READ(inumcoo) clfield, zdta
732      IF(clfield /= 'E1U  ') THEN
733         IF(lwp) THEN
734            WRITE(numout,cform_err)
735            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E1U  '
736         ENDIF
737         nstop = nstop + 1
738      ENDIF
739      DO jj = 1, nlcj
740         DO ji = 1, nlci
741            e1u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
742         END DO
743      END DO
744      READ(inumcoo) clfield, zdta
745      IF(clfield /= 'E1V  ') THEN
746         IF(lwp) THEN
747            WRITE(numout,cform_err)
748            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E1V  '
749         ENDIF
750         nstop = nstop + 1
751      ENDIF
752      DO jj = 1, nlcj
753         DO ji = 1, nlci
754            e1v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
755         END DO
756      END DO
757      READ(inumcoo) clfield, zdta
758      IF(clfield /= 'E1F  ') THEN
759         IF(lwp) THEN
760            WRITE(numout,cform_err)
761            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E1F  '
762         ENDIF
763         nstop = nstop + 1
764      ENDIF
765      DO jj = 1, nlcj
766         DO ji = 1, nlci
767            e1f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
768         END DO
769      END DO
770      READ(inumcoo) clfield, zdta
771      IF(clfield /= 'E2T  ') THEN
772         IF(lwp) THEN
773            WRITE(numout,cform_err)
774            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E2T  '
775         ENDIF
776         nstop = nstop + 1
777      ENDIF
778      DO jj = 1, nlcj
779         DO ji = 1, nlci
780            e2t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
781         END DO
782      END DO
783      READ(inumcoo) clfield, zdta
784      IF(clfield /= 'E2U  ') THEN
785         IF(lwp) THEN
786            WRITE(numout,cform_err)
787            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E2U  '
788         ENDIF
789         nstop = nstop + 1
790      ENDIF
791      DO jj = 1, nlcj
792         DO ji = 1, nlci
793            e2u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
794         END DO
795      END DO
796      READ(inumcoo) clfield, zdta
797      IF(clfield /= 'E2V  ') THEN
798         IF(lwp) THEN
799            WRITE(numout,cform_err)
800            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E2V  '
801         ENDIF
802         nstop = nstop + 1
803      ENDIF
804      DO jj = 1, nlcj
805         DO ji = 1, nlci
806            e2v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
807         END DO
808      END DO
809      READ(inumcoo) clfield, zdta
810      IF(clfield /= 'E2F  ') THEN
811         IF(lwp) THEN
812            WRITE(numout,cform_err)
813            WRITE(numout,*) 'hgrcoo: bad read',clfield,' E2F  '
814         ENDIF
815         nstop = nstop + 1
816      ENDIF
817      DO jj = 1, nlcj
818         DO ji = 1, nlci
819            e2f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
820         END DO
821      END DO
822
823      CLOSE( inumcoo )
824
825      ! set extra rows add in mpp to none zero values
826      DO jj = nlcj+1, jpj
827         DO ji = 1, nlci
828            glamt(ji,jj) = glamt(ji,1)   ;   gphit(ji,jj) = gphit(ji,1)
829            glamu(ji,jj) = glamu(ji,1)   ;   gphiu(ji,jj) = gphiu(ji,1)
830            glamv(ji,jj) = glamv(ji,1)   ;   gphiv(ji,jj) = gphiv(ji,1)
831            glamf(ji,jj) = glamf(ji,1)   ;   gphif(ji,jj) = gphif(ji,1)
832            e1t  (ji,jj) = e1t  (ji,1)   ;   e2t  (ji,jj) = e2t  (ji,1)
833            e1u  (ji,jj) = e1u  (ji,1)   ;   e2u  (ji,jj) = e2u  (ji,1)
834            e1v  (ji,jj) = e1v  (ji,1)   ;   e2v  (ji,jj) = e2v  (ji,1)
835            e1f  (ji,jj) = e1f  (ji,1)   ;   e2f  (ji,jj) = e2f  (ji,1)
836         END DO
837      END DO
838
839      ! set extra columns add in mpp to none zero values
840      DO ji = nlci+1, jpi
841         glamt(ji,:) = glamt(1,:)   ;   gphit(ji,:) = gphit(1,:)
842         glamu(ji,:) = glamu(1,:)   ;   gphiu(ji,:) = gphiu(1,:)
843         glamv(ji,:) = glamv(1,:)   ;   gphiv(ji,:) = gphiv(1,:)
844         glamf(ji,:) = glamf(1,:)   ;   gphif(ji,:) = gphif(1,:)
845         e1t  (ji,:) = e1t  (1,:)   ;   e2t  (ji,:) = e2t  (1,:)
846         e1u  (ji,:) = e1u  (1,:)   ;   e2u  (ji,:) = e2u  (1,:)
847         e1v  (ji,:) = e1v  (1,:)   ;   e2v  (ji,:) = e2v  (1,:)
848         e1f  (ji,:) = e1f  (1,:)   ;   e2f  (ji,:) = e2f  (1,:)
849      END DO
850
851   END SUBROUTINE hgr_read_fdir
852
853   !!======================================================================
854END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.