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

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package

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