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.
usrdef_zgr.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/USR – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_zgr.F90 @ 7702

Last change on this file since 7702 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
File size: 13.4 KB
Line 
1MODULE usrdef_zgr
2   !!======================================================================
3   !!                       ***  MODULE  usrdef_zgr  ***
4   !!
5   !!                       ===  GYRE configuration  ===
6   !!
7   !! User defined : vertical coordinate system of a user configuration
8   !!======================================================================
9   !! History :  4.0  ! 2016-06  (G. Madec)  Original code
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 depth_e3       ! depth <=> e3
21   !
22   USE in_out_manager ! I/O manager
23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp        ! distributed memory computing library
25   USE wrk_nemo       ! Memory allocation
26   USE timing         ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   usr_def_zgr        ! called by domzgr.F90
32
33  !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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 : GYRE 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    = .TRUE.         ! GYRE case:  z-coordinate without ocean cavities
76      ld_zps    = .FALSE.
77      ld_sco    = .FALSE.
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      REAL(wp) ::   zt, zw   ! local scalars
131      REAL(wp) ::   zsur, za0, za1, zkth, zacr   ! Values for the Madec & Imbard (1996) function 
132      !!----------------------------------------------------------------------
133      !
134      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
135      !
136      ! Set parameters of z(k) function
137      ! -------------------------------
138      zsur = -2033.194295283385_wp       
139      za0  =   155.8325369664153_wp 
140      za1  =   146.3615918601890_wp
141      zkth =    17.28520372419791_wp   
142      zacr =     5.0_wp       
143      !
144      IF(lwp) THEN            ! Parameter print
145         WRITE(numout,*)
146         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates '
147         WRITE(numout,*) '    ~~~~~~~'
148         WRITE(numout,*) '       GYRE case : MI96 function with the following coefficients :'
149         WRITE(numout,*) '                 zsur = ', zsur
150         WRITE(numout,*) '                 za0  = ', za0
151         WRITE(numout,*) '                 za1  = ', za1
152         WRITE(numout,*) '                 zkth = ', zkth
153         WRITE(numout,*) '                 zacr = ', zacr
154      ENDIF
155
156      !
157      ! 1D Reference z-coordinate    (using Madec & Imbard 1996 function)
158      ! -------------------------
159      !
160      DO jk = 1, jpk          ! depth at T and W-points
161         zw = REAL( jk , wp )
162         zt = REAL( jk , wp ) + 0.5_wp
163         pdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr *  LOG( COSH( (zw-zkth) / zacr ) )  )
164         pdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr *  LOG( COSH( (zt-zkth) / zacr ) )  )
165      END DO
166      !
167      !                       ! e3t and e3w from depth
168      CALL depth_to_e3( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d ) 
169      !
170      !                       ! recompute depths from SUM(e3)  <== needed
171      CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) 
172      !
173      IF(lwp) THEN                        ! control print
174         WRITE(numout,*)
175         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
176         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
177         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
178      ENDIF
179      !
180      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
181      !
182   END SUBROUTINE zgr_z
183
184
185   SUBROUTINE zgr_msk_top_bot( k_top , k_bot )
186      !!----------------------------------------------------------------------
187      !!                    ***  ROUTINE zgr_msk_top_bot  ***
188      !!
189      !! ** Purpose :   set the masked top and bottom ocean t-levels
190      !!
191      !! ** Method  :   GYRE case = closed flat box ocean without ocean cavities
192      !!                   k_top = 1     except along north, south, east and west boundaries
193      !!                   k_bot = jpk-1 except along north, south, east and west boundaries
194      !!
195      !! ** Action  : - k_top : first wet ocean level index
196      !!              - k_bot : last  wet ocean level index
197      !!----------------------------------------------------------------------
198      INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last wet ocean level
199      !
200      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace
201
202      INTEGER  ::   ji, jj
203      !!----------------------------------------------------------------------
204      !
205      IF(lwp) WRITE(numout,*)
206      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.'
207      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
208      IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities'
209      !
210!$OMP PARALLEL DO schedule(static) private(jj, ji)
211      DO jj = 1, jpj
212         DO ji = 1, jpi
213            z2d(ji,jj) = REAL( jpkm1 , wp )          ! flat bottom
214         END DO
215      END DO
216      !
217      CALL lbc_lnk( z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
218      !
219!$OMP PARALLEL DO schedule(static) private(jj, ji)
220      DO jj = 1, jpj
221         DO ji = 1, jpi
222            k_bot(ji,jj) = INT( z2d(ji,jj) )           ! =jpkm1 over the ocean point, =0 elsewhere
223            !
224            k_top(ji,jj) = MIN( 1 , k_bot(ji,jj) )     ! = 1    over the ocean point, =0 elsewhere
225         END DO
226      END DO
227      !
228   END SUBROUTINE zgr_msk_top_bot
229   
230
231   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
232      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
233      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
234      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
235      !!----------------------------------------------------------------------
236      !!                  ***  ROUTINE zgr_zco  ***
237      !!
238      !! ** Purpose :   define the reference z-coordinate system
239      !!
240      !! ** Method  :   set 3D coord. arrays to reference 1D array
241      !!----------------------------------------------------------------------
242      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
243      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
244      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
245      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
246      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
247      !
248      INTEGER  ::   ji, jj, jk
249      !!----------------------------------------------------------------------
250      !
251      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
252      !
253!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
254      DO jk = 1, jpk
255         DO jj = 1, jpj
256            DO ji = 1, jpi
257               pdept(ji,jj,jk) = pdept_1d(jk)
258               pdepw(ji,jj,jk) = pdepw_1d(jk)
259               pe3t (ji,jj,jk) = pe3t_1d (jk)
260               pe3u (ji,jj,jk) = pe3t_1d (jk)
261               pe3v (ji,jj,jk) = pe3t_1d (jk)
262               pe3f (ji,jj,jk) = pe3t_1d (jk)
263               pe3w (ji,jj,jk) = pe3w_1d (jk)
264               pe3uw(ji,jj,jk) = pe3w_1d (jk)
265               pe3vw(ji,jj,jk) = pe3w_1d (jk)
266            END DO
267         END DO
268      END DO
269      !
270      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
271      !
272   END SUBROUTINE zgr_zco
273
274   !!======================================================================
275END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.