source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DOM/domwri.F90 @ 12165

Last change on this file since 12165 was 11586, checked in by gsamson, 16 months ago

dev_r11265_ABL : see #2131

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