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

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 6593

Last change on this file since 6593 was 6593, checked in by flavoni, 8 years ago

#1692: simplified domhgr and usrdef modules
  • Property svn:keywords set to Id
File size: 10.7 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   !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse
17   !!                                       add optional read of e1e2u & e1e2v
18   !!             -   ! 2016-04  (S. Flavoni) change configuration's interface:
19   !!                            read file or CALL usr_def module to compute
20   !!                            horizontal grid (example given for GYRE)
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dom_hgr       : initialize the horizontal mesh
25   !!   hgr_read      : read "coordinate" NetCDF file
26   !!----------------------------------------------------------------------
27   USE dom_oce        ! ocean space and time domain
28   USE par_oce        ! ocean space and time domain
29   USE phycst         ! physical constants
30   USE usrdef         ! User defined routine
31!   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
32   !
33   USE in_out_manager ! I/O manager
34   USE lib_mpp        ! MPP library
35   USE timing         ! Timing
36
37   IMPLICIT NONE
38   PRIVATE
39
40   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in namelist !SF modify some routines????
41
42   PUBLIC   dom_hgr   ! called by domain.F90
43
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE dom_hgr
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE dom_hgr  ***
54      !!
55      !! ** Purpose :   Compute the geographical position (in degre) of the
56      !!      model grid-points,  the horizontal scale factors (in meters) and
57      !!      the Coriolis factor (in s-1).
58      !!
59      !! ** Method  :   Controlled by ln_read_cfg logical
60      !!              =T : all needed arrays are read in mesh_mask.nc file
61      !!              =F : user-defined configuration, all needed arrays
62      !!                   are computed in usr-def_hgr subroutine
63      !!
64      !!                If Coriolis factor is neither read nor computed (iff=0)
65      !!              it is computed from gphit assuming that the mesh is
66      !!              defined on the sphere :
67      !!                   ff = 2.*omega*sin(gphif)      (in s-1)
68      !!   
69      !!                If u- & v-surfaces are neither read nor computed (ie1e2u_v=0)
70      !!              (i.e. no use of reduced scale factors in some straits)
71      !!              they are computed from e1u, e2u, e1v and e2v as:
72      !!                   e1e2u = e1u*e2u   and   e1e2v = e1v*e2v 
73      !!   
74      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)
75      !!              - define Coriolis parameter at f-point                   (in 1/s)
76      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
77      !!              - define associated horizontal metrics at t-, u-, v- and f-points
78      !!                (inverse of scale factors 1/e1 & 1/e2, surface e1*e2, ratios e1/e2 & e2/e1)
79      !!----------------------------------------------------------------------
80      INTEGER  ::   ji, jj                   ! dummy loop indices
81      !INTEGER  ::   ii0, ii1, ij0, ij1, iff  ! temporary integers
82      !INTEGER  ::   ijeq                     ! index of equator T point (used in case 4)
83      !REAL(wp) ::   zti, zui, zvi, zfi       ! local scalars
84      !REAL(wp) ::   ztj, zuj, zvj, zfj       !   -      -
85      !REAL(wp) ::   zphi0, zbeta, znorme     !
86      !REAL(wp) ::   zarg, zf0, zminff, zmaxff
87      !REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
88      !REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
89      !INTEGER  ::   isrow                    ! index for ORCA1 starting row
90      INTEGER  ::   ie1e2u_v                 ! flag for u- & v-surfaces
91      INTEGER  ::   iff                      ! flag for Coriolis parameter
92      !!----------------------------------------------------------------------
93      !
94      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
95      !
96      IF(lwp) THEN
97         WRITE(numout,*)
98         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
99         WRITE(numout,*) '~~~~~~~   '
100      ENDIF
101      !
102      !
103      IF( ln_read_cfg ) THEN        !==  read in mesh_mask.nc file  ==!
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) '          read horizontal mesh in "mesh_mask" file'
106         !
107         CALL hgr_read( ie1e2u_v, iff )     
108         !
109      ELSE   
110         IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module'
111         CALL usr_def_hgr( nbench , jp_cfg, iff   , ff    ,    &
112         &                 glamt  , glamu , glamv , glamf ,    &
113         &                 gphit  , gphiu , gphiv , gphif ,    &
114         &                 e1t    , e1u   , e1v   , e1f   ,    &
115         &                 e2t    , e2u   , e2v   , e2f   ,    &
116         &                 ie1e2u_v  )
117         !
118      ENDIF
119      !
120      ! associated horizontal metrics
121      ! -----------------------------
122      !
123      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
124      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
125      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
126      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
127      !
128      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
129      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
130      IF( ie1e2u_v == 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
131         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
132         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
133      ENDIF
134      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
135      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
136      !   
137      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
138      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
139     
140      ! ================= !
141      !  Coriolis factor  !
142      ! ================= !
143      IF ( iff == 0 ) THEN       ! Coriolis parameter has not been defined: compute it on the sphere
144         !
145         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
146         !
147      ENDIF
148      !
149      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
150      !
151   END SUBROUTINE dom_hgr
152
153
154   SUBROUTINE hgr_read( ke1e2u_v, kff )
155      !!---------------------------------------------------------------------
156      !!              ***  ROUTINE hgr_read  ***
157      !!
158      !! ** Purpose :   Read a mesh_mask file in NetCDF format using IOM
159      !!
160      !!----------------------------------------------------------------------
161      USE iom
162      !!
163      INTEGER, INTENT( out ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in mesh_mask file (=1) or not (=0)
164      INTEGER, INTENT( out ) ::   kff        ! fag: kff read in mesh_mask file (=1) or not (=0)
165      !
166      INTEGER ::   inum   ! temporary logical unit
167      !!----------------------------------------------------------------------
168      !
169      IF(lwp) THEN
170         WRITE(numout,*)
171         WRITE(numout,*) 'hgr_read : read the horizontal coordinates in mesh_mask'
172         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
173      ENDIF
174      !
175      CALL iom_open( 'mesh_mask', inum )
176      !
177      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
178      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
179      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
180      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
181      !
182      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
183      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
184      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
185      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
186      !
187      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
188      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
189      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
190      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
191      !
192      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
193      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
194      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
195      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
196      !SF add read coriolis in mesh_mask file
197      CALL iom_get( inum, jpdom_data, 'ff'   , ff   , lrowattr=ln_use_jattr )
198      !
199      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
200         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in mesh_mask file'
201         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
202         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
203         ke1e2u_v = 1
204      ELSE
205         ke1e2u_v = 0
206      ENDIF
207      !
208      IF( iom_varid( inum, 'kff', ldstop = .FALSE. ) > 0 ) THEN
209         IF(lwp) WRITE(numout,*) 'hgr_read : ff read in mesh_mask file'
210         CALL iom_get( inum, jpdom_data, 'ff'  , ff  , lrowattr=ln_use_jattr )
211         kff = 1
212      ELSE
213         kff = 0
214      ENDIF
215      !
216      CALL iom_close( inum )
217     
218    END SUBROUTINE hgr_read
219   
220   !!======================================================================
221END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.