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 NEMO/branches/UKMO/NEMO_4.0.2_SEAMOUNT/tests/SEAMOUNT/MY_SRC – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.2_SEAMOUNT/tests/SEAMOUNT/MY_SRC/usrdef_zgr.F90 @ 12868

Last change on this file since 12868 was 12868, checked in by ayoung, 4 years ago

SEAMOUNT configuration files (EXPREF and MY_SRC)

File size: 15.0 KB
Line 
1MODULE usrdef_zgr
2   !!======================================================================
3   !!                   ***  MODULE  usrdef_zgr  ***
4   !!
5   !!                   ===  SEAMOUNT Test Case  ===
6   !!
7   !! Ocean domain : user defined vertical coordinate system
8   !!======================================================================
9   !! History :  4.0  ! 2016-06  (G. Madec)  Original code
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_zgr   : user defined vertical coordinate system (required)
14   !!       zgr_z     : reference 1D z-coordinate
15   !!---------------------------------------------------------------------
16   USE oce            ! ocean variables
17   USE dom_oce ,  ONLY: ht_0, mi0, mi1, nimpp, njmpp,  &
18                      & mj0, mj1, glamt, gphit         ! ocean space and time domain
19   USE usrdef_nam     ! User defined : namelist variables
20   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld  ! Wetting and drying
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: usrdef_zgr.F90 10425 2018-12-19 21:54:16Z smasson $
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(in) ::   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  ::   ji, jj, jk        ! dummy indices
62      INTEGER  ::   ik                ! local integers
63      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
64      REAL(wp) ::   ze3min            ! local scalar
65      REAL(wp) ::   zlam_mid, zphi_mid! local scalar
66      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, zhv, z2d   ! 2D workspace
67      !!----------------------------------------------------------------------
68      !
69      IF(lwp) WRITE(numout,*)
70      IF(lwp) WRITE(numout,*) 'usr_def_zgr : SEAMOUNT_TEST_CASE configuration (s-coordinate closed box ocean without cavities)'
71      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
72      !
73      !
74      ! type of vertical coordinate ==>>>   here SEAMOUNT_TEST_CASE : z, z-ps and s-coordinate available
75      ! ---------------------------
76      ld_isfcav = .FALSE.      ! ISF Ice Shelves Flag
77      !
78      !
79      ! Build the vertical coordinate system
80      ! ------------------------------------
81      !
82      !                       !==  UNmasked meter bathymetry  ==!
83      !
84      !
85      !
86      IF(lwp) WRITE(numout,*)
87      IF(lwp) WRITE(numout,*) 'usr_def_zgr (SEAMOUNT) : Isolated 2D Gaussian bump in E-W periodic channel'
88      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
89      !
90      zlam_mid = 0.5_wp * 1000._wp * rn_length
91      zphi_mid = 0.5_wp * 1000._wp * rn_width
92      zht = 0._wp
93      DO ji = 1, jpi
94         DO jj = 1, jpj
95            zht(ji,jj) = rn_bathy - rn_seamountheight * EXP( &
96                       & - ( ( 1000._wp * glamt(ji,jj) - zlam_mid) ** 2 + ( 1000._wp * gphit(ji,jj) - zphi_mid) ** 2) & 
97                       & / rn_l ** 2)
98         END DO
99      END DO
100      !
101      ! at u-point: averaging zht
102      DO ji = 1, jpim1
103         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
104      END DO
105 !     CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points
106      !                                ! ==>>>  set by hand non-zero value on first/last columns & rows
107      DO ji = mi0(1), mi1(1)              ! first row of global domain only
108         zhu(ji,:) = zht(1,:)
109      END DO
110       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only
111         zhu(ji,:) = zht(jpi,:)
112      END DO
113      ! at v-point: averaging zht
114      zhv = 0._wp
115      DO jj = 1, jpjm1
116         zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) )
117      END DO
118!      CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points
119      DO jj = mj0(1), mj1(1)   ! first  row of global domain only
120         zhv(:,jj) = zht(:,jj)
121      END DO
122      DO jj = mj0(jpjglo), mj1(jpjglo)   ! last  row of global domain only
123         zhv(:,jj) = zht(:,jj)
124      END DO
125      !     
126      CALL zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
127      !
128      !
129      !                       !==  top masked level bathymetry  ==!  (all coordinates)
130      !
131      ! no ocean cavities : top ocean level is ONE, except over land
132      ! the ocean basin surrounnded by land (1 grid-point) set through lbc_lnk call as jperio=0
133      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
134!      z2d(mi0(1):mi1(1),:) = 0._wp
135!      z2d(mi0(jpiglo):mi1(jpiglo),:) = 0._wp
136!      z2d(:,mj0(1):mj1(1)) = 0._wp
137!      z2d(:,mj0(jpjglo):mj1(jpjglo)) = 0._wp
138      !
139      !
140!      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
141      k_top(:,:) = NINT( z2d(:,:) )
142      !
143      !                             
144      !
145      IF ( ld_zco ) THEN      !==  z-coordinate  ==!   (step-like topography)
146         !
147         !                                !* bottom ocean compute from the depth of grid-points
148         k_bot(:,:) = jpkm1 * k_top(:,:)     ! here use k_top as a land mask
149         DO jk = 1, jpkm1
150            WHERE( pdept_1d(jk) < zht(:,:) .AND. zht(:,:) <= pdept_1d(jk+1) )   k_bot(:,:) = jk * k_top(:,:)
151         END DO
152         !                                !* horizontally uniform coordinate (reference z-co everywhere)
153         DO jk = 1, jpk
154            pdept(:,:,jk) = pdept_1d(jk)
155            pdepw(:,:,jk) = pdepw_1d(jk)
156            pe3t (:,:,jk) = pe3t_1d (jk)
157            pe3u (:,:,jk) = pe3t_1d (jk)
158            pe3v (:,:,jk) = pe3t_1d (jk)
159            pe3f (:,:,jk) = pe3t_1d (jk)
160            pe3w (:,:,jk) = pe3w_1d (jk)
161            pe3uw(:,:,jk) = pe3w_1d (jk)
162            pe3vw(:,:,jk) = pe3w_1d (jk)
163         END DO
164      ENDIF
165      IF ( ld_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps)
166         !
167         ze3min = 0.1_wp * rn_dz
168         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min
169         !
170         !
171         !                                !* bottom ocean compute from the depth of grid-points
172         k_bot(:,:) = jpkm1
173         DO jk = jpkm1, 1, -1
174            WHERE( zht(:,:) < pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
175         END DO
176         !
177         !                                !* vertical coordinate system
178         DO jk = 1, jpk                      ! initialization to the reference z-coordinate
179            pdept(:,:,jk) = pdept_1d(jk)
180            pdepw(:,:,jk) = pdepw_1d(jk)
181            pe3t (:,:,jk) = pe3t_1d (jk)
182            pe3u (:,:,jk) = pe3t_1d (jk)
183            pe3v (:,:,jk) = pe3t_1d (jk)
184            pe3f (:,:,jk) = pe3t_1d (jk)
185            pe3w (:,:,jk) = pe3w_1d (jk)
186            pe3uw(:,:,jk) = pe3w_1d (jk)
187            pe3vw(:,:,jk) = pe3w_1d (jk)
188         END DO
189         DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points
190            DO ji = 1, jpi
191               ik = k_bot(ji,jj)
192                  pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
193                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
194                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )
195                  !
196                  pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
197                  pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
198                  pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  )
199            END DO
200         END DO
201         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points
202         !                                   ! usually Computed as the minimum of neighbooring scale factors
203         pe3u (:,:,:) = pe3t(:,:,:)          ! HERE OVERFLOW configuration :
204         pe3v (:,:,:) = pe3t(:,:,:)          !    e3 increases with i-index and identical with j-index
205         pe3f (:,:,:) = pe3t(:,:,:)          !    so e3 minimum of (i,i+1) points is (i) point
206         pe3uw(:,:,:) = pe3w(:,:,:)          !    in j-direction e3v=e3t and e3f=e3v
207         pe3vw(:,:,:) = pe3w(:,:,:)          !    ==>>  no need of lbc_lnk calls
208         !     
209      ENDIF
210      !
211      IF ( ld_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
212         !
213         ht_0 = zht
214         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
215         !                                !* terrain-following coordinate with e3.(k)=cst)
216         !                                !  OVERFLOW case : identical with j-index (T=V, U=F)
217         DO jj = 1, jpjm1
218            DO ji = 1, jpim1
219              z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp)
220              DO jk = 1, jpk
221                  pdept(ji,jj,jk) = zht(ji,jj) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
222                  pdepw(ji,jj,jk) = zht(ji,jj) * z1_jpkm1 * ( REAL( jk-1 , wp )          )
223                  pe3t (ji,jj,jk) = zht(ji,jj) * z1_jpkm1
224                  pe3w (ji,jj,jk) = zht(ji,jj) * z1_jpkm1
225                  pe3u (ji,jj,jk) = zhu(ji,jj) * z1_jpkm1
226                  pe3uw(ji,jj,jk) = zhu(ji,jj) * z1_jpkm1
227                  pe3f (ji,jj,jk) = zhu(ji,jj) * z1_jpkm1
228                  pe3v (ji,jj,jk) = zhv(ji,jj) * z1_jpkm1
229                  pe3vw(ji,jj,jk) = zhv(ji,jj) * z1_jpkm1
230              END DO     
231           END DO     
232         END DO     
233         CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1. )
234         CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1. )
235         CALL lbc_lnk( 'usrdef_zgr', pe3t , 'T', 1. )
236         CALL lbc_lnk( 'usrdef_zgr', pe3w , 'T', 1. )
237         CALL lbc_lnk( 'usrdef_zgr', pe3u , 'U', 1. )
238         CALL lbc_lnk( 'usrdef_zgr', pe3uw, 'U', 1. )
239         CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1. )
240         CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1. )
241         CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1. )
242         WHERE( pe3t (:,:,:) == 0._wp )   pe3t (:,:,:) = 1._wp
243         WHERE( pe3u (:,:,:) == 0._wp )   pe3u (:,:,:) = 1._wp
244         WHERE( pe3v (:,:,:) == 0._wp )   pe3v (:,:,:) = 1._wp
245         WHERE( pe3f (:,:,:) == 0._wp )   pe3f (:,:,:) = 1._wp
246         WHERE( pe3w (:,:,:) == 0._wp )   pe3w (:,:,:) = 1._wp
247         WHERE( pe3uw(:,:,:) == 0._wp )   pe3uw(:,:,:) = 1._wp
248         WHERE( pe3vw(:,:,:) == 0._wp )   pe3vw(:,:,:) = 1._wp
249      ENDIF
250      !
251      !
252      !
253   END SUBROUTINE usr_def_zgr
254
255
256   SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
257      !!----------------------------------------------------------------------
258      !!                   ***  ROUTINE zgr_z  ***
259      !!
260      !! ** Purpose :   set the depth of model levels and the resulting
261      !!      vertical scale factors.
262      !!
263      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
264      !!       The depth of model levels is set from dep(k), an analytical function:
265      !!                   w-level: depw_1d  = dep(k)
266      !!                   t-level: dept_1d  = dep(k+0.5)
267      !!       The scale factors are the discrete derivative of the depth:
268      !!                   e3w_1d(jk) = dk[ dept_1d ]
269      !!                   e3t_1d(jk) = dk[ depw_1d ]
270      !!
271      !!            ===    Here constant vertical resolution   ===
272      !!
273      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
274      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
275      !!----------------------------------------------------------------------
276      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
277      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
278      !
279      INTEGER  ::   jk       ! dummy loop indices
280      REAL(wp) ::   zt, zw   ! local scalar
281      !!----------------------------------------------------------------------
282      !
283      IF(lwp) THEN                         ! Parameter print
284         WRITE(numout,*)
285         WRITE(numout,*) '    zgr_z : Reference vertical z-coordinates: uniform dz = ', rn_dz
286         WRITE(numout,*) '    ~~~~~~~'
287      ENDIF
288      !
289      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
290      ! ----------------------
291      DO jk = 1, jpk
292         zw = REAL( jk , wp )
293         zt = REAL( jk , wp ) + 0.5_wp
294         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
295         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
296         pe3w_1d (jk) =    rn_dz
297         pe3t_1d (jk) =    rn_dz
298      END DO
299      !
300      IF(lwp) THEN                        ! control print
301         WRITE(numout,*)
302         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
303         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
304         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
305      ENDIF
306      !
307   END SUBROUTINE zgr_z
308   
309   !!======================================================================
310END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.