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

source: utils/tools_ticket2457/DOMAINcfg/src/domwri.F90 @ 12871

Last change on this file since 12871 was 12414, checked in by smueller, 4 years ago

Reintegration of 2019 development branch /utils/tools_MERGE_2019 into the tools directory (/utils/tools)

File size: 15.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 !  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(:,:) = 1.                             ! default definition
218      ! fill only the inner part of the cpu with llbl converted into real
219      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
220      !
221   END SUBROUTINE dom_uniq
222
223
224   SUBROUTINE dom_stiff( px1 )
225      !!----------------------------------------------------------------------
226      !!                  ***  ROUTINE dom_stiff  ***
227      !!                     
228      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
229      !!
230      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
231      !!                Save the maximum in the vertical direction
232      !!                (this number is only relevant in s-coordinates)
233      !!
234      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
235      !!----------------------------------------------------------------------
236      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
237      !
238      INTEGER  ::   ji, jj, jk 
239      REAL(wp) ::   zrxmax
240      REAL(wp), DIMENSION(4) ::   zr1
241      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
242      !!----------------------------------------------------------------------
243      zx1(:,:) = 0._wp
244      zrxmax   = 0._wp
245      zr1(:)   = 0._wp
246      !
247      DO ji = 2, jpim1
248         DO jj = 2, jpjm1
249            DO jk = 1, jpkm1
250!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
251!!             especially since it is gde3w which is used to compute the pressure gradient
252!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
253!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
254               zr1(1) = ABS(  ( 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) )             &
256                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
257                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
258               zr1(2) = ABS(  ( 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) )             &
260                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
261                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
262               zr1(3) = ABS(  ( 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) )             &
264                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
265                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
266               zr1(4) = ABS(  ( 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) )             &
268                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
269                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
270               zrxmax = MAXVAL( zr1(1:4) )
271               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
272            END DO
273         END DO
274      END DO
275      CALL lbc_lnk( 'domwri', zx1, 'T', 1. )
276      !
277      IF( PRESENT( px1 ) )    px1 = zx1
278      !
279      zrxmax = MAXVAL( zx1 )
280      !
281      CALL mpp_max( 'domwri', zrxmax ) ! max over the global domain
282      !
283      IF(lwp) THEN
284         WRITE(numout,*)
285         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
286         WRITE(numout,*) '~~~~~~~~~'
287      ENDIF
288      !
289   END SUBROUTINE dom_stiff
290
291   !!======================================================================
292END MODULE domwri
Note: See TracBrowser for help on using the repository browser.