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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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