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/r8395_obs_oper_update/NEMOGCM/TOOLS/DOMAINcfg/src – NEMO

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/TOOLS/DOMAINcfg/src/domwri.f90 @ 11350

Last change on this file since 11350 was 11350, checked in by jcastill, 5 years ago

Clear svn keywords

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