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 branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 6667

Last change on this file since 6667 was 6667, checked in by gm, 8 years ago

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