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

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

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