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

Last change on this file since 9023 was 9023, checked in by timgraham, 6 years ago

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC
  • 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 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   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( ln_timing )   CALL timing_start('dom_wri')
81      !
82      IF(lwp) WRITE(numout,*)
83      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
84      IF(lwp) WRITE(numout,*) '~~~~~~~'
85     
86      clnam = 'mesh_mask'  ! filename (mesh and mask informations)
87     
88      !                                  ! ============================
89      !                                  !  create 'mesh_mask.nc' file
90      !                                  ! ============================
91      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
92      !
93      !                                                         ! global domain size
94      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
95      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
96      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
97
98      !                                                         ! domain characteristics
99      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
100      !                                                         ! type of vertical coordinate
101      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
102      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
103      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
104      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
105      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
106      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
107      !                                                         ! ocean cavities under iceshelves
108      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
109      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
110 
111      !                                                         ! masks
112      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
113      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
114      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
115      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
116     
117      CALL dom_uniq( zprw, 'T' )
118      DO jj = 1, jpj
119         DO ji = 1, jpi
120            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
121         END DO
122      END DO                             !    ! unique point mask
123      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 ) 
124      CALL dom_uniq( zprw, 'U' )
125      DO jj = 1, jpj
126         DO ji = 1, jpi
127            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
128         END DO
129      END DO
130      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 ) 
131      CALL dom_uniq( zprw, 'V' )
132      DO jj = 1, jpj
133         DO ji = 1, jpi
134            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
135         END DO
136      END DO
137      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 ) 
138!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
139!!    Here we just remove the output of fmaskutil.
140!      CALL dom_uniq( zprw, 'F' )
141!      DO jj = 1, jpj
142!         DO ji = 1, jpi
143!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
144!         END DO
145!      END DO
146!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 ) 
147!!gm
148
149      !                                                         ! horizontal mesh (inum3)
150      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
151      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
152      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
153      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
154     
155      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
156      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
157      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
158      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
159     
160      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
161      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
162      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
163      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
164     
165      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
166      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
167      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
168      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
169     
170      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor
171      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
172     
173      ! note that mbkt is set to 1 over land ==> use surface tmask
174      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
175      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
176      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
177      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
178      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
179      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
180      !                                                         ! vertical mesh
181      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
182      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
183      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
184      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
185      !
186      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
187      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
188      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
189      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
190      !
191      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
192         CALL dom_stiff( zprt )
193         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio
194      ENDIF
195      !
196      IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
197
198      !                                     ! ============================
199      CALL iom_close( inum )                !        close the files
200      !                                     ! ============================
201      !
202      IF( ln_timing )   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), DIMENSION(jpi,jpj) ::   ztstref
223      !!----------------------------------------------------------------------
224      !
225      IF( ln_timing )   CALL timing_start('dom_uniq')
226      !
227      ! build an array with different values for each element
228      ! in mpp: make sure that these values are different even between process
229      ! -> apply a shift value according to the process number
230      zshift = jpi * jpj * ( narea - 1 )
231      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
232      !
233      puniq(:,:) = ztstref(:,:)                   ! default definition
234      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
235      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
236      !
237      puniq(:,:) = 1.                             ! default definition
238      ! fill only the inner part of the cpu with llbl converted into real
239      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
240      !
241      IF( ln_timing )   CALL timing_stop('dom_uniq')
242      !
243   END SUBROUTINE dom_uniq
244
245
246   SUBROUTINE dom_stiff( px1 )
247      !!----------------------------------------------------------------------
248      !!                  ***  ROUTINE dom_stiff  ***
249      !!                     
250      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
251      !!
252      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
253      !!                Save the maximum in the vertical direction
254      !!                (this number is only relevant in s-coordinates)
255      !!
256      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
257      !!----------------------------------------------------------------------
258      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
259      !
260      INTEGER  ::   ji, jj, jk 
261      REAL(wp) ::   zrxmax
262      REAL(wp), DIMENSION(4) ::   zr1
263      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
264      !!----------------------------------------------------------------------
265      zx1(:,:) = 0._wp
266      zrxmax   = 0._wp
267      zr1(:)   = 0._wp
268      !
269      DO ji = 2, jpim1
270         DO jj = 2, jpjm1
271            DO jk = 1, jpkm1
272!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
273!!             especially since it is gde3w which is used to compute the pressure gradient
274!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
275!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
276               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
277                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
278                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
279                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
280               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
281                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
282                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
283                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
284               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
285                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
286                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
287                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
288               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
289                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
290                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
291                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
292               zrxmax = MAXVAL( zr1(1:4) )
293               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
294            END DO
295         END DO
296      END DO
297      CALL lbc_lnk( zx1, 'T', 1. )
298      !
299      IF( PRESENT( px1 ) )    px1 = zx1
300      !
301      zrxmax = MAXVAL( zx1 )
302      !
303      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
304      !
305      IF(lwp) THEN
306         WRITE(numout,*)
307         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
308         WRITE(numout,*) '~~~~~~~~~'
309      ENDIF
310      !
311   END SUBROUTINE dom_stiff
312
313   !!======================================================================
314END MODULE domwri
Note: See TracBrowser for help on using the repository browser.