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

source: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 7277

Last change on this file since 7277 was 7277, checked in by flavoni, 8 years ago

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