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/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src – NEMO

source: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domwri.f90 @ 15126

Last change on this file since 15126 was 15121, checked in by dbruciaferri, 3 years ago

adding code for MEs-coordinate system v1.0 (same of PhD)
File size: 24.7 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   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dom_wri        : create and write mesh and mask file(s)
14   !!   dom_uniq       : identify unique point of a grid (TUVF)
15   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
16   !!----------------------------------------------------------------------
17   USE dom_oce         ! ocean space and time domain
18   USE in_out_manager  ! I/O manager
19   USE iom             ! I/O library
20   USE lbclnk          ! lateral boundary conditions - mpp exchanges
21   USE lib_mpp         ! MPP library
22   USE wrk_nemo        ! Memory allocation
23   USE timing          ! Timing
24   USE phycst
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   dom_wri              ! routine called by inidom.F90
30   PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90
31   PUBLIC   dom_stiff            ! routine called by inidom.F90
32
33   !!----------------------------------------------------------------------
34   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
35   !! $Id: vectopt_loop_substitute.h90 4990 2014-12-15 16:42:49Z timgraham $
36   !! Software governed by the CeCILL licence (./LICENSE)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE dom_wri_coordinate
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE dom_wri_coordinate  ***
43      !!                   
44      !! ** Purpose :   Create the NetCDF file which contains all the
45      !!              standard coordinate information plus the surface,
46      !!              e1e2u and e1e2v. By doing so, those surface will
47      !!              not be changed by the reduction of e1u or e2v scale
48      !!              factors in some straits.
49      !!                 NB: call just after the read of standard coordinate
50      !!              and the reduction of scale factors in some straits
51      !!
52      !! ** output file :   coordinate_e1e2u_v.nc
53      !!----------------------------------------------------------------------
54      INTEGER           ::   inum0    ! temprary units for 'coordinate_e1e2u_v.nc' file
55      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
56      !                                   !  workspaces
57      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
58      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
59      !!----------------------------------------------------------------------
60      !
61      IF( nn_timing == 1 )  CALL timing_start('dom_wri_coordinate')
62      !
63      IF(lwp) WRITE(numout,*)
64      IF(lwp) WRITE(numout,*) 'dom_wri_coordinate : create NetCDF coordinate file'
65      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
66     
67      clnam0 = 'coordinate_e1e2u_v'  ! filename (mesh and mask informations)
68     
69      !  create 'coordinate_e1e2u_v.nc' file
70      ! ============================
71      !
72      CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
73      !
74      !                                                         ! horizontal mesh (inum3)
75      CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
76      CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r8 )
77      CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r8 )
78      CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r8 )
79     
80      CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
81      CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r8 )
82      CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r8 )
83      CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r8 )
84     
85      CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
86      CALL iom_rstput( 0, 0, inum0, 'e1u', e1u, ktype = jp_r8 )
87      CALL iom_rstput( 0, 0, inum0, 'e1v', e1v, ktype = jp_r8 )
88      CALL iom_rstput( 0, 0, inum0, 'e1f', e1f, ktype = jp_r8 )
89     
90      CALL iom_rstput( 0, 0, inum0, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
91      CALL iom_rstput( 0, 0, inum0, 'e2u', e2u, ktype = jp_r8 )
92      CALL iom_rstput( 0, 0, inum0, 'e2v', e2v, ktype = jp_r8 )
93      CALL iom_rstput( 0, 0, inum0, 'e2f', e2f, ktype = jp_r8 )
94     
95      CALL iom_rstput( 0, 0, inum0, 'e1e2u', e1e2u, ktype = jp_r8 )
96      CALL iom_rstput( 0, 0, inum0, 'e1e2v', e1e2v, ktype = jp_r8 )
97
98      CALL iom_close( inum0 )
99      !
100      IF( nn_timing == 1 )  CALL timing_stop('dom_wri_coordinate')
101      !
102   END SUBROUTINE dom_wri_coordinate
103
104
105   SUBROUTINE dom_wri
106      !!----------------------------------------------------------------------
107      !!                  ***  ROUTINE dom_wri  ***
108      !!                   
109      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
110      !!      ocean domain informations (mesh and mask arrays). This (these)
111      !!      file(s) is (are) used for visualisation (SAXO software) and
112      !!      diagnostic computation.
113      !!
114      !! ** Method  :   Write in a file all the arrays generated in routines
115      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
116      !!      the vertical coord. used (z-coord, partial steps, s-coord)
117      !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file
118      !!                         = 2  :   'mesh.nc' and mask.nc' files
119      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
120      !!                                  'mask.nc' files
121      !!      For huge size domain, use option 2 or 3 depending on your
122      !!      vertical coordinate.
123      !!
124      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
125      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
126      !!                        corresponding to the depth of the bottom t- and w-points
127      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
128      !!                        thickness (e3[tw]_ps) of the bottom points
129      !!
130      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
131      !!                                   masks, depth and vertical scale factors
132      !!----------------------------------------------------------------------
133      !!
134      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
135      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
136      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
137      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
138      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
139      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
140      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
141      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
142      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
143      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
144      INTEGER           ::   ji, jj, jk   ! dummy loop indices
145      !                                   !  workspaces
146      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
147      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
148      !!----------------------------------------------------------------------
149      !
150      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
151      !
152      CALL wrk_alloc( jpi, jpj, zprt, zprw )
153      CALL wrk_alloc( jpi, jpj, jpk, zdepu, zdepv )
154      !
155      IF(lwp) WRITE(numout,*)
156      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
157      IF(lwp) WRITE(numout,*) '~~~~~~~'
158     
159      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
160      clnam1 = 'mesh'       ! filename (mesh informations)
161      clnam2 = 'mask'       ! filename (mask informations)
162      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
163      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
164     
165      SELECT CASE ( MOD(nmsh, 3) )
166         !                                  ! ============================
167      CASE ( 1 )                            !  create 'mesh_mask.nc' file
168         !                                  ! ============================
169         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
170         inum2 = inum0                                            ! put all the informations
171         inum3 = inum0                                            ! in unit inum0
172         inum4 = inum0
173         
174         !                                  ! ============================
175      CASE ( 2 )                            !  create 'mesh.nc' and
176         !                                  !         'mask.nc' files
177         !                                  ! ============================
178         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
179         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
180         inum3 = inum1                                            ! put mesh informations
181         inum4 = inum1                                            ! in unit inum1
182         !                                  ! ============================
183      CASE ( 0 )                            !  create 'mesh_hgr.nc'
184         !                                  !         'mesh_zgr.nc' and
185         !                                  !         'mask.nc'     files
186         !                                  ! ============================
187         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
188         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
189         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
190         !
191      END SELECT
192     
193      !                                                         ! masks (inum2)
194      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
195      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
196      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
197      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
198     
199      CALL dom_uniq( zprw, 'T' )
200      DO jj = 1, jpj
201         DO ji = 1, jpi
202            jk=mikt(ji,jj) 
203            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
204         END DO
205      END DO                             !    ! unique point mask
206      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
207      CALL dom_uniq( zprw, 'U' )
208      DO jj = 1, jpj
209         DO ji = 1, jpi
210            jk=miku(ji,jj) 
211            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
212         END DO
213      END DO
214      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
215      CALL dom_uniq( zprw, 'V' )
216      DO jj = 1, jpj
217         DO ji = 1, jpi
218            jk=mikv(ji,jj) 
219            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
220         END DO
221      END DO
222      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
223      CALL dom_uniq( zprw, 'F' )
224      DO jj = 1, jpj
225         DO ji = 1, jpi
226            jk=mikf(ji,jj) 
227            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
228         END DO
229      END DO
230      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
231
232      !                                                         ! horizontal mesh (inum3)
233      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
234      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r8 )
235      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r8 )
236      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r8 )
237     
238      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
239      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r8 )
240      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r8 )
241      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r8 )
242     
243      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
244      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
245      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
246      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
247     
248      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
249      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
250      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
251      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
252     
253      CALL iom_rstput( 0, 0, inum3, 'ff_f', ff_f, ktype = jp_r8 )           !    ! coriolis factor
254      CALL iom_rstput( 0, 0, inum3, 'ff_t', ff_t, ktype = jp_r8 )           !    ! coriolis factor
255     
256      ! note that mbkt is set to 1 over land ==> use surface tmask
257      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
258      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
259      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
260      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points
261      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
262      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r8 )       !    ! nb of ocean T-points
263           
264      IF( ln_sco ) THEN                                         ! s-coordinate
265         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )
266         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
267         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
268         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
269         !
270         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
271         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
272         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
273         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
274         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
275      ENDIF
276
277      IF( ln_sco .OR. ln_mes ) THEN                             ! s- / MEs-coordinate
278         !
279         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
280         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
281         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
282         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
283         !
284         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
285         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
286         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )     
287         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
288         CALL dom_stiff( zprt )
289         CALL iom_rstput( 0, 0, inum4, 'stiffness', zprt )       !    ! Max. grid stiffness ratio
290      ENDIF
291     
292      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
293         !
294         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
295            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
296            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
297            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
298            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
299         ELSE                                                   !    ! 2D masked bottom ocean scale factors
300            DO jj = 1,jpj   
301               DO ji = 1,jpi
302                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
303                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
304               END DO
305            END DO
306            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
307            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
308         END IF
309         !
310         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
311            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )     
312            DO jk = 1,jpk   
313               DO jj = 1, jpjm1   
314                  DO ji = 1, jpim1   ! vector opt.
315                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
316                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
317                  END DO   
318               END DO   
319            END DO
320            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
321            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r8 )
322            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r8 )
323            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
324         ELSE                                                   !    ! 2D bottom depth
325            DO jj = 1,jpj   
326               DO ji = 1,jpi
327                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj)
328                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj)
329               END DO
330            END DO
331            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r8 )     
332            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r8 ) 
333         ENDIF
334         !
335         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord.
336         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
337         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
338         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
339      ENDIF
340     
341      IF( ln_zco ) THEN
342         !                                                      ! z-coordinate - full steps
343         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
344         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
345         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
346         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
347      ENDIF
348      !                                     ! ============================
349      !                                     !        close the files
350      !                                     ! ============================
351      SELECT CASE ( MOD(nmsh, 3) )
352      CASE ( 1 )               
353         CALL iom_close( inum0 )
354      CASE ( 2 )
355         CALL iom_close( inum1 )
356         CALL iom_close( inum2 )
357      CASE ( 0 )
358         CALL iom_close( inum2 )
359         CALL iom_close( inum3 )
360         CALL iom_close( inum4 )
361      END SELECT
362      !
363      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
364      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
365      !
366      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
367      !
368   END SUBROUTINE dom_wri
369
370
371   SUBROUTINE dom_uniq( puniq, cdgrd )
372      !!----------------------------------------------------------------------
373      !!                  ***  ROUTINE dom_uniq  ***
374      !!                   
375      !! ** Purpose :   identify unique point of a grid (TUVF)
376      !!
377      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
378      !!                2) check which elements have been changed
379      !!----------------------------------------------------------------------
380      !
381      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
382      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
383      !
384      REAL(wp) ::  zshift   ! shift value link to the process number
385      INTEGER  ::  ji       ! dummy loop indices
386      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
387      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
388      !!----------------------------------------------------------------------
389      !
390      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
391      !
392      CALL wrk_alloc( jpi, jpj, ztstref )
393      !
394      ! build an array with different values for each element
395      ! in mpp: make sure that these values are different even between process
396      ! -> apply a shift value according to the process number
397      zshift = jpi * jpj * ( narea - 1 )
398      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
399      !
400      puniq(:,:) = ztstref(:,:)                   ! default definition
401      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
402      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
403      !
404      puniq(:,:) = 1.                             ! default definition
405      ! fill only the inner part of the cpu with llbl converted into real
406      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
407      !
408      CALL wrk_dealloc( jpi, jpj, ztstref )
409      !
410      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
411      !
412   END SUBROUTINE dom_uniq
413
414
415   SUBROUTINE dom_stiff( px1 )
416      !!----------------------------------------------------------------------
417      !!                  ***  ROUTINE dom_stiff  ***
418      !!                     
419      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
420      !!
421      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
422      !!                Save the maximum in the vertical direction
423      !!                (this number is only relevant in s-coordinates)
424      !!
425      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
426      !!----------------------------------------------------------------------
427      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
428      !
429      INTEGER  ::   ji, jj, jk 
430      REAL(wp) ::   zrxmax
431      REAL(wp), DIMENSION(4) ::   zr1
432      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
433      !!----------------------------------------------------------------------
434      zx1(:,:) = 0._wp
435      zrxmax   = 0._wp
436      zr1(:)   = 0._wp
437      !
438      DO ji = 2, jpim1
439         DO jj = 2, jpjm1
440            DO jk = 1, jpkm1
441!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
442!!             especially since it is gde3w which is used to compute the pressure gradient
443!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
444!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
445               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
446                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
447                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
448                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
449               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
450                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
451                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
452                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
453               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
454                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
455                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
456                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
457               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
458                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
459                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
460                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
461               zrxmax = MAXVAL( zr1(1:4) )
462               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
463            END DO
464         END DO
465      END DO
466      CALL lbc_lnk( zx1, 'T', 1. )
467      !
468      IF( PRESENT( px1 ) )    px1 = zx1
469      !
470      zrxmax = MAXVAL( zx1 )
471      !
472      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
473      !
474      IF(lwp) THEN
475         WRITE(numout,*)
476         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
477         WRITE(numout,*) '~~~~~~~~~'
478      ENDIF
479      !
480   END SUBROUTINE dom_stiff
481
482   !!======================================================================
483END MODULE domwri
Note: See TracBrowser for help on using the repository browser.