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

source: branches/UKMO/dev_r5785_SSS_obsoper/trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 7740

Last change on this file since 7740 was 7740, checked in by mattmartin, 7 years ago

Removed svn keywords.

  • Property svn:keywords set to Id
File size: 32.8 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
474      ! Control printing : Grid informations (if not restart)
475      ! ----------------
476
477      IF( lwp .AND. .NOT.ln_rstart ) THEN
478         WRITE(numout,*)
479         WRITE(numout,*) '          longitude and e1 scale factors'
480         WRITE(numout,*) '          ------------------------------'
481         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
482            glamv(ji,1), glamf(ji,1),   &
483            e1t(ji,1), e1u(ji,1),   &
484            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4859300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
486            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
487         
488         WRITE(numout,*)
489         WRITE(numout,*) '          latitude and e2 scale factors'
490         WRITE(numout,*) '          -----------------------------'
491         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
492            &                     gphiv(1,jj), gphif(1,jj),   &
493            &                     e2t  (1,jj), e2u  (1,jj),   &
494            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
495      ENDIF
496
497     
498      IF( nprint == 1 .AND. lwp ) THEN
499         WRITE(numout,*) '          e1u e2u '
500         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
501         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
502         WRITE(numout,*) '          e1v e2v  '
503         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
504         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
505         WRITE(numout,*) '          e1f e2f  '
506         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
507         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
508      ENDIF
509
510
511      ! ================= !
512      !  Coriolis factor  !
513      ! ================= !
514
515      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
516
517      CASE ( 0, 1, 4 )               ! mesh on the sphere
518
519         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
520
521      CASE ( 2 )                     ! f-plane at ppgphi0
522
523         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
524
525         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
526
527      CASE ( 3 )                     ! beta-plane
528
529         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
530         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
531         
532#if defined key_agrif
533         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
534            IF( .NOT. Agrif_Root() ) THEN
535              zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
536            ENDIF
537         ENDIF
538#endif         
539         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
540
541         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
542         
543         IF(lwp) THEN
544            WRITE(numout,*) 
545            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
546            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
547         ENDIF
548         IF( lk_mpp ) THEN
549            zminff=ff(nldi,nldj)
550            zmaxff=ff(nldi,nlej)
551            CALL mpp_min( zminff )   ! min over the global domain
552            CALL mpp_max( zmaxff )   ! max over the global domain
553            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
554         END IF
555
556      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
557
558         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
559         zphi0 = 15.e0                                                      ! latitude of the first row F-points
560         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
561
562         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
563
564         IF(lwp) THEN
565            WRITE(numout,*) 
566            WRITE(numout,*) '          Beta-plane and rotated domain : '
567            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
568         ENDIF
569
570         IF( lk_mpp ) THEN
571            zminff=ff(nldi,nldj)
572            zmaxff=ff(nldi,nlej)
573            CALL mpp_min( zminff )   ! min over the global domain
574            CALL mpp_max( zmaxff )   ! max over the global domain
575            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
576         END IF
577
578      END SELECT
579
580
581      ! Control of domain for symetrical condition
582      ! ------------------------------------------
583      ! The equator line must be the latitude coordinate axe
584
585      IF( nperio == 2 ) THEN
586         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
587         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
588      ENDIF
589      !
590      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
591      !
592   END SUBROUTINE dom_hgr
593
594
595   SUBROUTINE hgr_read
596      !!---------------------------------------------------------------------
597      !!              ***  ROUTINE hgr_read  ***
598      !!
599      !! ** Purpose :   Read a coordinate file in NetCDF format
600      !!
601      !! ** Method  :   The mesh file has been defined trough a analytical
602      !!      or semi-analytical method. It is read in a NetCDF file.
603      !!     
604      !!----------------------------------------------------------------------
605      USE iom
606
607      INTEGER ::   inum   ! temporary logical unit
608      !!----------------------------------------------------------------------
609
610      IF(lwp) THEN
611         WRITE(numout,*)
612         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
613         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
614      ENDIF
615     
616      CALL iom_open( 'coordinates', inum )
617     
618      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
619      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
620      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
621      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
622     
623      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
624      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
625      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
626      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
627     
628      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
629      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
630      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
631      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
632     
633      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
634      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
635      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
636      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
637     
638      CALL iom_close( inum )
639     
640! need to be define for the extended grid south of -80S
641! some point are undefined but you need to have e1 and e2 .NE. 0
642      WHERE (e1t==0.0_wp)
643         e1t=1.0e2
644      END WHERE
645      WHERE (e1v==0.0_wp)
646         e1v=1.0e2
647      END WHERE
648      WHERE (e1u==0.0_wp)
649         e1u=1.0e2
650      END WHERE
651      WHERE (e1f==0.0_wp)
652         e1f=1.0e2
653      END WHERE
654      WHERE (e2t==0.0_wp)
655         e2t=1.0e2
656      END WHERE
657      WHERE (e2v==0.0_wp)
658         e2v=1.0e2
659      END WHERE
660      WHERE (e2u==0.0_wp)
661         e2u=1.0e2
662      END WHERE
663      WHERE (e2f==0.0_wp)
664         e2f=1.0e2
665      END WHERE
666       
667    END SUBROUTINE hgr_read
668   
669   !!======================================================================
670END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.