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

source: branches/2012/dev_LOCEAN_UKMO_CMCC_INGV_2012/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 3666

Last change on this file since 3666 was 3666, checked in by cetlod, 11 years ago

commit the changes resulting for the merged branches, see ticket #1025

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