source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/tests/ISOMIP/MY_SRC/usrdef_zgr.F90 @ 13138

Last change on this file since 13138 was 13138, checked in by smasson, 4 months ago

Extra_Halo: minor bugfixes and cleaning, see #2366

  • Property svn:keywords set to Id
File size: 11.7 KB
RevLine 
[7715]1MODULE usrdef_zgr
2   !!======================================================================
3   !!                     ***  MODULE usrdef_zgr  ***
4   !!
5   !!                       ===  ISOMIP 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   !!                 ! 2017-02  (P. Mathiot, S. Flavoni)  Adapt code to ISOMIP case
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   usr_def_zgr   : user defined vertical coordinate system (required)
15   !!       zgr_z1d   : reference 1D z-coordinate
16   !!---------------------------------------------------------------------
17   USE oce            ! ocean variables
[12939]18   USE dom_oce ,  ONLY: mj0   , mj1    ! ocean space and time domain
19   USE dom_oce ,  ONLY: glamt , gphit  ! ocean space and time domain
[7715]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 timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   usr_def_zgr   ! called by domzgr.F90
31
[12939]32   !! * Substitutions
33#  include "do_loop_substitute.h90"
[7715]34   !!----------------------------------------------------------------------
[10073]35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[7715]36   !! $Id$
[10073]37   !! Software governed by the CeCILL license (see ./LICENSE)
[7715]38   !!----------------------------------------------------------------------
39CONTAINS             
40
41   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
42      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
43      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
44      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
45      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
46      &                    k_top  , k_bot    )                             ! top & bottom ocean level
47      !!---------------------------------------------------------------------
48      !!              ***  ROUTINE usr_def_zgr  ***
49      !!
50      !! ** Purpose :   User defined the vertical coordinates
51      !!
52      !!----------------------------------------------------------------------
[9135]53      LOGICAL                   , INTENT(in   ) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags ( read in namusr_def )
54      LOGICAL                   , INTENT(  out) ::   ld_isfcav                   ! under iceshelf cavity flag
55      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
56      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
57      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth        [m]
58      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
59      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
60      INTEGER , DIMENSION(:,:)  , INTENT(  out) ::   k_top, k_bot                ! first & last ocean level
[7715]61      !
62      INTEGER  ::   ji , jj, jk       ! dummy indices
63      INTEGER  ::   ij0, ij1          ! dummy indices 
64      INTEGER  ::   ik                ! local integers
65      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
66      REAL(wp) ::   ze3min, zdepth    ! local scalar
67      REAL(wp), DIMENSION(jpi,jpj) ::   zht  , zhu         ! bottom depth
68      REAL(wp), DIMENSION(jpi,jpj) ::   zhisf, zhisfu      ! top depth
69      !!----------------------------------------------------------------------
70      !
71      IF(lwp) WRITE(numout,*)
72      IF(lwp) WRITE(numout,*) 'usr_def_zgr : ISOMIP configuration (z(ps)- or s-coordinate closed box ocean without cavities)'
73      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
74      !
75      !
76      ! type of vertical coordinate
77      ! ---------------------------
[8018]78      ! set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
79      ld_isfcav = .TRUE.       ! ISF Ice Shelves Flag
[7715]80      !
81      !
82      ! Build the vertical coordinate system
83      ! ------------------------------------
84      !
85      !                       !==  isfdraft  ==!
86      !
87      zht  (:,:) = rbathy 
88      zhisf(:,:) = 200._wp
[13138]89      ij0 = 1   ;   ij1 = 40+nn_hls
[7715]90      DO jj = mj0(ij0), mj1(ij1)
91         zhisf(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
92      END DO
93      !
94      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
95      !
96      !                       !==  top masked level bathymetry  ==!  (all coordinates)
97      !
[9135]98      IF ( ld_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps)
[7715]99         !
100         ze3min = 0.1_wp * rn_e3
101         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min
102         !
103         !                                !* bottom ocean compute from the depth of grid-points
104         k_bot(:,:) = jpkm1
105         DO jk = jpkm1, 1, -1
106            WHERE( zht(:,:) < pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
107         END DO
108         !                                !* top ocean compute from the depth of grid-points
109         k_top(:,:) = 1                   !
110         DO jk = 2, jpkm1
111            zdepth = pdepw_1d(jk+1) - ze3min
112            WHERE( zhisf(:,:) > 0.0 .AND. zhisf(:,:) >= zdepth )   k_top(:,:) = (jk + 1) 
113         END DO
114         !
115         !                                   !* vertical coordinate system
116         DO jk = 1, jpk                      ! initialization to the reference z-coordinate
117            pdept(:,:,jk) = pdept_1d(jk)
118            pdepw(:,:,jk) = pdepw_1d(jk)
119            pe3t (:,:,jk) = pe3t_1d (jk)
120            pe3u (:,:,jk) = pe3t_1d (jk)
121            pe3v (:,:,jk) = pe3t_1d (jk)
122            pe3f (:,:,jk) = pe3t_1d (jk)
123            pe3w (:,:,jk) = pe3w_1d (jk)
124            pe3uw(:,:,jk) = pe3w_1d (jk)
125            pe3vw(:,:,jk) = pe3w_1d (jk)
126         END DO
[12939]127         ! top scale factors and depth at T- and W-points
128         DO_2D_11_11
129            ik = k_top(ji,jj)
130            IF ( ik > 2 ) THEN
131               ! pdeptw at the interface
132               pdepw(ji,jj,ik  ) = MAX( zhisf(ji,jj) , pdepw(ji,jj,ik) )
133               ! e3t in both side of the interface
[7715]134               pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
[12939]135               ! pdept in both side of the interface (from previous e3t)
[7715]136               pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
[12939]137               pdept(ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pe3t (ji,jj,ik  ) * 0.5_wp
138               ! pe3w on both side of the interface
139               pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik  )
140               pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1)
141               ! e3t into the ice shelf
142               pe3t (ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pdepw(ji,jj,ik-1)
143               pe3w (ji,jj,ik-1) = pdept(ji,jj,ik-1) - pdept(ji,jj,ik-2)
144            END IF
145         END_2D
146         ! bottom scale factors and depth at T- and W-points
147         DO_2D_11_11
148            ik = k_bot(ji,jj)
149            pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
150            pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
151            pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
152            !
153            pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
154            pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
155            pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)
156         END_2D       
[7715]157         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points
[9089]158         pe3u (:,:,:) = pe3t(:,:,:)
159         pe3uw(:,:,:) = pe3w(:,:,:)
[12939]160         DO_3D_00_00( 1, jpk )
161         !                                   ! Computed as the minimum of neighbooring scale factors
162            pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) )
163            pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) )
164            pe3f (ji,jj,jk) = pe3v(ji,jj,jk)
165         END_3D
[10425]166         CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1._wp )   ;   CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1._wp )
167         CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1._wp )
[9089]168         DO jk = 1, jpk
[7715]169            ! set to z-scale factor if zero (i.e. along closed boundaries) because of lbclnk
170            WHERE( pe3u (:,:,jk) == 0._wp )   pe3u (:,:,jk) = pe3t_1d(jk)
171            WHERE( pe3v (:,:,jk) == 0._wp )   pe3v (:,:,jk) = pe3t_1d(jk)
[9089]172            WHERE( pe3f (:,:,jk) == 0._wp )   pe3f (:,:,jk) = pe3t_1d(jk)
[7715]173            WHERE( pe3uw(:,:,jk) == 0._wp )   pe3uw(:,:,jk) = pe3w_1d(jk)
174            WHERE( pe3vw(:,:,jk) == 0._wp )   pe3vw(:,:,jk) = pe3w_1d(jk)
175         END DO
176         !
177      ENDIF
178      !
179   END SUBROUTINE usr_def_zgr
180
181
182   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
183      !!----------------------------------------------------------------------
184      !!                   ***  ROUTINE zgr_z1d  ***
185      !!
186      !! ** Purpose :   set the depth of model levels and the resulting
187      !!      vertical scale factors.
188      !!
189      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
190      !!       The depth of model levels is set from dep(k), an analytical function:
191      !!                   w-level: depw_1d  = dep(k)
192      !!                   t-level: dept_1d  = dep(k+0.5)
193      !!       The scale factors are the discrete derivative of the depth:
194      !!                   e3w_1d(jk) = dk[ dept_1d ]
195      !!                   e3t_1d(jk) = dk[ depw_1d ]
196      !!
197      !!            ===    Here constant vertical resolution   ===
198      !!
199      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
200      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
201      !!----------------------------------------------------------------------
202      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
203      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
204      !
205      INTEGER  ::   jk       ! dummy loop indices
206      REAL(wp) ::   zt, zw   ! local scalar
207      !!----------------------------------------------------------------------
208      !
209      IF(lwp) THEN                         ! Parameter print
210         WRITE(numout,*)
211         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_e3
212         WRITE(numout,*) '    ~~~~~~~'
213      ENDIF
214      !
215      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
216      ! ----------------------
217      DO jk = 1, jpk
218         zw = REAL( jk , wp )
219         zt = REAL( jk , wp ) + 0.5_wp
220         pdepw_1d(jk) =    rn_e3 *   REAL( jk-1 , wp )
221         pdept_1d(jk) =    rn_e3 * ( REAL( jk-1 , wp ) + 0.5_wp )
222         pe3w_1d (jk) =    rn_e3
223         pe3t_1d (jk) =    rn_e3
224      END DO
225      !
226      IF(lwp) THEN                        ! control print
227         WRITE(numout,*)
228         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
229         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
230         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
231      ENDIF
232      !
233   END SUBROUTINE zgr_z1d
234   
235   !!======================================================================
236END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.