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

source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 8895

Last change on this file since 8895 was 8895, checked in by andmirek, 6 years ago

#1978 new timers for NEMO restart write
File size: 19.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   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      IF(nn_timing == 2)  CALL timing_start('rst_put')
129      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
130      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
131      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
132      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
133      IF(nn_timing == 2)  CALL timing_stop('rst_put')
134     
135      CALL dom_uniq( zprw, 'T' )
136      DO jj = 1, jpj
137         DO ji = 1, jpi
138            jk=mikt(ji,jj) 
139            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
140         END DO
141      END DO                             !    ! unique point mask
142      IF(nn_timing == 2)  CALL timing_start('rst_put')
143      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
144      IF(nn_timing == 2)  CALL timing_stop('rst_put')
145      CALL dom_uniq( zprw, 'U' )
146      DO jj = 1, jpj
147         DO ji = 1, jpi
148            jk=miku(ji,jj) 
149            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
150         END DO
151      END DO
152      IF(nn_timing == 2)  CALL timing_start('rst_put')
153      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
154      IF(nn_timing == 2)  CALL timing_stop('rst_put')
155      CALL dom_uniq( zprw, 'V' )
156      DO jj = 1, jpj
157         DO ji = 1, jpi
158            jk=mikv(ji,jj) 
159            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
160         END DO
161      END DO
162      IF(nn_timing == 2)  CALL timing_start('rst_put')
163      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
164      IF(nn_timing == 2)  CALL timing_stop('rst_put')
165      CALL dom_uniq( zprw, 'F' )
166      DO jj = 1, jpj
167         DO ji = 1, jpi
168            jk=mikf(ji,jj) 
169            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
170         END DO
171      END DO
172      IF(nn_timing == 2)  CALL timing_start('rst_put')
173      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
174
175      !                                                         ! horizontal mesh (inum3)
176      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
177      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
178      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
179      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
180     
181      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
182      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
183      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
184      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
185     
186      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
187      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
188      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
189      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
190     
191      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
192      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
193      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
194      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
195     
196      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
197      IF(nn_timing == 2)  CALL timing_stop('rst_put') 
198      ! note that mbkt is set to 1 over land ==> use surface tmask
199      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
200      IF(nn_timing == 2)  CALL timing_start('rst_put')
201      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
202      IF(nn_timing == 2)  CALL timing_stop('rst_put')
203      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
204      IF(nn_timing == 2)  CALL timing_start('rst_put')
205      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points
206      IF(nn_timing == 2)  CALL timing_stop('rst_put')
207      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
208      IF(nn_timing == 2)  CALL timing_start('rst_put')
209      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points
210      IF(nn_timing == 2)  CALL timing_stop('rst_put')
211           
212      IF( ln_sco ) THEN                                         ! s-coordinate
213         IF(nn_timing == 2)  CALL timing_start('rst_put')
214         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )
215         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
216         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
217         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
218         !
219         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
220         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
221         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
222         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
223         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
224         !
225         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
226         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
227         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
228         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
229         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio
230         !
231         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
232         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
233         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
234         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )     
235         IF(nn_timing == 2)  CALL timing_stop('rst_put')
236      ENDIF
237     
238      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
239         !
240         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
241            IF(nn_timing == 2)  CALL timing_start('rst_put')
242            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
243            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
244            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
245            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
246            IF(nn_timing == 2)  CALL timing_stop('rst_put')
247         ELSE                                                   !    ! 2D masked bottom ocean scale factors
248            DO jj = 1,jpj   
249               DO ji = 1,jpi
250                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
251                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
252               END DO
253            END DO
254            IF(nn_timing == 2)  CALL timing_start('rst_put')
255            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
256            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
257            IF(nn_timing == 2)  CALL timing_stop('rst_put')
258         END IF
259         !
260         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
261            IF(nn_timing == 2)  CALL timing_start('rst_put')
262            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
263            IF(nn_timing == 2)  CALL timing_stop('rst_put')
264            DO jk = 1,jpk   
265               DO jj = 1, jpjm1   
266                  DO ji = 1, fs_jpim1   ! vector opt.
267                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
268                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
269                  END DO   
270               END DO   
271            END DO
272            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
273            IF(nn_timing == 2)  CALL timing_start('rst_put')
274            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
275            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
276            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
277            IF(nn_timing == 2)  CALL timing_stop('rst_put')
278         ELSE                                                   !    ! 2D bottom depth
279            DO jj = 1,jpj   
280               DO ji = 1,jpi
281                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj)
282                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj)
283               END DO
284            END DO
285            IF(nn_timing == 2)  CALL timing_start('rst_put')
286            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
287            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
288            IF(nn_timing == 2)  CALL timing_stop('rst_put')
289         ENDIF
290         !
291         IF(nn_timing == 2)  CALL timing_start('rst_put')
292         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord.
293         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
294         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
295         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
296         IF(nn_timing == 2)  CALL timing_stop('rst_put')
297      ENDIF
298     
299      IF( ln_zco ) THEN
300         !                                                      ! z-coordinate - full steps
301         IF(nn_timing == 2)  CALL timing_start('rst_put')
302         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
303         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
304         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
305         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
306         IF(nn_timing == 2)  CALL timing_stop('rst_put')
307      ENDIF
308      !                                     ! ============================
309      !                                     !        close the files
310      !                                     ! ============================
311      SELECT CASE ( MOD(nmsh, 3) )
312      CASE ( 1 )               
313         CALL iom_close( inum0 )
314      CASE ( 2 )
315         CALL iom_close( inum1 )
316         CALL iom_close( inum2 )
317      CASE ( 0 )
318         CALL iom_close( inum2 )
319         CALL iom_close( inum3 )
320         CALL iom_close( inum4 )
321      END SELECT
322      !
323      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
324      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
325      !
326      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
327      !
328   END SUBROUTINE dom_wri
329
330
331   SUBROUTINE dom_uniq( puniq, cdgrd )
332      !!----------------------------------------------------------------------
333      !!                  ***  ROUTINE dom_uniq  ***
334      !!                   
335      !! ** Purpose :   identify unique point of a grid (TUVF)
336      !!
337      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
338      !!                2) check which elements have been changed
339      !!----------------------------------------------------------------------
340      !
341      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
342      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
343      !
344      REAL(wp) ::  zshift   ! shift value link to the process number
345      INTEGER  ::  ji       ! dummy loop indices
346      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
347      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
348      !!----------------------------------------------------------------------
349      !
350      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
351      !
352      CALL wrk_alloc( jpi, jpj, ztstref )
353      !
354      ! build an array with different values for each element
355      ! in mpp: make sure that these values are different even between process
356      ! -> apply a shift value according to the process number
357      zshift = jpi * jpj * ( narea - 1 )
358      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
359      !
360      puniq(:,:) = ztstref(:,:)                   ! default definition
361      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
362      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
363      !
364      puniq(:,:) = 1.                             ! default definition
365      ! fill only the inner part of the cpu with llbl converted into real
366      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
367      !
368      CALL wrk_dealloc( jpi, jpj, ztstref )
369      !
370      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
371      !
372   END SUBROUTINE dom_uniq
373
374   !!======================================================================
375END MODULE domwri
Note: See TracBrowser for help on using the repository browser.