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 @ 6923

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

#1692 - branch SIMPLIF_2_usrdef: update comments in usrdef modules
  • Property svn:keywords set to Id
File size: 16.5 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      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0 )         !    ! scale factors
185      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0 )
186      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0 )
187      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0 )
188      !
189      CALL dom_stiff( zprt )
190      CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )      !    ! Max. grid stiffness ratio
191      !
192      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d )  !    ! stretched system
193      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d )
194      CALL iom_rstput( 0, 0, inum, 'gdept_0', gdept_0, ktype = jp_r8 )
195      CALL iom_rstput( 0, 0, inum, 'gdepw_0', gdepw_0, ktype = jp_r8 )
196     
197      !                                     ! ============================
198      CALL iom_close( inum )                !        close the files
199      !                                     ! ============================
200      !
201      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
202      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
203      !
204      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
205      !
206   END SUBROUTINE dom_wri
207
208
209   SUBROUTINE dom_uniq( puniq, cdgrd )
210      !!----------------------------------------------------------------------
211      !!                  ***  ROUTINE dom_uniq  ***
212      !!                   
213      !! ** Purpose :   identify unique point of a grid (TUVF)
214      !!
215      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
216      !!                2) check which elements have been changed
217      !!----------------------------------------------------------------------
218      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
219      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
220      !
221      REAL(wp) ::  zshift   ! shift value link to the process number
222      INTEGER  ::  ji       ! dummy loop indices
223      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
224      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
225      !!----------------------------------------------------------------------
226      !
227      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
228      !
229      CALL wrk_alloc( jpi, jpj, ztstref )
230      !
231      ! build an array with different values for each element
232      ! in mpp: make sure that these values are different even between process
233      ! -> apply a shift value according to the process number
234      zshift = jpi * jpj * ( narea - 1 )
235      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
236      !
237      puniq(:,:) = ztstref(:,:)                   ! default definition
238      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
239      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
240      !
241      puniq(:,:) = 1.                             ! default definition
242      ! fill only the inner part of the cpu with llbl converted into real
243      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
244      !
245      CALL wrk_dealloc( jpi, jpj, ztstref )
246      !
247      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
248      !
249   END SUBROUTINE dom_uniq
250
251
252   SUBROUTINE dom_stiff( px1 )
253      !!----------------------------------------------------------------------
254      !!                  ***  ROUTINE dom_stiff  ***
255      !!                     
256      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
257      !!
258      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
259      !!                Save the maximum in the vertical direction
260      !!                (this number is only relevant in s-coordinates)
261      !!
262      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
263      !!----------------------------------------------------------------------
264      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
265      !
266      INTEGER  ::   ji, jj, jk 
267      REAL(wp) ::   zrxmax
268      REAL(wp), DIMENSION(4) ::   zr1
269      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
270      !!----------------------------------------------------------------------
271      zx1(:,:) = 0._wp
272      zrxmax   = 0._wp
273      zr1(:)   = 0._wp
274      !
275      DO ji = 2, jpim1
276         DO jj = 2, jpjm1
277            DO jk = 1, jpkm1
278!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
279!!             especially since it is gde3w which is used to compute the pressure gradient
280!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
281!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
282               zr1(1) = ABS(  ( 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) )             &
284                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
285                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
286               zr1(2) = ABS(  ( 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) )             &
288                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
289                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
290               zr1(3) = ABS(  ( 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) )             &
292                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
293                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
294               zr1(4) = ABS(  ( 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) )             &
296                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
297                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
298               zrxmax = MAXVAL( zr1(1:4) )
299               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
300            END DO
301         END DO
302      END DO
303      CALL lbc_lnk( zx1, 'T', 1. )
304      !
305      IF( PRESENT( px1 ) )    px1 = zx1
306      !
307      zrxmax = MAXVAL( zx1 )
308      !
309      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
310      !
311      IF(lwp) THEN
312         WRITE(numout,*)
313         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
314         WRITE(numout,*) '~~~~~~~~~'
315      ENDIF
316      !
317   END SUBROUTINE dom_stiff
318
319   !!======================================================================
320END MODULE domwri
Note: See TracBrowser for help on using the repository browser.