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

source: NEMO/trunk/src/OCE/DOM/domwri.F90 @ 14433

Last change on this file since 14433 was 14433, checked in by smasson, 3 years ago

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