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 @ 6923

Last change on this file since 6923 was 6923, checked in by gm, 8 years ago

#1692 - branch SIMPLIF_2_usrdef: update comments in usrdef modules

File size: 13.1 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  :   z-coordinate system (use in all type of coordinate)
221      !!      The depth of model levels is defined from an analytical
222      !!      function the derivative of which gives the scale factors.
223      !!      both depth and scale factors only depend on k (1d arrays).
224      !!              w-level: pdepw_1d  = pdep(k)
225      !!                       pe3w_1d(k) = dk(pdep)(k)     = e3(k)
226      !!              t-level: pdept_1d  = pdep(k+0.5)
227      !!                       pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5)
228      !!
229      !!            ===    Here constant vertical resolution   ===
230      !!
231      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
232      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
233      !!----------------------------------------------------------------------
234      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
235      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
236      !
237      INTEGER  ::   jk       ! dummy loop indices
238      REAL(wp) ::   zt, zw   ! local scalar
239      !!----------------------------------------------------------------------
240      !
241      IF(lwp) THEN                         ! Parameter print
242         WRITE(numout,*)
243         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
244         WRITE(numout,*) '    ~~~~~~~'
245      ENDIF
246      !
247      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
248      ! ----------------------
249      DO jk = 1, jpk
250         zw = REAL( jk , wp )
251         zt = REAL( jk , wp ) + 0.5_wp
252         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
253         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
254         pe3w_1d (jk) =    rn_dz
255         pe3t_1d (jk) =    rn_dz
256      END DO
257      !
258      IF(lwp) THEN                        ! control print
259         WRITE(numout,*)
260         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
261         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
262         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
263      ENDIF
264      !
265   END SUBROUTINE zgr_z1d
266   
267   !!======================================================================
268END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.