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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 10 years ago

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