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_r12558_HPC-08_epico_Extra_Halo/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DOM/domwri.F90 @ 12978

Last change on this file since 12978 was 12807, checked in by smasson, 4 years ago

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