source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 34.5 KB
Line 
1MODULE domhgr
2   !!==============================================================================
3   !!                       ***  MODULE domhgr   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  ! 1988-03  (G. Madec) Original code
7   !!            7.0  ! 1996-01  (G. Madec)  terrain following coordinates
8   !!            8.0  ! 1997-02  (G. Madec)  print mesh informations
9   !!            8.1  ! 1999-11  (M. Imbard) NetCDF format with IO-IPSL
10   !!            8.2  ! 2000-08  (D. Ludicone) Reduced section at Bab el Mandeb
11   !!             -   ! 2001-09  (M. Levy)  eel config: grid in km, beta-plane
12   !!  NEMO      1.0  ! 2002-08  (G. Madec)  F90: Free form and module, namelist
13   !!             -   ! 2004-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
14   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist
15   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration
16   !!            4.0  ! 2011-02  (G. Madec) add cell surface (e1e2t)
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   dom_hgr       : initialize the horizontal mesh
21   !!   hgr_read      : read "coordinate" NetCDF file
22   !!----------------------------------------------------------------------
23   USE dom_oce        ! ocean space and time domain
24   USE phycst         ! physical constants
25   USE in_out_manager ! I/O manager
26   USE lib_mpp        ! MPP library
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
33
34   PUBLIC   dom_hgr   ! called by domain.F90
35
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dom_hgr
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE dom_hgr  ***
46      !!
47      !! ** Purpose :   Compute the geographical position (in degre) of the
48      !!      model grid-points,  the horizontal scale factors (in meters) and
49      !!      the Coriolis factor (in s-1).
50      !!
51      !! ** Method  :   The geographical position of the model grid-points is
52      !!      defined from analytical functions, fslam and fsphi, the deriva-
53      !!      tives of which gives the horizontal scale factors e1,e2.
54      !!      Defining two function fslam and fsphi and their derivatives in
55      !!      the two horizontal directions (fse1 and fse2), the model grid-
56      !!      point position and scale factors are given by:
57      !!         t-point:
58      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
59      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
60      !!         u-point:
61      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
62      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
63      !!         v-point:
64      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
65      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
66      !!            f-point:
67      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
68      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
69      !!      Where fse1 and fse2 are defined by:
70      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
71      !!                                     +          di(fsphi) **2 )(i,j)
72      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
73      !!                                     +          dj(fsphi) **2 )(i,j)
74      !!
75      !!        The coriolis factor is given at z-point by:
76      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
77      !!
78      !!        This routine is given as an example, it must be modified
79      !!      following the user s desiderata. nevertheless, the output as
80      !!      well as the way to compute the model grid-point position and
81      !!      horizontal scale factors must be respected in order to insure
82      !!      second order accuracy schemes.
83      !!
84      !! N.B. If the domain is periodic, verify that scale factors are also
85      !!      periodic, and the coriolis term again.
86      !!
87      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
88      !!                u-, v- and f-points (in degre)
89      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
90      !!               u-, v-  and f-points (in degre)
91      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
92      !!      scale factors (in meters) at t-, u-, v-, and f-points.
93      !!        define ff: coriolis factor at f-point
94      !!
95      !! References :   Marti, Madec and Delecluse, 1992, JGR
96      !!                Madec, Imbard, 1996, Clim. Dyn.
97      !!----------------------------------------------------------------------
98      INTEGER  ::   ji, jj               ! dummy loop indices
99      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
100      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
101      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
102      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
103      REAL(wp) ::   zphi0, zbeta, znorme !
104      REAL(wp) ::   zarg, zf0, zminff, zmaxff
105      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
106      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
107      INTEGER  ::   isrow                ! index for ORCA1 starting row
108
109      !!----------------------------------------------------------------------
110      !
111      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
112      !
113      IF(lwp) THEN
114         WRITE(numout,*)
115         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
116         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
117         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
118         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
119         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
120         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
121         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
122         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
123      ENDIF
124
125
126      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
127
128      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
129
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
132
133         CALL hgr_read           ! Defaultl option  :   NetCDF file
134
135         !                                                ! =====================
136         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
137            !                                             ! =====================
138            IF( nn_cla == 0 ) THEN
139               !
140               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
141               ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
142               IF(lwp) WRITE(numout,*)
143               IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km'
144               !
145               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
146               ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
147                                               e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3
148               IF(lwp) WRITE(numout,*)
149               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km'
150               IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km'
151            ENDIF
152
153            ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km)
154            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
155            IF(lwp) WRITE(numout,*)
156            IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km'
157            !
158         ENDIF
159
160            !                                             ! =====================
161         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
162            !                                             ! =====================
163            ! This dirty section will be suppressed by simplification process: all this will come back in input files
164            ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration
165            ! which had a grid-size of 362x292.
166            ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6.
167            ! The original domain can still be used optionally if the isf code is not activated.
168            ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
169            !
170            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
171            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain.
172            ENDIF     
173            !
174            ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km)
175            ij0 = 201 + isrow   ;   ij1 = 201 + isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
176            IF(lwp) WRITE(numout,*)
177            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km'
178
179            ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km)
180            ij0 = 208 + isrow   ;   ij1 = 208 + isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
181            IF(lwp) WRITE(numout,*)
182            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km'
183
184            ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km)
185            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
186            IF(lwp) WRITE(numout,*)
187            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km'
188
189            ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
190            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
191            IF(lwp) WRITE(numout,*)
192            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km'
193
194            ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km)
195            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
196            IF(lwp) WRITE(numout,*)
197            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km'
198
199            ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km)
200            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
201            IF(lwp) WRITE(numout,*)
202            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km'
203
204            ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km)
205            ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
206            IF(lwp) WRITE(numout,*)
207            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km'
208
209            ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km)
210            ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
211            IF(lwp) WRITE(numout,*)
212            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km'
213            !
214            !
215         ENDIF
216
217         !                                                ! ======================
218         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
219            !                                             ! ======================
220            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
221            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
222            IF(lwp) WRITE(numout,*)
223            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
224            !
225            ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km)
226            ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
227            IF(lwp) WRITE(numout,*)
228            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait'
229            !
230            ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km)
231            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3
232            IF(lwp) WRITE(numout,*)
233            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait'
234            !
235            ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km)
236            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
237            IF(lwp) WRITE(numout,*)
238            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait'
239            !
240            ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km)
241            ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
242            IF(lwp) WRITE(numout,*)
243            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait'
244            !
245            ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km)
246            ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
247            IF(lwp) WRITE(numout,*)
248            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait'
249            !
250            !
251            ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km)
252            ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3
253            IF(lwp) WRITE(numout,*)
254            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb'
255            !
256         ENDIF
257
258
259         ! N.B. :  General case, lat and long function of both i and j indices:
260         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
261         !                                  + (                           fsdiph( zti, ztj ) )**2  )
262         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
263         !                                  + (                           fsdiph( zui, zuj ) )**2  )
264         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
265         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
266         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
267         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
268         !
269         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
270         !                                  + (                           fsdjph( zti, ztj ) )**2  )
271         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
272         !                                  + (                           fsdjph( zui, zuj ) )**2  )
273         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
274         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
275         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
276         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
277
278
279      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
280
281         IF(lwp) WRITE(numout,*)
282         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
283         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
284
285         DO jj = 1, jpj
286            DO ji = 1, jpi
287               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
288               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
289               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
290               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
291         ! Longitude
292               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
293               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
294               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
295               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
296         ! Latitude
297               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
298               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
299               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
300               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
301         ! e1
302               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
303               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
304               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
305               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
306         ! e2
307               e2t(ji,jj) = ra * rad * ppe2_deg
308               e2u(ji,jj) = ra * rad * ppe2_deg
309               e2v(ji,jj) = ra * rad * ppe2_deg
310               e2f(ji,jj) = ra * rad * ppe2_deg
311            END DO
312         END DO
313
314
315      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
316
317         IF(lwp) WRITE(numout,*)
318         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
319         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
320
321         ! Position coordinates (in kilometers)
322         !                          ==========
323         glam0 = 0.e0
324         gphi0 = - ppe2_m * 1.e-3
325         
326#if defined key_agrif 
327         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
328            IF( .NOT. Agrif_Root() ) THEN
329              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
330              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
331              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
332              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
333            ENDIF
334         ENDIF
335#endif         
336         DO jj = 1, jpj
337            DO ji = 1, jpi
338               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
339               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
340               glamv(ji,jj) = glamt(ji,jj)
341               glamf(ji,jj) = glamu(ji,jj)
342   
343               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
344               gphiu(ji,jj) = gphit(ji,jj)
345               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
346               gphif(ji,jj) = gphiv(ji,jj)
347            END DO
348         END DO
349
350         ! Horizontal scale factors (in meters)
351         !                              ======
352         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
353         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
354         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
355         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
356
357      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
358
359         IF(lwp) WRITE(numout,*)
360         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
361         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
362         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
363
364         !  Find index corresponding to the equator, given the grid spacing e1_deg
365         !  and the (approximate) southern latitude ppgphi0.
366         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
367         !  The formula should work even if the equator is outside the domain.
368         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
369         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
370         IF(  ppgphi0 > 0 )  ijeq = -ijeq
371
372         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
373
374         DO jj = 1, jpj
375            DO ji = 1, jpi
376               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
377               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
378               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
379               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
380         ! Longitude
381               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
382               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
383               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
384               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
385         ! Latitude
386               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
387               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
388               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
389               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
390         ! e1
391               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
392               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
393               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
394               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
395         ! e2
396               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
397               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
398               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
399               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
400            END DO
401         END DO
402
403      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
404
405         IF(lwp) WRITE(numout,*)
406         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
407         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
408
409         ! Position coordinates (in kilometers)
410         !                          ==========
411
412         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
413         zlam1 = -85
414         zphi1 = 29
415         ! resolution in meters
416         ze1 = 106000. / FLOAT(jp_cfg)           
417         ! benchmark: forced the resolution to be about 100 km
418         IF( nbench /= 0 )   ze1 = 106000.e0     
419         zsin_alpha = - SQRT( 2. ) / 2.
420         zcos_alpha =   SQRT( 2. ) / 2.
421         ze1deg = ze1 / (ra * rad)
422         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
423         !                                                          ! at the right jp_cfg resolution
424         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
425         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
426
427         IF( nprint==1 .AND. lwp )   THEN
428            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
429            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
430         ENDIF
431
432         DO jj = 1, jpj
433           DO ji = 1, jpi
434             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
435             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
436
437             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
438             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
439
440             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
441             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
442
443             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
444             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
445
446             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
447             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
448           END DO
449          END DO
450
451         ! Horizontal scale factors (in meters)
452         !                              ======
453         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
454         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
455         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
456         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
457
458      CASE DEFAULT
459         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
460         CALL ctl_stop( ctmp1 )
461
462      END SELECT
463     
464      ! T-cell surface
465      ! --------------
466      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
467   
468      ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning)
469      ! ---------------------------------------------------------------------------
470      e12t    (:,:) = e1t(:,:) * e2t(:,:)
471      e12u    (:,:) = e1u(:,:) * e2u(:,:)
472      e12v    (:,:) = e1v(:,:) * e2v(:,:)
473      e12f    (:,:) = e1f(:,:) * e2f(:,:)
474      r1_e12t (:,:) = 1._wp    / e12t(:,:)
475      r1_e12u (:,:) = 1._wp    / e12u(:,:)
476      r1_e12v (:,:) = 1._wp    / e12v(:,:)
477      r1_e12f (:,:) = 1._wp    / e12f(:,:)
478      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
479      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
480      r1_e1t  (:,:) = 1._wp    / e1t(:,:)
481      r1_e1u  (:,:) = 1._wp    / e1u(:,:)
482      r1_e1v  (:,:) = 1._wp    / e1v(:,:)
483      r1_e1f  (:,:) = 1._wp    / e1f(:,:)
484      r1_e2t  (:,:) = 1._wp    / e2t(:,:)
485      r1_e2u  (:,:) = 1._wp    / e2u(:,:)
486      r1_e2v  (:,:) = 1._wp    / e2v(:,:)
487      r1_e2f  (:,:) = 1._wp    / e2f(:,:)
488
489      ! Control printing : Grid informations (if not restart)
490      ! ----------------
491
492      IF( lwp .AND. .NOT.ln_rstart ) THEN
493         WRITE(numout,*)
494         WRITE(numout,*) '          longitude and e1 scale factors'
495         WRITE(numout,*) '          ------------------------------'
496         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
497            glamv(ji,1), glamf(ji,1),   &
498            e1t(ji,1), e1u(ji,1),   &
499            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
5009300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
501            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
502         
503         WRITE(numout,*)
504         WRITE(numout,*) '          latitude and e2 scale factors'
505         WRITE(numout,*) '          -----------------------------'
506         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
507            &                     gphiv(1,jj), gphif(1,jj),   &
508            &                     e2t  (1,jj), e2u  (1,jj),   &
509            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
510      ENDIF
511
512     
513      IF( nprint == 1 .AND. lwp ) THEN
514         WRITE(numout,*) '          e1u e2u '
515         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
516         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
517         WRITE(numout,*) '          e1v e2v  '
518         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
519         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
520         WRITE(numout,*) '          e1f e2f  '
521         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
522         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
523      ENDIF
524
525
526      ! ================= !
527      !  Coriolis factor  !
528      ! ================= !
529
530      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
531
532      CASE ( 0, 1, 4 )               ! mesh on the sphere
533
534         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
535
536      CASE ( 2 )                     ! f-plane at ppgphi0
537
538         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
539
540         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
541
542      CASE ( 3 )                     ! beta-plane
543
544         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
545         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
546         
547#if defined key_agrif
548         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
549            IF( .NOT. Agrif_Root() ) THEN
550              zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
551            ENDIF
552         ENDIF
553#endif         
554         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
555
556         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
557         
558         IF(lwp) THEN
559            WRITE(numout,*) 
560            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
561            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
562         ENDIF
563         IF( lk_mpp ) THEN
564            zminff=ff(nldi,nldj)
565            zmaxff=ff(nldi,nlej)
566            CALL mpp_min( zminff )   ! min over the global domain
567            CALL mpp_max( zmaxff )   ! max over the global domain
568            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
569         END IF
570
571      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
572
573         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
574         zphi0 = 15.e0                                                      ! latitude of the first row F-points
575         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
576
577         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
578
579         IF(lwp) THEN
580            WRITE(numout,*) 
581            WRITE(numout,*) '          Beta-plane and rotated domain : '
582            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
583         ENDIF
584
585         IF( lk_mpp ) THEN
586            zminff=ff(nldi,nldj)
587            zmaxff=ff(nldi,nlej)
588            CALL mpp_min( zminff )   ! min over the global domain
589            CALL mpp_max( zmaxff )   ! max over the global domain
590            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
591         END IF
592
593      END SELECT
594
595
596      ! Control of domain for symetrical condition
597      ! ------------------------------------------
598      ! The equator line must be the latitude coordinate axe
599
600      IF( nperio == 2 ) THEN
601         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
602         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
603      ENDIF
604      !
605      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
606      !
607   END SUBROUTINE dom_hgr
608
609
610   SUBROUTINE hgr_read
611      !!---------------------------------------------------------------------
612      !!              ***  ROUTINE hgr_read  ***
613      !!
614      !! ** Purpose :   Read a coordinate file in NetCDF format
615      !!
616      !! ** Method  :   The mesh file has been defined trough a analytical
617      !!      or semi-analytical method. It is read in a NetCDF file.
618      !!     
619      !!----------------------------------------------------------------------
620      USE iom
621
622      INTEGER ::   inum   ! temporary logical unit
623      !!----------------------------------------------------------------------
624
625      IF(lwp) THEN
626         WRITE(numout,*)
627         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
628         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
629      ENDIF
630     
631      CALL iom_open( 'coordinates', inum )
632     
633      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
634      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
635      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
636      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
637     
638      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
639      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
640      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
641      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
642     
643      CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr )
644      CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr )
645      CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr )
646      CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr )
647     
648      CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr )
649      CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr )
650      CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr )
651      CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr )
652     
653      CALL iom_close( inum )
654     
655! need to be define for the extended grid south of -80S
656! some point are undefined but you need to have e1 and e2 .NE. 0
657      WHERE (e1t==0.0_wp)
658         e1t=1.0e2
659      END WHERE
660      WHERE (e1v==0.0_wp)
661         e1v=1.0e2
662      END WHERE
663      WHERE (e1u==0.0_wp)
664         e1u=1.0e2
665      END WHERE
666      WHERE (e1f==0.0_wp)
667         e1f=1.0e2
668      END WHERE
669      WHERE (e2t==0.0_wp)
670         e2t=1.0e2
671      END WHERE
672      WHERE (e2v==0.0_wp)
673         e2v=1.0e2
674      END WHERE
675      WHERE (e2u==0.0_wp)
676         e2u=1.0e2
677      END WHERE
678      WHERE (e2f==0.0_wp)
679         e2f=1.0e2
680      END WHERE
681       
682    END SUBROUTINE hgr_read
683   
684   !!======================================================================
685END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.