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 utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/domwri.F90 @ 14931

Last change on this file since 14931 was 14931, checked in by jchanut, 3 years ago

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