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/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 9732

Last change on this file since 9732 was 9732, checked in by mathiot, 6 years ago

replace the specific treatment of shallow ice shelves by a more general case
File size: 17.7 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   USE wrk_nemo        ! Memory allocation
23   USE timing          ! Timing
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC dom_wri        ! routine called by inidom.F90
29
30   !! * Substitutions
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
34   !! $Id$
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE dom_wri
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE dom_wri  ***
42      !!                   
43      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
44      !!      ocean domain informations (mesh and mask arrays). This (these)
45      !!      file(s) is (are) used for visualisation (SAXO software) and
46      !!      diagnostic computation.
47      !!
48      !! ** Method  :   Write in a file all the arrays generated in routines
49      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
50      !!      the vertical coord. used (z-coord, partial steps, s-coord)
51      !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file
52      !!                         = 2  :   'mesh.nc' and mask.nc' files
53      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
54      !!                                  'mask.nc' files
55      !!      For huge size domain, use option 2 or 3 depending on your
56      !!      vertical coordinate.
57      !!
58      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
59      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
60      !!                        corresponding to the depth of the bottom t- and w-points
61      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
62      !!                        thickness (e3[tw]_ps) of the bottom points
63      !!
64      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
65      !!                                   masks, depth and vertical scale factors
66      !!----------------------------------------------------------------------
67      !!
68      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
69      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
70      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
71      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
72      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
73      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
74      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
75      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
76      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
77      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
78      INTEGER           ::   ji, jj, jk   ! dummy loop indices
79      !                                   !  workspaces
80      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
81      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
85      !
86      CALL wrk_alloc( jpi, jpj, zprt, zprw )
87      CALL wrk_alloc( jpi, jpj, jpk, zdepu, zdepv )
88      !
89      IF(lwp) WRITE(numout,*)
90      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
91      IF(lwp) WRITE(numout,*) '~~~~~~~'
92     
93      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
94      clnam1 = 'mesh'       ! filename (mesh informations)
95      clnam2 = 'mask'       ! filename (mask informations)
96      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
97      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
98     
99      SELECT CASE ( MOD(nmsh, 3) )
100         !                                  ! ============================
101      CASE ( 1 )                            !  create 'mesh_mask.nc' file
102         !                                  ! ============================
103         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
104         inum2 = inum0                                            ! put all the informations
105         inum3 = inum0                                            ! in unit inum0
106         inum4 = inum0
107         
108         !                                  ! ============================
109      CASE ( 2 )                            !  create 'mesh.nc' and
110         !                                  !         'mask.nc' files
111         !                                  ! ============================
112         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
113         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
114         inum3 = inum1                                            ! put mesh informations
115         inum4 = inum1                                            ! in unit inum1
116         !                                  ! ============================
117      CASE ( 0 )                            !  create 'mesh_hgr.nc'
118         !                                  !         'mesh_zgr.nc' and
119         !                                  !         'mask.nc'     files
120         !                                  ! ============================
121         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
122         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
123         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
124         !
125      END SELECT
126     
127      !                                                         ! masks (inum2)
128      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
129      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
130      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
131      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
132     
133      CALL dom_uniq( zprw, 'T' )
134      DO jj = 1, jpj
135         DO ji = 1, jpi
136            jk=mikt(ji,jj) 
137            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
138         END DO
139      END DO                             !    ! unique point mask
140      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
141      CALL dom_uniq( zprw, 'U' )
142      DO jj = 1, jpj
143         DO ji = 1, jpi
144            jk=miku(ji,jj) 
145            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
146         END DO
147      END DO
148      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
149      CALL dom_uniq( zprw, 'V' )
150      DO jj = 1, jpj
151         DO ji = 1, jpi
152            jk=mikv(ji,jj) 
153            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
154         END DO
155      END DO
156      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
157      CALL dom_uniq( zprw, 'F' )
158      DO jj = 1, jpj
159         DO ji = 1, jpi
160            jk=mikf(ji,jj) 
161            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
162         END DO
163      END DO
164      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
165
166      !                                                         ! horizontal mesh (inum3)
167      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
168      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
169      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
170      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
171     
172      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
173      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
174      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
175      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
176     
177      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
178      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
179      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
180      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
181     
182      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
183      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
184      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
185      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
186     
187      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
188     
189      ! note that mbkt is set to 1 over land ==> use surface tmask
190      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
191      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
192      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
193      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points
194      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
195      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points
196      zprt(:,:) = ssmask(:,:) * REAL( bathy(:,:) , wp )
197      CALL iom_rstput( 0, 0, inum4, 'bathy', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points
198           
199           
200      IF( ln_sco ) THEN                                         ! s-coordinate
201         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )
202         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
203         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
204         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
205         !
206         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
207         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
208         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
209         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
210         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
211         !
212         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
213         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
214         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
215         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
216         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio
217         !
218         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
219         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
220         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
221         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )     
222      ENDIF
223     
224      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
225         !
226         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
227            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
228            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
229            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
230            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
231         ELSE                                                   !    ! 2D masked bottom ocean scale factors
232            DO jj = 1,jpj   
233               DO ji = 1,jpi
234                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
235                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
236               END DO
237            END DO
238            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
239            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
240         END IF
241         !
242         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
243            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
244            DO jk = 1,jpk   
245               DO jj = 1, jpjm1   
246                  DO ji = 1, fs_jpim1   ! vector opt.
247                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
248                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
249                  END DO   
250               END DO   
251            END DO
252            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
253            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
254            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
255            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
256         ELSE                                                   !    ! 2D bottom depth
257            DO jj = 1,jpj   
258               DO ji = 1,jpi
259                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj)
260                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj)
261               END DO
262            END DO
263            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
264            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
265         ENDIF
266         !
267         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord.
268         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
269         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
270         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
271      ENDIF
272     
273      IF( ln_zco ) THEN
274         !                                                      ! z-coordinate - full steps
275         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
276         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
277         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
278         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
279      ENDIF
280      !                                     ! ============================
281      !                                     !        close the files
282      !                                     ! ============================
283      SELECT CASE ( MOD(nmsh, 3) )
284      CASE ( 1 )               
285         CALL iom_close( inum0 )
286      CASE ( 2 )
287         CALL iom_close( inum1 )
288         CALL iom_close( inum2 )
289      CASE ( 0 )
290         CALL iom_close( inum2 )
291         CALL iom_close( inum3 )
292         CALL iom_close( inum4 )
293      END SELECT
294      !
295      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
296      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
297      !
298      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
299      !
300   END SUBROUTINE dom_wri
301
302
303   SUBROUTINE dom_uniq( puniq, cdgrd )
304      !!----------------------------------------------------------------------
305      !!                  ***  ROUTINE dom_uniq  ***
306      !!                   
307      !! ** Purpose :   identify unique point of a grid (TUVF)
308      !!
309      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
310      !!                2) check which elements have been changed
311      !!----------------------------------------------------------------------
312      !
313      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
314      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
315      !
316      REAL(wp) ::  zshift   ! shift value link to the process number
317      INTEGER  ::  ji       ! dummy loop indices
318      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
319      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
320      !!----------------------------------------------------------------------
321      !
322      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
323      !
324      CALL wrk_alloc( jpi, jpj, ztstref )
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      CALL wrk_dealloc( jpi, jpj, ztstref )
341      !
342      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
343      !
344   END SUBROUTINE dom_uniq
345
346   !!======================================================================
347END MODULE domwri
Note: See TracBrowser for help on using the repository browser.