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/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 2648

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

Changed OFF_SRC component to use dynamic memory
  • Property svn:keywords set to Id
File size: 17.7 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 wrk_nemo, ONLY: wrk_in_use, wrk_not_released
20
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC   dom_rea    ! routine called by inidom.F90
25   !!----------------------------------------------------------------------
26   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
27   !! $Id$
28   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
29   !!----------------------------------------------------------------------
30CONTAINS
31
32   SUBROUTINE dom_rea
33      !!----------------------------------------------------------------------
34      !!                  ***  ROUTINE dom_rea  ***
35      !!                   
36      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
37      !!      ocean domain informations (mesh and mask arrays). This (these)
38      !!      file(s) is (are) used for visualisation (SAXO software) and
39      !!      diagnostic computation.
40      !!
41      !! ** Method  :   Read in a file all the arrays generated in routines
42      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
43      !!      the vertical coord. used (z-coord, partial steps, s-coord)
44      !!                    nmsh = 1  :   'mesh_mask.nc' file
45      !!                         = 2  :   'mesh.nc' and mask.nc' files
46      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
47      !!                                  'mask.nc' files
48      !!      For huge size domain, use option 2 or 3 depending on your
49      !!      vertical coordinate.
50      !!
51      !! ** input file :
52      !!      meshmask.nc  : domain size, horizontal grid-point position,
53      !!                     masks, depth and vertical scale factors
54      !!----------------------------------------------------------------------
55      USE iom
56      USE wrk_nemo, ONLY: zmbk => wrk_2d_1, zprt => wrk_2d_2, zprw => wrk_2d_3
57      !!
58      INTEGER  ::   ji, jj, jk   ! dummy loop indices
59      INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
60      REAL(wp) ::   zrefdep         ! local real
61      !!----------------------------------------------------------------------
62
63      IF(lwp) WRITE(numout,*)
64      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
65      IF(lwp) WRITE(numout,*) '~~~~~~~'
66
67      IF( wrk_in_use(2, 1,2,3)  ) THEN
68         CALL ctl_stop('dom_rea: ERROR: requested workspace arrays unavailable.') ; RETURN
69      END IF
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         ENDIF
256
257!!gm BUG in s-coordinate this does not work!
258      ! deepest/shallowest W level Above/Below ~10m
259      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_0) )                  ! ref. depth with tolerance (10% of minimum layer thickness)
260      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below ~10m
261      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
262!!gm end bug
263
264      ! Control printing : Grid informations (if not restart)
265      ! ----------------
266
267      IF(lwp .AND. .NOT.ln_rstart ) THEN
268         WRITE(numout,*)
269         WRITE(numout,*) '          longitude and e1 scale factors'
270         WRITE(numout,*) '          ------------------------------'
271         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
272            glamv(ji,1), glamf(ji,1),   &
273            e1t(ji,1), e1u(ji,1),   &
274            e1v(ji,1), ji = 1, jpi,10)
2759300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
276            f19.10, 1x, f19.10, 1x, f19.10 )
277
278         WRITE(numout,*)
279         WRITE(numout,*) '          latitude and e2 scale factors'
280         WRITE(numout,*) '          -----------------------------'
281         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
282            &                     gphiv(1,jj), gphif(1,jj),   &
283            &                     e2t  (1,jj), e2u  (1,jj),   &
284            &                     e2v  (1,jj), jj = 1, jpj, 10 )
285      ENDIF
286
287
288      IF( nprint == 1 .AND. lwp ) THEN
289         WRITE(numout,*) '          e1u e2u '
290         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
291         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
292         WRITE(numout,*) '          e1v e2v  '
293         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
294         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
295      ENDIF
296
297      IF(lwp) THEN
298         WRITE(numout,*)
299         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
300         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
301         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
302      ENDIF
303
304      DO jk = 1, jpk
305         IF( e3w_0  (jk) <= 0._wp .OR. e3t_0  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_0 or e3t_0 =< 0 ' )
306         IF( gdepw_0(jk) <  0._wp .OR. gdept_0(jk) <  0._wp )   CALL ctl_stop( ' gdepw_0 or gdept_0 < 0 ' )
307      END DO
308      !                                     ! ============================
309      !                                     !        close the files
310      !                                     ! ============================
311      SELECT CASE ( nmsh )
312         CASE ( 1 )               
313            CALL iom_close( inum0 )
314         CASE ( 2 )
315            CALL iom_close( inum1 )
316            CALL iom_close( inum2 )
317         CASE ( 3 )
318            CALL iom_close( inum2 )
319            CALL iom_close( inum3 )
320            CALL iom_close( inum4 )
321      END SELECT
322      !
323      IF( wrk_not_released(2, 1,2,3)  ) CALL ctl_stop('dom_rea:failed to release workspace arrays.')
324      !
325   END SUBROUTINE dom_rea
326
327
328   SUBROUTINE zgr_bot_level
329      !!----------------------------------------------------------------------
330      !!                    ***  ROUTINE zgr_bot_level  ***
331      !!
332      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
333      !!
334      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
335      !!
336      !! ** Action  : - update mbathy: level bathymetry (in level index)
337      !!----------------------------------------------------------------------
338      USE wrk_nemo, ONLY: zmbk => wrk_2d_4
339      !
340      INTEGER ::   ji, jj   ! dummy loop indices
341      !!----------------------------------------------------------------------
342
343      !
344      IF(lwp) WRITE(numout,*)
345      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
346      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
347      !
348      IF( wrk_in_use(2, 4) ) THEN
349         CALL ctl_stop('dom_rea: ERROR: requested workspace arrays unavailable.')  ;  RETURN
350      END IF
351      !
352      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
353      !
354      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
355         DO ji = 1, jpim1
356            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
357            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
358         END DO
359      END DO
360      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
361      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
362      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
363      !
364      IF( wrk_not_released(2, 4) ) CALL ctl_stop('dom_rea:failed to release workspace arrays.')
365      !
366   END SUBROUTINE zgr_bot_level
367
368   !!======================================================================
369END MODULE domrea
Note: See TracBrowser for help on using the repository browser.