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

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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