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.
domwri.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 16.3 KB
Line 
1MODULE domwri
2   !!======================================================================
3   !!                       ***  MODULE domwri  ***
4   !! Ocean initialization : write the ocean domain mesh file(s)
5   !!======================================================================
6   !! History :  OPA  ! 1997-02  (G. Madec)  Original code
7   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL
8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file
9   !!            3.0  ! 2008-01  (S. Masson) add dom_uniq
10   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dom_wri        : create and write mesh and mask file(s)
15   !!   dom_uniq       :
16   !!----------------------------------------------------------------------
17   USE dom_oce         ! ocean space and time domain
18   USE in_out_manager  ! I/O manager
19   USE iom             ! I/O library
20   USE lbclnk          ! lateral boundary conditions - mpp exchanges
21   USE lib_mpp         ! MPP library
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC dom_wri        ! routine called by inidom.F90
27
28   !! * Substitutions
29#  include "vectopt_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
32   !! $Id$
33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE dom_wri
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE dom_wri  ***
40      !!                   
41      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
42      !!      ocean domain informations (mesh and mask arrays). This (these)
43      !!      file(s) is (are) used for visualisation (SAXO software) and
44      !!      diagnostic computation.
45      !!
46      !! ** Method  :   Write in a file all the arrays generated in routines
47      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
48      !!      the vertical coord. used (z-coord, partial steps, s-coord)
49      !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file
50      !!                         = 2  :   'mesh.nc' and mask.nc' files
51      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
52      !!                                  'mask.nc' files
53      !!      For huge size domain, use option 2 or 3 depending on your
54      !!      vertical coordinate.
55      !!
56      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
57      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
58      !!                        corresponding to the depth of the bottom t- and w-points
59      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
60      !!                        thickness (e3[tw]_ps) of the bottom points
61      !!
62      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
63      !!                                   masks, depth and vertical scale factors
64      !!----------------------------------------------------------------------
65      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
66      USE wrk_nemo, ONLY:   zprt  => wrk_2d_1 , zprw  => wrk_2d_2    ! 2D workspace
67      USE wrk_nemo, ONLY:   zdepu => wrk_3d_1 , zdepv => wrk_3d_2    ! 3D     -
68      !!
69      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
70      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
71      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
72      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
73      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
74      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
75      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
76      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
77      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
78      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
79      INTEGER           ::   ji, jj, jk   ! dummy loop indices
80      !!----------------------------------------------------------------------
81
82      IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 1,2) )THEN
83         CALL ctl_stop('dom_wri: requested workspace arrays unavailable')   ;   RETURN
84      END IF
85
86      IF(lwp) WRITE(numout,*)
87      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
88      IF(lwp) WRITE(numout,*) '~~~~~~~'
89     
90      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
91      clnam1 = 'mesh'       ! filename (mesh informations)
92      clnam2 = 'mask'       ! filename (mask informations)
93      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
94      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
95     
96      SELECT CASE ( MOD(nmsh, 3) )
97         !                                  ! ============================
98      CASE ( 1 )                            !  create 'mesh_mask.nc' file
99         !                                  ! ============================
100         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
101         inum2 = inum0                                            ! put all the informations
102         inum3 = inum0                                            ! in unit inum0
103         inum4 = inum0
104         
105         !                                  ! ============================
106      CASE ( 2 )                            !  create 'mesh.nc' and
107         !                                  !         'mask.nc' files
108         !                                  ! ============================
109         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
110         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
111         inum3 = inum1                                            ! put mesh informations
112         inum4 = inum1                                            ! in unit inum1
113         !                                  ! ============================
114      CASE ( 0 )                            !  create 'mesh_hgr.nc'
115         !                                  !         'mesh_zgr.nc' and
116         !                                  !         'mask.nc'     files
117         !                                  ! ============================
118         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
119         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
120         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
121         !
122      END SELECT
123     
124      !                                                         ! masks (inum2)
125      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
126      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
127      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
128      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
129     
130      CALL dom_uniq( zprw, 'T' )
131      zprt = tmask(:,:,1) * zprw                               !    ! unique point mask
132      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
133      CALL dom_uniq( zprw, 'U' )
134      zprt = umask(:,:,1) * zprw
135      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
136      CALL dom_uniq( zprw, 'V' )
137      zprt = vmask(:,:,1) * zprw
138      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
139      CALL dom_uniq( zprw, 'F' )
140      zprt = fmask(:,:,1) * zprw
141      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
142
143      !                                                         ! horizontal mesh (inum3)
144      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
145      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
146      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
147      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
148     
149      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
150      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
151      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
152      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
153     
154      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
155      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
156      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
157      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
158     
159      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
160      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
161      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
162      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
163     
164      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
165     
166      ! note that mbkt is set to 1 over land ==> use surface tmask
167      zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp )
168      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
169           
170      IF( ln_sco ) THEN                                         ! s-coordinate
171         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth
172         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
173         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
174         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
175         !
176         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
177         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
178         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
179         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
180         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
181         !
182         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors
183         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
184         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
185         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
186         !
187         CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system
188         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 )
189      ENDIF
190     
191      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
192         !
193         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
194            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )         
195            CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
196            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
197            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
198         ELSE                                                   !    ! 2D masked bottom ocean scale factors
199            DO jj = 1,jpj   
200               DO ji = 1,jpi
201                  e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
202                  e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
203               END DO
204            END DO
205            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
206            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
207         END IF
208         !
209         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
210            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )     
211            DO jk = 1,jpk   
212               DO jj = 1, jpjm1   
213                  DO ji = 1, fs_jpim1   ! vector opt.
214                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) )
215                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) )
216                  END DO   
217               END DO   
218            END DO
219            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
220            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
221            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
222            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
223         ELSE                                                   !    ! 2D bottom depth
224            DO jj = 1,jpj   
225               DO ji = 1,jpi
226                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
227                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
228               END DO
229            END DO
230            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
231            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
232         ENDIF
233         !
234         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! reference z-coord.
235         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
236         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )
237         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
238      ENDIF
239     
240      IF( ln_zco ) THEN
241         !                                                      ! z-coordinate - full steps
242         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! depth
243         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
244         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )     !    ! scale factors
245         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
246      ENDIF
247      !                                     ! ============================
248      !                                     !        close the files
249      !                                     ! ============================
250      SELECT CASE ( MOD(nmsh, 3) )
251      CASE ( 1 )               
252         CALL iom_close( inum0 )
253      CASE ( 2 )
254         CALL iom_close( inum1 )
255         CALL iom_close( inum2 )
256      CASE ( 0 )
257         CALL iom_close( inum2 )
258         CALL iom_close( inum3 )
259         CALL iom_close( inum4 )
260      END SELECT
261      !
262      IF( wrk_not_released(2, 1,2)  .OR.   &
263          wrk_not_released(3, 1,2)  )   CALL ctl_stop('dom_wri: failed to release workspace arrays')
264      !
265   END SUBROUTINE dom_wri
266
267
268   SUBROUTINE dom_uniq( puniq, cdgrd )
269      !!----------------------------------------------------------------------
270      !!                  ***  ROUTINE dom_uniq  ***
271      !!                   
272      !! ** Purpose :   identify unique point of a grid (TUVF)
273      !!
274      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
275      !!                2) check which elements have been changed
276      !!----------------------------------------------------------------------
277      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
278      USE wrk_nemo, ONLY:   ztstref => wrk_2d_3      ! array with different values for each element
279      !
280      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
281      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
282      !
283      REAL(wp) ::  zshift   ! shift value link to the process number
284      INTEGER  ::  ji       ! dummy loop indices
285      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
286      !!----------------------------------------------------------------------
287
288      IF( wrk_in_use(2, 3) ) THEN
289         CALL ctl_stop('dom_uniq: requested workspace array unavailable')   ;   RETURN
290      ENDIF
291
292      ! build an array with different values for each element
293      ! in mpp: make sure that these values are different even between process
294      ! -> apply a shift value according to the process number
295      zshift = jpi * jpj * ( narea - 1 )
296      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
297      !
298      puniq(:,:) = ztstref(:,:)                   ! default definition
299      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
300      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
301      !
302      puniq(:,:) = 1.                             ! default definition
303      ! fill only the inner part of the cpu with llbl converted into real
304      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
305      !
306      IF( wrk_not_released(2, 3) )   CALL ctl_stop('dom_uniq: failed to release workspace array')
307      !
308   END SUBROUTINE dom_uniq
309
310   !!======================================================================
311END MODULE domwri
Note: See TracBrowser for help on using the repository browser.