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 branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR – NEMO

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

Last change on this file since 7164 was 7164, checked in by gm, 7 years ago

#1692 - branch SIMPLIF_2_usrdef: comment update only
File size: 13.5 KB
RevLine 
[6667]1MODULE usrdef_zgr
[6923]2   !!======================================================================
3   !!                       ***  MODULE  usrdef_zgr  ***
[6667]4   !!
[6923]5   !!                       ===  GYRE configuration  ===
[6667]6   !!
[6923]7   !! User defined : vertical coordinate system of a user configuration
8   !!======================================================================
[6667]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   !
21   USE in_out_manager    ! I/O manager
22   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp           ! distributed memory computing library
24   USE wrk_nemo          ! Memory allocation
25   USE timing            ! Timing
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)
[6923]36   !! $Id$
[6667]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 : GYRE configuration (z-coordinate closed flat box ocean without cavities)'
69      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
70      !
71      !
72      ! type of vertical coordinate
73      ! ---------------------------
74      ld_zco    = .TRUE.         ! GYRE case:  z-coordinate & no 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      !
[6904]84      CALL zgr_msk_top_bot( k_top , k_bot )                  ! masked top and bottom ocean t-level indices
[6667]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      !!
[7164]99      !! ** Purpose :   set the 1D depth of model levels and the resulting
[6904]100      !!              vertical scale factors.
[6667]101      !!
102      !! ** Method  :   z-coordinate system (use in all type of coordinate)
[6904]103      !!              The depth of model levels is defined from an analytical
104      !!              function the derivative of which gives the scale factors.
105      !!              both depth and scale factors only depend on k (1d arrays).
106      !!                 w-level: pdepw_1d  = pdep(k)
107      !!                          pe3w_1d(k) = dk(pdep)(k)     = e3(k)
108      !!                 t-level: pdept_1d  = pdep(k+0.5)
109      !!                          pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5)
[6667]110      !!
111      !!      Here the Madec & Imbard (1996) function is used
112      !!
113      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
114      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
115      !!
116      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
117      !!             Madec and Imbard, 1996, Clim. Dyn.
118      !!----------------------------------------------------------------------
119      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
120      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
121      !
122      INTEGER  ::   jk       ! dummy loop indices
123      REAL(wp) ::   zt, zw   ! local scalars
124      REAL(wp) ::   zsur, za0, za1, zkth, zacr   ! Values for the Madec & Imbard (1996) function 
125      !!----------------------------------------------------------------------
126      !
127      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
128      !
[7164]129      ! Set parameters of z(k) function
130      ! -------------------------------
[6667]131      zsur = -2033.194295283385_wp       
132      za0  =   155.8325369664153_wp 
133      za1  =   146.3615918601890_wp
134      zkth =    17.28520372419791_wp   
135      zacr =     5.0_wp       
136      !
[7164]137      IF(lwp) THEN            ! Parameter print
[6667]138         WRITE(numout,*)
139         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates '
140         WRITE(numout,*) '    ~~~~~~~'
141         WRITE(numout,*) '       GYRE case : MI96 function with the following coefficients :'
142         WRITE(numout,*) '                 zsur = ', zsur
143         WRITE(numout,*) '                 za0  = ', za0
144         WRITE(numout,*) '                 za1  = ', za1
145         WRITE(numout,*) '                 zkth = ', zkth
146         WRITE(numout,*) '                 zacr = ', zacr
147      ENDIF
148
[7164]149      !
150      ! 1D Reference z-coordinate    (using Madec & Imbard 1996 function)
151      ! -------------------------
152      !
153      DO jk = 1, jpk          ! depth at T and W-points
[6667]154         zw = REAL( jk , wp )
155         zt = REAL( jk , wp ) + 0.5_wp
156         pdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr *  LOG( COSH( (zw-zkth) / zacr ) )  )
157         pdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr *  LOG( COSH( (zt-zkth) / zacr ) )  )
158         pe3w_1d (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
159         pe3t_1d (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
160      END DO
161      pdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
162
163
164!!gm   This should become the reference !
165!      IF ( ln_isfcav ) THEN
166! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
167! define pe3t_0 and pe3w_0 as the differences between pdept and pdepw respectively
168!         DO jk = 1, jpkm1
169!            pe3t_1d(jk) = pdepw_1d(jk+1)-pdepw_1d(jk)
170!         END DO
171!         pe3t_1d(jpk) = pe3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
172!
173!         DO jk = 2, jpk
174!            pe3w_1d(jk) = pdept_1d(jk) - pdept_1d(jk-1)
175!         END DO
176!         pe3w_1d(1  ) = 2._wp * (pdept_1d(1) - pdepw_1d(1))
177!      END IF
178!!gm end
179
180      IF(lwp) THEN                        ! control print
181         WRITE(numout,*)
182         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
183         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
184         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
185      ENDIF
186      DO jk = 1, jpk                      ! control positivity
187         IF( pe3w_1d (jk) <= 0._wp .OR. pe3t_1d (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
188         IF( pdepw_1d(jk) <  0._wp .OR. pdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
189      END DO
190      !
191      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
192      !
193   END SUBROUTINE zgr_z
194
195
[6904]196   SUBROUTINE zgr_msk_top_bot( k_top , k_bot )
[6667]197      !!----------------------------------------------------------------------
[6904]198      !!                    ***  ROUTINE zgr_msk_top_bot  ***
[6667]199      !!
[6904]200      !! ** Purpose :   set the masked top and bottom ocean t-levels
[6667]201      !!
202      !! ** Method  :   GYRE case = closed flat box ocean without ocean cavities
203      !!                   k_top = 1     except along north, south, east and west boundaries
204      !!                   k_bot = jpk-1 except along north, south, east and west boundaries
205      !!
[7164]206      !! ** Action  : - k_top : first wet ocean level index
207      !!              - k_bot : last  wet ocean level index
[6667]208      !!----------------------------------------------------------------------
[7164]209      INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last wet ocean level
[6667]210      !
211      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace
212      !!----------------------------------------------------------------------
213      !
214      IF(lwp) WRITE(numout,*)
[7164]215      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.'
[6667]216      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
217      IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities'
218      !
219      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
220      !
[6904]221      CALL lbc_lnk( z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
[6667]222      !
223      k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
224      !
[6904]225      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere
[6667]226      !
[6904]227   END SUBROUTINE zgr_msk_top_bot
[6667]228   
229
230   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
231      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
232      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
233      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
234      !!----------------------------------------------------------------------
235      !!                  ***  ROUTINE zgr_zco  ***
236      !!
237      !! ** Purpose :   define the reference z-coordinate system
238      !!
239      !! ** Method  :   set 3D coord. arrays to reference 1D array
240      !!----------------------------------------------------------------------
241      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
242      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
243      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
244      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
245      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
246      !
247      INTEGER  ::   jk
248      !!----------------------------------------------------------------------
249      !
250      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
251      !
252      DO jk = 1, jpk
253         pdept(:,:,jk) = pdept_1d(jk)
254         pdepw(:,:,jk) = pdepw_1d(jk)
255         pe3t (:,:,jk) = pe3t_1d (jk)
256         pe3u (:,:,jk) = pe3t_1d (jk)
257         pe3v (:,:,jk) = pe3t_1d (jk)
258         pe3f (:,:,jk) = pe3t_1d (jk)
259         pe3w (:,:,jk) = pe3w_1d (jk)
260         pe3uw(:,:,jk) = pe3w_1d (jk)
261         pe3vw(:,:,jk) = pe3w_1d (jk)
262      END DO
263      !
264      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
265      !
266   END SUBROUTINE zgr_zco
267
268   !!======================================================================
269END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.