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

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

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