source: branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 5605

Last change on this file since 5605 was 5605, checked in by mathiot, 6 years ago

Marine glacier termini: initial commit

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