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 utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src – NEMO

source: utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src/domwri.F90 @ 10727

Last change on this file since 10727 was 10727, checked in by rblod, 6 years ago

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