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

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

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