source: NEMO/trunk/tests/SWG/MY_SRC/usrdef_zgr.F90 @ 15033

Last change on this file since 15033 was 15033, checked in by smasson, 6 months ago

trunk: suppress jpim1 et jpjm1, #2699

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