source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/tests/BENCH/MY_SRC/usrdef_zgr.F90 @ 11364

Last change on this file since 11364 was 11364, checked in by smasson, 2 years ago

dev_r10984_HPC-13 : update BENCH to be more robust, see #2285

File size: 12.2 KB
Line 
1MODULE usrdef_zgr
2   !!======================================================================
3   !!                       ***  MODULE  usrdef_zgr  ***
4   !!
5   !!                      ===  BENCH configuration  ===
6   !!
7   !! User defined : vertical coordinate system of a user configuration
8   !!======================================================================
9   !! History :  4.0  !
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_zgr   : user defined vertical coordinate system
14   !!      zgr_z      : reference 1D z-coordinate
15   !!      zgr_top_bot: ocean top and bottom level indices
16   !!      zgr_zco    : 3D verticl coordinate in pure z-coordinate case
17   !!---------------------------------------------------------------------
18   USE oce            ! ocean variables
19   USE dom_oce        ! ocean domain
20   USE phycst         ! physical constants
21   USE depth_e3       ! depth <=> e3
22   !
23   USE in_out_manager ! I/O manager
24   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
25   USE lib_mpp        ! distributed memory computing library
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   usr_def_zgr        ! called by domzgr.F90
31
32  !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS             
40
41   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
42      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
43      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
44      &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors
45      &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      -
46      &                    k_top  , k_bot    )                             ! top & bottom ocean level
47      !!---------------------------------------------------------------------
48      !!              ***  ROUTINE usr_def_zgr  ***
49      !!
50      !! ** Purpose :   User defined the vertical coordinates
51      !!
52      !!----------------------------------------------------------------------
53      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
54      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
55      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
56      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
57      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
58      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
59      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
60      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
61      !
62      INTEGER  ::   inum   ! local logical unit
63      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
64      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
65      !!----------------------------------------------------------------------
66      !
67      IF(lwp) WRITE(numout,*)
68      IF(lwp) WRITE(numout,*) 'usr_def_zgr : BENCH configuration (z-coordinate closed flat box ocean)'
69      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
70      !
71      !
72      ! type of vertical coordinate
73      ! ---------------------------
74      ld_zco    = .TRUE.         ! BENCH case:  z-coordinate without ocean cavities
75      ld_zps    = .FALSE.
76      ld_sco    = .FALSE.
77      ld_isfcav = .FALSE.
78      !
79      !
80      ! Build the vertical coordinate system
81      ! ------------------------------------
82      CALL zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
83      !
84      CALL zgr_msk_top_bot( k_top , k_bot )                 ! masked top and bottom ocean t-level indices
85      !
86      !                                                     ! z-coordinate (3D arrays) from the 1D z-coord.
87      CALL zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in  : 1D reference vertical coordinate
88         &          pdept   , pdepw   ,                     &   ! out : 3D t & w-points depth
89         &          pe3t    , pe3u    , pe3v   , pe3f   ,   &   !       vertical scale factors
90         &          pe3w    , pe3uw   , pe3vw             )     !           -      -      -
91      !
92   END SUBROUTINE usr_def_zgr
93
94
95   SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
96      !!----------------------------------------------------------------------
97      !!                   ***  ROUTINE zgr_z  ***
98      !!
99      !! ** Purpose :   set the 1D depth of model levels and the resulting
100      !!              vertical scale factors.
101      !!
102      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
103      !!       The depth of model levels is set from dep(k), an analytical function:
104      !!                   w-level: depw_1d  = dep(k)
105      !!                   t-level: dept_1d  = dep(k+0.5)
106      !!       The scale factors are the discrete derivative of the depth:
107      !!                   e3w_1d(jk) = dk[ dept_1d ]
108      !!                   e3t_1d(jk) = dk[ depw_1d ]
109      !!           with at top and bottom :
110      !!                   e3w_1d( 1 ) = 2 * ( dept_1d( 1 ) - depw_1d( 1 ) )
111      !!                   e3t_1d(jpk) = 2 * ( dept_1d(jpk) - depw_1d(jpk) )
112      !!       The depth are then re-computed from the sum of e3. This ensures
113      !!    that depths are identical when reading domain configuration file.
114      !!    Indeed, only e3. are saved in this file, depth are compute by a call
115      !!    to the e3_to_depth subroutine.
116      !!
117      !!       Here the Madec & Imbard (1996) function is used.
118      !!
119      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
120      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
121      !!
122      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
123      !!             Madec and Imbard, 1996, Clim. Dyn.
124      !!----------------------------------------------------------------------
125      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
126      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
127      !
128      INTEGER  ::   jk       ! dummy loop indices
129      REAL(wp) ::   zd       ! local scalar
130      !!----------------------------------------------------------------------
131      !
132      zd = 5000./FLOAT(jpkm1)
133      !
134      IF(lwp) THEN            ! Parameter print
135         WRITE(numout,*)
136         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates '
137         WRITE(numout,*) '    ~~~~~~~'
138         WRITE(numout,*) '       BENCH case : uniform vertical grid :'
139         WRITE(numout,*) '                     with thickness = ', zd
140      ENDIF
141
142      !
143      ! 1D Reference z-coordinate    (using Madec & Imbard 1996 function)
144      ! -------------------------
145      !
146      pdepw_1d(1) = 0._wp
147      pdept_1d(1) = 0.5_wp * zd
148      !
149      DO jk = 2, jpk          ! depth at T and W-points
150         pdepw_1d(jk) = pdepw_1d(jk-1) + zd 
151         pdept_1d(jk) = pdept_1d(jk-1) + zd 
152      END DO
153      !
154      !                       ! e3t and e3w from depth
155      CALL depth_to_e3( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d ) 
156      !
157      !                       ! recompute depths from SUM(e3)  <== needed
158      CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) 
159      !
160      IF(lwp) THEN                        ! control print
161         WRITE(numout,*)
162         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
163         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
164         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
165      ENDIF
166      !
167   END SUBROUTINE zgr_z
168
169
170   SUBROUTINE zgr_msk_top_bot( k_top , k_bot )
171      !!----------------------------------------------------------------------
172      !!                    ***  ROUTINE zgr_msk_top_bot  ***
173      !!
174      !! ** Purpose :   set the masked top and bottom ocean t-levels
175      !!
176      !! ** Method  :   BENCH case = closed flat box ocean without ocean cavities
177      !!                   k_top = 1     except along north, south, east and west boundaries
178      !!                   k_bot = jpk-1 except along north, south, east and west boundaries
179      !!
180      !! ** Action  : - k_top : first wet ocean level index
181      !!              - k_bot : last  wet ocean level index
182      !!----------------------------------------------------------------------
183      INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last wet ocean level
184      !
185      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace
186      REAL(wp)                     ::   zmaxlam, zscl
187      !!----------------------------------------------------------------------
188      !
189      IF(lwp) WRITE(numout,*)
190      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.'
191      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
192      IF(lwp) WRITE(numout,*) '       BENCH case : closed flat box ocean without ocean cavities'
193      !
194      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
195      !
196      IF( jperio == 3 .OR. jperio ==4 ) THEN   ! add a small island in the upper corners to avoid model instabilities...
197         z2d(mi0(       1):mi1(     3),mj0(jpjglo-2):mj1(jpjglo)) = 0.
198         z2d(mi0(jpiglo-2):mi1(jpiglo),mj0(jpjglo-2):mj1(jpjglo)) = 0.
199      ENDIF
200      !
201      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
202      !
203      k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
204      !
205      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere
206      !
207   END SUBROUTINE zgr_msk_top_bot
208   
209
210   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
211      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
212      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
213      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
214      !!----------------------------------------------------------------------
215      !!                  ***  ROUTINE zgr_zco  ***
216      !!
217      !! ** Purpose :   define the reference z-coordinate system
218      !!
219      !! ** Method  :   set 3D coord. arrays to reference 1D array
220      !!----------------------------------------------------------------------
221      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
222      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
223      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
224      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
225      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
226      !
227      INTEGER  ::   jk
228      !!----------------------------------------------------------------------
229      !
230      DO jk = 1, jpk
231         pdept(:,:,jk) = pdept_1d(jk)
232         pdepw(:,:,jk) = pdepw_1d(jk)
233         pe3t (:,:,jk) = pe3t_1d (jk)
234         pe3u (:,:,jk) = pe3t_1d (jk)
235         pe3v (:,:,jk) = pe3t_1d (jk)
236         pe3f (:,:,jk) = pe3t_1d (jk)
237         pe3w (:,:,jk) = pe3w_1d (jk)
238         pe3uw(:,:,jk) = pe3w_1d (jk)
239         pe3vw(:,:,jk) = pe3w_1d (jk)
240      END DO
241      !
242   END SUBROUTINE zgr_zco
243
244   !!======================================================================
245END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.