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

Last change on this file since 15153 was 15153, checked in by dbruciaferri, 18 months ago

introducing v1 of saw-tooth patterns diagnostic

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