New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domhgr.F90 in branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago

Dev. branch CMCC4_simplification ticket #1456

  • Property svn:keywords set to Id
File size: 32.7 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      !!----------------------------------------------------------------------
108      !
109      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
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" file'
130
131         CALL hgr_read           ! Defaultl option  :   NetCDF file
132
133         !                                                ! =====================
134         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
135            !                                             ! =====================
136               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
137               ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
138               IF(lwp) WRITE(numout,*)
139               IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km'
140               !
141               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
142               ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
143                                               e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3
144               IF(lwp) WRITE(numout,*)
145               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km'
146               IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km'
147
148            ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km)
149            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
150            IF(lwp) WRITE(numout,*)
151            IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km'
152            !
153         ENDIF
154
155            !                                             ! =====================
156         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
157            !                                             ! =====================
158
159            ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u = 20 km)
160            ij0 = 200   ;   ij1 = 200   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
161            IF(lwp) WRITE(numout,*)
162            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km'
163
164            ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km)
165            ij0 = 208   ;   ij1 = 208   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
166            IF(lwp) WRITE(numout,*)
167            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km'
168
169            ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v = 13 km)
170            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
171            IF(lwp) WRITE(numout,*)
172            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km'
173
174            ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
175            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
176            IF(lwp) WRITE(numout,*)
177            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km'
178
179            ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v = 13 km)
180            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
181            IF(lwp) WRITE(numout,*)
182            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km'
183
184            ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v = 20 km)
185            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
186            IF(lwp) WRITE(numout,*)
187            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km'
188
189            ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km)
190            ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
191            IF(lwp) WRITE(numout,*)
192            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km'
193
194            ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km)
195            ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
196            IF(lwp) WRITE(numout,*)
197            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km'
198
199            !
200
201            !
202            !
203            !
204            !
205         ENDIF
206
207         !                                                ! ======================
208         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
209            !                                             ! ======================
210            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
211            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
212            IF(lwp) WRITE(numout,*)
213            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
214            !
215            ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km)
216            ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
217            IF(lwp) WRITE(numout,*)
218            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait'
219            !
220            ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km)
221            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3
222            IF(lwp) WRITE(numout,*)
223            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait'
224            !
225            ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km)
226            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
227            IF(lwp) WRITE(numout,*)
228            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait'
229            !
230            ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km)
231            ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
232            IF(lwp) WRITE(numout,*)
233            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait'
234            !
235            ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km)
236            ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
237            IF(lwp) WRITE(numout,*)
238            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait'
239            !
240            !
241            ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km)
242            ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3
243            IF(lwp) WRITE(numout,*)
244            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb'
245            !
246         ENDIF
247
248
249         ! N.B. :  General case, lat and long function of both i and j indices:
250         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
251         !                                  + (                           fsdiph( zti, ztj ) )**2  )
252         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
253         !                                  + (                           fsdiph( zui, zuj ) )**2  )
254         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
255         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
256         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
257         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
258         !
259         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
260         !                                  + (                           fsdjph( zti, ztj ) )**2  )
261         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
262         !                                  + (                           fsdjph( zui, zuj ) )**2  )
263         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
264         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
265         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
266         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
267
268
269      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
270
271         IF(lwp) WRITE(numout,*)
272         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
273         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
274
275         DO jj = 1, jpj
276            DO ji = 1, jpi
277               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
278               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
279               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
280               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
281         ! Longitude
282               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
283               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
284               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
285               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
286         ! Latitude
287               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
288               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
289               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
290               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
291         ! e1
292               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
293               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
294               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
295               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
296         ! e2
297               e2t(ji,jj) = ra * rad * ppe2_deg
298               e2u(ji,jj) = ra * rad * ppe2_deg
299               e2v(ji,jj) = ra * rad * ppe2_deg
300               e2f(ji,jj) = ra * rad * ppe2_deg
301            END DO
302         END DO
303
304
305      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
306
307         IF(lwp) WRITE(numout,*)
308         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
309         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
310
311         ! Position coordinates (in kilometers)
312         !                          ==========
313         glam0 = 0.e0
314         gphi0 = - ppe2_m * 1.e-3
315         
316#if defined key_agrif 
317         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
318            IF( .NOT. Agrif_Root() ) THEN
319              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
320              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
321              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
322              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
323            ENDIF
324         ENDIF
325#endif         
326         DO jj = 1, jpj
327            DO ji = 1, jpi
328               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
329               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
330               glamv(ji,jj) = glamt(ji,jj)
331               glamf(ji,jj) = glamu(ji,jj)
332   
333               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
334               gphiu(ji,jj) = gphit(ji,jj)
335               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
336               gphif(ji,jj) = gphiv(ji,jj)
337            END DO
338         END DO
339
340         ! Horizontal scale factors (in meters)
341         !                              ======
342         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
343         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
344         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
345         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
346
347      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
348
349         IF(lwp) WRITE(numout,*)
350         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
351         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
352         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
353
354         !  Find index corresponding to the equator, given the grid spacing e1_deg
355         !  and the (approximate) southern latitude ppgphi0.
356         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
357         !  The formula should work even if the equator is outside the domain.
358         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
359         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
360         IF(  ppgphi0 > 0 )  ijeq = -ijeq
361
362         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
363
364         DO jj = 1, jpj
365            DO ji = 1, jpi
366               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
367               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
368               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
369               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
370         ! Longitude
371               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
372               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
373               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
374               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
375         ! Latitude
376               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
377               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
378               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
379               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
380         ! e1
381               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
382               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
383               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
384               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
385         ! e2
386               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
387               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
388               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
389               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
390            END DO
391         END DO
392
393      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
394
395         IF(lwp) WRITE(numout,*)
396         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
397         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
398
399         ! Position coordinates (in kilometers)
400         !                          ==========
401
402         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
403         zlam1 = -85
404         zphi1 = 29
405         ! resolution in meters
406         ze1 = 106000. / FLOAT(jp_cfg)           
407         ! benchmark: forced the resolution to be about 100 km
408         IF( nbench /= 0 )   ze1 = 106000.e0     
409         zsin_alpha = - SQRT( 2. ) / 2.
410         zcos_alpha =   SQRT( 2. ) / 2.
411         ze1deg = ze1 / (ra * rad)
412         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
413         !                                                          ! at the right jp_cfg resolution
414         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
415         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
416
417         IF( nprint==1 .AND. lwp )   THEN
418            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
419            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
420         ENDIF
421
422         DO jj = 1, jpj
423           DO ji = 1, jpi
424             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
425             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
426
427             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
428             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
429
430             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
431             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
432
433             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
434             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
435
436             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
437             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
438           END DO
439          END DO
440
441         ! Horizontal scale factors (in meters)
442         !                              ======
443         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
444         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
445         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
446         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
447
448      CASE DEFAULT
449         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
450         CALL ctl_stop( ctmp1 )
451
452      END SELECT
453     
454      ! T-cell surface
455      ! --------------
456      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
457   
458      ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning)
459      ! ---------------------------------------------------------------------------
460      e12t    (:,:) = e1t(:,:) * e2t(:,:)
461      e12u    (:,:) = e1u(:,:) * e2u(:,:)
462      e12v    (:,:) = e1v(:,:) * e2v(:,:)
463      e12f    (:,:) = e1f(:,:) * e2f(:,:)
464      r1_e12t (:,:) = 1._wp    / e12t(:,:)
465      r1_e12u (:,:) = 1._wp    / e12u(:,:)
466      r1_e12v (:,:) = 1._wp    / e12v(:,:)
467      r1_e12f (:,:) = 1._wp    / e12f(:,:)
468      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
469      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
470
471      ! Control printing : Grid informations (if not restart)
472      ! ----------------
473
474      IF( lwp .AND. .NOT.ln_rstart ) THEN
475         WRITE(numout,*)
476         WRITE(numout,*) '          longitude and e1 scale factors'
477         WRITE(numout,*) '          ------------------------------'
478         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
479            glamv(ji,1), glamf(ji,1),   &
480            e1t(ji,1), e1u(ji,1),   &
481            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4829300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
483            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
484         
485         WRITE(numout,*)
486         WRITE(numout,*) '          latitude and e2 scale factors'
487         WRITE(numout,*) '          -----------------------------'
488         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
489            &                     gphiv(1,jj), gphif(1,jj),   &
490            &                     e2t  (1,jj), e2u  (1,jj),   &
491            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
492      ENDIF
493
494     
495      IF( nprint == 1 .AND. lwp ) THEN
496         WRITE(numout,*) '          e1u e2u '
497         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
498         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
499         WRITE(numout,*) '          e1v e2v  '
500         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
501         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
502         WRITE(numout,*) '          e1f e2f  '
503         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
504         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
505      ENDIF
506
507
508      ! ================= !
509      !  Coriolis factor  !
510      ! ================= !
511
512      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
513
514      CASE ( 0, 1, 4 )               ! mesh on the sphere
515
516         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
517
518      CASE ( 2 )                     ! f-plane at ppgphi0
519
520         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
521
522         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
523
524      CASE ( 3 )                     ! beta-plane
525
526         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
527         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
528         
529#if defined key_agrif
530         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
531            IF( .NOT. Agrif_Root() ) THEN
532              zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
533            ENDIF
534         ENDIF
535#endif         
536         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
537
538         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
539         
540         IF(lwp) THEN
541            WRITE(numout,*) 
542            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
543            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
544         ENDIF
545         IF( lk_mpp ) THEN
546            zminff=ff(nldi,nldj)
547            zmaxff=ff(nldi,nlej)
548            CALL mpp_min( zminff )   ! min over the global domain
549            CALL mpp_max( zmaxff )   ! max over the global domain
550            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
551         END IF
552
553      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
554
555         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
556         zphi0 = 15.e0                                                      ! latitude of the first row F-points
557         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
558
559         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
560
561         IF(lwp) THEN
562            WRITE(numout,*) 
563            WRITE(numout,*) '          Beta-plane and rotated domain : '
564            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
565         ENDIF
566
567         IF( lk_mpp ) THEN
568            zminff=ff(nldi,nldj)
569            zmaxff=ff(nldi,nlej)
570            CALL mpp_min( zminff )   ! min over the global domain
571            CALL mpp_max( zmaxff )   ! max over the global domain
572            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
573         END IF
574
575      END SELECT
576
577
578      ! Control of domain for symetrical condition
579      ! ------------------------------------------
580      ! The equator line must be the latitude coordinate axe
581
582      IF( nperio == 2 ) THEN
583         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
584         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
585      ENDIF
586      !
587      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
588      !
589   END SUBROUTINE dom_hgr
590
591
592   SUBROUTINE hgr_read
593      !!---------------------------------------------------------------------
594      !!              ***  ROUTINE hgr_read  ***
595      !!
596      !! ** Purpose :   Read a coordinate file in NetCDF format
597      !!
598      !! ** Method  :   The mesh file has been defined trough a analytical
599      !!      or semi-analytical method. It is read in a NetCDF file.
600      !!     
601      !!----------------------------------------------------------------------
602      USE iom
603
604      INTEGER ::   inum   ! temporary logical unit
605      !!----------------------------------------------------------------------
606
607      IF(lwp) THEN
608         WRITE(numout,*)
609         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
610         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
611      ENDIF
612     
613      CALL iom_open( 'coordinates', inum )
614     
615      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
616      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
617      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
618      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
619     
620      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
621      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
622      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
623      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
624     
625      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
626      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
627      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
628      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
629     
630      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
631      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
632      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
633      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
634     
635      CALL iom_close( inum )
636     
637! need to be define for the extended grid south of -80S
638! some point are undefined but you need to have e1 and e2 .NE. 0
639      WHERE (e1t==0.0_wp)
640         e1t=1.0e2
641      END WHERE
642      WHERE (e1v==0.0_wp)
643         e1v=1.0e2
644      END WHERE
645      WHERE (e1u==0.0_wp)
646         e1u=1.0e2
647      END WHERE
648      WHERE (e1f==0.0_wp)
649         e1f=1.0e2
650      END WHERE
651      WHERE (e2t==0.0_wp)
652         e2t=1.0e2
653      END WHERE
654      WHERE (e2v==0.0_wp)
655         e2v=1.0e2
656      END WHERE
657      WHERE (e2u==0.0_wp)
658         e2u=1.0e2
659      END WHERE
660      WHERE (e2f==0.0_wp)
661         e2f=1.0e2
662      END WHERE
663       
664    END SUBROUTINE hgr_read
665   
666   !!======================================================================
667END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.