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/CONFIG/OVERFLOW/MY_SRC – NEMO

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_zgr.F90 @ 7188

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

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

File size: 13.0 KB
Line 
1MODULE usrdef_zgr
2   !!======================================================================
3   !!                     ***  MODULE usrdef_zgr  ***
4   !!
5   !!                       ===  OVERFLOW case  ===
6   !!
7   !! user defined :  vertical coordinate system of a user configuration
8   !!======================================================================
9   !! History :  4.0  ! 2016-08  (G. Madec, S. Flavoni)  Original code
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_zgr   : user defined vertical coordinate system (required)
14   !!       zgr_z1d   : reference 1D z-coordinate
15   !!---------------------------------------------------------------------
16   USE oce            ! ocean variables
17   USE dom_oce ,  ONLY: ln_zco, ln_zps, ln_sco   ! ocean space and time domain
18   USE dom_oce ,  ONLY: mi0, mi1, nimpp, njmpp   ! ocean space and time domain
19   USE dom_oce ,  ONLY: glamt                    ! ocean space and time domain
20   USE usrdef_nam     ! User defined : namelist variables
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  ::   ji, jj, jk        ! dummy indices
64      INTEGER  ::   ik                ! local integers
65      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
66      REAL(wp) ::   ze3min            ! local scalar
67      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, z2d   ! 2D workspace
68      !!----------------------------------------------------------------------
69      !
70      IF(lwp) WRITE(numout,*)
71      IF(lwp) WRITE(numout,*) 'usr_def_zgr : OVERFLOW configuration (z(ps)- or s-coordinate closed box ocean without cavities)'
72      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
73      !
74      !
75      ! type of vertical coordinate
76      ! ---------------------------
77      ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
78      ld_isfcav = .FALSE.
79      !
80      !
81      ! Build the vertical coordinate system
82      ! ------------------------------------
83      !
84      !                       !==  UNmasked meter bathymetry  ==!
85      !
86      ! western continental shelf (500m deep) and eastern deep ocean (2000m deep)
87      ! (set through the jpk and jpi (see userdef_nam.F90))
88      ! with a hyperbolic tangent transition centered at 40km
89      ! NB: here glamt is in kilometers
90      !
91      zht(:,:) = + (  500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,:) - 40.) / 7. ) )  )
92      !
93      ! at u-point: averaging zht
94      DO ji = 1, jpim1
95         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
96      END DO
97      CALL lbc_lnk( zhu, 'U', 1. )     ! boundary condition: this mask the surrouding grid-points
98      !                                ! ==>>>  set by hand non-zero value on first/last columns & rows
99      DO ji = mi0(1), mi1(1)              ! first row of global domain only
100         zhu(ji,2) = zht(1,2)
101      END DO
102       DO ji = mi0(jpi), mi1(jpi)         ! last  row of global domain only
103         zhu(ji,2) = zht(jpi,2)
104      END DO
105      zhu(:,1) = zhu(:,2)
106      zhu(:,3) = zhu(:,2)
107      !     
108      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
109      !
110      !
111      !                       !==  top masked level bathymetry  ==!  (all coordinates)
112      !
113      ! no ocean cavities : top ocean level is ONE, except over land
114      ! the ocean basin surrounded by land (1 grid-point) set through lbc_lnk call as jperio=0
115      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
116      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
117      k_top(:,:) = NINT( z2d(:,:) )
118      !
119      !                             
120      !
121      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
122         !
123         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
124         !
125         !                                !* terrain-following coordinate with e3.(k)=cst)
126         !                                !  OVERFLOW case : identical with j-index (T=V, U=F)
127         z1_jpkm1 = 1._wp / REAL( jpkm1 , wp)
128         DO jk = 1, jpk
129            pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
130            pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp )          )
131            pe3t (:,:,jk) = zht(:,:) * z1_jpkm1
132            pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1
133            pe3v (:,:,jk) = zht(:,:) * z1_jpkm1
134            pe3f (:,:,jk) = zhu(:,:) * z1_jpkm1
135            pe3w (:,:,jk) = zht(:,:) * z1_jpkm1
136            pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1
137            pe3vw(:,:,jk) = zht(:,:) * z1_jpkm1
138         END DO     
139      ENDIF
140      !
141      !
142      IF ( ln_zco ) THEN      !==  z-coordinate  ==!   (step-like topography)
143         !
144         !                                !* bottom ocean compute from the depth of grid-points
145         k_bot(:,:) = jpkm1 * k_top(:,:)     ! here use k_top as a land mask
146         DO jk = 1, jpkm1
147            WHERE( pdept_1d(jk) < zht(:,:) .AND. zht(:,:) <= pdept_1d(jk+1) )   k_bot(:,:) = jk * k_top(:,:)
148         END DO
149         !                                !* horizontally uniform coordinate (reference z-co everywhere)
150         DO jk = 1, jpk
151            pdept(:,:,jk) = pdept_1d(jk)
152            pdepw(:,:,jk) = pdepw_1d(jk)
153            pe3t (:,:,jk) = pe3t_1d (jk)
154            pe3u (:,:,jk) = pe3t_1d (jk)
155            pe3v (:,:,jk) = pe3t_1d (jk)
156            pe3f (:,:,jk) = pe3t_1d (jk)
157            pe3w (:,:,jk) = pe3w_1d (jk)
158            pe3uw(:,:,jk) = pe3w_1d (jk)
159            pe3vw(:,:,jk) = pe3w_1d (jk)
160         END DO
161      ENDIF
162      !
163      !
164      IF ( ln_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps)
165         !
166         ze3min = 0.1_wp * rn_dz
167         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min
168         !
169         !
170         !                                !* bottom ocean compute from the depth of grid-points
171         k_bot(:,:) = jpkm1
172         DO jk = jpkm1, 1, -1
173            WHERE( zht(:,:) < pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
174         END DO
175         !
176         !                                !* vertical coordinate system
177         DO jk = 1, jpk                      ! initialization to the reference z-coordinate
178            pdept(:,:,jk) = pdept_1d(jk)
179            pdepw(:,:,jk) = pdepw_1d(jk)
180            pe3t (:,:,jk) = pe3t_1d (jk)
181            pe3u (:,:,jk) = pe3t_1d (jk)
182            pe3v (:,:,jk) = pe3t_1d (jk)
183            pe3f (:,:,jk) = pe3t_1d (jk)
184            pe3w (:,:,jk) = pe3w_1d (jk)
185            pe3uw(:,:,jk) = pe3w_1d (jk)
186            pe3vw(:,:,jk) = pe3w_1d (jk)
187         END DO
188         DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points
189            DO ji = 1, jpi
190               ik = k_bot(ji,jj)
191                  pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
192                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
193                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
194                  !
195                  pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
196                  pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
197                  pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  )
198            END DO
199         END DO         
200         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points
201         !                                   ! usually Computed as the minimum of neighbooring scale factors
202         pe3u (:,:,:) = pe3t(:,:,:)          ! HERE OVERFLOW configuration :
203         pe3v (:,:,:) = pe3t(:,:,:)          !    e3 increases with i-index and identical with j-index
204         pe3f (:,:,:) = pe3t(:,:,:)          !    so e3 minimum of (i,i+1) points is (i) point
205         pe3uw(:,:,:) = pe3w(:,:,:)          !    in j-direction e3v=e3t and e3f=e3v
206         pe3vw(:,:,:) = pe3w(:,:,:)          !    ==>>  no need of lbc_lnk calls
207         !     
208      ENDIF
209      !
210   END SUBROUTINE usr_def_zgr
211
212
213   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
214      !!----------------------------------------------------------------------
215      !!                   ***  ROUTINE zgr_z1d  ***
216      !!
217      !! ** Purpose :   set the depth of model levels and the resulting
218      !!      vertical scale factors.
219      !!
220      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
221      !!       The depth of model levels is set from dep(k), an analytical function:
222      !!                   w-level: depw_1d  = dep(k)
223      !!                   t-level: dept_1d  = dep(k+0.5)
224      !!       The scale factors are the discrete derivative of the depth:
225      !!                   e3w_1d(jk) = dk[ dept_1d ]
226      !!                   e3t_1d(jk) = dk[ depw_1d ]
227      !!
228      !!            ===    Here constant vertical resolution   ===
229      !!
230      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
231      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
232      !!----------------------------------------------------------------------
233      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
234      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
235      !
236      INTEGER  ::   jk       ! dummy loop indices
237      REAL(wp) ::   zt, zw   ! local scalar
238      !!----------------------------------------------------------------------
239      !
240      IF(lwp) THEN                         ! Parameter print
241         WRITE(numout,*)
242         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
243         WRITE(numout,*) '    ~~~~~~~'
244      ENDIF
245      !
246      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
247      ! ----------------------
248      DO jk = 1, jpk
249         zw = REAL( jk , wp )
250         zt = REAL( jk , wp ) + 0.5_wp
251         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
252         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
253         pe3w_1d (jk) =    rn_dz
254         pe3t_1d (jk) =    rn_dz
255      END DO
256      !
257      IF(lwp) THEN                        ! control print
258         WRITE(numout,*)
259         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
260         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
261         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
262      ENDIF
263      !
264   END SUBROUTINE zgr_z1d
265   
266   !!======================================================================
267END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.