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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 31.2 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
28   IMPLICIT NONE
29   PRIVATE
30
31   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
32
33   PUBLIC   dom_hgr   ! called by domain.F90
34
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dom_hgr
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dom_hgr  ***
45      !!
46      !! ** Purpose :   Compute the geographical position (in degre) of the
47      !!      model grid-points,  the horizontal scale factors (in meters) and
48      !!      the Coriolis factor (in s-1).
49      !!
50      !! ** Method  :   The geographical position of the model grid-points is
51      !!      defined from analytical functions, fslam and fsphi, the deriva-
52      !!      tives of which gives the horizontal scale factors e1,e2.
53      !!      Defining two function fslam and fsphi and their derivatives in
54      !!      the two horizontal directions (fse1 and fse2), the model grid-
55      !!      point position and scale factors are given by:
56      !!         t-point:
57      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
58      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
59      !!         u-point:
60      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
61      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
62      !!         v-point:
63      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
64      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
65      !!            f-point:
66      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
67      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
68      !!      Where fse1 and fse2 are defined by:
69      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
70      !!                                     +          di(fsphi) **2 )(i,j)
71      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
72      !!                                     +          dj(fsphi) **2 )(i,j)
73      !!
74      !!        The coriolis factor is given at z-point by:
75      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
76      !!
77      !!        This routine is given as an example, it must be modified
78      !!      following the user s desiderata. nevertheless, the output as
79      !!      well as the way to compute the model grid-point position and
80      !!      horizontal scale factors must be respected in order to insure
81      !!      second order accuracy schemes.
82      !!
83      !! N.B. If the domain is periodic, verify that scale factors are also
84      !!      periodic, and the coriolis term again.
85      !!
86      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,
87      !!                u-, v- and f-points (in degre)
88      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
89      !!               u-, v-  and f-points (in degre)
90      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
91      !!      scale factors (in meters) at t-, u-, v-, and f-points.
92      !!        define ff: coriolis factor at f-point
93      !!
94      !! References :   Marti, Madec and Delecluse, 1992, JGR
95      !!                Madec, Imbard, 1996, Clim. Dyn.
96      !!----------------------------------------------------------------------
97      INTEGER  ::   ji, jj               ! dummy loop indices
98      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
99      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
100      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
101      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
102      REAL(wp) ::   zphi0, zbeta, znorme !
103      REAL(wp) ::   zarg, zf0, zminff, zmaxff
104      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
105      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
106      !!----------------------------------------------------------------------
107
108      IF(lwp) THEN
109         WRITE(numout,*)
110         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
111         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
112         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
113         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
114         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
115         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
116         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m 
117         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m 
118      ENDIF
119
120
121      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
122
123      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
124
125         IF(lwp) WRITE(numout,*)
126         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
127
128         CALL hgr_read           ! Defaultl option  :   NetCDF file
129
130         !                                                ! =====================
131         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
132            !                                             ! =====================
133            IF( nn_cla == 0 ) THEN
134               !
135               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
136               ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
137               IF(lwp) WRITE(numout,*)
138               IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km'
139               !
140               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
141               ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
142                                               e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3
143               IF(lwp) WRITE(numout,*)
144               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km'
145               IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km'
146            ENDIF
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 && defined key_eel_r6
317         IF( .NOT. Agrif_Root() ) THEN
318           glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
319           gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
320           ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
321           ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
322         ENDIF
323#endif         
324         DO jj = 1, jpj
325            DO ji = 1, jpi
326               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
327               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
328               glamv(ji,jj) = glamt(ji,jj)
329               glamf(ji,jj) = glamu(ji,jj)
330   
331               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
332               gphiu(ji,jj) = gphit(ji,jj)
333               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
334               gphif(ji,jj) = gphiv(ji,jj)
335            END DO
336         END DO
337
338         ! Horizontal scale factors (in meters)
339         !                              ======
340         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
341         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
342         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
343         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
344
345      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
346
347         IF(lwp) WRITE(numout,*)
348         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
349         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
350         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
351
352         !  Find index corresponding to the equator, given the grid spacing e1_deg
353         !  and the (approximate) southern latitude ppgphi0.
354         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
355         !  The formula should work even if the equator is outside the domain.
356         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
357         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
358         IF(  ppgphi0 > 0 )  ijeq = -ijeq
359
360         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
361
362         DO jj = 1, jpj
363            DO ji = 1, jpi
364               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
365               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
366               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
367               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
368         ! Longitude
369               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
370               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
371               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
372               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
373         ! Latitude
374               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
375               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
376               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
377               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
378         ! e1
379               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
380               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
381               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
382               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
383         ! e2
384               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
385               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
386               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
387               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
388            END DO
389         END DO
390
391      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
392
393         IF(lwp) WRITE(numout,*)
394         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
395         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
396
397         ! Position coordinates (in kilometers)
398         !                          ==========
399
400         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
401         zlam1 = -85
402         zphi1 = 29
403         ! resolution in meters
404         ze1 = 106000. / FLOAT(jp_cfg)           
405         ! benchmark: forced the resolution to be about 100 km
406         IF( nbench /= 0 )   ze1 = 106000.e0     
407         zsin_alpha = - SQRT( 2. ) / 2.
408         zcos_alpha =   SQRT( 2. ) / 2.
409         ze1deg = ze1 / (ra * rad)
410         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
411         !                                                          ! at the right jp_cfg resolution
412         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
413         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
414
415         IF( nprint==1 .AND. lwp )   THEN
416            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
417            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
418         ENDIF
419
420         DO jj = 1, jpj
421           DO ji = 1, jpi
422             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
423             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
424
425             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
426             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
427
428             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
429             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
430
431             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
432             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
433
434             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
435             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
436           END DO
437          END DO
438
439         ! Horizontal scale factors (in meters)
440         !                              ======
441         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
442         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
443         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
444         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
445
446      CASE DEFAULT
447         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
448         CALL ctl_stop( ctmp1 )
449
450      END SELECT
451     
452      ! T-cell surface
453      ! --------------
454      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
455
456
457      ! Control printing : Grid informations (if not restart)
458      ! ----------------
459
460      IF( lwp .AND. .NOT.ln_rstart ) THEN
461         WRITE(numout,*)
462         WRITE(numout,*) '          longitude and e1 scale factors'
463         WRITE(numout,*) '          ------------------------------'
464         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
465            glamv(ji,1), glamf(ji,1),   &
466            e1t(ji,1), e1u(ji,1),   &
467            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4689300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
469            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
470         
471         WRITE(numout,*)
472         WRITE(numout,*) '          latitude and e2 scale factors'
473         WRITE(numout,*) '          -----------------------------'
474         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
475            &                     gphiv(1,jj), gphif(1,jj),   &
476            &                     e2t  (1,jj), e2u  (1,jj),   &
477            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
478      ENDIF
479
480     
481      IF( nprint == 1 .AND. lwp ) THEN
482         WRITE(numout,*) '          e1u e2u '
483         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
484         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
485         WRITE(numout,*) '          e1v e2v  '
486         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
487         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
488         WRITE(numout,*) '          e1f e2f  '
489         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
490         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
491      ENDIF
492
493
494      ! ================= !
495      !  Coriolis factor  !
496      ! ================= !
497
498      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
499
500      CASE ( 0, 1, 4 )               ! mesh on the sphere
501
502         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
503
504      CASE ( 2 )                     ! f-plane at ppgphi0
505
506         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
507
508         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
509
510      CASE ( 3 )                     ! beta-plane
511
512         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
513         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
514         
515#if defined key_agrif && defined key_eel_r6
516         IF( .NOT. Agrif_Root() ) THEN
517           zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
518         ENDIF
519#endif         
520         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
521
522         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
523         
524         IF(lwp) THEN
525            WRITE(numout,*) 
526            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
527            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
528         ENDIF
529         IF( lk_mpp ) THEN
530            zminff=ff(nldi,nldj)
531            zmaxff=ff(nldi,nlej)
532            CALL mpp_min( zminff )   ! min over the global domain
533            CALL mpp_max( zmaxff )   ! max over the global domain
534            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
535         END IF
536
537      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
538
539         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
540         zphi0 = 15.e0                                                      ! latitude of the first row F-points
541         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
542
543         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
544
545         IF(lwp) THEN
546            WRITE(numout,*) 
547            WRITE(numout,*) '          Beta-plane and rotated domain : '
548            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
549         ENDIF
550
551         IF( lk_mpp ) THEN
552            zminff=ff(nldi,nldj)
553            zmaxff=ff(nldi,nlej)
554            CALL mpp_min( zminff )   ! min over the global domain
555            CALL mpp_max( zmaxff )   ! max over the global domain
556            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
557         END IF
558
559      END SELECT
560
561
562      ! Control of domain for symetrical condition
563      ! ------------------------------------------
564      ! The equator line must be the latitude coordinate axe
565
566      IF( nperio == 2 ) THEN
567         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
568         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
569      ENDIF
570
571   END SUBROUTINE dom_hgr
572
573
574   SUBROUTINE hgr_read
575      !!---------------------------------------------------------------------
576      !!              ***  ROUTINE hgr_read  ***
577      !!
578      !! ** Purpose :   Read a coordinate file in NetCDF format
579      !!
580      !! ** Method  :   The mesh file has been defined trough a analytical
581      !!      or semi-analytical method. It is read in a NetCDF file.
582      !!     
583      !!----------------------------------------------------------------------
584      USE iom
585
586      INTEGER ::   inum   ! temporary logical unit
587      !!----------------------------------------------------------------------
588
589      IF(lwp) THEN
590         WRITE(numout,*)
591         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
592         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
593      ENDIF
594     
595      CALL iom_open( 'coordinates', inum )
596     
597      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
598      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
599      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
600      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
601     
602      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
603      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
604      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
605      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
606     
607      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
608      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
609      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
610      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
611     
612      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
613      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
614      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
615      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
616     
617      CALL iom_close( inum )
618     
619    END SUBROUTINE hgr_read
620   
621   !!======================================================================
622END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.