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.
zgrmes.F90 in NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src – NEMO

source: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90 @ 15129

Last change on this file since 15129 was 15129, checked in by dbruciaferri, 3 years ago

remove rn_ebot_{max,min} from namelist
File size: 13.7 KB
Line 
1MODULE zgrmes
2   !!==============================================================================
3   !!                       ***  MODULE zgrmes   ***
4   !! Ocean initialization : Multiple Enveloped s coordinate (MES)
5   !!==============================================================================
6   !!  NEMO      4.0  ! 2021-07  (D. Bruciaferri)   
7   !!----------------------------------------------------------------------
8   !
9   USE oce               ! ocean variables
10   USE dom_oce           ! ocean domain
11   USE depth_e3          ! depth <=> e3
12   USE closea            ! closed seas
13   USE mes               ! MEs-coordinates
14   !
15   USE in_out_manager    ! I/O manager
16   USE iom               ! I/O library
17   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
18   USE lib_mpp           ! distributed memory computing library
19   USE wrk_nemo          ! Memory allocation
20   USE timing            ! Timing
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   zgr_mes        ! called by domzgr.F90
26   ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.)
27   !
28   REAL(wp), POINTER, DIMENSION(:)     :: e3t_1d_in  , e3w_1d_in   !: ref. scale factors (m)
29   REAL(wp), POINTER, DIMENSION(:)     :: gdept_1d_in, gdepw_1d_in !: ref. depth (m)
30   !
31   REAL(wp), POINTER, DIMENSION(:,:)   :: mbathy_in 
32   !
33   REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m]
34   REAL(wp), POINTER, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in         !:   -      -
35   !
36   REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_in, gdepw_in               !: depths [m]
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40
41CONTAINS
42
43! =====================================================================================================
44
45   SUBROUTINE zgr_mes
46      !!---------------------------------------------------------------------
47      !!              ***  ROUTINE zgr_mes  ***
48      !!
49      !! ** Purpose : Wrap the steps to generate gloabal or localised
50      !!              MEs-coordinates systems
51      !!----------------------------------------------------------------------
52      !
53      ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.)
54      !
55      !-----------------------------------------------------------------------
56      !
57      ! Generating a global MEs vertical grid
58      CALL mes_build
59
60      ! Local MEs
61      IF( ln_loc_mes ) THEN
62         !
63         CALL wrk_alloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in)
64         CALL wrk_alloc( jpi, jpj, mbathy_in)
65         CALL wrk_alloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in)
66         CALL wrk_alloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in)
67         !
68         ! Reading the input vertical grid that will be used globally
69         CALL zgr_dom_read
70         ! Creating the local MEs within the global input vertical grid
71         CALL zgr_mes_local
72         !
73         CALL wrk_dealloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in)
74         CALL wrk_dealloc( jpi, jpj, mbathy_in)
75         CALL wrk_dealloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in)
76         CALL wrk_dealloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in)
77         !
78      END IF
79
80   END SUBROUTINE zgr_mes
81
82! =====================================================================================================
83
84   SUBROUTINE zgr_dom_read
85     !!---------------------------------------------------------------------
86     !!              ***  ROUTINE zgr_dom_read  ***
87     !!
88     !! ** Purpose :   Read the vertical information in the domain configuration file
89     !!
90     !!----------------------------------------------------------------------
91     INTEGER  ::   jk     ! dummy loop index
92     INTEGER  ::   inum   ! local logical unit
93     REAL(WP) ::   z_cav
94     REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
95     !!----------------------------------------------------------------------
96     !
97     IF(lwp) THEN
98       WRITE(numout,*)
99       WRITE(numout,*) '   zgr_dom_read : read the vertical coordinates in domain_cfg_global.nc file'
100       WRITE(numout,*) '   ~~~~~~~~'
101     ENDIF
102     !
103     CALL iom_open( 'domain_cfg_global.nc', inum )
104     !
105     !                          !* ocean cavities under iceshelves
106     CALL iom_get( inum, 'ln_isfcav', z_cav )
107     IF( z_cav == 0._wp ) THEN   ;   ln_isfcav = .false.   ;   ELSE   ;   ln_isfcav = .true.   ;   ENDIF
108     !
109     !                          !* vertical scale factors
110     CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d_in ) ! 1D reference coordinate
111     CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d_in )
112     !
113     CALL iom_get( inum, jpdom_global, 'bottom_level' , mbathy_in )     ! 2D mbathy
114     !
115     CALL iom_get( inum, jpdom_global, 'e3t_0'  , e3t_in  )     ! 3D coordinate
116     CALL iom_get( inum, jpdom_global, 'e3u_0'  , e3u_in  )
117     CALL iom_get( inum, jpdom_global, 'e3v_0'  , e3v_in  )
118     CALL iom_get( inum, jpdom_global, 'e3f_0'  , e3f_in  )
119     CALL iom_get( inum, jpdom_global, 'e3w_0'  , e3w_in  )
120     CALL iom_get( inum, jpdom_global, 'e3uw_0' , e3uw_in )
121     CALL iom_get( inum, jpdom_global, 'e3vw_0' , e3vw_in )
122     !
123     !                          !* depths
124     !                                   !- old depth definition (obsolescent feature)
125     IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
126        & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
127        & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
128        & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
129        CALL ctl_warn( 'zgr_dom_read : old definition of depths and scale factors used ', & 
130           &           '           depths at t- and w-points read in the domain configuration file')
131        CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d_in )   
132        CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d_in )
133        CALL iom_get( inum, jpdom_global , 'gdept_0' , gdept_in )
134        CALL iom_get( inum, jpdom_global , 'gdepw_0' , gdepw_in )
135        !
136     ELSE                                !- depths computed from e3. scale factors
137        CALL e3_to_depth( e3t_1d_in, e3w_1d_in, gdept_1d_in, gdepw_1d_in )    ! 1D reference depth
138        CALL e3_to_depth( e3t_in   , e3w_in   , gdept_in   , gdepw_in    )    ! 3D depths
139        IF(lwp) THEN
140          WRITE(numout,*)
141          WRITE(numout,*) '              GLOBAL reference 1D z-coordinate depth and scale factors:'
142          WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
143          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d_in(jk), gdepw_1d_in(jk), e3t_1d_in(jk), e3w_1d_in(jk), jk = 1, jpk )
144        ENDIF
145     ENDIF
146     !
147     CALL iom_close( inum )
148     !
149   END SUBROUTINE zgr_dom_read
150
151! =====================================================================================================
152
153   SUBROUTINE zgr_mes_local
154     !!---------------------------------------------------------------------
155     !!              ***  ROUTINE zgr_dom_merge  ***
156     !!
157     !! ** Purpose :   Create a local MEs grid within a gloabal grid
158     !!                using different vertical coordinates.
159     !!
160     !!----------------------------------------------------------------------
161     INTEGER                      ::   ji, jj, jk ! dummy loop index
162     INTEGER                      ::   inum       ! local logical unit
163     REAL(wp), DIMENSION(jpi,jpj) ::   s2z_msk    ! mask for MEs-area (=2),
164                                                  !          transition area (=1),
165                                                  !          z-area (=0)
166     REAL(wp), DIMENSION(jpi,jpj) ::   s2z_wgt    ! weigths for computing model
167                                                  ! levels in transition area
168     !!----------------------------------------------------------------------
169
170     IF(lwp) THEN
171       WRITE(numout,*)
172       WRITE(numout,*) '   zgr_mes_merge : read the vertical coordinates in domain_cfg_global.nc file'
173       WRITE(numout,*) '   ~~~~~~~~'
174     ENDIF
175     !
176     CALL iom_open( 'bathy_meter.nc', inum )
177     CALL iom_get( inum, jpdom_data, 's2z_msk', s2z_msk)
178     CALL iom_get( inum, jpdom_data, 's2z_wgt', s2z_wgt)
179
180     DO jj = 1,jpj
181        DO ji = 1,jpi
182           SELECT CASE (INT(s2z_msk(ji,jj)))
183             CASE (0) ! global zps area
184                  mbathy (ji,jj  ) = mbathy_in(ji,jj)
185                  gdept_0(ji,jj,:) = gdept_in(ji,jj,:)
186                  gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:)
187                  e3t_0  (ji,jj,:) = e3t_in (ji,jj,:)
188                  e3w_0  (ji,jj,:) = e3w_in (ji,jj,:)
189                  e3u_0  (ji,jj,:) = e3u_in (ji,jj,:)
190                  e3v_0  (ji,jj,:) = e3v_in (ji,jj,:)
191                  e3f_0  (ji,jj,:) = e3f_in (ji,jj,:)
192                  e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:)
193                  e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:)
194             CASE (1) ! MEs to zps transition area
195                  gdept_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdept_0(ji,jj,:) + &
196                    &                ( 1._wp - s2z_wgt(ji,jj) ) * gdept_in(ji,jj,:)
197                  gdepw_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdepw_0(ji,jj,:) + &
198                    &                ( 1._wp - s2z_wgt(ji,jj) ) * gdepw_in(ji,jj,:)
199             CASE (2) ! MEs area
200                  CYCLE     
201           END SELECT
202        END DO
203     END DO
204     !
205     ! e3t, e3w for transition zone
206     ! as finite differences
207     DO jj = 1, jpj
208        DO ji = 1, jpi
209           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
210              DO jk = 1,jpkm1
211                 e3t_0(ji,jj,jk)   = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk)
212                 e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk)
213              ENDDO
214              ! Surface
215              jk = 1
216              e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1))
217              !
218              ! Bottom
219              jk = jpk
220              e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk))
221           END IF
222        END DO
223     END DO
224     !
225     ! MBATHY transition zone
226     DO jj = 1, jpj
227        DO ji = 1, jpi
228           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
229              DO jk = 1, jpkm1
230                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
231                 IF( scobot(ji,jj) == 0.e0              ) mbathy(ji,jj) = 0
232                 IF( scobot(ji,jj) < 0.e0               ) mbathy(ji,jj) = INT( scobot(ji,jj)) ! do we need it?
233              END DO
234           END IF
235        END DO
236     END DO
237     !
238     ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0
239     ! for transition zone
240     !
241     DO jj = 1, jpjm1
242        DO ji = 1, jpim1
243           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
244              DO jk = 1, jpk
245                 e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
246                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) /         &
247                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
248
249                 e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
250                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) /         &
251                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
252
253                 e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
254                                   MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) /         &
255                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
256
257                 e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
258                                   MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) /         &
259                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
260
261                 e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)         +       &
262                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk)      +       &
263                                  MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) +       &
264                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) )    /       &
265                                  MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  &
266                                +   MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))  )
267              END DO
268           END IF
269        END DO
270     END DO
271     ! Computing gdep3w
272     gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
273     DO jk = 2, jpk
274        gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
275     END DO
276     !
277     ! From here equal to sco code - domzgr.F90 line 2079
278     ! Lateral B.C. we apply only for transition zone
279     ! since for s- and zps- areas has already been applied
280
281     CALL lbc_lnk( e3t_0 , 'T', 1._wp )
282     CALL lbc_lnk( e3u_0 , 'U', 1._wp )
283     CALL lbc_lnk( e3v_0 , 'V', 1._wp )
284     CALL lbc_lnk( e3f_0 , 'F', 1._wp )
285     CALL lbc_lnk( e3w_0 , 'W', 1._wp )
286     CALL lbc_lnk( e3uw_0, 'U', 1._wp )
287     CALL lbc_lnk( e3vw_0, 'V', 1._wp )
288
289     WHERE (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
290     WHERE (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
291     WHERE (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
292     WHERE (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
293     WHERE (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
294     WHERE (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
295     WHERE (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
296
297     IF( lwp ) WRITE(numout,*) 'Refine mbathy'
298     DO jj = 1, jpj
299        DO ji = 1, jpi
300           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
301              DO jk = 1, jpkm1
302                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
303              END DO
304              IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
305           END IF
306        END DO
307     END DO
308
309   END SUBROUTINE zgr_mes_local
310
311END MODULE zgrmes
Note: See TracBrowser for help on using the repository browser.