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/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DOM/domwri.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496
  • Property svn:keywords set to Id
File size: 13.1 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_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
16   !!----------------------------------------------------------------------
17   !
18   USE dom_oce         ! ocean space and time domain
19   USE domutl          !
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 "do_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      !                                                         ! domain characteristics
77      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
78      !                                                         ! type of vertical coordinate
79      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
80      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
81      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
82      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
83      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
84      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
85      !                                                         ! ocean cavities under iceshelves
86      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
87      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
88 
89      !                                                         ! masks
90      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
91      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
92      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
93      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
94     
95      CALL dom_uniq( zprw, 'T' )
96      DO_2D( 1, 1, 1, 1 )
97         zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
98      END_2D
99      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 ) 
100      CALL dom_uniq( zprw, 'U' )
101      DO_2D( 1, 1, 1, 1 )
102         zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
103      END_2D
104      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 ) 
105      CALL dom_uniq( zprw, 'V' )
106      DO_2D( 1, 1, 1, 1 )
107         zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
108      END_2D
109      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 ) 
110!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
111!!    Here we just remove the output of fmaskutil.
112!      CALL dom_uniq( zprw, 'F' )
113!      DO jj = 1, jpj
114!         DO ji = 1, jpi
115!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
116!         END DO
117!      END DO
118!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 ) 
119!!gm
120
121      !                                                         ! horizontal mesh (inum3)
122      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
123      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
124      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
125      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
126     
127      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
128      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
129      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
130      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
131     
132      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
133      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
134      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
135      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
136     
137      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
138      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
139      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
140      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
141     
142      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor
143      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
144     
145      ! note that mbkt is set to 1 over land ==> use surface tmask
146      zprt(:,:) = REAL( mbkt(:,:) , wp )
147      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
148      zprt(:,:) = REAL( mikt(:,:) , wp )
149      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
150      !                                                         ! vertical mesh
151      CALL iom_rstput( 0, 0, inum, 'e3t_1d', e3t_1d, ktype = jp_r8  )    !    ! scale factors
152      CALL iom_rstput( 0, 0, inum, 'e3w_1d', e3w_1d, ktype = jp_r8  )
153     
154      CALL iom_rstput( 0, 0, inum, 'e3t_0' , e3t_0 , ktype = jp_r8  )
155      CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8  )
156      CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8  )
157      CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8  )
158      CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8  )
159      CALL iom_rstput( 0, 0, inum, 'e3uw_0', e3uw_0, ktype = jp_r8  )
160      CALL iom_rstput( 0, 0, inum, 'e3vw_0', e3vw_0, ktype = jp_r8  )
161      !
162      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
163      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
164      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
165      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
166      !
167      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
168         CALL dom_stiff( zprt )
169         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio
170      ENDIF
171      !
172      IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
173
174      !                                     ! ============================
175      CALL iom_close( inum )                !        close the files
176      !                                     ! ============================
177   END SUBROUTINE dom_wri
178
179
180   SUBROUTINE dom_stiff( px1 )
181      !!----------------------------------------------------------------------
182      !!                  ***  ROUTINE dom_stiff  ***
183      !!                     
184      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
185      !!
186      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
187      !!                Save the maximum in the vertical direction
188      !!                (this number is only relevant in s-coordinates)
189      !!
190      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
191      !!----------------------------------------------------------------------
192      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
193      !
194      INTEGER  ::   ji, jj, jk 
195      REAL(wp) ::   zrxmax
196      REAL(wp), DIMENSION(4) ::   zr1
197      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
198      !!----------------------------------------------------------------------
199      zx1(:,:) = 0._wp
200      zrxmax   = 0._wp
201      zr1(:)   = 0._wp
202      !
203      DO ji = 2, jpim1
204         DO jj = 2, jpjm1
205            DO jk = 1, jpkm1
206!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
207!!             especially since it is gde3w which is used to compute the pressure gradient
208!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
209!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
210               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
211                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
212                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
213                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
214               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
215                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
216                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
217                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
218               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
219                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
220                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
221                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
222               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
223                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
224                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
225                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
226               zrxmax = MAXVAL( zr1(1:4) )
227               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
228            END DO
229         END DO
230      END DO
231#if defined key_mpi3
232      CALL lbc_lnk_nc_multi( 'domwri', zx1, 'T', 1.0_wp )
233#else
234      CALL lbc_lnk( 'domwri', zx1, 'T', 1.0_wp )
235#endif
236      !
237      IF( PRESENT( px1 ) )    px1 = zx1
238      !
239      zrxmax = MAXVAL( zx1 )
240      !
241      CALL mpp_max( 'domwri', zrxmax ) ! max over the global domain
242      !
243      IF(lwp) THEN
244         WRITE(numout,*)
245         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
246         WRITE(numout,*) '~~~~~~~~~'
247      ENDIF
248      !
249   END SUBROUTINE dom_stiff
250
251   !!======================================================================
252END MODULE domwri
Note: See TracBrowser for help on using the repository browser.