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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 9124

Last change on this file since 9124 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

  • Property svn:keywords set to Id
File size: 16.2 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  ! 2016-01  (G. Madec)  simplified mesh_mask.nc file
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dom_wri        : create and write mesh and mask file(s)
15   !!   dom_uniq       : identify unique point of a grid (TUVF)
16   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
17   !!----------------------------------------------------------------------
18   USE dom_oce         ! ocean space and time domain
19   USE phycst ,   ONLY :   rsmall
20   USE wet_dry,   ONLY :   ll_wd  ! Wetting and drying
21   !
22   USE in_out_manager  ! I/O manager
23   USE iom             ! I/O library
24   USE lbclnk          ! lateral boundary conditions - mpp exchanges
25   USE lib_mpp         ! MPP library
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dom_wri              ! routine called by inidom.F90
31   PUBLIC   dom_stiff            ! routine called by inidom.F90
32
33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dom_wri
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dom_wri  ***
45      !!                   
46      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
47      !!      ocean domain informations (mesh and mask arrays). This (these)
48      !!      file(s) is (are) used for visualisation (SAXO software) and
49      !!      diagnostic computation.
50      !!
51      !! ** Method  :   Write in a file all the arrays generated in routines
52      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
53      !!      the vertical coord. used (z-coord, partial steps, s-coord)
54      !!            MOD(nn_msh, 3) = 1  :   'mesh_mask.nc' file
55      !!                         = 2  :   'mesh.nc' and mask.nc' files
56      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
57      !!                                  'mask.nc' files
58      !!      For huge size domain, use option 2 or 3 depending on your
59      !!      vertical coordinate.
60      !!
61      !!      if     nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
62      !!      if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
63      !!                        corresponding to the depth of the bottom t- and w-points
64      !!      if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the
65      !!                        thickness (e3[tw]_ps) of the bottom points
66      !!
67      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
68      !!                                   masks, depth and vertical scale factors
69      !!----------------------------------------------------------------------
70      INTEGER           ::   inum    ! temprary units for 'mesh_mask.nc' file
71      CHARACTER(len=21) ::   clnam   ! filename (mesh and mask informations)
72      INTEGER           ::   ji, jj, jk   ! dummy loop indices
73      INTEGER           ::   izco, izps, isco, icav
74      !                               
75      REAL(wp), DIMENSION(jpi,jpj)     ::   zprt, zprw     ! 2D workspace
76      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu, zdepv   ! 3D workspace
77      !!----------------------------------------------------------------------
78      !
79      IF(lwp) WRITE(numout,*)
80      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
81      IF(lwp) WRITE(numout,*) '~~~~~~~'
82     
83      clnam = 'mesh_mask'  ! filename (mesh and mask informations)
84     
85      !                                  ! ============================
86      !                                  !  create 'mesh_mask.nc' file
87      !                                  ! ============================
88      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
89      !
90      !                                                         ! global domain size
91      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
92      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
93      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
94
95      !                                                         ! domain characteristics
96      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
97      !                                                         ! type of vertical coordinate
98      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
99      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
100      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
101      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
102      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
103      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
104      !                                                         ! ocean cavities under iceshelves
105      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
106      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
107 
108      !                                                         ! masks
109      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
110      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
111      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
112      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
113     
114      CALL dom_uniq( zprw, 'T' )
115      DO jj = 1, jpj
116         DO ji = 1, jpi
117            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
118         END DO
119      END DO                             !    ! unique point mask
120      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 ) 
121      CALL dom_uniq( zprw, 'U' )
122      DO jj = 1, jpj
123         DO ji = 1, jpi
124            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
125         END DO
126      END DO
127      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 ) 
128      CALL dom_uniq( zprw, 'V' )
129      DO jj = 1, jpj
130         DO ji = 1, jpi
131            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
132         END DO
133      END DO
134      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 ) 
135!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
136!!    Here we just remove the output of fmaskutil.
137!      CALL dom_uniq( zprw, 'F' )
138!      DO jj = 1, jpj
139!         DO ji = 1, jpi
140!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
141!         END DO
142!      END DO
143!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 ) 
144!!gm
145
146      !                                                         ! horizontal mesh (inum3)
147      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
148      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
149      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
150      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
151     
152      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
153      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
154      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
155      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
156     
157      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
158      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
159      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
160      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
161     
162      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
163      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
164      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
165      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
166     
167      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor
168      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
169     
170      ! note that mbkt is set to 1 over land ==> use surface tmask
171      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
172      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
173      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
174      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
175      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
176      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
177      !                                                         ! vertical mesh
178      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
179      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
180      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
181      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
182      !
183      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
184      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
185      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
186      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
187      !
188      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
189         CALL dom_stiff( zprt )
190         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio
191      ENDIF
192      !
193      IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
194
195      !                                     ! ============================
196      CALL iom_close( inum )                !        close the files
197      !                                     ! ============================
198      !
199   END SUBROUTINE dom_wri
200
201
202   SUBROUTINE dom_uniq( puniq, cdgrd )
203      !!----------------------------------------------------------------------
204      !!                  ***  ROUTINE dom_uniq  ***
205      !!                   
206      !! ** Purpose :   identify unique point of a grid (TUVF)
207      !!
208      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
209      !!                2) check which elements have been changed
210      !!----------------------------------------------------------------------
211      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
212      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
213      !
214      REAL(wp) ::  zshift   ! shift value link to the process number
215      INTEGER  ::  ji       ! dummy loop indices
216      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
217      REAL(wp), DIMENSION(jpi,jpj) ::   ztstref
218      !!----------------------------------------------------------------------
219      !
220      ! build an array with different values for each element
221      ! in mpp: make sure that these values are different even between process
222      ! -> apply a shift value according to the process number
223      zshift = jpi * jpj * ( narea - 1 )
224      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
225      !
226      puniq(:,:) = ztstref(:,:)                   ! default definition
227      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
228      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
229      !
230      puniq(:,:) = 1.                             ! default definition
231      ! fill only the inner part of the cpu with llbl converted into real
232      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
233      !
234   END SUBROUTINE dom_uniq
235
236
237   SUBROUTINE dom_stiff( px1 )
238      !!----------------------------------------------------------------------
239      !!                  ***  ROUTINE dom_stiff  ***
240      !!                     
241      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
242      !!
243      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
244      !!                Save the maximum in the vertical direction
245      !!                (this number is only relevant in s-coordinates)
246      !!
247      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
248      !!----------------------------------------------------------------------
249      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
250      !
251      INTEGER  ::   ji, jj, jk 
252      REAL(wp) ::   zrxmax
253      REAL(wp), DIMENSION(4) ::   zr1
254      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
255      !!----------------------------------------------------------------------
256      zx1(:,:) = 0._wp
257      zrxmax   = 0._wp
258      zr1(:)   = 0._wp
259      !
260      DO ji = 2, jpim1
261         DO jj = 2, jpjm1
262            DO jk = 1, jpkm1
263!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
264!!             especially since it is gde3w which is used to compute the pressure gradient
265!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
266!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
267               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
268                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
269                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
270                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
271               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
272                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
273                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
274                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
275               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
276                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
277                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
278                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
279               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
280                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
281                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
282                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
283               zrxmax = MAXVAL( zr1(1:4) )
284               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
285            END DO
286         END DO
287      END DO
288      CALL lbc_lnk( zx1, 'T', 1. )
289      !
290      IF( PRESENT( px1 ) )    px1 = zx1
291      !
292      zrxmax = MAXVAL( zx1 )
293      !
294      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
295      !
296      IF(lwp) THEN
297         WRITE(numout,*)
298         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
299         WRITE(numout,*) '~~~~~~~~~'
300      ENDIF
301      !
302   END SUBROUTINE dom_stiff
303
304   !!======================================================================
305END MODULE domwri
Note: See TracBrowser for help on using the repository browser.