source: NEMO/trunk/src/OCE/USR/usrdef_zgr.F90 @ 10425

Last change on this file since 10425 was 10425, checked in by smasson, 21 months ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

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