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/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications
  • Property svn:keywords set to Id
File size: 17.7 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  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dom_wri        : create and write mesh and mask file(s)
15   !!   dom_uniq       :
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
30   !! * Substitutions
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
34   !! $Id$
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE dom_wri
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE dom_wri  ***
42      !!                   
43      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
44      !!      ocean domain informations (mesh and mask arrays). This (these)
45      !!      file(s) is (are) used for visualisation (SAXO software) and
46      !!      diagnostic computation.
47      !!
48      !! ** Method  :   Write in a file all the arrays generated in routines
49      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
50      !!      the vertical coord. used (z-coord, partial steps, s-coord)
51      !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file
52      !!                         = 2  :   'mesh.nc' and mask.nc' files
53      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
54      !!                                  'mask.nc' files
55      !!      For huge size domain, use option 2 or 3 depending on your
56      !!      vertical coordinate.
57      !!
58      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
59      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
60      !!                        corresponding to the depth of the bottom t- and w-points
61      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
62      !!                        thickness (e3[tw]_ps) of the bottom points
63      !!
64      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
65      !!                                   masks, depth and vertical scale factors
66      !!----------------------------------------------------------------------
67      !!
68      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
69      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
70      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
71      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
72      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
73      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
74      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
75      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
76      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
77      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
78      INTEGER           ::   ji, jj, jk   ! dummy loop indices
79      !                                   !  workspaces
80      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
81      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
85      !
86      CALL wrk_alloc( jpi, jpj, zprt, zprw )
87      CALL wrk_alloc( jpi, jpj, jpk, zdepu, zdepv )
88      !
89      IF(lwp) WRITE(numout,*)
90      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
91      IF(lwp) WRITE(numout,*) '~~~~~~~'
92     
93      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
94      clnam1 = 'mesh'       ! filename (mesh informations)
95      clnam2 = 'mask'       ! filename (mask informations)
96      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
97      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
98     
99      SELECT CASE ( MOD(nmsh, 3) )
100         !                                  ! ============================
101      CASE ( 1 )                            !  create 'mesh_mask.nc' file
102         !                                  ! ============================
103         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
104         inum2 = inum0                                            ! put all the informations
105         inum3 = inum0                                            ! in unit inum0
106         inum4 = inum0
107         
108         !                                  ! ============================
109      CASE ( 2 )                            !  create 'mesh.nc' and
110         !                                  !         'mask.nc' files
111         !                                  ! ============================
112         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
113         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
114         inum3 = inum1                                            ! put mesh informations
115         inum4 = inum1                                            ! in unit inum1
116         !                                  ! ============================
117      CASE ( 0 )                            !  create 'mesh_hgr.nc'
118         !                                  !         'mesh_zgr.nc' and
119         !                                  !         'mask.nc'     files
120         !                                  ! ============================
121         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
122         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
123         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
124         !
125      END SELECT
126     
127      !                                                         ! masks (inum2)
128      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
129      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
130      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
131      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
132     
133      CALL dom_uniq( zprw, 'T' )
134      zprt = tmask(:,:,1) * zprw                               !    ! unique point mask
135      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
136      CALL dom_uniq( zprw, 'U' )
137      zprt = umask(:,:,1) * zprw
138      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
139      CALL dom_uniq( zprw, 'V' )
140      zprt = vmask(:,:,1) * zprw
141      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
142      CALL dom_uniq( zprw, 'F' )
143      zprt = fmask(:,:,1) * zprw
144      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
145
146      !                                                         ! horizontal mesh (inum3)
147      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
148      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
149      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
150      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
151     
152      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
153      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
154      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
155      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
156     
157      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
158      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
159      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
160      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
161     
162      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
163      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
164      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
165      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
166     
167      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
168     
169      ! note that mbkt is set to 1 over land ==> use surface tmask
170      zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp )
171      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
172           
173      IF( ln_sco ) THEN                                         ! s-coordinate
174         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth
175         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
176         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
177         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
178         !
179#if defined key_smsh
180          CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt3 )         !    ! scaling coef.
181          CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw3 )
182          CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w3 )
183          CALL iom_rstput( 0, 0, inum4, 'esigt', esigt3 )
184          CALL iom_rstput( 0, 0, inum4, 'esigw', esigw3 )
185#else
186          CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
187          CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )
188          CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
189          CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
190          CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
191#endif
192
193         !
194         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors
195         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
196         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
197         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
198         !
199#if defined key_smsh
200            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )     
201            DO jk = 1,jpk   
202               DO jj = 1, jpjm1   
203                  DO ji = 1, fs_jpim1   ! vector opt.
204                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) )
205                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) )
206                  END DO   
207               END DO   
208            END DO
209            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
210            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
211            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
212            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
213            DO jj = 1,jpj   
214               DO ji = 1,jpi
215                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
216                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
217               END DO
218            END DO
219            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
220            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
221#endif
222         CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system
223         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 )
224      ENDIF
225     
226      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
227         !
228         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
229            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )         
230            CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
231            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
232            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
233         ELSE                                                   !    ! 2D masked bottom ocean scale factors
234            DO jj = 1,jpj   
235               DO ji = 1,jpi
236                  e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
237                  e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
238               END DO
239            END DO
240            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
241            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
242         END IF
243         !
244         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
245            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )     
246            DO jk = 1,jpk   
247               DO jj = 1, jpjm1   
248                  DO ji = 1, fs_jpim1   ! vector opt.
249                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) )
250                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) )
251                  END DO   
252               END DO   
253            END DO
254            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
255            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
256            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
257            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
258         ELSE                                                   !    ! 2D bottom depth
259            DO jj = 1,jpj   
260               DO ji = 1,jpi
261                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
262                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
263               END DO
264            END DO
265            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
266            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
267         ENDIF
268         !
269         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! reference z-coord.
270         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
271         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )
272         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
273      ENDIF
274     
275      IF( ln_zco ) THEN
276         !                                                      ! z-coordinate - full steps
277         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! depth
278         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
279         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )     !    ! scale factors
280         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
281      ENDIF
282      !                                     ! ============================
283      !                                     !        close the files
284      !                                     ! ============================
285      SELECT CASE ( MOD(nmsh, 3) )
286      CASE ( 1 )               
287         CALL iom_close( inum0 )
288      CASE ( 2 )
289         CALL iom_close( inum1 )
290         CALL iom_close( inum2 )
291      CASE ( 0 )
292         CALL iom_close( inum2 )
293         CALL iom_close( inum3 )
294         CALL iom_close( inum4 )
295      END SELECT
296      !
297      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
298      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
299      !
300      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
301      !
302   END SUBROUTINE dom_wri
303
304
305   SUBROUTINE dom_uniq( puniq, cdgrd )
306      !!----------------------------------------------------------------------
307      !!                  ***  ROUTINE dom_uniq  ***
308      !!                   
309      !! ** Purpose :   identify unique point of a grid (TUVF)
310      !!
311      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
312      !!                2) check which elements have been changed
313      !!----------------------------------------------------------------------
314      !
315      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
316      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
317      !
318      REAL(wp) ::  zshift   ! shift value link to the process number
319      INTEGER  ::  ji       ! dummy loop indices
320      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
321      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
322      !!----------------------------------------------------------------------
323      !
324      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
325      !
326      CALL wrk_alloc( jpi, jpj, ztstref )
327      !
328      ! build an array with different values for each element
329      ! in mpp: make sure that these values are different even between process
330      ! -> apply a shift value according to the process number
331      zshift = jpi * jpj * ( narea - 1 )
332      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
333      !
334      puniq(:,:) = ztstref(:,:)                   ! default definition
335      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
336      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
337      !
338      puniq(:,:) = 1.                             ! default definition
339      ! fill only the inner part of the cpu with llbl converted into real
340      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
341      !
342      CALL wrk_dealloc( jpi, jpj, ztstref )
343      !
344      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
345      !
346   END SUBROUTINE dom_uniq
347
348   !!======================================================================
349END MODULE domwri
Note: See TracBrowser for help on using the repository browser.