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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 5067

Last change on this file since 5067 was 5067, checked in by clem, 9 years ago

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