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/2017/dev_r7663_ISOMIP_TEST_CASE/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/MY_SRC – NEMO

source: branches/2017/dev_r7663_ISOMIP_TEST_CASE/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/MY_SRC/usrdef_zgr.F90 @ 7668

Last change on this file since 7668 was 7668, checked in by mathiot, 7 years ago

Change location of ISOMIP, add plot routine and expected melt/moc/psi figure and README
File size: 12.7 KB
Line 
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
18   USE dom_oce ,  ONLY: ln_zco, ln_zps, ln_sco          ! ocean space and time domain
19   USE dom_oce ,  ONLY: mj0   , mj1   , nimpp , njmpp   ! ocean space and time domain
20   USE dom_oce ,  ONLY: glamt , gphit                   ! 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$
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      INTEGER  ::   ij0, ij1          ! dummy indices 
66      INTEGER  ::   ik                ! local integers
67      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
68      REAL(wp) ::   ze3min, zdepth    ! local scalar
69      REAL(wp), DIMENSION(jpi,jpj) ::   zht  , zhu         ! bottom depth
70      REAL(wp), DIMENSION(jpi,jpj) ::   zhisf, zhisfu      ! top depth
71      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk 
72      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                ! 2d workspace
73      !!----------------------------------------------------------------------
74      !
75      IF(lwp) WRITE(numout,*)
76      IF(lwp) WRITE(numout,*) 'usr_def_zgr : ISOMIP configuration (z(ps)- or s-coordinate closed box ocean without cavities)'
77      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
78      !
79      !
80      ! type of vertical coordinate
81      ! ---------------------------
82      ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
83      ld_isfcav = .TRUE.
84      !
85      !
86      ! Build the vertical coordinate system
87      ! ------------------------------------
88      !
89      !                       !==  isfdraft  ==!
90      !
91      ! the ocean basin surrounded by land (1 grid-point) set through lbc_lnk call as jperio=0
92      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
93      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
94      zmsk(:,:) = NINT( z2d(:,:) )
95      !
96      !
97      zht  (:,:) = rbathy 
98      zhisf(:,:) = 200._wp
99      ij0 = 1 ; ij1 = 40
100      DO jj = mj0(ij0), mj1(ij1)
101         zhisf(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
102      END DO
103      zhisf(:,:) = zhisf(:,:) * zmsk(:,:)
104      !
105      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
106      !
107      !                       !==  top masked level bathymetry  ==!  (all coordinates)
108      !
109      IF ( ln_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps)
110         !
111         ze3min = 0.1_wp * rn_e3
112         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min
113         !
114         !                                !* bottom ocean compute from the depth of grid-points
115         k_bot(:,:) = jpkm1
116         DO jk = jpkm1, 1, -1
117            WHERE( zht(:,:) < pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
118         END DO
119         !                                !* top ocean compute from the depth of grid-points
120         k_top(:,:) = 1                   !
121         DO jk = 2, jpkm1
122            zdepth = pdepw_1d(jk+1) - ze3min
123            WHERE( zhisf(:,:) > 0.0 .AND. zhisf(:,:) >= zdepth )   k_top(:,:) = (jk + 1) 
124         END DO
125         !
126         !                                   !* vertical coordinate system
127         DO jk = 1, jpk                      ! initialization to the reference z-coordinate
128            pdept(:,:,jk) = pdept_1d(jk)
129            pdepw(:,:,jk) = pdepw_1d(jk)
130            pe3t (:,:,jk) = pe3t_1d (jk)
131            pe3u (:,:,jk) = pe3t_1d (jk)
132            pe3v (:,:,jk) = pe3t_1d (jk)
133            pe3f (:,:,jk) = pe3t_1d (jk)
134            pe3w (:,:,jk) = pe3w_1d (jk)
135            pe3uw(:,:,jk) = pe3w_1d (jk)
136            pe3vw(:,:,jk) = pe3w_1d (jk)
137         END DO
138         DO jj = 1, jpj                      ! top scale factors and depth at T- and W-points
139            DO ji = 1, jpi
140               ik = k_top(ji,jj)
141               IF ( ik > 2 ) THEN
142                  ! pdeptw at the interface
143                  pdepw(ji,jj,ik  ) = MAX( zhisf(ji,jj) , pdepw(ji,jj,ik) )
144                  ! e3t in both side of the interface
145                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
146                  ! pdept in both side of the interface (from previous e3t)
147                  pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
148                  pdept(ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pe3t (ji,jj,ik  ) * 0.5_wp
149                  ! pe3w on both side of the interface
150                  pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik  )
151                  pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1)
152                  ! e3t into the ice shelf
153                  pe3t (ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pdepw(ji,jj,ik-1)
154               END IF
155            END DO
156         END DO         
157         DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points
158            DO ji = 1, jpi
159               ik = k_bot(ji,jj)
160               pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
161               pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
162               pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
163               !
164               pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
165               pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
166               pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)
167            END DO
168         END DO         
169         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points
170         DO jk = 1,jpkm1                     ! Computed as the minimum of neighbooring scale factors
171            DO jj = 1, jpjm1
172               DO ji = 1, fs_jpim1   ! vector opt.
173                  pe3u (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji+1,jj,jk) )
174                  pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) )
175                  pe3uw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji+1,jj,jk) )
176                  pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) )
177               END DO
178            END DO
179            DO jj = 1, jpjm1
180               DO ji = 1, fs_jpim1   ! vector opt.
181                 pe3f(ji,jj,jk) = MIN( pe3v(ji,jj,jk), pe3v(ji+1,jj,jk) )
182               END DO
183            END DO
184         END DO
185         CALL lbc_lnk( pe3u , 'U', 1._wp )   ;   CALL lbc_lnk( pe3uw, 'U', 1._wp )   ! lateral boundary conditions
186         CALL lbc_lnk( pe3v , 'V', 1._wp )   ;   CALL lbc_lnk( pe3vw, 'V', 1._wp )
187         CALL lbc_lnk( pe3f , 'F', 1._wp )
188         DO jk = 1,jpkm1
189            ! set to z-scale factor if zero (i.e. along closed boundaries) because of lbclnk
190            WHERE( pe3u (:,:,jk) == 0._wp )   pe3u (:,:,jk) = pe3t_1d(jk)
191            WHERE( pe3v (:,:,jk) == 0._wp )   pe3v (:,:,jk) = pe3t_1d(jk)
192            WHERE( pe3uw(:,:,jk) == 0._wp )   pe3uw(:,:,jk) = pe3w_1d(jk)
193            WHERE( pe3vw(:,:,jk) == 0._wp )   pe3vw(:,:,jk) = pe3w_1d(jk)
194         END DO
195         !
196      ENDIF
197      !
198   END SUBROUTINE usr_def_zgr
199
200
201   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
202      !!----------------------------------------------------------------------
203      !!                   ***  ROUTINE zgr_z1d  ***
204      !!
205      !! ** Purpose :   set the depth of model levels and the resulting
206      !!      vertical scale factors.
207      !!
208      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
209      !!       The depth of model levels is set from dep(k), an analytical function:
210      !!                   w-level: depw_1d  = dep(k)
211      !!                   t-level: dept_1d  = dep(k+0.5)
212      !!       The scale factors are the discrete derivative of the depth:
213      !!                   e3w_1d(jk) = dk[ dept_1d ]
214      !!                   e3t_1d(jk) = dk[ depw_1d ]
215      !!
216      !!            ===    Here constant vertical resolution   ===
217      !!
218      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
219      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
220      !!----------------------------------------------------------------------
221      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
222      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
223      !
224      INTEGER  ::   jk       ! dummy loop indices
225      REAL(wp) ::   zt, zw   ! local scalar
226      !!----------------------------------------------------------------------
227      !
228      IF(lwp) THEN                         ! Parameter print
229         WRITE(numout,*)
230         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_e3
231         WRITE(numout,*) '    ~~~~~~~'
232      ENDIF
233      !
234      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
235      ! ----------------------
236      DO jk = 1, jpk
237         zw = REAL( jk , wp )
238         zt = REAL( jk , wp ) + 0.5_wp
239         pdepw_1d(jk) =    rn_e3 *   REAL( jk-1 , wp )
240         pdept_1d(jk) =    rn_e3 * ( REAL( jk-1 , wp ) + 0.5_wp )
241         pe3w_1d (jk) =    rn_e3
242         pe3t_1d (jk) =    rn_e3
243      END DO
244      !
245      IF(lwp) THEN                        ! control print
246         WRITE(numout,*)
247         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
248         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
249         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
250      ENDIF
251      !
252   END SUBROUTINE zgr_z1d
253   
254   !!======================================================================
255END MODULE usrdef_zgr
Note: See TracBrowser for help on using the repository browser.