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 branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 2683

Last change on this file since 2683 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • 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
67      USE wrk_nemo, ONLY: zdepu => wrk_3d_1, zdepv => wrk_3d_2
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: ERROR - requested workspace arrays unavailable.')
84         RETURN
85      END IF
86
87      IF(lwp) WRITE(numout,*)
88      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
89      IF(lwp) WRITE(numout,*) '~~~~~~~'
90     
91      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
92      clnam1 = 'mesh'       ! filename (mesh informations)
93      clnam2 = 'mask'       ! filename (mask informations)
94      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
95      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
96     
97      SELECT CASE ( MOD(nmsh, 3) )
98         !                                  ! ============================
99      CASE ( 1 )                            !  create 'mesh_mask.nc' file
100         !                                  ! ============================
101         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
102         inum2 = inum0                                            ! put all the informations
103         inum3 = inum0                                            ! in unit inum0
104         inum4 = inum0
105         
106         !                                  ! ============================
107      CASE ( 2 )                            !  create 'mesh.nc' and
108         !                                  !         'mask.nc' files
109         !                                  ! ============================
110         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
111         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
112         inum3 = inum1                                            ! put mesh informations
113         inum4 = inum1                                            ! in unit inum1
114         !                                  ! ============================
115      CASE ( 0 )                            !  create 'mesh_hgr.nc'
116         !                                  !         'mesh_zgr.nc' and
117         !                                  !         'mask.nc'     files
118         !                                  ! ============================
119         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
120         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
121         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
122         !
123      END SELECT
124     
125      !                                                         ! masks (inum2)
126      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
127      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
128      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
129      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
130     
131      CALL dom_uniq( zprw, 'T' )
132      zprt = tmask(:,:,1) * zprw                               !    ! unique point mask
133      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
134      CALL dom_uniq( zprw, 'U' )
135      zprt = umask(:,:,1) * zprw
136      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
137      CALL dom_uniq( zprw, 'V' )
138      zprt = vmask(:,:,1) * zprw
139      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
140      CALL dom_uniq( zprw, 'F' )
141      zprt = fmask(:,:,1) * zprw
142      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
143
144      !                                                         ! horizontal mesh (inum3)
145      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
146      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
147      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
148      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
149     
150      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
151      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
152      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
153      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
154     
155      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
156      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
157      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
158      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
159     
160      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
161      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
162      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
163      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
164     
165      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
166     
167      ! note that mbkt is set to 1 over land ==> use surface tmask
168      zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp )
169      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
170           
171      IF( ln_sco ) THEN                                         ! s-coordinate
172         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth
173         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
174         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
175         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
176         !
177         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
178         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
179         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
180         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
181         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
182         !
183         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors
184         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
185         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
186         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
187         !
188         CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system
189         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 )
190      ENDIF
191     
192      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
193         !
194         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
195            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )         
196            CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
197            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
198            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
199         ELSE                                                   !    ! 2D masked bottom ocean scale factors
200            DO jj = 1,jpj   
201               DO ji = 1,jpi
202                  e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
203                  e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
204               END DO
205            END DO
206            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
207            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
208         END IF
209         !
210         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
211            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )     
212            DO jk = 1,jpk   
213               DO jj = 1, jpjm1   
214                  DO ji = 1, fs_jpim1   ! vector opt.
215                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) )
216                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) )
217                  END DO   
218               END DO   
219            END DO
220            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
221            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
222            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
223            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
224         ELSE                                                   !    ! 2D bottom depth
225            DO jj = 1,jpj   
226               DO ji = 1,jpi
227                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
228                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
229               END DO
230            END DO
231            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
232            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
233         ENDIF
234         !
235         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! reference z-coord.
236         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
237         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )
238         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
239      ENDIF
240     
241      IF( ln_zco ) THEN
242         !                                                      ! z-coordinate - full steps
243         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! depth
244         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
245         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )     !    ! scale factors
246         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
247      ENDIF
248      !                                     ! ============================
249      !                                     !        close the files
250      !                                     ! ============================
251      SELECT CASE ( MOD(nmsh, 3) )
252      CASE ( 1 )               
253         CALL iom_close( inum0 )
254      CASE ( 2 )
255         CALL iom_close( inum1 )
256         CALL iom_close( inum2 )
257      CASE ( 0 )
258         CALL iom_close( inum2 )
259         CALL iom_close( inum3 )
260         CALL iom_close( inum4 )
261      END SELECT
262      !
263      IF( wrk_not_released(2, 1,2) .OR. wrk_not_released(3, 1,2) )THEN
264         CALL ctl_stop('dom_wri: ERROR - failed to release workspace arrays.')
265      END IF
266      !
267   END SUBROUTINE dom_wri
268
269
270   SUBROUTINE dom_uniq( puniq, cdgrd )
271      !!----------------------------------------------------------------------
272      !!                  ***  ROUTINE dom_uniq  ***
273      !!                   
274      !! ** Purpose :   identify unique point of a grid (TUVF)
275      !!
276      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
277      !!                2) check which elements have been changed
278      !!----------------------------------------------------------------------
279      !!
280      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
281      USE wrk_nemo, ONLY: ztstref => wrk_2d_1      ! array with different values for each element
282     !!
283      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
284      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
285      !
286      REAL(wp) ::  zshift   ! shift value link to the process number
287      INTEGER  ::  ji       ! dummy loop indices
288      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
289      !!----------------------------------------------------------------------
290
291      IF( wrk_in_use(2, 1) ) THEN
292         CALL ctl_stop('dom_uniq : requested workspace array unavailable.')   ;   RETURN
293         RETURN
294      END IF
295
296      ! build an array with different values for each element
297      ! in mpp: make sure that these values are different even between process
298      ! -> apply a shift value according to the process number
299      zshift = jpi * jpj * ( narea - 1 )
300      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
301      !
302      puniq(:,:) = ztstref(:,:)                   ! default definition
303      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
304      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
305      !
306      puniq(:,:) = 1.                             ! default definition
307      ! fill only the inner part of the cpu with llbl converted into real
308      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
309      !
310      IF( wrk_not_released(2, 1) ) THEN
311         CALL ctl_stop('dom_uniq : failed to release workspace array.')   ;   RETURN
312      END IF
313      !
314   END SUBROUTINE dom_uniq
315
316   !!======================================================================
317END MODULE domwri
Note: See TracBrowser for help on using the repository browser.