source: NEMO/trunk/src/OCE/DOM/domzgr.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 9 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 18.8 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 information in the domain configuration file
26   !!   zgr_top_bot   : 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   USE closea         ! closed seas
32   USE depth_e3       ! depth <=> e3
33   USE wet_dry,   ONLY: ll_wd, ssh_ref  ! Wetting and drying
34   !
35   USE in_out_manager ! I/O manager
36   USE iom            ! I/O library
37   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
38   USE lib_mpp        ! distributed memory computing library
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   dom_zgr        ! called by dom_init.F90
44
45  !! * Substitutions
46#  include "do_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS       
53
54   SUBROUTINE dom_zgr( k_top, k_bot )
55      !!----------------------------------------------------------------------
56      !!                ***  ROUTINE dom_zgr  ***
57      !!                   
58      !! ** Purpose :   set the depth of model levels and the resulting
59      !!              vertical scale factors.
60      !!
61      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
62      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
63      !!              - vertical coordinate (gdep., e3.) depending on the
64      !!                coordinate chosen :
65      !!                   ln_zco=T   z-coordinate   
66      !!                   ln_zps=T   z-coordinate with partial steps
67      !!                   ln_zco=T   s-coordinate
68      !!
69      !! ** Action  :   define gdep., e3., mbathy and bathy
70      !!----------------------------------------------------------------------
71      INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices
72      !
73      INTEGER  ::   ji,jj,jk            ! dummy loop index
74      INTEGER  ::   ikt, ikb            ! top/bot index
75      INTEGER  ::   ioptio, ibat, ios   ! local integer
76      REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
77      !!----------------------------------------------------------------------
78      !
79      IF(lwp) THEN                     ! Control print
80         WRITE(numout,*)
81         WRITE(numout,*) 'dom_zgr : vertical coordinate'
82         WRITE(numout,*) '~~~~~~~'
83      ENDIF
84
85      IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time'
86
87
88      IF( ln_read_cfg ) THEN        !==  read in mesh_mask.nc file  ==!
89         IF(lwp) WRITE(numout,*)
90         IF(lwp) WRITE(numout,*) '   ==>>>   Read vertical mesh in ', TRIM( cn_domcfg ), ' file'
91         !
92         CALL zgr_read   ( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   & 
93            &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
94            &              gdept_0 , gdepw_0                    ,   &    ! gridpoints depth
95            &              e3t_0   , e3u_0   , e3v_0 , e3f_0    ,   &    ! vertical scale factors
96            &              e3w_0   , e3uw_0  , e3vw_0           ,   &    ! vertical scale factors
97            &              k_top   , k_bot            )                  ! 1st & last ocean level
98            !
99      ELSE                          !==  User defined configuration  ==!
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) '          User defined vertical mesh (usr_def_zgr)'
102         !
103         CALL usr_def_zgr( 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      ENDIF
111      !
112!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
113      ! Compute gde3w_0 (vertical sum of e3w)
114      gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
115      DO jk = 2, jpk
116         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
117      END DO
118      !
119      ! Any closed seas (defined by closea_mask > 0 in domain_cfg file) to be filled
120      ! in at runtime if ln_closea=.false.
121      IF( ln_closea ) THEN
122         IF ( ln_maskcs ) THEN
123            ! mask all the closed sea
124            CALL clo_msk( k_top, k_bot, mask_opnsea, 'mask_opensea' )
125         ELSE IF ( ln_mask_csundef ) THEN
126            ! defined closed sea are kept
127            ! mask all the undefined closed sea
128            CALL clo_msk( k_top, k_bot, mask_csundef, 'mask_csundef' )
129         END IF
130      END IF
131      !
132      IF(lwp) THEN                     ! Control print
133         WRITE(numout,*)
134         WRITE(numout,*) '   Type of vertical coordinate (read in ', TRIM( cn_domcfg ), ' file or set in userdef_nam) :'
135         WRITE(numout,*) '      z-coordinate - full steps      ln_zco    = ', ln_zco
136         WRITE(numout,*) '      z-coordinate - partial steps   ln_zps    = ', ln_zps
137         WRITE(numout,*) '      s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
138         WRITE(numout,*) '      ice shelf cavities             ln_isfcav = ', ln_isfcav
139      ENDIF
140
141      ioptio = 0                       ! Check Vertical coordinate options
142      IF( ln_zco      )   ioptio = ioptio + 1
143      IF( ln_zps      )   ioptio = ioptio + 1
144      IF( ln_sco      )   ioptio = ioptio + 1
145      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
146
147
148      !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
149      CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1
150      !
151      !                                ! ice shelf draft and bathymetry
152      DO_2D_11_11
153         ikt = mikt(ji,jj)
154         ikb = mbkt(ji,jj)
155         bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1)
156         risfdep(ji,jj) = gdepw_0(ji,jj,ikt  )
157      END_2D
158      !
159      !                                ! deepest/shallowest W level Above/Below ~10m
160!!gm BUG in s-coordinate this does not work!
161      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
162      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
163      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
164!!gm end bug
165      !
166      IF( nprint == 1 .AND. lwp )   THEN
167         WRITE(numout,*) ' MIN val k_top   ', MINVAL(   k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) )
168         WRITE(numout,*) ' MIN val k_bot   ', MINVAL(   k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) )
169         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
170            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
171         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  &
172            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  &
173            &                          ' uw', MINVAL(  e3uw_0(:,:,:) ), ' vw', MINVAL(  e3vw_0(:,:,:)),   &
174            &                          ' w ', MINVAL(   e3w_0(:,:,:) )
175
176         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
177            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
178         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  &
179            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  &
180            &                          ' uw', MAXVAL(  e3uw_0(:,:,:) ), ' vw', MAXVAL(  e3vw_0(:,:,:) ),  &
181            &                          ' w ', MAXVAL(   e3w_0(:,:,:) )
182      ENDIF
183      !
184   END SUBROUTINE dom_zgr
185
186
187   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate
188      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate
189      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth
190      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors
191      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      -
192      &                 k_top  , k_bot    )                            ! top & bottom ocean level
193      !!---------------------------------------------------------------------
194      !!              ***  ROUTINE zgr_read  ***
195      !!
196      !! ** Purpose :   Read the vertical information in the domain configuration file
197      !!
198      !!----------------------------------------------------------------------
199      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
200      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
201      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
202      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
203      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m]
204      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
205      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
206      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
207      !
208      INTEGER  ::   jk     ! dummy loop index
209      INTEGER  ::   inum   ! local logical unit
210      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
211      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
212      !!----------------------------------------------------------------------
213      !
214      IF(lwp) THEN
215         WRITE(numout,*)
216         WRITE(numout,*) '   zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
217         WRITE(numout,*) '   ~~~~~~~~'
218      ENDIF
219      !
220      CALL iom_open( cn_domcfg, inum )
221      !
222      !                          !* type of vertical coordinate
223      CALL iom_get( inum, 'ln_zco'   , z_zco )
224      CALL iom_get( inum, 'ln_zps'   , z_zps )
225      CALL iom_get( inum, 'ln_sco'   , z_sco )
226      IF( z_zco == 0._wp ) THEN   ;   ld_zco = .false.   ;   ELSE   ;   ld_zco = .true.   ;   ENDIF
227      IF( z_zps == 0._wp ) THEN   ;   ld_zps = .false.   ;   ELSE   ;   ld_zps = .true.   ;   ENDIF
228      IF( z_sco == 0._wp ) THEN   ;   ld_sco = .false.   ;   ELSE   ;   ld_sco = .true.   ;   ENDIF
229      !
230      !                          !* ocean cavities under iceshelves
231      CALL iom_get( inum, 'ln_isfcav', z_cav )
232      IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF
233      !
234      !                          !* vertical scale factors
235      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )                     ! 1D reference coordinate
236      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  )
237      !
238      CALL iom_get( inum, jpdom_data, 'e3t_0'  , pe3t  , lrowattr=ln_use_jattr )    ! 3D coordinate
239      CALL iom_get( inum, jpdom_data, 'e3u_0'  , pe3u  , lrowattr=ln_use_jattr )
240      CALL iom_get( inum, jpdom_data, 'e3v_0'  , pe3v  , lrowattr=ln_use_jattr )
241      CALL iom_get( inum, jpdom_data, 'e3f_0'  , pe3f  , lrowattr=ln_use_jattr )
242      CALL iom_get( inum, jpdom_data, 'e3w_0'  , pe3w  , lrowattr=ln_use_jattr )
243      CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr )
244      CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr )
245      !
246      !                          !* depths
247      !                                   !- old depth definition (obsolescent feature)
248      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
249         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
250         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
251         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
252         CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 
253            &           '           depths at t- and w-points read in the domain configuration file')
254         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )   
255         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
256         CALL iom_get( inum, jpdom_data   , 'gdept_0' , pdept , lrowattr=ln_use_jattr )
257         CALL iom_get( inum, jpdom_data   , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr )
258         !
259      ELSE                                !- depths computed from e3. scale factors
260         CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth
261         CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths
262         IF(lwp) THEN
263            WRITE(numout,*)
264            WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
265            WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
266            WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
267         ENDIF
268      ENDIF
269      !
270      !                          !* ocean top and bottom level
271      CALL iom_get( inum, jpdom_data, 'top_level'    , z2d  , lrowattr=ln_use_jattr )   ! 1st wet T-points (ISF)
272      k_top(:,:) = NINT( z2d(:,:) )
273      CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d  , lrowattr=ln_use_jattr )   ! last wet T-points
274      k_bot(:,:) = NINT( z2d(:,:) )
275      !
276      ! reference depth for negative bathy (wetting and drying only)
277      IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   )
278      !
279      CALL iom_close( inum )
280      !
281   END SUBROUTINE zgr_read
282
283
284   SUBROUTINE zgr_top_bot( k_top, k_bot )
285      !!----------------------------------------------------------------------
286      !!                    ***  ROUTINE zgr_top_bot  ***
287      !!
288      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
289      !!
290      !! ** Method  :   computes from k_top and k_bot with a minimum value of 1 over land
291      !!
292      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest
293      !!                                     ocean level at t-, u- & v-points
294      !!                                     (min value = 1)
295      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
296      !!                                     ocean level at t-, u- & v-points
297      !!                                     (min value = 1 over land)
298      !!----------------------------------------------------------------------
299      INTEGER , DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! top & bottom ocean level indices
300      !
301      INTEGER ::   ji, jj   ! dummy loop indices
302      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace
303      !!----------------------------------------------------------------------
304      !
305      IF(lwp) WRITE(numout,*)
306      IF(lwp) WRITE(numout,*) '    zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels '
307      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
308      !
309      mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land)
310      !
311      mbkt(:,:) = MAX( k_bot(:,:) , 1 )    ! bottom ocean k-index of T-level (=1 over land)
312 
313      !                                    ! N.B.  top     k-index of W-level = mikt
314      !                                    !       bottom  k-index of W-level = mbkt+1
315      DO_2D_10_10
316         miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
317         mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
318         mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
319         !
320         mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
321         mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
322      END_2D
323      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
324      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
325      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mikv(:,:) = MAX( NINT( zk(:,:) ), 1 )
326      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'F', 1. )   ;   mikf(:,:) = MAX( NINT( zk(:,:) ), 1 )
327      !
328      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   mbku(:,:) = MAX( NINT( zk(:,:) ), 1 )
329      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 )
330      !
331   END SUBROUTINE zgr_top_bot
332
333   !!======================================================================
334END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.