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

Last change on this file since 3319 was 3294, checked in by rblod, 12 years ago

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