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

source: NEMO/trunk/tools/DOMAINcfg/src/domwri.f90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 7 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks
File size: 24.6 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         !
276         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
277         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
278         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
279         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
280         !
281         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
282         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
283         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )     
284         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
285         CALL dom_stiff( zprt )
286         CALL iom_rstput( 0, 0, inum4, 'stiffness', zprt )       !    ! Max. grid stiffness ratio
287      ENDIF
288     
289      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
290         !
291         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
292            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
293            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
294            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
295            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
296         ELSE                                                   !    ! 2D masked bottom ocean scale factors
297            DO jj = 1,jpj   
298               DO ji = 1,jpi
299                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
300                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
301               END DO
302            END DO
303            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
304            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
305         END IF
306         !
307         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
308            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )     
309            DO jk = 1,jpk   
310               DO jj = 1, jpjm1   
311                  DO ji = 1, jpim1   ! vector opt.
312                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
313                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
314                  END DO   
315               END DO   
316            END DO
317            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
318            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r8 )
319            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r8 )
320            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
321         ELSE                                                   !    ! 2D bottom depth
322            DO jj = 1,jpj   
323               DO ji = 1,jpi
324                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj)
325                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj)
326               END DO
327            END DO
328            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r8 )     
329            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r8 ) 
330         ENDIF
331         !
332         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord.
333         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
334         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
335         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
336      ENDIF
337     
338      IF( ln_zco ) THEN
339         !                                                      ! z-coordinate - full steps
340         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
341         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
342         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
343         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
344      ENDIF
345      !                                     ! ============================
346      !                                     !        close the files
347      !                                     ! ============================
348      SELECT CASE ( MOD(nmsh, 3) )
349      CASE ( 1 )               
350         CALL iom_close( inum0 )
351      CASE ( 2 )
352         CALL iom_close( inum1 )
353         CALL iom_close( inum2 )
354      CASE ( 0 )
355         CALL iom_close( inum2 )
356         CALL iom_close( inum3 )
357         CALL iom_close( inum4 )
358      END SELECT
359      !
360      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
361      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
362      !
363      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
364      !
365   END SUBROUTINE dom_wri
366
367
368   SUBROUTINE dom_uniq( puniq, cdgrd )
369      !!----------------------------------------------------------------------
370      !!                  ***  ROUTINE dom_uniq  ***
371      !!                   
372      !! ** Purpose :   identify unique point of a grid (TUVF)
373      !!
374      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
375      !!                2) check which elements have been changed
376      !!----------------------------------------------------------------------
377      !
378      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
379      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
380      !
381      REAL(wp) ::  zshift   ! shift value link to the process number
382      INTEGER  ::  ji       ! dummy loop indices
383      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
384      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
385      !!----------------------------------------------------------------------
386      !
387      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
388      !
389      CALL wrk_alloc( jpi, jpj, ztstref )
390      !
391      ! build an array with different values for each element
392      ! in mpp: make sure that these values are different even between process
393      ! -> apply a shift value according to the process number
394      zshift = jpi * jpj * ( narea - 1 )
395      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
396      !
397      puniq(:,:) = ztstref(:,:)                   ! default definition
398      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
399      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
400      !
401      puniq(:,:) = 1.                             ! default definition
402      ! fill only the inner part of the cpu with llbl converted into real
403      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
404      !
405      CALL wrk_dealloc( jpi, jpj, ztstref )
406      !
407      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
408      !
409   END SUBROUTINE dom_uniq
410
411
412   SUBROUTINE dom_stiff( px1 )
413      !!----------------------------------------------------------------------
414      !!                  ***  ROUTINE dom_stiff  ***
415      !!                     
416      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
417      !!
418      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
419      !!                Save the maximum in the vertical direction
420      !!                (this number is only relevant in s-coordinates)
421      !!
422      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
423      !!----------------------------------------------------------------------
424      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
425      !
426      INTEGER  ::   ji, jj, jk 
427      REAL(wp) ::   zrxmax
428      REAL(wp), DIMENSION(4) ::   zr1
429      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
430      !!----------------------------------------------------------------------
431      zx1(:,:) = 0._wp
432      zrxmax   = 0._wp
433      zr1(:)   = 0._wp
434      !
435      DO ji = 2, jpim1
436         DO jj = 2, jpjm1
437            DO jk = 1, jpkm1
438!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
439!!             especially since it is gde3w which is used to compute the pressure gradient
440!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
441!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
442               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
443                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
444                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
445                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
446               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
447                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
448                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
449                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
450               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
451                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
452                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
453                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
454               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
455                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
456                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
457                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
458               zrxmax = MAXVAL( zr1(1:4) )
459               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
460            END DO
461         END DO
462      END DO
463      CALL lbc_lnk( zx1, 'T', 1. )
464      !
465      IF( PRESENT( px1 ) )    px1 = zx1
466      !
467      zrxmax = MAXVAL( zx1 )
468      !
469      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
470      !
471      IF(lwp) THEN
472         WRITE(numout,*)
473         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
474         WRITE(numout,*) '~~~~~~~~~'
475      ENDIF
476      !
477   END SUBROUTINE dom_stiff
478
479   !!======================================================================
480END MODULE domwri
Note: See TracBrowser for help on using the repository browser.