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 branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 2788

Last change on this file since 2788 was 2787, checked in by cetlod, 13 years ago

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