source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_zgr.F90 @ 7188

Last change on this file since 7188 was 7188, checked in by gm, 4 years ago

#1692 - branch SIMPLIF_2_usrdef: e3.=dk[dep.] (discret derivative)

File size: 12.9 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_cfg.nc file. Indeed,
115      !!    Only e3. are saved in this file, depth are compute by a call to
116      !!    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      !
203      IF(lwp) WRITE(numout,*)
204      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.'
205      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
206      IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities'
207      !
208      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
209      !
210      CALL lbc_lnk( z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
211      !
212      k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
213      !
214      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere
215      !
216   END SUBROUTINE zgr_msk_top_bot
217   
218
219   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
220      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
221      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
222      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
223      !!----------------------------------------------------------------------
224      !!                  ***  ROUTINE zgr_zco  ***
225      !!
226      !! ** Purpose :   define the reference z-coordinate system
227      !!
228      !! ** Method  :   set 3D coord. arrays to reference 1D array
229      !!----------------------------------------------------------------------
230      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
231      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
232      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
233      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
234      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
235      !
236      INTEGER  ::   jk
237      !!----------------------------------------------------------------------
238      !
239      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
240      !
241      DO jk = 1, jpk
242         pdept(:,:,jk) = pdept_1d(jk)
243         pdepw(:,:,jk) = pdepw_1d(jk)
244         pe3t (:,:,jk) = pe3t_1d (jk)
245         pe3u (:,:,jk) = pe3t_1d (jk)
246         pe3v (:,:,jk) = pe3t_1d (jk)
247         pe3f (:,:,jk) = pe3t_1d (jk)
248         pe3w (:,:,jk) = pe3w_1d (jk)
249         pe3uw(:,:,jk) = pe3w_1d (jk)
250         pe3vw(:,:,jk) = pe3w_1d (jk)
251      END DO
252      !
253      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
254      !
255   END SUBROUTINE zgr_zco
256
257   !!======================================================================
258END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.