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.
domzgr.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/domzgr.F90 @ 6667

Last change on this file since 6667 was 6667, checked in by gm, 9 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: 17.4 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean domain : definition of the vertical coordinate system
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
20   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dom_zgr       : read or set the ocean vertical coordinate system
25   !!   zgr_read      : read the vertical domain coordinate and mask in domain_cfg file
26   !!   zgr_tb_level  : ocean top and bottom level for t-, u, and v-points with 1 as minimum value
27   !!---------------------------------------------------------------------
28   USE oce            ! ocean variables
29   USE dom_oce        ! ocean domain
30   USE usrdef_zgr     ! user defined vertical coordinate system
31   !
32   USE in_out_manager ! I/O manager
33   USE iom            ! I/O library
34   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
35   USE lib_mpp        ! distributed memory computing library
36   USE wrk_nemo       ! Memory allocation
37   USE timing         ! Timing
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   dom_zgr        ! called by dom_init.F90
43
44  !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS       
52
53   SUBROUTINE dom_zgr( k_top, k_bot )
54      !!----------------------------------------------------------------------
55      !!                ***  ROUTINE dom_zgr  ***
56      !!                   
57      !! ** Purpose :   set the depth of model levels and the resulting
58      !!              vertical scale factors.
59      !!
60      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
61      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
62      !!              - vertical coordinate (gdep., e3.) depending on the
63      !!                coordinate chosen :
64      !!                   ln_zco=T   z-coordinate   
65      !!                   ln_zps=T   z-coordinate with partial steps
66      !!                   ln_zco=T   s-coordinate
67      !!
68      !! ** Action  :   define gdep., e3., mbathy and bathy
69      !!----------------------------------------------------------------------
70      INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices
71      !
72      INTEGER  ::   jk                  ! dummy loop index
73      INTEGER  ::   ioptio, ibat, ios   ! local integer
74      REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
75      !!
76!      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
77      !!----------------------------------------------------------------------
78      !
79      IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
80      !
81!      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
82!      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
83!901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
84!
85!      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
86!      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
87!902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
88!      IF(lwm) WRITE ( numond, namzgr )
89
90      IF(lwp) THEN                     ! Control print
91         WRITE(numout,*)
92         WRITE(numout,*) 'dom_zgr : vertical coordinate'
93         WRITE(numout,*) '~~~~~~~'
94      ENDIF
95
96      IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time'
97
98
99      IF( ln_read_cfg ) THEN        !==  read in mesh_mask.nc file  ==!
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) '          Read vertical mesh in "domain_cfg" file'
102         !
103         CALL zgr_read   ( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   & 
104            &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
105            &              gdept_0 , gdepw_0                    ,   &    ! gridpoints depth
106            &              e3t_0   , e3u_0   , e3v_0 , e3f_0    ,   &    ! vertical scale factors
107            &              e3w_0   , e3uw_0  , e3vw_0           ,   &    ! vertical scale factors
108            &              k_top   , k_bot            )                  ! 1st & last ocean level
109         !
110      ELSE                          !==  User defined configuration  ==!
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) '          User defined horizontal mesh (usr_def_hgr)'
113         !
114         CALL usr_def_zgr( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   & 
115            &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
116            &              gdept_0 , gdepw_0                    ,   &    ! gridpoints depth
117            &              e3t_0   , e3u_0   , e3v_0 , e3f_0    ,   &    ! vertical scale factors
118            &              e3w_0   , e3uw_0  , e3vw_0           ,   &    ! vertical scale factors
119            &              k_top   , k_bot            )                  ! 1st & last ocean level
120         !
121      ENDIF
122      !
123!!gm to be remove when changing the definition of e3 scale factors so that gde3w disappears
124      ! Compute gde3w_0 (vertical sum of e3w)
125      gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
126      DO jk = 2, jpk
127         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
128      END DO
129      !
130      IF(lwp) THEN                     ! Control print
131         WRITE(numout,*)
132         WRITE(numout,*) '   Read in domain_cfg.nc or user defined type of vertical coordinate:'
133         WRITE(numout,*) '      z-coordinate - full steps      ln_zco    = ', ln_zco
134         WRITE(numout,*) '      z-coordinate - partial steps   ln_zps    = ', ln_zps
135         WRITE(numout,*) '      s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
136         WRITE(numout,*) '      ice shelf cavities             ln_isfcav = ', ln_isfcav
137      ENDIF
138
139      ioptio = 0                       ! Check Vertical coordinate options
140      IF( ln_zco      )   ioptio = ioptio + 1
141      IF( ln_zps      )   ioptio = ioptio + 1
142      IF( ln_sco      )   ioptio = ioptio + 1
143      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
144
145
146      !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
147      CALL zgr_tb_level( k_top, k_bot )      ! with a minimum value set to 1
148     
149
150      !                                ! deepest/shallowest W level Above/Below ~10m
151!!gm BUG in s-coordinate this does not work!
152      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
153      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
154      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
155!!gm end bug
156      !
157      !
158      IF( nprint == 1 .AND. lwp )   THEN
159         WRITE(numout,*) ' MIN val k_top   ', MINVAL(   k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) )
160         WRITE(numout,*) ' MIN val k_bot   ', MINVAL(   k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) )
161         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
162            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
163         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  &
164            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  &
165            &                          ' uw', MINVAL(  e3uw_0(:,:,:) ), ' vw', MINVAL(  e3vw_0(:,:,:)),   &
166            &                          ' w ', MINVAL(   e3w_0(:,:,:) )
167
168         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
169            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
170         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  &
171            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  &
172            &                          ' uw', MAXVAL(  e3uw_0(:,:,:) ), ' vw', MAXVAL(  e3vw_0(:,:,:) ),  &
173            &                          ' w ', MAXVAL(   e3w_0(:,:,:) )
174      ENDIF
175      !
176      IF( nn_timing == 1 )  CALL timing_stop('dom_zgr')
177      !
178   END SUBROUTINE dom_zgr
179
180
181   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate
182      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate
183      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth
184      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors
185      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      -
186      &                 k_top  , k_bot    )                            ! top & bottom ocean level
187      !!---------------------------------------------------------------------
188      !!              ***  ROUTINE zgr_read  ***
189      !!
190      !! ** Purpose :   Read the vertical information in a domain_cfg.nc file
191      !!
192      !!----------------------------------------------------------------------
193      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
194      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
195      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
196      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
197      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m]
198      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
199      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
200      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
201      !
202      INTEGER  ::   inum   ! local logical unit
203      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
204      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
205      !!----------------------------------------------------------------------
206      !
207      IF(lwp) THEN
208         WRITE(numout,*)
209         WRITE(numout,*) 'hgr_read : read the vertical coordinates in "domain_cfg.nc" file'
210         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo
211      ENDIF
212      !
213      CALL iom_open( 'domain_cfg', inum )
214      !
215      !                                   ! type of vertical coordinate
216      CALL iom_get( inum, 'ln_zco'   , z_zco )
217      CALL iom_get( inum, 'ln_zps'   , z_zps )
218      CALL iom_get( inum, 'ln_sco'   , z_sco )
219      IF( z_zco == 0._wp ) THEN   ;   ld_zco = .false.   ;   ELSE   ;   ld_zco = .true.   ;   ENDIF
220      IF( z_zps == 0._wp ) THEN   ;   ld_zps = .false.   ;   ELSE   ;   ld_zps = .true.   ;   ENDIF
221      IF( z_sco == 0._wp ) THEN   ;   ld_sco = .false.   ;   ELSE   ;   ld_sco = .true.   ;   ENDIF
222      !
223      !                                   ! ocean cavities under iceshelves
224      CALL iom_get( inum, 'ln_isfcav', z_cav )
225      IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF
226      !
227      !                                   ! reference 1D-coordinate
228      CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )   
229      CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
230      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )
231      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  )
232      !
233      !                                   ! depth
234      CALL iom_get( inum, jpdom_data, 'gdept_0', pdept  )
235      CALL iom_get( inum, jpdom_data, 'gdepw_0', pdepw  )
236!      CALL iom_get( inum, jpdom_data, 'gde3w_0', pde3w )
237      !
238      !                                   ! vertical scale factors
239      CALL iom_get( inum, jpdom_data, 'e3t_0'  , pe3t  )
240      CALL iom_get( inum, jpdom_data, 'e3u_0'  , pe3u  )
241      CALL iom_get( inum, jpdom_data, 'e3v_0'  , pe3v  )
242      CALL iom_get( inum, jpdom_data, 'e3f_0'  , pe3f  )
243      CALL iom_get( inum, jpdom_data, 'e3w_0'  , pe3w  )
244      CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw )
245      CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw )
246
247      !                                   ! ocean top and bottom level
248      CALL iom_get( inum, jpdom_data, 'bottom level' , z2d    )  ! nb of ocean T-points
249      k_bot(:,:) = INT( z2d(:,:) )
250      CALL iom_get( inum, jpdom_data, 'top    level' , z2d     )   ! nb of ocean T-points (ISF)
251      k_top(:,:) = INT( z2d(:,:) )
252      !
253      CALL iom_close( inum )
254      !
255   END SUBROUTINE zgr_read
256
257
258   SUBROUTINE zgr_tb_level( k_top, k_bot )
259      !!----------------------------------------------------------------------
260      !!                    ***  ROUTINE zgr_tb_level  ***
261      !!
262      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
263      !!
264      !! ** Method  :   computes from k_top and k_bot with a minimum value of 1 over land
265      !!
266      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
267      !!                                     ocean level at t-, u- & v-points
268      !!                                     (min value = 1)
269      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
270      !!                                     ocean level at t-, u- & v-points
271      !!                                     (min value = 1 over land)
272      !!----------------------------------------------------------------------
273      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot   ! top & bottom ocean level indices
274      !
275      INTEGER ::   ji, jj   ! dummy loop indices
276      REAL(wp), POINTER, DIMENSION(:,:) ::  zk
277      !!----------------------------------------------------------------------
278      !
279      IF( nn_timing == 1 )  CALL timing_start('zgr_bot_level')
280      !
281      CALL wrk_alloc( jpi,jpj,   zk )
282      !
283      IF(lwp) WRITE(numout,*)
284      IF(lwp) WRITE(numout,*) '    zgr_tb_level : ocean top and bottom k-index of T-, U-, V- and W-levels '
285      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
286      !
287      mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land)
288      !
289      mbkt(:,:) = MAX( k_bot(:,:) , 1 )    ! bottom ocean k-index of T-level (=1 over land)
290 
291      !                                    ! N.B.  top     k-index of W-level = mikt
292      !                                    !       bottom  k-index of W-level = mbkt+1
293      DO jj = 1, jpjm1
294         DO ji = 1, jpim1
295            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
296            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
297            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
298            !
299            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
300            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
301         END DO
302      END DO
303      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
304      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( zk, 'U', 1. )   ;   miku(:,:) = MAX( INT( zk(:,:) ), 1 )
305      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( zk, 'V', 1. )   ;   mikv(:,:) = MAX( INT( zk(:,:) ), 1 )
306      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( zk, 'F', 1. )   ;   mikf(:,:) = MAX( INT( zk(:,:) ), 1 )
307      !
308      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( zk, 'U', 1. )   ;   mbku(:,:) = MAX( INT( zk(:,:) ), 1 )
309      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( zk, 'V', 1. )   ;   mbkv(:,:) = MAX( INT( zk(:,:) ), 1 )
310      !
311      CALL wrk_dealloc( jpi,jpj,   zk )
312      !
313      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level')
314      !
315   END SUBROUTINE zgr_tb_level
316
317   !!======================================================================
318END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.