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 @ 6667

Last change on this file since 6667 was 6667, checked in by gm, 8 years ago

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