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 NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domwri.F90 @ 11201

Last change on this file since 11201 was 11201, checked in by mathiot, 5 years ago

ENHANCE-03_domcfg : add management of closed seas in domain cfg by flood filling and lat/lon seed instead of i/j box definition (ticket #2143)

File size: 15.2 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 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      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
160      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
161      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
162      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
163      !                                                         ! vertical mesh
164      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
165      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
166      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
167      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
168      !
169      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
170      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
171      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
172      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
173      !
174      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
175         CALL dom_stiff( zprt )
176         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio
177      ENDIF
178      !
179      !                                     ! ============================
180      CALL iom_close( inum )                !        close the files
181      !                                     ! ============================
182   END SUBROUTINE dom_wri
183
184
185   SUBROUTINE dom_uniq( puniq, cdgrd )
186      !!----------------------------------------------------------------------
187      !!                  ***  ROUTINE dom_uniq  ***
188      !!                   
189      !! ** Purpose :   identify unique point of a grid (TUVF)
190      !!
191      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
192      !!                2) check which elements have been changed
193      !!----------------------------------------------------------------------
194      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
195      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
196      !
197      REAL(wp) ::  zshift   ! shift value link to the process number
198      INTEGER  ::  ji       ! dummy loop indices
199      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
200      REAL(wp), DIMENSION(jpi,jpj) ::   ztstref
201      !!----------------------------------------------------------------------
202      !
203      ! build an array with different values for each element
204      ! in mpp: make sure that these values are different even between process
205      ! -> apply a shift value according to the process number
206      zshift = jpi * jpj * ( narea - 1 )
207      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
208      !
209      puniq(:,:) = ztstref(:,:)                   ! default definition
210      CALL lbc_lnk( 'domwri', puniq, cdgrd, 1. )            ! apply boundary conditions
211      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
212      !
213      puniq(:,:) = 1.                             ! default definition
214      ! fill only the inner part of the cpu with llbl converted into real
215      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
216      !
217   END SUBROUTINE dom_uniq
218
219
220   SUBROUTINE dom_stiff( px1 )
221      !!----------------------------------------------------------------------
222      !!                  ***  ROUTINE dom_stiff  ***
223      !!                     
224      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
225      !!
226      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
227      !!                Save the maximum in the vertical direction
228      !!                (this number is only relevant in s-coordinates)
229      !!
230      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
231      !!----------------------------------------------------------------------
232      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
233      !
234      INTEGER  ::   ji, jj, jk 
235      REAL(wp) ::   zrxmax
236      REAL(wp), DIMENSION(4) ::   zr1
237      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
238      !!----------------------------------------------------------------------
239      zx1(:,:) = 0._wp
240      zrxmax   = 0._wp
241      zr1(:)   = 0._wp
242      !
243      DO ji = 2, jpim1
244         DO jj = 2, jpjm1
245            DO jk = 1, jpkm1
246!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
247!!             especially since it is gde3w which is used to compute the pressure gradient
248!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
249!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
250               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
251                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
252                    &       / ( 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) + rsmall )  ) * umask(ji-1,jj,jk)
254               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
255                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
256                    &       / ( 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) + rsmall )  ) * umask(ji  ,jj,jk)
258               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
259                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
260                    &       / ( 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) + rsmall )  ) * vmask(ji,jj  ,jk)
262               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
263                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
264                    &       / ( 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) + rsmall )  ) * vmask(ji,jj-1,jk)
266               zrxmax = MAXVAL( zr1(1:4) )
267               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
268            END DO
269         END DO
270      END DO
271      CALL lbc_lnk( 'domwri', zx1, 'T', 1. )
272      !
273      IF( PRESENT( px1 ) )    px1 = zx1
274      !
275      zrxmax = MAXVAL( zx1 )
276      !
277      CALL mpp_max( 'domwri', zrxmax ) ! max over the global domain
278      !
279      IF(lwp) THEN
280         WRITE(numout,*)
281         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
282         WRITE(numout,*) '~~~~~~~~~'
283      ENDIF
284      !
285   END SUBROUTINE dom_stiff
286
287   !!======================================================================
288END MODULE domwri
Note: See TracBrowser for help on using the repository browser.