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/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC/domwri.F90 @ 8777

Last change on this file since 8777 was 8777, checked in by deazer, 6 years ago

Includes Bug fix for asseling filter and 3d rivers
Includes limiter in surface fluxesin shallow water near WAD points as a temporary measure.
File size: 16.9 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, ln_rwd, 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 7646 2017-02-06 09:25:03Z timgraham $
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      IF( ln_wd  ) CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 )
201      IF( ln_wd .OR. ln_rwd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
202
203      !                                     ! ============================
204      CALL iom_close( inum )                !        close the files
205      !                                     ! ============================
206      !
207      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
208      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
209      !
210      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
211      !
212   END SUBROUTINE dom_wri
213
214
215   SUBROUTINE dom_uniq( puniq, cdgrd )
216      !!----------------------------------------------------------------------
217      !!                  ***  ROUTINE dom_uniq  ***
218      !!                   
219      !! ** Purpose :   identify unique point of a grid (TUVF)
220      !!
221      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
222      !!                2) check which elements have been changed
223      !!----------------------------------------------------------------------
224      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
225      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
226      !
227      REAL(wp) ::  zshift   ! shift value link to the process number
228      INTEGER  ::  ji       ! dummy loop indices
229      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
230      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
231      !!----------------------------------------------------------------------
232      !
233      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
234      !
235      CALL wrk_alloc( jpi, jpj, ztstref )
236      !
237      ! build an array with different values for each element
238      ! in mpp: make sure that these values are different even between process
239      ! -> apply a shift value according to the process number
240      zshift = jpi * jpj * ( narea - 1 )
241      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
242      !
243      puniq(:,:) = ztstref(:,:)                   ! default definition
244      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
245      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
246      !
247      puniq(:,:) = 1.                             ! default definition
248      ! fill only the inner part of the cpu with llbl converted into real
249      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
250      !
251      CALL wrk_dealloc( jpi, jpj, ztstref )
252      !
253      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
254      !
255   END SUBROUTINE dom_uniq
256
257
258   SUBROUTINE dom_stiff( px1 )
259      !!----------------------------------------------------------------------
260      !!                  ***  ROUTINE dom_stiff  ***
261      !!                     
262      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
263      !!
264      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
265      !!                Save the maximum in the vertical direction
266      !!                (this number is only relevant in s-coordinates)
267      !!
268      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
269      !!----------------------------------------------------------------------
270      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
271      !
272      INTEGER  ::   ji, jj, jk 
273      REAL(wp) ::   zrxmax
274      REAL(wp), DIMENSION(4) ::   zr1
275      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
276      !!----------------------------------------------------------------------
277      zx1(:,:) = 0._wp
278      zrxmax   = 0._wp
279      zr1(:)   = 0._wp
280      !
281      DO ji = 2, jpim1
282         DO jj = 2, jpjm1
283            DO jk = 1, jpkm1
284!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
285!!             especially since it is gde3w which is used to compute the pressure gradient
286!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
287!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
288               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
289                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
290                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
291                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
292               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
293                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
294                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
295                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
296               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
297                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
298                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
299                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
300               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
301                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
302                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
303                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
304               zrxmax = MAXVAL( zr1(1:4) )
305               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
306            END DO
307         END DO
308      END DO
309      CALL lbc_lnk( zx1, 'T', 1. )
310      !
311      IF( PRESENT( px1 ) )    px1 = zx1
312      !
313      zrxmax = MAXVAL( zx1 )
314      !
315      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
316      !
317      IF(lwp) THEN
318         WRITE(numout,*)
319         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
320         WRITE(numout,*) '~~~~~~~~~'
321      ENDIF
322      !
323   END SUBROUTINE dom_stiff
324
325   !!======================================================================
326END MODULE domwri
Note: See TracBrowser for help on using the repository browser.