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

Last change on this file since 15271 was 15271, checked in by dbruciaferri, 17 months ago

rewriting saw tooth dignostic

File size: 13.8 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     mes_msk(:,:) = 1.
181     WHERE (s2z_msk(:,:).eq.0.0)  mes_msk(:,:) = 0.0
182
183     DO jj = 1,jpj
184        DO ji = 1,jpi
185           SELECT CASE (INT(s2z_msk(ji,jj)))
186             CASE (0) ! global zps area
187                  mbathy (ji,jj  ) = mbathy_in(ji,jj)
188                  gdept_0(ji,jj,:) = gdept_in(ji,jj,:)
189                  gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:)
190                  e3t_0  (ji,jj,:) = e3t_in (ji,jj,:)
191                  e3w_0  (ji,jj,:) = e3w_in (ji,jj,:)
192                  e3u_0  (ji,jj,:) = e3u_in (ji,jj,:)
193                  e3v_0  (ji,jj,:) = e3v_in (ji,jj,:)
194                  e3f_0  (ji,jj,:) = e3f_in (ji,jj,:)
195                  e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:)
196                  e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:)
197             CASE (1) ! MEs to zps transition area
198                  gdept_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdept_0(ji,jj,:) + &
199                    &                ( 1._wp - s2z_wgt(ji,jj) ) * gdept_in(ji,jj,:)
200                  gdepw_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdepw_0(ji,jj,:) + &
201                    &                ( 1._wp - s2z_wgt(ji,jj) ) * gdepw_in(ji,jj,:)
202             CASE (2) ! MEs area
203                  CYCLE     
204           END SELECT
205        END DO
206     END DO
207     !
208     ! e3t, e3w for transition zone
209     ! as finite differences
210     DO jj = 1, jpj
211        DO ji = 1, jpi
212           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
213              DO jk = 1,jpkm1
214                 e3t_0(ji,jj,jk)   = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk)
215                 e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk)
216              ENDDO
217              ! Surface
218              jk = 1
219              e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1))
220              !
221              ! Bottom
222              jk = jpk
223              e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk))
224           END IF
225        END DO
226     END DO
227     !
228     ! MBATHY transition zone
229     DO jj = 1, jpj
230        DO ji = 1, jpi
231           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
232              DO jk = 1, jpkm1
233                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
234                 IF( scobot(ji,jj) == 0.e0              ) mbathy(ji,jj) = 0
235                 IF( scobot(ji,jj) < 0.e0               ) mbathy(ji,jj) = INT( scobot(ji,jj)) ! do we need it?
236              END DO
237           END IF
238        END DO
239     END DO
240     !
241     ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0
242     ! for transition zone
243     !
244     DO jj = 1, jpjm1
245        DO ji = 1, jpim1
246           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
247              DO jk = 1, jpk
248                 e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
249                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) /         &
250                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
251
252                 e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
253                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) /         &
254                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
255
256                 e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
257                                   MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) /         &
258                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
259
260                 e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
261                                   MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) /         &
262                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
263
264                 e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)         +       &
265                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk)      +       &
266                                  MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) +       &
267                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) )    /       &
268                                  MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  &
269                                +   MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))  )
270              END DO
271           END IF
272        END DO
273     END DO
274     ! Computing gdep3w
275     gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
276     DO jk = 2, jpk
277        gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
278     END DO
279     !
280     ! From here equal to sco code - domzgr.F90 line 2079
281     ! Lateral B.C. we apply only for transition zone
282     ! since for s- and zps- areas has already been applied
283
284     CALL lbc_lnk( e3t_0 , 'T', 1._wp )
285     CALL lbc_lnk( e3u_0 , 'U', 1._wp )
286     CALL lbc_lnk( e3v_0 , 'V', 1._wp )
287     CALL lbc_lnk( e3f_0 , 'F', 1._wp )
288     CALL lbc_lnk( e3w_0 , 'W', 1._wp )
289     CALL lbc_lnk( e3uw_0, 'U', 1._wp )
290     CALL lbc_lnk( e3vw_0, 'V', 1._wp )
291
292     WHERE (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
293     WHERE (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
294     WHERE (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
295     WHERE (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
296     WHERE (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
297     WHERE (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
298     WHERE (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
299
300     IF( lwp ) WRITE(numout,*) 'Refine mbathy'
301     DO jj = 1, jpj
302        DO ji = 1, jpi
303           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
304              DO jk = 1, jpkm1
305                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
306              END DO
307              IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
308           END IF
309        END DO
310     END DO
311
312   END SUBROUTINE zgr_mes_local
313
314END MODULE zgrmes
Note: See TracBrowser for help on using the repository browser.