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/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 5075

Last change on this file since 5075 was 5075, checked in by timgraham, 9 years ago

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