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 branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 8809

Last change on this file since 8809 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

  • Property svn:keywords set to Id
File size: 17.0 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 :   ln_wd, ht_wd
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   USE wrk_nemo        ! Memory allocation
27   USE timing          ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dom_wri              ! routine called by inidom.F90
33   PUBLIC   dom_stiff            ! routine called by inidom.F90
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE dom_wri
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE dom_wri  ***
47      !!                   
48      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
49      !!      ocean domain informations (mesh and mask arrays). This (these)
50      !!      file(s) is (are) used for visualisation (SAXO software) and
51      !!      diagnostic computation.
52      !!
53      !! ** Method  :   Write in a file all the arrays generated in routines
54      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
55      !!      the vertical coord. used (z-coord, partial steps, s-coord)
56      !!            MOD(nn_msh, 3) = 1  :   'mesh_mask.nc' file
57      !!                         = 2  :   'mesh.nc' and mask.nc' files
58      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
59      !!                                  'mask.nc' files
60      !!      For huge size domain, use option 2 or 3 depending on your
61      !!      vertical coordinate.
62      !!
63      !!      if     nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
64      !!      if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
65      !!                        corresponding to the depth of the bottom t- and w-points
66      !!      if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the
67      !!                        thickness (e3[tw]_ps) of the bottom points
68      !!
69      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
70      !!                                   masks, depth and vertical scale factors
71      !!----------------------------------------------------------------------
72      INTEGER           ::   inum    ! temprary units for 'mesh_mask.nc' file
73      CHARACTER(len=21) ::   clnam   ! filename (mesh and mask informations)
74      INTEGER           ::   ji, jj, jk   ! dummy loop indices
75      INTEGER           ::   izco, izps, isco, icav
76      !                               
77      REAL(wp), POINTER, DIMENSION(:,:)   ::   zprt, zprw     ! 2D workspace
78      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdepu, zdepv   ! 3D workspace
79      !!----------------------------------------------------------------------
80      !
81      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
82      !
83      CALL wrk_alloc( jpi,jpj,       zprt , zprw  )
84      CALL wrk_alloc( jpi,jpj,jpk,   zdepu, zdepv )
85      !
86      IF(lwp) WRITE(numout,*)
87      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
88      IF(lwp) WRITE(numout,*) '~~~~~~~'
89     
90      clnam = 'mesh_mask'  ! filename (mesh and mask informations)
91     
92      !                                  ! ============================
93      !                                  !  create 'mesh_mask.nc' file
94      !                                  ! ============================
95      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
96      !
97      !                                                         ! global domain size
98      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
99      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
100      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
101
102      !                                                         ! domain characteristics
103      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
104      !                                                         ! type of vertical coordinate
105      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
106      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
107      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
108      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
109      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
110      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
111      !                                                         ! ocean cavities under iceshelves
112      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
113      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
114 
115      !                                                         ! masks
116      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
117      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
118      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
119      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
120     
121      CALL dom_uniq( zprw, 'T' )
122      DO jj = 1, jpj
123         DO ji = 1, jpi
124            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
125         END DO
126      END DO                             !    ! unique point mask
127      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 ) 
128      CALL dom_uniq( zprw, 'U' )
129      DO jj = 1, jpj
130         DO ji = 1, jpi
131            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
132         END DO
133      END DO
134      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 ) 
135      CALL dom_uniq( zprw, 'V' )
136      DO jj = 1, jpj
137         DO ji = 1, jpi
138            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
139         END DO
140      END DO
141      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 ) 
142!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
143!!    Here we just remove the output of fmaskutil.
144!      CALL dom_uniq( zprw, 'F' )
145!      DO jj = 1, jpj
146!         DO ji = 1, jpi
147!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
148!         END DO
149!      END DO
150!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 ) 
151!!gm
152
153      !                                                         ! horizontal mesh (inum3)
154      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
155      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
156      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
157      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
158     
159      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
160      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
161      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
162      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
163     
164      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
165      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
166      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
167      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
168     
169      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
170      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
171      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
172      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
173     
174      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor
175      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
176     
177      ! note that mbkt is set to 1 over land ==> use surface tmask
178      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
179      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
180      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
181      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
182      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
183      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
184      !                                                         ! vertical mesh
185      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
186      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
187      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
188      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
189      !
190      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
191      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
192      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
193      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
194      !
195      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
196         CALL dom_stiff( zprt )
197         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio
198      ENDIF
199      !
200      IF( ln_wd ) THEN                                          ! wetting and drying domain
201         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
202         CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 )
203      ENDIF
204      !                                     ! ============================
205      CALL iom_close( inum )                !        close the files
206      !                                     ! ============================
207      !
208      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
209      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
210      !
211      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
212      !
213   END SUBROUTINE dom_wri
214
215
216   SUBROUTINE dom_uniq( puniq, cdgrd )
217      !!----------------------------------------------------------------------
218      !!                  ***  ROUTINE dom_uniq  ***
219      !!                   
220      !! ** Purpose :   identify unique point of a grid (TUVF)
221      !!
222      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
223      !!                2) check which elements have been changed
224      !!----------------------------------------------------------------------
225      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
226      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
227      !
228      REAL(wp) ::  zshift   ! shift value link to the process number
229      INTEGER  ::  ji       ! dummy loop indices
230      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
231      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
232      !!----------------------------------------------------------------------
233      !
234      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
235      !
236      CALL wrk_alloc( jpi, jpj, ztstref )
237      !
238      ! build an array with different values for each element
239      ! in mpp: make sure that these values are different even between process
240      ! -> apply a shift value according to the process number
241      zshift = jpi * jpj * ( narea - 1 )
242      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
243      !
244      puniq(:,:) = ztstref(:,:)                   ! default definition
245      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
246      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
247      !
248      puniq(:,:) = 1.                             ! default definition
249      ! fill only the inner part of the cpu with llbl converted into real
250      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
251      !
252      CALL wrk_dealloc( jpi, jpj, ztstref )
253      !
254      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
255      !
256   END SUBROUTINE dom_uniq
257
258
259   SUBROUTINE dom_stiff( px1 )
260      !!----------------------------------------------------------------------
261      !!                  ***  ROUTINE dom_stiff  ***
262      !!                     
263      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
264      !!
265      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
266      !!                Save the maximum in the vertical direction
267      !!                (this number is only relevant in s-coordinates)
268      !!
269      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
270      !!----------------------------------------------------------------------
271      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
272      !
273      INTEGER  ::   ji, jj, jk 
274      REAL(wp) ::   zrxmax
275      REAL(wp), DIMENSION(4) ::   zr1
276      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
277      !!----------------------------------------------------------------------
278      zx1(:,:) = 0._wp
279      zrxmax   = 0._wp
280      zr1(:)   = 0._wp
281      !
282      DO ji = 2, jpim1
283         DO jj = 2, jpjm1
284            DO jk = 1, jpkm1
285!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
286!!             especially since it is gde3w which is used to compute the pressure gradient
287!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
288!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
289               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
290                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
291                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
292                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
293               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
294                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
295                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
296                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
297               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
298                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
299                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
300                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
301               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
302                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
303                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
304                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
305               zrxmax = MAXVAL( zr1(1:4) )
306               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
307            END DO
308         END DO
309      END DO
310      CALL lbc_lnk( zx1, 'T', 1. )
311      !
312      IF( PRESENT( px1 ) )    px1 = zx1
313      !
314      zrxmax = MAXVAL( zx1 )
315      !
316      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
317      !
318      IF(lwp) THEN
319         WRITE(numout,*)
320         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
321         WRITE(numout,*) '~~~~~~~~~'
322      ENDIF
323      !
324   END SUBROUTINE dom_stiff
325
326   !!======================================================================
327END MODULE domwri
Note: See TracBrowser for help on using the repository browser.