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 @ 15271

Last change on this file since 15271 was 15271, checked in by dbruciaferri, 17 months ago

rewriting saw tooth dignostic

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