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
RevLine 
[6895]1MODULE usrdef_zgr
[6923]2   !!======================================================================
3   !!                     ***  MODULE usrdef_zgr  ***
[6895]4   !!
[6923]5   !!                       ===  OVERFLOW case  ===
[6895]6   !!
[6923]7   !! user defined :  vertical coordinate system of a user configuration
8   !!======================================================================
[6914]9   !! History :  4.0  ! 2016-08  (G. Madec, S. Flavoni)  Original code
[6895]10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
[6914]13   !!   usr_def_zgr   : user defined vertical coordinate system (required)
14   !!       zgr_z1d   : reference 1D z-coordinate
[6895]15   !!---------------------------------------------------------------------
[6914]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
[6895]21   !
[6914]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
[6895]27
28   IMPLICIT NONE
29   PRIVATE
30
[6923]31   PUBLIC   usr_def_zgr   ! called by domzgr.F90
[6895]32
33  !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
[6923]37   !! $Id$
[6895]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
[6904]45      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
46      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
[6895]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      !
[6904]63      INTEGER  ::   ji, jj, jk        ! dummy indices
[6905]64      INTEGER  ::   ik                ! local integers
[6904]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
[6895]68      !!----------------------------------------------------------------------
69      !
70      IF(lwp) WRITE(numout,*)
[6914]71      IF(lwp) WRITE(numout,*) 'usr_def_zgr : OVERFLOW configuration (z(ps)- or s-coordinate closed box ocean without cavities)'
[6895]72      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
73      !
74      !
75      ! type of vertical coordinate
76      ! ---------------------------
[6904]77      ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
[6895]78      ld_isfcav = .FALSE.
79      !
80      !
81      ! Build the vertical coordinate system
82      ! ------------------------------------
83      !
[6904]84      !                       !==  UNmasked meter bathymetry  ==!
[6895]85      !
[6904]86      ! western continental shelf (500m deep) and eastern deep ocean (2000m deep)
[6914]87      ! (set through the jpk and jpi (see userdef_nam.F90))
[6904]88      ! with a hyperbolic tangent transition centered at 40km
89      ! NB: here glamt is in kilometers
[6895]90      !
[6904]91      zht(:,:) = + (  500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,:) - 40.) / 7. ) )  )
92      !
[6914]93      ! at u-point: averaging zht
94      DO ji = 1, jpim1
95         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
[6904]96      END DO
[6914]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      !     
[6904]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
[6914]116      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
117      k_top(:,:) = NINT( z2d(:,:) )
[6904]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)
[6905]126         !                                !  OVERFLOW case : identical with j-index (T=V, U=F)
[6904]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
[6914]173            WHERE( zht(:,:) < pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
[6904]174         END DO
[6905]175         !
[6904]176         !                                !* vertical coordinate system
[6905]177         DO jk = 1, jpk                      ! initialization to the reference z-coordinate
[6904]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
[6905]188         DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points
[6904]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) )
[6905]192                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
[6914]193                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
[6905]194                  !
[6914]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  )
[6904]198            END DO
199         END DO         
[6905]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
[6914]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
[6905]206         pe3vw(:,:,:) = pe3w(:,:,:)          !    ==>>  no need of lbc_lnk calls
207         !     
[6904]208      ENDIF
209      !
[6895]210   END SUBROUTINE usr_def_zgr
211
212
[6904]213   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
[6895]214      !!----------------------------------------------------------------------
[6904]215      !!                   ***  ROUTINE zgr_z1d  ***
[6895]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      !!
[6904]229      !!            ===    Here constant vertical resolution   ===
[6895]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      !!----------------------------------------------------------------------
[6904]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]
[6895]236      !
237      INTEGER  ::   jk       ! dummy loop indices
[6904]238      REAL(wp) ::   zt, zw   ! local scalar
[6895]239      !!----------------------------------------------------------------------
240      !
241      IF(lwp) THEN                         ! Parameter print
242         WRITE(numout,*)
[6904]243         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
[6895]244         WRITE(numout,*) '    ~~~~~~~'
245      ENDIF
[6904]246      !
[6895]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
[6904]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
[6895]256      END DO
[6904]257      !
[6895]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      !
[6904]265   END SUBROUTINE zgr_z1d
[6895]266   
267   !!======================================================================
268END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.