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/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domwri.F90 @ 12172

Last change on this file since 12172 was 12166, checked in by cetlod, 4 years ago

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