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

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications

  • Property svn:keywords set to Id
File size: 31.4 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   !! test - delete this line
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 && defined key_eel_r6
320         IF( .NOT. Agrif_Root() ) THEN
321           glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
322           gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
323           ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
324           ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()         
325         ENDIF
326#endif         
327         DO jj = 1, jpj
328            DO ji = 1, jpi
329               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
330               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
331               glamv(ji,jj) = glamt(ji,jj)
332               glamf(ji,jj) = glamu(ji,jj)
333   
334               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
335               gphiu(ji,jj) = gphit(ji,jj)
336               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
337               gphif(ji,jj) = gphiv(ji,jj)
338            END DO
339         END DO
340
341         ! Horizontal scale factors (in meters)
342         !                              ======
343         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
344         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
345         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
346         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
347
348      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
349
350         IF(lwp) WRITE(numout,*)
351         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
352         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
353         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
354
355         !  Find index corresponding to the equator, given the grid spacing e1_deg
356         !  and the (approximate) southern latitude ppgphi0.
357         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
358         !  The formula should work even if the equator is outside the domain.
359         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
360         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
361         IF(  ppgphi0 > 0 )  ijeq = -ijeq
362
363         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
364
365         DO jj = 1, jpj
366            DO ji = 1, jpi
367               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
368               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
369               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
370               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
371         ! Longitude
372               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
373               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
374               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
375               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
376         ! Latitude
377               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
378               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
379               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
380               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
381         ! e1
382               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
383               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
384               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
385               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
386         ! e2
387               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
388               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
389               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
390               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
391            END DO
392         END DO
393
394      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
395
396         IF(lwp) WRITE(numout,*)
397         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
398         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
399
400         ! Position coordinates (in kilometers)
401         !                          ==========
402
403         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
404         zlam1 = -85
405         zphi1 = 29
406         ! resolution in meters
407         ze1 = 106000. / FLOAT(jp_cfg)           
408         ! benchmark: forced the resolution to be about 100 km
409         IF( nbench /= 0 )   ze1 = 106000.e0     
410         zsin_alpha = - SQRT( 2. ) / 2.
411         zcos_alpha =   SQRT( 2. ) / 2.
412         ze1deg = ze1 / (ra * rad)
413         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
414         !                                                          ! at the right jp_cfg resolution
415         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
416         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
417
418         IF( nprint==1 .AND. lwp )   THEN
419            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
420            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
421         ENDIF
422
423         DO jj = 1, jpj
424           DO ji = 1, jpi
425             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
426             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
427
428             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
429             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
430
431             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
432             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
433
434             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
435             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
436
437             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
438             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
439           END DO
440          END DO
441
442         ! Horizontal scale factors (in meters)
443         !                              ======
444         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
445         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
446         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
447         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
448
449      CASE DEFAULT
450         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
451         CALL ctl_stop( ctmp1 )
452
453      END SELECT
454     
455      ! T-cell surface
456      ! --------------
457      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
458
459
460      ! Control printing : Grid informations (if not restart)
461      ! ----------------
462
463      IF( lwp .AND. .NOT.ln_rstart ) THEN
464         WRITE(numout,*)
465         WRITE(numout,*) '          longitude and e1 scale factors'
466         WRITE(numout,*) '          ------------------------------'
467         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
468            glamv(ji,1), glamf(ji,1),   &
469            e1t(ji,1), e1u(ji,1),   &
470            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
4719300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
472            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
473         
474         WRITE(numout,*)
475         WRITE(numout,*) '          latitude and e2 scale factors'
476         WRITE(numout,*) '          -----------------------------'
477         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
478            &                     gphiv(1,jj), gphif(1,jj),   &
479            &                     e2t  (1,jj), e2u  (1,jj),   &
480            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
481      ENDIF
482
483     
484      IF( nprint == 1 .AND. lwp ) THEN
485         WRITE(numout,*) '          e1u e2u '
486         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
487         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
488         WRITE(numout,*) '          e1v e2v  '
489         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
490         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
491         WRITE(numout,*) '          e1f e2f  '
492         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
493         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
494      ENDIF
495
496
497      ! ================= !
498      !  Coriolis factor  !
499      ! ================= !
500
501      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
502
503      CASE ( 0, 1, 4 )               ! mesh on the sphere
504
505         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
506
507      CASE ( 2 )                     ! f-plane at ppgphi0
508
509         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
510
511         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
512
513      CASE ( 3 )                     ! beta-plane
514
515         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
516         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
517         
518#if defined key_agrif && defined key_eel_r6
519         IF( .NOT. Agrif_Root() ) THEN
520           zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
521         ENDIF
522#endif         
523         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
524
525         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
526         
527         IF(lwp) THEN
528            WRITE(numout,*) 
529            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
530            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
531         ENDIF
532         IF( lk_mpp ) THEN
533            zminff=ff(nldi,nldj)
534            zmaxff=ff(nldi,nlej)
535            CALL mpp_min( zminff )   ! min over the global domain
536            CALL mpp_max( zmaxff )   ! max over the global domain
537            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
538         END IF
539
540      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
541
542         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
543         zphi0 = 15.e0                                                      ! latitude of the first row F-points
544         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
545
546         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
547
548         IF(lwp) THEN
549            WRITE(numout,*) 
550            WRITE(numout,*) '          Beta-plane and rotated domain : '
551            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
552         ENDIF
553
554         IF( lk_mpp ) THEN
555            zminff=ff(nldi,nldj)
556            zmaxff=ff(nldi,nlej)
557            CALL mpp_min( zminff )   ! min over the global domain
558            CALL mpp_max( zmaxff )   ! max over the global domain
559            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
560         END IF
561
562      END SELECT
563
564
565      ! Control of domain for symetrical condition
566      ! ------------------------------------------
567      ! The equator line must be the latitude coordinate axe
568
569      IF( nperio == 2 ) THEN
570         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
571         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
572      ENDIF
573      !
574      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
575      !
576   END SUBROUTINE dom_hgr
577
578
579   SUBROUTINE hgr_read
580      !!---------------------------------------------------------------------
581      !!              ***  ROUTINE hgr_read  ***
582      !!
583      !! ** Purpose :   Read a coordinate file in NetCDF format
584      !!
585      !! ** Method  :   The mesh file has been defined trough a analytical
586      !!      or semi-analytical method. It is read in a NetCDF file.
587      !!     
588      !!----------------------------------------------------------------------
589      USE iom
590
591      INTEGER ::   inum   ! temporary logical unit
592      !!----------------------------------------------------------------------
593
594      IF(lwp) THEN
595         WRITE(numout,*)
596         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
597         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
598      ENDIF
599     
600      CALL iom_open( 'coordinates', inum )
601     
602      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
603      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
604      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
605      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
606     
607      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
608      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
609      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
610      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
611     
612      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
613      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
614      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
615      CALL iom_get( inum, jpdom_data, 'e1f', e1f )
616     
617      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
618      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
619      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
620      CALL iom_get( inum, jpdom_data, 'e2f', e2f )
621     
622      CALL iom_close( inum )
623     
624    END SUBROUTINE hgr_read
625   
626   !!======================================================================
627END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.