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

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

v3.6_alpha : defines the surface cells in offline, see ticket #1191
  • Property svn:keywords set to Id
File size: 18.2 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         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', fse3t_n(:,:,:) ) ! scale factors
176            CALL iom_get( inum4, jpdom_data, 'e3u', fse3u_n(:,:,:) )
177            CALL iom_get( inum4, jpdom_data, 'e3v', fse3v_n(:,:,:) )
178            CALL iom_get( inum4, jpdom_data, 'e3w', fse3w_n(:,:,:) )
179
180            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
181            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
182         ENDIF
183
184 
185         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps
186            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth
187            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
188            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
189            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
190            !
191            IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors
192               CALL iom_get( inum4, jpdom_data, 'e3t', fse3t_n(:,:,:) )
193               CALL iom_get( inum4, jpdom_data, 'e3u', fse3u_n(:,:,:) )
194               CALL iom_get( inum4, jpdom_data, 'e3v', fse3v_n(:,:,:) )
195               CALL iom_get( inum4, jpdom_data, 'e3w', fse3w_n(:,:,:) )
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                  fse3t_n(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
202                  fse3u_n(:,:,jk) = e3t_1d(jk)
203                  fse3v_n(:,:,jk) = e3t_1d(jk)
204                  fse3w_n(:,:,jk) = e3w_1d(jk)
205               END DO
206               DO jj = 1,jpj                                                  ! adjust the deepest values
207                  DO ji = 1,jpi
208                     ik = mbkt(ji,jj)
209                     fse3t_n(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
210                     fse3w_n(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(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                        fse3u_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji+1,jj,jk) )
217                        fse3v_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji,jj+1,jk) )
218                     END DO
219                  END DO
220               END DO
221               CALL lbc_lnk( fse3u_n(:,:,:) , 'U', 1._wp )   ;   CALL lbc_lnk( fse3uw_n(:,:,:), 'U', 1._wp )   ! lateral boundary conditions
222               CALL lbc_lnk( fse3v_n(:,:,:) , 'V', 1._wp )   ;   CALL lbc_lnk( fse3vw_n(:,:,:), 'V', 1._wp )
223               !
224               DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
225                  WHERE( fse3u_n(:,:,jk) == 0._wp )   fse3u_n(:,:,jk) = e3t_1d(jk)
226                  WHERE( fse3v_n(:,:,jk) == 0._wp )   fse3v_n(:,:,jk) = e3t_1d(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', fsdept_n(:,:,:) )
232               CALL iom_get( inum4, jpdom_data, 'gdepw', fsdepw_n(:,:,:) )
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                  fsdept_n(:,:,jk) = gdept_1d(jk)
239                  fsdepw_n(:,:,jk) = gdepw_1d(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                        fsdepw_n(ji,jj,ik+1) = zprw(ji,jj)
246                        fsdept_n(ji,jj,ik  ) = zprt(ji,jj)
247                        fsdept_n(ji,jj,ik+1) = fsdept_n(ji,jj,ik) + fse3t_n(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_1d', gdept_1d ) ! depth
257            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
258            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )
259            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
260            DO jk = 1, jpk
261               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
262               fse3u_n(:,:,jk) = e3t_1d(jk)
263               fse3v_n(:,:,jk) = e3t_1d(jk)
264               fse3w_n(:,:,jk) = e3w_1d(jk)
265               fsdept_n(:,:,jk) = gdept_1d(jk)
266               fsdepw_n(:,:,jk) = gdepw_1d(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_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
273      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > 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_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
315      ENDIF
316
317      DO jk = 1, jpk
318         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
319         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 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.