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

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

#1692 - branch SIMPLIF_2_usrdef: OVERFLOW configuration (zco & sco) + small bug corrections
File size: 18.3 KB
Line 
1MODULE usrdef_zgr
2   !!==============================================================================
3   !!                       ***  MODULE usrdef_zgr  ***
4   !! Ocean domain : user defined vertical coordinate system
5   !!
6   !!                       ===      OVERFLOW case      ===
7   !!
8   !!==============================================================================
9   !! History :  4.0  ! 2016-08  (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   !!       zgr_zps      : 3D vertical coordinate in z-partial cell coordinate
16   !!---------------------------------------------------------------------
17   USE oce               ! ocean variables
18   USE dom_oce  ,  ONLY: ln_zco, ln_zps, ln_sco   ! ocean space and time domain
19   USE dom_oce  ,  ONLY: nimpp, njmpp             ! ocean space and time domain
20   USE dom_oce  ,  ONLY: glamt                    ! ocean space and time domain
21   USE usrdef_nam        ! User defined : namelist variables
22   !
23   USE in_out_manager    ! I/O manager
24   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
25   USE lib_mpp           ! distributed memory computing library
26   USE wrk_nemo          ! Memory allocation
27   USE timing            ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   usr_def_zgr        ! called by domzgr.F90
33
34  !! * Substitutions
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
38   !! $Id: domzgr.F90 6624 2016-05-26 08:59:48Z gm $
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS             
42
43   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
44      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
45      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
46      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
47      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
48      &                    k_top  , k_bot    )                             ! top & bottom ocean level
49      !!---------------------------------------------------------------------
50      !!              ***  ROUTINE usr_def_zgr  ***
51      !!
52      !! ** Purpose :   User defined the vertical coordinates
53      !!
54      !!----------------------------------------------------------------------
55      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
56      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
57      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
58      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
59      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
60      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
61      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
62      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
63      !
64      INTEGER  ::   ji, jj, jk        ! dummy indices
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      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
71      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
72      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t
73      REAL(wp) ::   zdiff            ! temporary scalar
74      REAL(wp) ::   zmax             ! temporary scalar
75
76
77      !!----------------------------------------------------------------------
78      !
79      IF(lwp) WRITE(numout,*)
80      IF(lwp) WRITE(numout,*) 'usr_def_zgr : OVERFLOW configuration (z-coordinate closed box ocean without cavities)'
81      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
82      !
83      !
84      ! type of vertical coordinate
85      ! ---------------------------
86      ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
87      ld_isfcav = .FALSE.
88      !
89      !
90      ! Build the vertical coordinate system
91      ! ------------------------------------
92      !
93      !                       !==  UNmasked meter bathymetry  ==!
94      !
95      ! western continental shelf (500m deep) and eastern deep ocean (2000m deep)
96      ! with a hyperbolic tangent transition centered at 40km
97      ! NB: here glamt is in kilometers
98      !
99      zht(:,:) = + (  500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,:) - 40.) / 7. ) )  )
100      !
101      ! at u-point: recompute from glamu (see usrdef_hgr.F90) to avoid the masking when using lbc_lnk
102      zfact = rn_dx * 1.e-3         ! conversion in km
103      DO ji = 1, jpi
104         zhu(ji,:) = + (  500. + 0.5 * 1500. * ( 1.0 + tanh( ( zfact * REAL( ji-1 + nimpp-1 , wp ) - 40.) / 7. ) )  )
105      END DO
106      !
107      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
108      !
109      !
110      !                       !==  top masked level bathymetry  ==!  (all coordinates)
111      !
112      ! no ocean cavities : top ocean level is ONE, except over land
113      ! the ocean basin surrounded by land (1 grid-point) set through lbc_lnk call as jperio=0
114      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
115      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin
116      k_top(:,:) = z2d(:,:)
117      !
118      !                             
119      !
120      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
121         !
122         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
123         !
124         !                                !* terrain-following coordinate with e3.(k)=cst)
125         z1_jpkm1 = 1._wp / REAL( jpkm1 , wp)
126         DO jk = 1, jpk
127            pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
128            pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp )          )
129            pe3t (:,:,jk) = zht(:,:) * z1_jpkm1
130            pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1
131            pe3v (:,:,jk) = zht(:,:) * z1_jpkm1
132            pe3f (:,:,jk) = zhu(:,:) * z1_jpkm1
133            pe3w (:,:,jk) = zht(:,:) * z1_jpkm1
134            pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1
135            pe3vw(:,:,jk) = zht(:,:) * z1_jpkm1
136         END DO     
137      ENDIF
138      !
139      !
140      IF ( ln_zco ) THEN      !==  z-coordinate  ==!   (step-like topography)
141         !
142         !                                !* bottom ocean compute from the depth of grid-points
143         k_bot(:,:) = jpkm1 * k_top(:,:)     ! here use k_top as a land mask
144         DO jk = 1, jpkm1
145            WHERE( pdept_1d(jk) < zht(:,:) .AND. zht(:,:) <= pdept_1d(jk+1) )   k_bot(:,:) = jk * k_top(:,:)
146         END DO
147         !                                !* horizontally uniform coordinate (reference z-co everywhere)
148         DO jk = 1, jpk
149            pdept(:,:,jk) = pdept_1d(jk)
150            pdepw(:,:,jk) = pdepw_1d(jk)
151            pe3t (:,:,jk) = pe3t_1d (jk)
152            pe3u (:,:,jk) = pe3t_1d (jk)
153            pe3v (:,:,jk) = pe3t_1d (jk)
154            pe3f (:,:,jk) = pe3t_1d (jk)
155            pe3w (:,:,jk) = pe3w_1d (jk)
156            pe3uw(:,:,jk) = pe3w_1d (jk)
157            pe3vw(:,:,jk) = pe3w_1d (jk)
158         END DO
159      ENDIF
160      !
161      !
162      IF ( ln_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps)
163     
164     
165     
166         CALL ctl_stop( 'STOP', ' zps coordinate not yet coded' )
167
168
169
170     
171      !!----------------------------------------------------------------------
172      !!                  ***  ROUTINE zgr_zps  ***
173      !!                     
174      !! ** Purpose :   the depth and vertical scale factor in partial step
175      !!              reference z-coordinate case
176      !!
177      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
178      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
179      !!      a partial step representation of bottom topography.
180      !!
181      !!     
182      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
183      !!----------------------------------------------------------------------
184      !!---------------------------------------------------------------------
185         !
186         ze3min = 0.1_wp * rn_dz
187         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min
188         !
189         !
190         !                                !* bottom ocean compute from the depth of grid-points
191         k_bot(:,:) = jpkm1
192         DO jk = jpkm1, 1, -1
193            WHERE( zht(:,:) <= pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
194         END DO
195
196         !                                !* vertical coordinate system
197         !
198         DO jk = 1, jpk                         ! initialization to the reference z-coordinate
199            pdept(:,:,jk) = pdept_1d(jk)
200            pdepw(:,:,jk) = pdepw_1d(jk)
201            pe3t (:,:,jk) = pe3t_1d (jk)
202            pe3u (:,:,jk) = pe3t_1d (jk)
203            pe3v (:,:,jk) = pe3t_1d (jk)
204            pe3f (:,:,jk) = pe3t_1d (jk)
205            pe3w (:,:,jk) = pe3w_1d (jk)
206            pe3uw(:,:,jk) = pe3w_1d (jk)
207            pe3vw(:,:,jk) = pe3w_1d (jk)
208         END DO
209         !
210         DO jj = 1, jpj                         ! bottom scale factors and depth at T- and W-points
211            DO ji = 1, jpi
212               ik = k_bot(ji,jj)
213               IF( ik /= jpkm1 ) THEN                 ! last level ==> ref 1d z-coordinate
214                  pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
215                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik) - pdepw(ji,jj,ik+1)
216                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik) 
217               
218                  pdept(ji,jj,ik  ) = pdept(ji,jj,ik-1) + pe3t(ji,jj,ik) * 0.5_wp
219                  pe3w (ji,jj,ik+1) = pdepw(ji,jj,ik) - pdepw(ji,jj,ik+1)
220               ENDIF
221            END DO
222         END DO         
223         
224         
225         !
226         DO jj = 1, jpj                         ! bottom scale factors and depth at T- and W-points
227            DO ji = 1, jpi
228               ik = k_bot(ji,jj)
229                  !
230                  IF( zht(ji,jj) <= pdepw_1d(ik+1) ) THEN  ;   pdepw(ji,jj,ik+1) = zht(ji,jj)
231                  ELSE                                     ;   pdepw(ji,jj,ik+1) = pdepw_1d(ik+1)
232                  ENDIF
233   !gm Bug?  check the gdepw_1d
234                  !       ... on ik
235                  pdept(ji,jj,ik) = pdepw_1d(ik) + ( pdepw   (ji,jj,ik+1) - pdepw_1d(ik) )   &
236                     &                           * ( pdept_1d(      ik  ) - pdepw_1d(ik) )   &
237                     &                           / ( pdepw_1d(      ik+1) - pdepw_1d(ik) )
238                  pe3t (ji,jj,ik) = pe3t_1d (ik) * ( pdepw   (ji,jj,ik+1) - pdepw_1d(ik) )   & 
239                     &                           / ( pdepw_1d(      ik+1) - pdepw_1d(ik) ) 
240                  pe3w (ji,jj,ik) = 0.5_wp * ( pdepw(ji,jj,ik+1) + pdepw_1d(ik+1) - 2._wp * pdepw_1d(ik) )   &
241                     &                     * ( pe3w_1d(ik) / ( pdepw_1d(ik+1) - pdepw_1d(ik) ) )
242                  !       ... on ik+1
243                  pe3w (ji,jj,ik+1) = pe3t (ji,jj,ik)
244                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik)
245                  pdept(ji,jj,ik+1) = pdept(ji,jj,ik) + pe3t(ji,jj,ik)
246            END DO
247         END DO
248         !
249      !
250      !
251      !                                      ! bottom scale factors and depth at  U-, V-, UW and VW-points
252      !
253      DO jk = 1,jpk                          ! Computed as the minimum of neighbooring scale factors
254         DO jj = 1, jpjm1
255            DO ji = 1, jpim1
256               pe3u (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji+1,jj,jk) )
257               pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) )
258               pe3uw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji+1,jj,jk) )
259               pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) )
260            END DO
261         END DO
262      END DO
263      !                                      ! lateral boundary conditions
264      CALL lbc_lnk( pe3u , 'U', 1._wp )   ;   CALL lbc_lnk( pe3uw, 'U', 1._wp )   
265      CALL lbc_lnk( pe3v , 'V', 1._wp )   ;   CALL lbc_lnk( pe3vw, 'V', 1._wp )
266      !
267
268      DO jk = 1, jpk                         ! set to z-scale factor if zero (i.e. along closed boundaries)
269         WHERE( pe3u (:,:,jk) == 0._wp )   pe3u (:,:,jk) = pe3t_1d(jk)
270         WHERE( pe3v (:,:,jk) == 0._wp )   pe3v (:,:,jk) = pe3t_1d(jk)
271         WHERE( pe3uw(:,:,jk) == 0._wp )   pe3uw(:,:,jk) = pe3w_1d(jk)
272         WHERE( pe3vw(:,:,jk) == 0._wp )   pe3vw(:,:,jk) = pe3w_1d(jk)
273      END DO
274     
275      ! Scale factor at F-point
276      DO jk = 1, jpk                        ! initialisation to z-scale factors
277         pe3f(:,:,jk) = pe3t_1d(jk)
278      END DO
279      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
280         DO jj = 1, jpjm1
281            DO ji = 1, fs_jpim1   ! vector opt.
282               pe3f(ji,jj,jk) = MIN( pe3v(ji,jj,jk), pe3v(ji+1,jj,jk) )
283            END DO
284         END DO
285      END DO
286      CALL lbc_lnk( pe3f, 'F', 1._wp )       ! Lateral boundary conditions
287      !
288      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
289         WHERE( pe3f(:,:,jk) == 0._wp )   pe3f(:,:,jk) = pe3t_1d(jk)
290      END DO
291!!gm  bug ? :  must be a do loop with mj0,mj1
292      !
293     
294!!gm  DO it differently !
295!      pe3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
296!      pe3w(:,mj0(1),:) = e3w(:,mj0(2),:)
297!      pe3u(:,mj0(1),:) = e3u(:,mj0(2),:)
298!      pe3v(:,mj0(1),:) = e3v(:,mj0(2),:)
299!      pe3f(:,mj0(1),:) = e3f(:,mj0(2),:)
300
301      ! Control of the sign
302      IF( MINVAL( pe3t (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
303      IF( MINVAL( pe3w (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
304      IF( MINVAL( pdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
305      IF( MINVAL( pdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
306     
307      ! Compute gde3w_0 (vertical sum of e3w)
308!      IF ( ln_isfcav ) THEN ! if cavity
309!         WHERE( misfdep == 0 )   misfdep = 1
310!         DO jj = 1,jpj
311!            DO ji = 1,jpi
312!               gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
313!               DO jk = 2, misfdep(ji,jj)
314!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
315!               END DO
316!               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
317!               DO jk = misfdep(ji,jj) + 1, jpk
318!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
319!               END DO
320!            END DO
321!         END DO
322!      ELSE ! no cavity
323!         gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
324!         DO jk = 2, jpk
325!            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
326!         END DO
327!      END IF
328      !
329      !
330     
331     
332     
333     
334      ENDIF
335      !
336   END SUBROUTINE usr_def_zgr
337
338
339   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
340      !!----------------------------------------------------------------------
341      !!                   ***  ROUTINE zgr_z1d  ***
342      !!
343      !! ** Purpose :   set the depth of model levels and the resulting
344      !!      vertical scale factors.
345      !!
346      !! ** Method  :   z-coordinate system (use in all type of coordinate)
347      !!      The depth of model levels is defined from an analytical
348      !!      function the derivative of which gives the scale factors.
349      !!      both depth and scale factors only depend on k (1d arrays).
350      !!              w-level: pdepw_1d  = pdep(k)
351      !!                       pe3w_1d(k) = dk(pdep)(k)     = e3(k)
352      !!              t-level: pdept_1d  = pdep(k+0.5)
353      !!                       pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5)
354      !!
355      !!            ===    Here constant vertical resolution   ===
356      !!
357      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
358      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
359      !!----------------------------------------------------------------------
360      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
361      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
362      !
363      INTEGER  ::   jk       ! dummy loop indices
364      REAL(wp) ::   zt, zw   ! local scalar
365      !!----------------------------------------------------------------------
366      !
367      IF(lwp) THEN                         ! Parameter print
368         WRITE(numout,*)
369         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
370         WRITE(numout,*) '    ~~~~~~~'
371      ENDIF
372      !
373      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
374      ! ----------------------
375      DO jk = 1, jpk
376         zw = REAL( jk , wp )
377         zt = REAL( jk , wp ) + 0.5_wp
378         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
379         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
380         pe3w_1d (jk) =    rn_dz
381         pe3t_1d (jk) =    rn_dz
382      END DO
383      !
384      IF(lwp) THEN                        ! control print
385         WRITE(numout,*)
386         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
387         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
388         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
389      ENDIF
390      !
391   END SUBROUTINE zgr_z1d
392   
393   !!======================================================================
394END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.