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_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC – NEMO

source: branches/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC/domwri.F90 @ 7485

Last change on this file since 7485 was 7467, checked in by acc, 8 years ago

Changes to WAD_TEST_CASES/MY_SRC files and namelist_cfg to successfully run with the merge branch. Some changes will be moved into the main source once verified
File size: 16.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   !!            4.0  ! 2016-01  (G. Madec)  simplified mesh_mask.nc file
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dom_wri        : create and write mesh and mask file(s)
15   !!   dom_uniq       : identify unique point of a grid (TUVF)
16   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
17   !!----------------------------------------------------------------------
18   USE dom_oce         ! ocean space and time domain
19   USE phycst ,   ONLY :   rsmall
20   USE wet_dry,   ONLY :   ln_wd, ht_wd
21   !
22   USE in_out_manager  ! I/O manager
23   USE iom             ! I/O library
24   USE lbclnk          ! lateral boundary conditions - mpp exchanges
25   USE lib_mpp         ! MPP library
26   USE wrk_nemo        ! Memory allocation
27   USE timing          ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dom_wri              ! routine called by inidom.F90
33   PUBLIC   dom_stiff            ! routine called by inidom.F90
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
39   !! $Id: domwri.F90 7421 2016-12-01 17:10:41Z flavoni $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE dom_wri
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE dom_wri  ***
47      !!                   
48      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
49      !!      ocean domain informations (mesh and mask arrays). This (these)
50      !!      file(s) is (are) used for visualisation (SAXO software) and
51      !!      diagnostic computation.
52      !!
53      !! ** Method  :   Write in a file all the arrays generated in routines
54      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
55      !!      the vertical coord. used (z-coord, partial steps, s-coord)
56      !!            MOD(nn_msh, 3) = 1  :   'mesh_mask.nc' file
57      !!                         = 2  :   'mesh.nc' and mask.nc' files
58      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
59      !!                                  'mask.nc' files
60      !!      For huge size domain, use option 2 or 3 depending on your
61      !!      vertical coordinate.
62      !!
63      !!      if     nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
64      !!      if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
65      !!                        corresponding to the depth of the bottom t- and w-points
66      !!      if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the
67      !!                        thickness (e3[tw]_ps) of the bottom points
68      !!
69      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
70      !!                                   masks, depth and vertical scale factors
71      !!----------------------------------------------------------------------
72      INTEGER           ::   inum    ! temprary units for 'mesh_mask.nc' file
73      CHARACTER(len=21) ::   clnam   ! filename (mesh and mask informations)
74      INTEGER           ::   ji, jj, jk   ! dummy loop indices
75      INTEGER           ::   izco, izps, isco, icav
76      !                               
77      REAL(wp), POINTER, DIMENSION(:,:)   ::   zprt, zprw     ! 2D workspace
78      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdepu, zdepv   ! 3D workspace
79      !!----------------------------------------------------------------------
80      !
81      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
82      !
83      CALL wrk_alloc( jpi,jpj,       zprt , zprw  )
84      CALL wrk_alloc( jpi,jpj,jpk,   zdepu, zdepv )
85      !
86      IF(lwp) WRITE(numout,*)
87      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
88      IF(lwp) WRITE(numout,*) '~~~~~~~'
89     
90      clnam = 'mesh_mask'  ! filename (mesh and mask informations)
91     
92      !                                  ! ============================
93      !                                  !  create 'mesh_mask.nc' file
94      !                                  ! ============================
95      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
96      !
97      !                                                         ! global domain size
98      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
99      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
100      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
101
102      !                                                         ! domain characteristics
103      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
104      !                                                         ! type of vertical coordinate
105      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
106      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
107      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
108      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
109      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
110      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
111      !                                                         ! ocean cavities under iceshelves
112      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
113      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
114 
115      !                                                         ! masks
116      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
117      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
118      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
119      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
120     
121      CALL dom_uniq( zprw, 'T' )
122      DO jj = 1, jpj
123         DO ji = 1, jpi
124            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
125         END DO
126      END DO                             !    ! unique point mask
127      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 ) 
128      CALL dom_uniq( zprw, 'U' )
129      DO jj = 1, jpj
130         DO ji = 1, jpi
131            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
132         END DO
133      END DO
134      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 ) 
135      CALL dom_uniq( zprw, 'V' )
136      DO jj = 1, jpj
137         DO ji = 1, jpi
138            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
139         END DO
140      END DO
141      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 ) 
142!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
143!!    Here we just remove the output of fmaskutil.
144!      CALL dom_uniq( zprw, 'F' )
145!      DO jj = 1, jpj
146!         DO ji = 1, jpi
147!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
148!         END DO
149!      END DO
150!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 ) 
151!!gm
152
153      !                                                         ! horizontal mesh (inum3)
154      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
155      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
156      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
157      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
158     
159      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
160      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
161      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
162      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
163     
164      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
165      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
166      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
167      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
168     
169      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
170      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
171      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
172      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
173     
174      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor
175      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
176     
177      ! note that mbkt is set to 1 over land ==> use surface tmask
178      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
179      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
180      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
181      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
182      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
183      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
184      !                                                         ! vertical mesh
185      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
186      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
187      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
188      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
189      !
190      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
191      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
192      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
193      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
194      !
195      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
196         CALL dom_stiff( zprt )
197         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )      !    ! Max. grid stiffness ratio
198      ENDIF
199      !                                     ! ============================
200      CALL iom_close( inum )                !        close the files
201      !                                     ! ============================
202      !
203      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
204      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
205      !
206      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
207      !
208   END SUBROUTINE dom_wri
209
210
211   SUBROUTINE dom_uniq( puniq, cdgrd )
212      !!----------------------------------------------------------------------
213      !!                  ***  ROUTINE dom_uniq  ***
214      !!                   
215      !! ** Purpose :   identify unique point of a grid (TUVF)
216      !!
217      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
218      !!                2) check which elements have been changed
219      !!----------------------------------------------------------------------
220      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
221      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
222      !
223      REAL(wp) ::  zshift   ! shift value link to the process number
224      INTEGER  ::  ji       ! dummy loop indices
225      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
226      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
227      !!----------------------------------------------------------------------
228      !
229      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
230      !
231      CALL wrk_alloc( jpi, jpj, ztstref )
232      !
233      ! build an array with different values for each element
234      ! in mpp: make sure that these values are different even between process
235      ! -> apply a shift value according to the process number
236      zshift = jpi * jpj * ( narea - 1 )
237      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
238      !
239      puniq(:,:) = ztstref(:,:)                   ! default definition
240      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
241      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
242      !
243      puniq(:,:) = 1.                             ! default definition
244      ! fill only the inner part of the cpu with llbl converted into real
245      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
246      !
247      CALL wrk_dealloc( jpi, jpj, ztstref )
248      !
249      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
250      !
251   END SUBROUTINE dom_uniq
252
253
254   SUBROUTINE dom_stiff( px1 )
255      !!----------------------------------------------------------------------
256      !!                  ***  ROUTINE dom_stiff  ***
257      !!                     
258      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
259      !!
260      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
261      !!                Save the maximum in the vertical direction
262      !!                (this number is only relevant in s-coordinates)
263      !!
264      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
265      !!----------------------------------------------------------------------
266      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
267      !
268      INTEGER  ::   ji, jj, jk 
269      REAL(wp) ::   zrxmax
270      REAL(wp), DIMENSION(4) ::   zr1
271      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
272      !!----------------------------------------------------------------------
273      zx1(:,:) = 0._wp
274      zrxmax   = 0._wp
275      zr1(:)   = 0._wp
276      !
277      DO ji = 2, jpim1
278         DO jj = 2, jpjm1
279            DO jk = 1, jpkm1
280!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
281!!             especially since it is gde3w which is used to compute the pressure gradient
282!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
283!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
284               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
285                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
286                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
287                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
288               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
289                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
290                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
291                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
292               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
293                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
294                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
295                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
296               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
297                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
298                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
299                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
300               zrxmax = MAXVAL( zr1(1:4) )
301               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
302            END DO
303         END DO
304      END DO
305      CALL lbc_lnk( zx1, 'T', 1. )
306      !
307      IF( PRESENT( px1 ) )    px1 = zx1
308      !
309      zrxmax = MAXVAL( zx1 )
310      !
311      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
312      !
313      IF(lwp) THEN
314         WRITE(numout,*)
315         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
316         WRITE(numout,*) '~~~~~~~~~'
317      ENDIF
318      !
319   END SUBROUTINE dom_stiff
320
321   !!======================================================================
322END MODULE domwri
Note: See TracBrowser for help on using the repository browser.