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 NEMO/trunk/src/OCE/DOM – NEMO

source: NEMO/trunk/src/OCE/DOM/domhgr.F90

Last change on this file was 15056, checked in by smasson, 3 years ago

trunk: correct [15052], #2701

  • Property svn:keywords set to Id
File size: 14.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   !!            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, G. Madec) new configuration interface: read or usrdef.F90
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dom_hgr       : initialize the horizontal mesh
23   !!   hgr_read      : read horizontal information in the domain configuration file
24   !!----------------------------------------------------------------------
25   USE dom_oce        ! ocean space and time domain
26   USE par_oce        ! ocean space and time domain
27   USE phycst         ! physical constants
28   USE usrdef_hgr     ! User defined routine
29   !
30   USE in_out_manager ! I/O manager
31   USE iom            ! I/O library
32   USE lib_mpp        ! MPP library
33   USE lbclnk         ! lateal boundary condition / mpp exchanges
34   USE timing         ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   dom_hgr   ! called by domain.F90
40
41   !!----------------------------------------------------------------------
42   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE dom_hgr
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dom_hgr  ***
51      !!
52      !! ** Purpose :   Read or compute the geographical position (in degrees) 
53      !!      of the model grid-points, the horizontal scale factors (in meters),
54      !!      the associated horizontal metrics, and the Coriolis factor (in s-1).
55      !!
56      !! ** Method  :   Controlled by ln_read_cfg logical
57      !!              =T : all needed arrays are read in mesh_mask.nc file
58      !!              =F : user-defined configuration, all needed arrays
59      !!                   are computed in usr-def_hgr subroutine
60      !!
61      !!                If Coriolis factor is neither read nor computed (iff=0)
62      !!              it is computed from gphit assuming that the mesh is
63      !!              defined on the sphere :
64      !!                   ff = 2.*omega*sin(gphif)      (in s-1)
65      !!   
66      !!                If u- & v-surfaces are neither read nor computed (ie1e2u_v=0)
67      !!              (i.e. no use of reduced scale factors in some straits)
68      !!              they are computed from e1u, e2u, e1v and e2v as:
69      !!                   e1e2u = e1u*e2u   and   e1e2v = e1v*e2v 
70      !!   
71      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)
72      !!              - define Coriolis parameter at f-point                   (in 1/s)
73      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
74      !!              - define associated horizontal metrics at t-, u-, v- and f-points
75      !!                (inverse of scale factors 1/e1 & 1/e2, surface e1*e2, ratios e1/e2 & e2/e1)
76      !!----------------------------------------------------------------------
77      INTEGER ::   ji, jj     ! dummy loop indices
78      INTEGER ::   ie1e2u_v   ! flag for u- & v-surfaces
79      INTEGER ::   iff        ! flag for Coriolis parameter
80      !!----------------------------------------------------------------------
81      !
82      IF( ln_timing )   CALL timing_start('dom_hgr')
83      !
84      IF(lwp) THEN
85         WRITE(numout,*)
86         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
87         WRITE(numout,*) '~~~~~~~   '
88         WRITE(numout,*) '   namcfg : read (=T) or user defined (=F) configuration    ln_read_cfg  = ', ln_read_cfg
89      ENDIF
90      !
91      IF( ln_read_cfg ) THEN        !==  read in mesh_mask.nc file  ==!
92         !
93         IF(lwp) WRITE(numout,*)
94         IF(lwp) WRITE(numout,*) '   ==>>>   read horizontal mesh in ', TRIM( cn_domcfg ), ' file'
95         !
96         CALL hgr_read   ( glamt , glamu , glamv , glamf ,   &    ! geographic position (required)
97            &              gphit , gphiu , gphiv , gphif ,   &    !     -        -
98            &              iff   , ff_f  , ff_t  ,           &    ! Coriolis parameter (if not on the sphere)
99            &              e1t   , e1u   , e1v   , e1f   ,   &    ! scale factors (required)
100            &              e2t   , e2u   , e2v   , e2f   ,   &    !    -     -        -
101            &              ie1e2u_v      , e1e2u , e1e2v     )    ! u- & v-surfaces (if gridsize reduction in some straits)
102         !
103      ELSE                          !==  User defined configuration  ==!
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) '          User defined horizontal mesh (usr_def_hgr)'
106         !
107         CALL usr_def_hgr( glamt , glamu , glamv , glamf ,   &    ! geographic position (required)
108            &              gphit , gphiu , gphiv , gphif ,   &    !
109            &              iff   , ff_f  , ff_t  ,           &    ! Coriolis parameter  (if domain not on the sphere)
110            &              e1t   , e1u   , e1v   , e1f   ,   &    ! scale factors       (required)
111            &              e2t   , e2u   , e2v   , e2f   ,   &    !
112            &              ie1e2u_v      , e1e2u , e1e2v     )    ! u- & v-surfaces (if gridsize reduction is used in strait(s))
113         !
114         ! make sure that periodicities are properly applied
115         CALL lbc_lnk( 'dom_hgr', glamt, 'T', 1._wp, glamu, 'U', 1._wp, glamv, 'V', 1._wp, glamf, 'F', 1._wp,   &
116            &                     gphit, 'T', 1._wp, gphiu, 'U', 1._wp, gphiv, 'V', 1._wp, gphif, 'F', 1._wp,   &
117            &                       e1t, 'T', 1._wp,   e1u, 'U', 1._wp,   e1v, 'V', 1._wp,   e1f, 'F', 1._wp,   &   
118            &                       e2t, 'T', 1._wp,   e2u, 'U', 1._wp,   e2v, 'V', 1._wp,   e2f, 'F', 1._wp,   &
119            &                     kfillmode = jpfillcopy )   ! do not put 0 over closed boundaries
120         !
121      ENDIF
122      !
123      !                             !==  Coriolis parameter  ==!   (if necessary)
124      !
125      IF( iff == 0 ) THEN                 ! Coriolis parameter has not been defined
126         IF(lwp) WRITE(numout,*) '          Coriolis parameter calculated on the sphere from gphif & gphit'
127         ff_f(:,:) = 2._wp * omega * SIN( rad * gphif(:,:) )     ! compute it on the sphere at f-point
128         ff_t(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )     !    -        -       -    at t-point
129      ELSE
130         IF( ln_read_cfg ) THEN
131            IF(lwp) WRITE(numout,*) '          Coriolis parameter have been read in ', TRIM( cn_domcfg ), ' file'
132         ELSE
133            CALL lbc_lnk( 'dom_hgr', ff_t, 'T', 1._wp, ff_f, 'F', 1._wp, kfillmode = jpfillcopy )   ! do not put 0 if closed
134            IF(lwp) WRITE(numout,*) '          Coriolis parameter have been set in usr_def_hgr routine'
135         ENDIF
136      ENDIF
137      !
138      !                             !==  associated horizontal metrics  ==!
139      !
140      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
141      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
142      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
143      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
144      !
145      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
146      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
147      IF( ie1e2u_v == 0 ) THEN               ! u- & v-surfaces have not been defined
148         IF(lwp) WRITE(numout,*) '          u- & v-surfaces calculated as e1 e2 product'
149         e1e2u (:,:) = e1u(:,:) * e2u(:,:)         ! compute them
150         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
151      ELSE
152         IF( ln_read_cfg ) THEN
153            IF(lwp) WRITE(numout,*) '          u- & v-surfaces have been read in ', TRIM( cn_domcfg ), ' file:'
154            IF(lwp) WRITE(numout,*) '                     grid size reduction in strait(s) is used'
155         ELSE
156            CALL lbc_lnk( 'dom_hgr', e1e2u, 'U', 1._wp, e1e2v, 'V', 1._wp, kfillmode = jpfillcopy )   ! do not put 0 if closed
157            IF(lwp) WRITE(numout,*) '          u- & v-surfaces have been have been set in usr_def_hgr routine'
158         ENDIF
159      ENDIF
160      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in any cases
161      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
162      !   
163      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
164      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
165      !
166      IF( ln_timing )   CALL timing_stop('dom_hgr')
167      !
168   END SUBROUTINE dom_hgr
169
170
171   SUBROUTINE hgr_read( plamt , plamu , plamv  , plamf  ,   &    ! gridpoints position (required)
172      &                 pphit , pphiu , pphiv  , pphif  ,   &     
173      &                 kff   , pff_f , pff_t  ,            &    ! Coriolis parameter  (if not on the sphere)
174      &                 pe1t  , pe1u  , pe1v   , pe1f   ,   &    ! scale factors       (required)
175      &                 pe2t  , pe2u  , pe2v   , pe2f   ,   &
176      &                 ke1e2u_v      , pe1e2u , pe1e2v     )    ! u- & v-surfaces (if gridsize reduction in some straits)
177      !!---------------------------------------------------------------------
178      !!              ***  ROUTINE hgr_read  ***
179      !!
180      !! ** Purpose :   Read a mesh_mask file in NetCDF format using IOM
181      !!
182      !!----------------------------------------------------------------------
183      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs
184      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs
185      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter read here, =0 otherwise
186      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point (if found in file)
187      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors
188      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors
189      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces read here, =0 otherwise
190      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v              ! u- & v-surfaces (if found in file)
191      !
192      INTEGER  ::   inum                  ! logical unit
193      !!----------------------------------------------------------------------
194      !
195      IF(lwp) THEN
196         WRITE(numout,*)
197         WRITE(numout,*) '   hgr_read : read the horizontal coordinates in mesh_mask'
198         WRITE(numout,*) '   ~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
199      ENDIF
200      !
201      CALL iom_open( cn_domcfg, inum )
202      !
203      CALL iom_get( inum, jpdom_global, 'glamt', plamt, cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
204      CALL iom_get( inum, jpdom_global, 'glamu', plamu, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
205      CALL iom_get( inum, jpdom_global, 'glamv', plamv, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
206      CALL iom_get( inum, jpdom_global, 'glamf', plamf, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
207      !
208      CALL iom_get( inum, jpdom_global, 'gphit', pphit, cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
209      CALL iom_get( inum, jpdom_global, 'gphiu', pphiu, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
210      CALL iom_get( inum, jpdom_global, 'gphiv', pphiv, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
211      CALL iom_get( inum, jpdom_global, 'gphif', pphif, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
212      !
213      CALL iom_get( inum, jpdom_global, 'e1t'  , pe1t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
214      CALL iom_get( inum, jpdom_global, 'e1u'  , pe1u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
215      CALL iom_get( inum, jpdom_global, 'e1v'  , pe1v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
216      CALL iom_get( inum, jpdom_global, 'e1f'  , pe1f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
217      !
218      CALL iom_get( inum, jpdom_global, 'e2t'  , pe2t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
219      CALL iom_get( inum, jpdom_global, 'e2u'  , pe2u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
220      CALL iom_get( inum, jpdom_global, 'e2v'  , pe2v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
221      CALL iom_get( inum, jpdom_global, 'e2f'  , pe2f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
222      !
223      IF(  iom_varid( inum, 'ff_f', ldstop = .FALSE. ) > 0  .AND.  &
224         & iom_varid( inum, 'ff_t', ldstop = .FALSE. ) > 0    ) THEN
225         IF(lwp) WRITE(numout,*) '           Coriolis factor at f- and t-points read in ', TRIM( cn_domcfg ), ' file'
226         CALL iom_get( inum, jpdom_global, 'ff_f', pff_f, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
227         CALL iom_get( inum, jpdom_global, 'ff_t', pff_t, cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
228         kff = 1
229      ELSE
230         kff = 0
231      ENDIF
232      !
233      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
234         IF(lwp) WRITE(numout,*) '           e1e2u & e1e2v read in ', TRIM( cn_domcfg ), ' file'
235         CALL iom_get( inum, jpdom_global, 'e1e2u', pe1e2u, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
236         CALL iom_get( inum, jpdom_global, 'e1e2v', pe1e2v, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
237         ke1e2u_v = 1
238      ELSE
239         ke1e2u_v = 0
240      ENDIF
241      !
242      CALL iom_close( inum )
243      !
244   END SUBROUTINE hgr_read
245   
246   !!======================================================================
247END MODULE domhgr
Note: See TracBrowser for help on using the repository browser.