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/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 8865

Last change on this file since 8865 was 8865, checked in by deazer, 6 years ago

Moving Changes from CS15mini config into NEMO main src
Updating TEST configs to run wit this version of the code
all sette tests pass at this revision other than AGRIF
Includes updates to dynnxt and tranxt required for 3D rives in WAD case to be conservative.

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