source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 16.6 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 :   ln_wd, ht_wd
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   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), DIMENSION(jpi,jpj)   ::   zprt, zprw     ! 2D workspace
77      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu, zdepv   ! 3D workspace
78      !!----------------------------------------------------------------------
79      !
80      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
81      !
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      !                                                         ! vertical mesh
182      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
183      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
184      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
185      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
186      !
187      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
188      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
189      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
190      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
191      !
192      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
193         CALL dom_stiff( zprt )
194         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio
195      ENDIF
196      !
197      IF( ln_wd ) THEN                                          ! wetting and drying domain
198         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
199         CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 )
200      ENDIF
201      !                                     ! ============================
202      CALL iom_close( inum )                !        close the files
203      !                                     ! ============================
204      !
205      !
206      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
207      !
208   END SUBROUTINE dom_wri
209
210
211   SUBROUTINE dom_uniq( puniq, cdgrd )
212      !!----------------------------------------------------------------------
213      !!                  ***  ROUTINE dom_uniq  ***
214      !!                   
215      !! ** Purpose :   identify unique point of a grid (TUVF)
216      !!
217      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
218      !!                2) check which elements have been changed
219      !!----------------------------------------------------------------------
220      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
221      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
222      !
223      REAL(wp) ::  zshift   ! shift value link to the process number
224      INTEGER  ::  ji       ! dummy loop indices
225      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
226      REAL(wp), DIMENSION(jpi,jpj) :: ztstref
227      !!----------------------------------------------------------------------
228      !
229      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
230      !
231      !
232      ! build an array with different values for each element
233      ! in mpp: make sure that these values are different even between process
234      ! -> apply a shift value according to the process number
235      zshift = jpi * jpj * ( narea - 1 )
236      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
237      !
238      puniq(:,:) = ztstref(:,:)                   ! default definition
239      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
240      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
241      !
242      puniq(:,:) = 1.                             ! default definition
243      ! fill only the inner part of the cpu with llbl converted into real
244      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
245      !
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.