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.
domrea.F90 in trunk/NEMOGCM/NEMO/OFF_SRC – NEMO

source: trunk/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 4569

Last change on this file since 4569 was 4569, checked in by vichi, 10 years ago

Correct problem with reading mesh_mask in OFF_SRC, Ticket #1253

The reading of mesh mask information is now compliant with the
grid variables written in mesh_mask.nc.
Be aware that if one is using a mesh_mask generated with a previous
version of NEMO (<3.6) variable e3t must be renamed to e3t_0.
This commit also adds information on the date which may be needed
by other biogeochemical models.
  • Property svn:keywords set to Id
File size: 18.9 KB
Line 
1MODULE domrea
2   !!======================================================================
3   !!                       ***  MODULE domrea  ***
4   !! Ocean initialization : read the ocean domain meshmask file(s)
5   !!======================================================================
6   !! History :  3.3  ! 2010-05  (C. Ethe)  Full reorganization of the off-line
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dom_rea        : read mesh and mask file(s)
11   !!                    nmsh = 1  :   mesh_mask file
12   !!                         = 2  :   mesh and mask file
13   !!                         = 3  :   mesh_hgr, mesh_zgr and mask
14   !!----------------------------------------------------------------------
15   USE dom_oce         ! ocean space and time domain
16   USE dommsk          ! domain: masks
17   USE lbclnk          ! lateral boundary condition - MPP exchanges
18   USE trc_oce         ! shared ocean/biogeochemical variables
19   USE lib_mpp 
20   USE in_out_manager
21   USE wrk_nemo 
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   dom_rea    ! routine called by inidom.F90
27  !! * Substitutions
28#  include "domzgr_substitute.h90"
29   !!----------------------------------------------------------------------
30   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE dom_rea
37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE dom_rea  ***
39      !!                   
40      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
41      !!      ocean domain informations (mesh and mask arrays). This (these)
42      !!      file(s) is (are) used for visualisation (SAXO software) and
43      !!      diagnostic computation.
44      !!
45      !! ** Method  :   Read in a file all the arrays generated in routines
46      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
47      !!      the vertical coord. used (z-coord, partial steps, s-coord)
48      !!                    nmsh = 1  :   'mesh_mask.nc' file
49      !!                         = 2  :   'mesh.nc' and mask.nc' files
50      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
51      !!                                  'mask.nc' files
52      !!      For huge size domain, use option 2 or 3 depending on your
53      !!      vertical coordinate.
54      !!
55      !! ** input file :
56      !!      meshmask.nc  : domain size, horizontal grid-point position,
57      !!                     masks, depth and vertical scale factors
58      !!----------------------------------------------------------------------
59      USE iom
60      !!
61      INTEGER  ::   ji, jj, jk   ! dummy loop indices
62      INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
63      REAL(wp) ::   zrefdep         ! local real
64      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw
65      !!----------------------------------------------------------------------
66
67      IF(lwp) WRITE(numout,*)
68      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
69      IF(lwp) WRITE(numout,*) '~~~~~~~'
70
71      CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw )
72
73      zmbk(:,:) = 0._wp
74
75      SELECT CASE (nmsh)
76         !                                     ! ============================
77         CASE ( 1 )                            !  create 'mesh_mask.nc' file
78            !                                  ! ============================
79
80            IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
81            CALL iom_open( 'mesh_mask', inum0 )
82
83            inum2 = inum0                                            ! put all the informations
84            inum3 = inum0                                            ! in unit inum0
85            inum4 = inum0
86
87            !                                  ! ============================
88         CASE ( 2 )                            !  create 'mesh.nc' and
89            !                                  !         'mask.nc' files
90            !                                  ! ============================
91
92            IF(lwp) WRITE(numout,*) '          two files in "mesh.nc" and "mask.nc" '
93            CALL iom_open( 'mesh', inum1 )
94            CALL iom_open( 'mask', inum2 )
95
96            inum3 = inum1                                            ! put mesh informations
97            inum4 = inum1                                            ! in unit inum1
98
99            !                                  ! ============================
100         CASE ( 3 )                            !  create 'mesh_hgr.nc'
101            !                                  !         'mesh_zgr.nc' and
102            !                                  !         'mask.nc'     files
103            !                                  ! ============================
104
105            IF(lwp) WRITE(numout,*) '          three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" '
106            CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc'
107            CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc'
108            CALL iom_open( 'mask'    , inum2 ) ! create 'mask.nc'
109
110            !                                  ! ===========================
111         CASE DEFAULT                          ! return error
112            !                                  ! mesh has to be provided
113            !                                  ! ===========================
114            CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ',   &
115            &                                 ' Invalid nn_msh value in the namelist (0 is not allowed)' )
116
117         END SELECT
118
119         !                                                         ! masks (inum2)
120         CALL iom_get( inum2, jpdom_data, 'tmask', tmask )
121         CALL iom_get( inum2, jpdom_data, 'umask', umask )
122         CALL iom_get( inum2, jpdom_data, 'vmask', vmask )
123         CALL iom_get( inum2, jpdom_data, 'fmask', fmask )
124
125         CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
126         CALL lbc_lnk( umask, 'U', 1._wp )     
127         CALL lbc_lnk( vmask, 'V', 1._wp )
128         CALL lbc_lnk( fmask, 'F', 1._wp )
129
130#if defined key_c1d
131         ! set umask and vmask equal tmask in 1D configuration
132         IF(lwp) WRITE(numout,*)
133         IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
134         IF(lwp) WRITE(numout,*) '**********                                                     ********'
135
136         umask(:,:,:) = tmask(:,:,:)
137         vmask(:,:,:) = tmask(:,:,:)
138#endif
139
140#if defined key_degrad
141         CALL iom_get( inum2, jpdom_data, 'facvolt', facvol )
142#endif
143
144         !                                                         ! horizontal mesh (inum3)
145         CALL iom_get( inum3, jpdom_data, 'glamt', glamt )
146         CALL iom_get( inum3, jpdom_data, 'glamu', glamu )
147         CALL iom_get( inum3, jpdom_data, 'glamv', glamv )
148         CALL iom_get( inum3, jpdom_data, 'glamf', glamf )
149
150         CALL iom_get( inum3, jpdom_data, 'gphit', gphit )
151         CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu )
152         CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv )
153         CALL iom_get( inum3, jpdom_data, 'gphif', gphif )
154
155         CALL iom_get( inum3, jpdom_data, 'e1t', e1t )
156         CALL iom_get( inum3, jpdom_data, 'e1u', e1u )
157         CALL iom_get( inum3, jpdom_data, 'e1v', e1v )
158         
159         CALL iom_get( inum3, jpdom_data, 'e2t', e2t )
160         CALL iom_get( inum3, jpdom_data, 'e2u', e2u )
161         CALL iom_get( inum3, jpdom_data, 'e2v', e2v )
162
163         CALL iom_get( inum3, jpdom_data, 'ff', ff )
164
165         CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
166         mbathy(:,:) = INT( zmbk(:,:) )
167         
168         CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
169         !
170         IF( ln_sco ) THEN                                         ! s-coordinate
171            CALL iom_get( inum4, jpdom_data, 'hbatt', hbatt )
172            CALL iom_get( inum4, jpdom_data, 'hbatu', hbatu )
173            CALL iom_get( inum4, jpdom_data, 'hbatv', hbatv )
174            CALL iom_get( inum4, jpdom_data, 'hbatf', hbatf )
175           
176            CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
177            CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw )
178            CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w ) 
179            CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt )
180            CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw )
181
182            CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors
183            CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
184            CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
185            CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
186
187            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
188            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
189         ENDIF
190
191 
192         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps
193            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth
194            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
195            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
196            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
197            !
198            IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors
199               CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )
200               CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
201               CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
202               CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
203            ELSE                                                        ! 2D bottom scale factors
204               CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp )
205               CALL iom_get( inum4, jpdom_data, 'e3w_ps', e3wp )
206               !                                                        ! deduces the 3D scale factors
207               DO jk = 1, jpk
208                  fse3t_n(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
209                  fse3u_n(:,:,jk) = e3t_1d(jk)
210                  fse3v_n(:,:,jk) = e3t_1d(jk)
211                  fse3w_n(:,:,jk) = e3w_1d(jk)
212               END DO
213               DO jj = 1,jpj                                                  ! adjust the deepest values
214                  DO ji = 1,jpi
215                     ik = mbkt(ji,jj)
216                     fse3t_n(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
217                     fse3w_n(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(1) * ( 1._wp - tmask(ji,jj,1) )
218                  END DO
219               END DO
220               DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
221                  DO jj = 1, jpjm1
222                     DO ji = 1, jpim1
223                        fse3u_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji+1,jj,jk) )
224                        fse3v_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji,jj+1,jk) )
225                     END DO
226                  END DO
227               END DO
228               CALL lbc_lnk( fse3u_n(:,:,:) , 'U', 1._wp )   ;   CALL lbc_lnk( fse3uw_n(:,:,:), 'U', 1._wp )   ! lateral boundary conditions
229               CALL lbc_lnk( fse3v_n(:,:,:) , 'V', 1._wp )   ;   CALL lbc_lnk( fse3vw_n(:,:,:), 'V', 1._wp )
230               !
231               DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
232                  WHERE( fse3u_n(:,:,jk) == 0._wp )   fse3u_n(:,:,jk) = e3t_1d(jk)
233                  WHERE( fse3v_n(:,:,jk) == 0._wp )   fse3v_n(:,:,jk) = e3t_1d(jk)
234               END DO
235            END IF
236
237            IF( iom_varid( inum4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
238               CALL iom_get( inum4, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )
239               CALL iom_get( inum4, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )
240            ELSE                                                           ! 2D bottom depth
241               CALL iom_get( inum4, jpdom_data, 'hdept', zprt )
242               CALL iom_get( inum4, jpdom_data, 'hdepw', zprw )
243               !
244               DO jk = 1, jpk                                              ! deduces the 3D depth
245                  fsdept_n(:,:,jk) = gdept_1d(jk)
246                  fsdepw_n(:,:,jk) = gdepw_1d(jk)
247               END DO
248               DO jj = 1, jpj
249                  DO ji = 1, jpi
250                     ik = mbkt(ji,jj)
251                     IF( ik > 0 ) THEN
252                        fsdepw_n(ji,jj,ik+1) = zprw(ji,jj)
253                        fsdept_n(ji,jj,ik  ) = zprt(ji,jj)
254                        fsdept_n(ji,jj,ik+1) = fsdept_n(ji,jj,ik) + fse3t_n(ji,jj,ik)
255                     ENDIF
256                  END DO
257               END DO
258            ENDIF
259            !
260         ENDIF
261
262         IF( ln_zco ) THEN           ! Vertical coordinates and scales factors
263            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
264            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
265            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )
266            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
267            DO jk = 1, jpk
268               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
269               fse3u_n(:,:,jk) = e3t_1d(jk)
270               fse3v_n(:,:,jk) = e3t_1d(jk)
271               fse3w_n(:,:,jk) = e3w_1d(jk)
272               fsdept_n(:,:,jk) = gdept_1d(jk)
273               fsdepw_n(:,:,jk) = gdepw_1d(jk)
274            END DO
275         ENDIF
276
277!!gm BUG in s-coordinate this does not work!
278      ! deepest/shallowest W level Above/Below ~10m
279      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
280      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
281      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
282!!gm end bug
283
284      ! Control printing : Grid informations (if not restart)
285      ! ----------------
286
287      IF(lwp .AND. .NOT.ln_rstart ) THEN
288         WRITE(numout,*)
289         WRITE(numout,*) '          longitude and e1 scale factors'
290         WRITE(numout,*) '          ------------------------------'
291         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
292            glamv(ji,1), glamf(ji,1),   &
293            e1t(ji,1), e1u(ji,1),   &
294            e1v(ji,1), ji = 1, jpi,10)
2959300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
296            f19.10, 1x, f19.10, 1x, f19.10 )
297
298         WRITE(numout,*)
299         WRITE(numout,*) '          latitude and e2 scale factors'
300         WRITE(numout,*) '          -----------------------------'
301         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
302            &                     gphiv(1,jj), gphif(1,jj),   &
303            &                     e2t  (1,jj), e2u  (1,jj),   &
304            &                     e2v  (1,jj), jj = 1, jpj, 10 )
305      ENDIF
306
307
308      IF( nprint == 1 .AND. lwp ) THEN
309         WRITE(numout,*) '          e1u e2u '
310         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
311         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
312         WRITE(numout,*) '          e1v e2v  '
313         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
314         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
315      ENDIF
316
317      IF(lwp) THEN
318         WRITE(numout,*)
319         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
320         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
321         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
322      ENDIF
323
324      DO jk = 1, jpk
325         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
326         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
327      END DO
328      !                                     ! ============================
329      !                                     !        close the files
330      !                                     ! ============================
331      SELECT CASE ( nmsh )
332         CASE ( 1 )               
333            CALL iom_close( inum0 )
334         CASE ( 2 )
335            CALL iom_close( inum1 )
336            CALL iom_close( inum2 )
337         CASE ( 3 )
338            CALL iom_close( inum2 )
339            CALL iom_close( inum3 )
340            CALL iom_close( inum4 )
341      END SELECT
342      !
343      CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw )
344      !
345   END SUBROUTINE dom_rea
346
347
348   SUBROUTINE zgr_bot_level
349      !!----------------------------------------------------------------------
350      !!                    ***  ROUTINE zgr_bot_level  ***
351      !!
352      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
353      !!
354      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
355      !!
356      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
357      !!                                     ocean level at t-, u- & v-points
358      !!                                     (min value = 1 over land)
359      !!----------------------------------------------------------------------
360      !
361      INTEGER ::   ji, jj   ! dummy loop indices
362      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
363      !!----------------------------------------------------------------------
364
365      !
366      IF(lwp) WRITE(numout,*)
367      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
368      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
369      !
370      CALL wrk_alloc( jpi, jpj, zmbk )
371      !
372      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
373      !                                     ! bottom k-index of W-level = mbkt+1
374      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
375         DO ji = 1, jpim1
376            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
377            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
378         END DO
379      END DO
380      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
381      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
382      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
383      !
384      CALL wrk_dealloc( jpi, jpj, zmbk )
385      !
386   END SUBROUTINE zgr_bot_level
387
388   !!======================================================================
389END MODULE domrea
Note: See TracBrowser for help on using the repository browser.