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

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

First attempt to put dynamic allocation on the trunk

  • 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 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         ENDIF
257
258!!gm BUG in s-coordinate this does not work!
259      ! deepest/shallowest W level Above/Below ~10m
260      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_0) )                  ! ref. depth with tolerance (10% of minimum layer thickness)
261      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below ~10m
262      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
263!!gm end bug
264
265      ! Control printing : Grid informations (if not restart)
266      ! ----------------
267
268      IF(lwp .AND. .NOT.ln_rstart ) THEN
269         WRITE(numout,*)
270         WRITE(numout,*) '          longitude and e1 scale factors'
271         WRITE(numout,*) '          ------------------------------'
272         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
273            glamv(ji,1), glamf(ji,1),   &
274            e1t(ji,1), e1u(ji,1),   &
275            e1v(ji,1), ji = 1, jpi,10)
2769300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
277            f19.10, 1x, f19.10, 1x, f19.10 )
278
279         WRITE(numout,*)
280         WRITE(numout,*) '          latitude and e2 scale factors'
281         WRITE(numout,*) '          -----------------------------'
282         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
283            &                     gphiv(1,jj), gphif(1,jj),   &
284            &                     e2t  (1,jj), e2u  (1,jj),   &
285            &                     e2v  (1,jj), jj = 1, jpj, 10 )
286      ENDIF
287
288
289      IF( nprint == 1 .AND. lwp ) THEN
290         WRITE(numout,*) '          e1u e2u '
291         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
292         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
293         WRITE(numout,*) '          e1v e2v  '
294         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
295         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
296      ENDIF
297
298      IF(lwp) THEN
299         WRITE(numout,*)
300         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
301         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
302         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
303      ENDIF
304
305      DO jk = 1, jpk
306         IF( e3w_0  (jk) <= 0._wp .OR. e3t_0  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_0 or e3t_0 =< 0 ' )
307         IF( gdepw_0(jk) <  0._wp .OR. gdept_0(jk) <  0._wp )   CALL ctl_stop( ' gdepw_0 or gdept_0 < 0 ' )
308      END DO
309      !                                     ! ============================
310      !                                     !        close the files
311      !                                     ! ============================
312      SELECT CASE ( nmsh )
313         CASE ( 1 )               
314            CALL iom_close( inum0 )
315         CASE ( 2 )
316            CALL iom_close( inum1 )
317            CALL iom_close( inum2 )
318         CASE ( 3 )
319            CALL iom_close( inum2 )
320            CALL iom_close( inum3 )
321            CALL iom_close( inum4 )
322      END SELECT
323      !
324      IF( wrk_not_released(2, 1,2,3)  ) CALL ctl_stop('dom_rea:failed to release workspace arrays.')
325      !
326   END SUBROUTINE dom_rea
327
328
329   SUBROUTINE zgr_bot_level
330      !!----------------------------------------------------------------------
331      !!                    ***  ROUTINE zgr_bot_level  ***
332      !!
333      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
334      !!
335      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
336      !!
337      !! ** Action  : - update mbathy: level bathymetry (in level index)
338      !!----------------------------------------------------------------------
339      USE wrk_nemo, ONLY: zmbk => wrk_2d_4
340      !
341      INTEGER ::   ji, jj   ! dummy loop indices
342      !!----------------------------------------------------------------------
343
344      !
345      IF(lwp) WRITE(numout,*)
346      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
347      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
348      !
349      IF( wrk_in_use(2, 4) ) THEN
350         CALL ctl_stop('dom_rea: ERROR: requested workspace arrays unavailable.')  ;  RETURN
351      END IF
352      !
353      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
354      !
355      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
356         DO ji = 1, jpim1
357            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
358            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
359         END DO
360      END DO
361      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
362      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
363      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
364      !
365      IF( wrk_not_released(2, 4) ) CALL ctl_stop('dom_rea:failed to release workspace arrays.')
366      !
367   END SUBROUTINE zgr_bot_level
368
369   !!======================================================================
370END MODULE domrea
Note: See TracBrowser for help on using the repository browser.