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

source: trunk/NEMO/OPA_SRC/DOM/domhgr.F90 @ 3

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

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.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   !!   dom_hgr_coo    : read "coordinate" file (except for EEL config.)
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE dom_oce         ! ocean space and time domain
13   USE phycst          ! physical constants
14   USE in_out_manager  ! I/O manager
15
16   IMPLICIT NONE
17   PRIVATE
18
19   !! * Module variables
20   REAL (wp)  :: glam0, gphi0        ! variables corresponding to parameters
21      !                              ! ppglam0 ppgphi0 set in par_oce
22
23   !! * Routine accessibility
24   PUBLIC dom_hgr        ! called by domain.F90
25   !!----------------------------------------------------------------------
26   !!   OPA 9.0 , LODYC-IPSL   (2003)
27   !!----------------------------------------------------------------------
28
29CONTAINS
30
31   SUBROUTINE dom_hgr
32      !!----------------------------------------------------------------------
33      !!                  ***  ROUTINE dom_hgr  ***
34      !!
35      !! ** Purpose :   Compute the geographical position (in degre) of the
36      !!      model grid-points,  the horizontal scale factors (in meters) and
37      !!      the Coriolis factor (in s-1).
38      !!
39      !! ** Method  :   The geographical position of the model grid-points is
40      !!      defined from analytical functions, fslam and fsphi, the deriva-
41      !!      tives of which gives the horizontal scale factors e1,e2.
42      !!      Defining two function fslam and fsphi and their derivatives in
43      !!      the two horizontal directions (fse1 and fse2), the model grid-
44      !!      point position and scale factors are given by:
45      !!            t-point: glamt(i,j) = fslam(i,j)
46      !!                     gphit(i,j) = fsphi(i,j)
47      !!                     e1t  (i,j) = fse1 (i,j)
48      !!                     e2t  (i,j) = fse2 (i,j)
49      !!            u-point: glamu(i,j) = fslam(i+0.5,j)
50      !!                     gphiu(i,j) = fsphi(i+0.5,j)
51      !!                     e1u  (i,j) = fse1 (i+0.5,j)
52      !!                     e2u  (i,j) = fse2 (i+0.5,j)
53      !!            v-point: glamv(i,j) = fslam(i,j+0.5)
54      !!                     gphiv(i,j) = fsphi(i,j+0.5)
55      !!                     e1v  (i,j) = fse1 (i,j+0.5)
56      !!                     e2v  (i,j) = fse2 (i,j+0.5)
57      !!            f-point: glamf(i,j) = fslam(i+0.5,j+0.5)
58      !!                     gphif(i,j) = fsphi(i+0.5,j+0.5)
59      !!                     e1f  (i,j) = fse1 (i+0.5,j+0.5)
60      !!                     e2f  (i,j) = fse2 (i+0.5,j+0.5)
61      !!      Where fse1 and fse2 are defined by:
62      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
63      !!                                     +          di(fsphi) **2 )(i,j)
64      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
65      !!                                     +          dj(fsphi) **2 )(i,j)
66      !!
67      !!        The coriolis factor is given at z-point by:
68      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
69      !!
70      !!        This routine is given as an example, it must be modified
71      !!      following the user s desiderata. nevertheless, the output as
72      !!      well as the way to compute the model grid-point position and
73      !!      horizontal scale factors must be respected in order to insure
74      !!      second order accuracy schemes.
75      !!
76      !! N.B. If the domain is periodic, verify that scale factors are also
77      !!      periodic, and the coriolis term again.
78      !!
79      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
80      !!                u-, v- and f-points (in degre)
81      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
82      !!               u-, v-  and f-points (in degre)
83      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
84      !!      scale factors (in meters) at t-, u-, v-, and f-points.
85      !!        define ff: coriolis factor at f-point
86      !!
87      !! References :
88      !!      Marti, Madec and Delecluse, 1992, j. geophys. res., in press.
89      !!
90      !! History :
91      !!        !  88-03  (G. Madec)
92      !!        !  91-11  (G. Madec)
93      !!        !  92-06  (M. Imbard)
94      !!        !  96-01  (G. Madec)  terrain following coordinates
95      !!        !  97-02  (G. Madec)  print mesh informations
96      !!        !  01-09  (M. Levy)  eel config: grid in km, beta-plane
97      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, namelist
98      !!   9.0  !  04-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
99      !!                  use of parameters in par_CONFIG-Rxx.h90, not in namelist
100      !!----------------------------------------------------------------------
101      !! * local declarations
102      INTEGER  ::   ji, jj          ! dummy loop indices
103      INTEGER  ::   jeq             ! index of equator T point (computed for case 4)
104      REAL(wp) ::   &
105         zti, zui, zvi, zfi,     &  ! temporary scalars
106         ztj, zuj, zvj, zfj,     &  !
107         zphi0, zbeta, znorme,   &  !
108         zarg, zf0
109      !!----------------------------------------------------------------------
110
111      IF(lwp) THEN
112         WRITE(numout,*)
113         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
114         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
115         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
116         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
117         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
118         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
119         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
120         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
121      ENDIF
122
123
124      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
125
126      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
127
128         IF(lwp) WRITE(numout,*)
129         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in coordinate.nc file'
130
131         CALL dom_hgr_coo
132
133         ! N.B. :  General case, lat and long function of both i and j indices:
134         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
135         !                                  + (                           fsdiph( zti, ztj ) )**2  )
136         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
137         !                                  + (                           fsdiph( zui, zuj ) )**2  )
138         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
139         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
140         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
141         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
142         !
143         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
144         !                                  + (                           fsdjph( zti, ztj ) )**2  )
145         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
146         !                                  + (                           fsdjph( zui, zuj ) )**2  )
147         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
148         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
149         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
150         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
151
152
153      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
154
155         IF(lwp) WRITE(numout,*)
156         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
157         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
158
159         DO jj = 1, jpj
160            DO ji = 1, jpi
161               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
162               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
163               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
164               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
165         ! Longitude
166               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
167               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
168               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
169               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
170         ! Latitude
171               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
172               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
173               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
174               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
175         ! e1
176               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
177               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
178               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
179               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
180         ! e2
181               e2t(ji,jj) = ra * rad * ppe2_deg
182               e2u(ji,jj) = ra * rad * ppe2_deg
183               e2v(ji,jj) = ra * rad * ppe2_deg
184               e2f(ji,jj) = ra * rad * ppe2_deg
185            END DO
186         END DO
187
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.e0
198         gphi0 = - ppe2_m * 1.e-3
199         DO jj = 1, jpj
200            DO ji = 1, jpi
201               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
202               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
203               glamv(ji,jj) = glamt(ji,jj)
204               glamf(ji,jj) = glamu(ji,jj)
205   
206               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
207               gphiu(ji,jj) = gphit(ji,jj)
208               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
209               gphif(ji,jj) = gphiv(ji,jj)
210            END DO
211         END DO
212
213         ! Horizontal scale factors (in meters)
214         !                              ======
215         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
216         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
217         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
218         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
219
220      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
221
222         IF(lwp) WRITE(numout,*)
223         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
224         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
225         IF ( ppgphi0 == -90 ) THEN
226                IF(lwp) WRITE(numout,*) ' Mercator grid cannot start at south pole !!!! '
227                IF(lwp) WRITE(numout,*) ' We stop '
228                STOP
229         ENDIF
230
231         !  Find index corresponding to the equator, given the grid spacing e1_deg
232         !  and the (approximate) southern latitude ppgphi0.
233         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
234         !  The formula should work even if the equator is outside the domain.
235         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
236         jeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
237
238         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', jeq
239
240         DO jj = 1, jpj
241            DO ji = 1, jpi
242               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - jeq + njmpp - 1 )
243               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - jeq + njmpp - 1 )
244               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - jeq + njmpp - 1 ) + 0.5
245               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - jeq + njmpp - 1 ) + 0.5
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 ) )
253               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
254               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
255               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
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
268
269
270
271      CASE DEFAULT
272         IF(lwp) WRITE(numout,cform_err)
273         IF(lwp) WRITE(numout,*) '          bad flag value for jphgr_msh = ', jphgr_msh
274         nstop = nstop + 1
275
276      END SELECT
277
278
279      ! Control printing : Grid informations (if not restart)
280      ! ----------------
281
282      IF(lwp .AND. .NOT.ln_rstart ) THEN
283         WRITE(numout,*)
284         WRITE(numout,*) '          longitude and e1 scale factors'
285         WRITE(numout,*) '          ------------------------------'
286         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
287            glamv(ji,1), glamf(ji,1),   &
288            e1t(ji,1), e1u(ji,1),   &
289            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
2909300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
291            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
292         
293         WRITE(numout,*)
294         WRITE(numout,*) '          latitude and e2 scale factors'
295         WRITE(numout,*) '          -----------------------------'
296         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
297            &                     gphiv(1,jj), gphif(1,jj),   &
298            &                     e2t  (1,jj), e2u  (1,jj),   &
299            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
300      ENDIF
301
302     
303      IF( nprint == 1 .AND. lwp ) THEN
304         WRITE(numout,*) '          e1u e2u '
305         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
306         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
307         WRITE(numout,*) '          e1v e2v  '
308         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
309         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
310         WRITE(numout,*) '          e1f e2f  '
311         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
312         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
313      ENDIF
314
315
316      ! ================= !
317      !  Coriolis factor  !
318      ! ================= !
319
320      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
321
322      CASE ( 0, 1, 4 )               ! mesh on the sphere
323
324         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
325
326      CASE ( 2 )                     ! f-plane at ppgphi0
327
328         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
329
330         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
331
332      CASE ( 3 )                     ! beta-plane
333
334         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
335         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m *1.e-3  / ( ra * rad ) ! latitude of the first row F-points
336         zf0     = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
337
338         ff(:,:) = ( zf0  + zbeta * gphif(:,:) )                          ! f = f0 +beta* y ( y=0 at south)
339
340         IF(lwp) WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(1,1)
341         IF(lwp) WRITE(numout,*) '                      Coriolis parameter varies from ', ff(1,1),' to ', ff(1,jpj)
342
343      END SELECT
344
345
346      ! Control of domain for symetrical condition
347      ! ------------------------------------------
348      ! The equator line must be the latitude coordinate axe
349
350      IF( nperio == 2 ) THEN
351         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
352         IF( znorme > 1.e-13 ) THEN
353            IF(lwp) WRITE(numout,cform_err)
354            IF(lwp) WRITE(numout,*) ' '
355            IF(lwp) WRITE(numout,*) ' ===>>>> : symmetrical condition'
356            IF(lwp) WRITE(numout,*) 'stop rerun with good equator line'
357            IF(lwp) WRITE(numout,*) '           ----'
358            nstop = nstop + 1
359         ENDIF
360      ENDIF
361
362   END SUBROUTINE dom_hgr
363
364#if defined key_fdir
365   !!---------------------------------------------------------------------
366   !!   'key-fdir                                       direct access file
367   !!---------------------------------------------------------------------
368#  include "domhgr_coo_fdir.h90"
369
370#else
371   !!---------------------------------------------------------------------
372   !!   Defaultl option  :                                     NetCDF file
373   !!---------------------------------------------------------------------
374
375   SUBROUTINE dom_hgr_coo
376      !!---------------------------------------------------------------------
377      !!              ***  ROUTINE dom_hgr_coo  ***
378      !!
379      !! ** Purpose :   Read a coordinate file in NetCDF format
380      !!
381      !! ** Method  :   The mesh file has been defined trough a analytical
382      !!      or semi-analytical method. It is read in a NetCDF file.
383      !!     
384      !! References :
385      !!      Marti, Madec and Delecluse, 1992, JGR, 97, 12,763-12,766.
386      !!      Madec, Imbard, 1996, Clim. Dyn., 12, 381-388.
387      !!
388      !! History :
389      !!        !         (O. Marti)  Original code
390      !!        !  91-03  (G. Madec)
391      !!        !  92-07  (M. Imbard)
392      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
393      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
394      !!   8.5  !  02-06  (G. Madec)  F90: Free form
395      !!----------------------------------------------------------------------
396      !! * Modules used
397      USE ioipsl
398
399      !! * Local declarations
400      LOGICAL ::   llog
401      CHARACTER(len=21) ::   clname
402      INTEGER  ::   ji, jj              ! dummy loop indices
403      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
404      INTEGER  ::   inum                ! temporary logical unit
405      INTEGER  ::   ilev, itime         ! temporary integers
406      REAL(wp) ::   zdt, zdate0         ! temporary scalars
407      REAL(wp) ::   zdept(1)            ! temporary workspace
408      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
409         zlamt, zphit, zdta             ! temporary workspace
410      !!----------------------------------------------------------------------
411
412
413      ! 1. Read of the grid coordinates and scale factors
414      ! -------------------------------------------------
415
416      IF(lwp) THEN
417         WRITE(numout,*)
418         WRITE(numout,*) 'dom_hgr_coo : read the horizontal coordinates'
419         WRITE(numout,*) '~~~~~~~~~~~'
420         WRITE(numout,*) '         jpiglo jpjglo jpk : ', jpiglo, jpjglo, jpk
421      ENDIF
422
423      ! read the file
424      itime = 0
425      clname = 'coordinates'
426      llog = .FALSE.
427      ilev = 1   
428      zlamt(:,:) = 0.e0
429      zphit(:,:) = 0.e0
430      CALL restini( clname,jpidta,jpjdta,zlamt,zphit,ilev,zdept,clname   &
431         &          ,itime,zdate0,zdt,inum)
432
433      CALL restget(inum,'glamt',jpidta,jpjdta,1,0,llog,zdta)
434      DO jj = 1, nlcj
435         DO ji = 1, nlci
436            glamt(ji,jj) = zdta(mig(ji),mjg(jj))
437         END DO
438      END DO
439      CALL restget(inum,'glamu',jpidta,jpjdta,1,0,llog,zdta)
440      DO jj = 1, nlcj
441         DO ji = 1, nlci
442            glamu(ji,jj) = zdta(mig(ji),mjg(jj))                   
443         END DO
444      END DO
445      CALL restget(inum,'glamv',jpidta,jpjdta,1,0,llog,zdta)
446      DO jj = 1, nlcj
447         DO ji = 1, nlci
448            glamv(ji,jj) = zdta(mig(ji),mjg(jj))                   
449         END DO
450      END DO
451      CALL restget(inum,'glamf',jpidta,jpjdta,1,0,llog,zdta)
452      DO jj = 1, nlcj
453         DO ji = 1, nlci
454            glamf(ji,jj) = zdta(mig(ji),mjg(jj))                   
455         END DO
456      END DO
457      CALL restget(inum,'gphit',jpidta,jpjdta,1,0,llog,zdta)
458      DO jj = 1, nlcj
459         DO ji = 1, nlci
460            gphit(ji,jj) = zdta(mig(ji),mjg(jj))                   
461         END DO
462      END DO
463      CALL restget(inum,'gphiu',jpidta,jpjdta,1,0,llog,zdta)
464      DO jj = 1, nlcj
465         DO ji = 1, nlci
466            gphiu(ji,jj) = zdta(mig(ji),mjg(jj))                   
467         END DO
468      END DO
469      CALL restget(inum,'gphiv',jpidta,jpjdta,1,0,llog,zdta)
470      DO jj = 1, nlcj
471         DO ji = 1, nlci
472            gphiv(ji,jj) = zdta(mig(ji),mjg(jj))                   
473         END DO
474      END DO
475      CALL restget(inum,'gphif',jpidta,jpjdta,1,0,llog,zdta)
476      DO jj = 1, nlcj
477         DO ji = 1, nlci
478            gphif(ji,jj) = zdta(mig(ji),mjg(jj))                   
479         END DO
480      END DO
481      CALL restget(inum,'e1t',jpidta,jpjdta,1,0,llog,zdta)
482      DO jj = 1, nlcj
483         DO ji = 1, nlci
484            e1t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
485         END DO
486      END DO
487      CALL restget(inum,'e1u',jpidta,jpjdta,1,0,llog,zdta)
488      DO jj = 1, nlcj
489         DO ji = 1, nlci
490            e1u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
491         END DO
492      END DO
493      CALL restget(inum,'e1v',jpidta,jpjdta,1,0,llog,zdta)
494      DO jj = 1, nlcj
495         DO ji = 1, nlci
496            e1v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
497         END DO
498      END DO
499      CALL restget(inum,'e1f',jpidta,jpjdta,1,0,llog,zdta)
500      DO jj = 1, nlcj
501         DO ji = 1, nlci
502            e1f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
503         END DO
504      END DO
505      CALL restget(inum,'e2t',jpidta,jpjdta,1,0,llog,zdta)
506      DO jj = 1, nlcj
507         DO ji = 1, nlci
508            e2t  (ji,jj) = zdta(mig(ji),mjg(jj))                   
509         END DO
510      END DO
511      CALL restget(inum,'e2u',jpidta,jpjdta,1,0,llog,zdta)
512      DO jj = 1, nlcj
513         DO ji = 1, nlci
514            e2u  (ji,jj) = zdta(mig(ji),mjg(jj))                   
515         END DO
516      END DO
517      CALL restget(inum,'e2v',jpidta,jpjdta,1,0,llog,zdta)
518      DO jj = 1, nlcj
519         DO ji = 1, nlci
520            e2v  (ji,jj) = zdta(mig(ji),mjg(jj))                   
521         END DO
522      END DO
523      CALL restget(inum,'e2f',jpidta,jpjdta,1,0,llog,zdta)
524      DO jj = 1, nlcj
525         DO ji = 1, nlci
526            e2f  (ji,jj) = zdta(mig(ji),mjg(jj))                   
527         END DO
528      END DO
529
530      CALL restclo(inum)
531
532      ! set extra rows add in mpp to none zero values
533      DO jj = nlcj+1, jpj
534         DO ji = 1, nlci
535            glamt(ji,jj) = glamt(ji,1)   ;   gphit(ji,jj) = gphit(ji,1)
536            glamu(ji,jj) = glamu(ji,1)   ;   gphiu(ji,jj) = gphiu(ji,1)
537            glamv(ji,jj) = glamv(ji,1)   ;   gphiv(ji,jj) = gphiv(ji,1)
538            glamf(ji,jj) = glamf(ji,1)   ;   gphif(ji,jj) = gphif(ji,1)
539            e1t  (ji,jj) = e1t  (ji,1)   ;   e2t  (ji,jj) = e2t  (ji,1)
540            e1u  (ji,jj) = e1u  (ji,1)   ;   e2u  (ji,jj) = e2u  (ji,1)
541            e1v  (ji,jj) = e1v  (ji,1)   ;   e2v  (ji,jj) = e2v  (ji,1)
542            e1f  (ji,jj) = e1f  (ji,1)   ;   e2f  (ji,jj) = e2f  (ji,1)
543         END DO
544      END DO
545
546      ! set extra columns add in mpp to none zero values
547      DO ji = nlci+1, jpi
548         glamt(ji,:) = glamt(1,:)   ;   gphit(ji,:) = gphit(1,:)
549         glamu(ji,:) = glamu(1,:)   ;   gphiu(ji,:) = gphiu(1,:)
550         glamv(ji,:) = glamv(1,:)   ;   gphiv(ji,:) = gphiv(1,:)
551         glamf(ji,:) = glamf(1,:)   ;   gphif(ji,:) = gphif(1,:)
552         e1t  (ji,:) = e1t  (1,:)   ;   e2t  (ji,:) = e2t  (1,:)
553         e1u  (ji,:) = e1u  (1,:)   ;   e2u  (ji,:) = e2u  (1,:)
554         e1v  (ji,:) = e1v  (1,:)   ;   e2v  (ji,:) = e2v  (1,:)
555         e1f  (ji,:) = e1f  (1,:)   ;   e2f  (ji,:) = e2f  (1,:)
556      END DO
557
558      !                                                ! =====================
559      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
560         !                                             ! =====================
561         IF( n_cla == 0 ) THEN
562            e2u( mi0(160):mi1(161) , mj0(88):mj1(88) ) =  18.e3   ! Bab el Mandeb (e2u = 18 km)
563            IF(lwp) WRITE(numout,*)
564            IF(lwp) WRITE(numout,*) '          Bab el Mandeb: e2u reduced to 18 km'
565         ENDIF
566 
567         e2u( mi0(145):mi1(146) , mj0(116):mj1(116) ) =  15.e3    ! Sound Strait (e2u = 15 km)
568         IF(lwp) WRITE(numout,*)
569         IF(lwp) WRITE(numout,*) '        : Reduced e2u at the Sound Strait'
570         !
571      ENDIF
572
573      !                                                ! ======================
574      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
575         !                                             ! ======================
576         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
577         ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij0) ) =  20.e3   
578         IF(lwp) WRITE(numout,*)
579         IF(lwp) WRITE(numout,*) '        : Reduced e2u at the Gibraltar Strait'
580         !
581      ENDIF
582
583   END SUBROUTINE dom_hgr_coo
584
585#endif
586   !!======================================================================
587END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.