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

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 6717

Last change on this file since 6717 was 6717, checked in by gm, 8 years ago

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