source: branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 8758

Last change on this file since 8758 was 8758, checked in by acc, 3 years ago

Branch 2017/dev_r8126_ROBUST08_no_ghost. Changes to eliminate ghost rows and columns. Currently the halo width is still fixed as 1 but a single variable (nn_hls) has been introduced for the halo-size in preparation for allowing this to vary. nn_hls replaces jpreci and jprecj. These changes have passed full SETTE tests but iceberg exchanges across the north-fold remain untested (SETTE tests only release bergs in the SO) and will require further attention. Note layout.dat now reports the jpi and jpj values for the reporting processor only.

  • Property svn:keywords set to Id
File size: 17.3 KB
Line 
1MODULE crsdomwri
2   !!======================================================================
3   !! Coarse Ocean initialization : write the coarse ocean domain mesh and mask files
4   !!======================================================================
5   !! History :  OPA  ! 1997-02  (G. Madec)  Original code
6   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL
7   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file
8   !!            3.0  ! 2008-01  (S. Masson) add dom_uniq_crs
9   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
10   !!                 ! 2012-06  (J. Simeon, C. Calone, C Ethe )  Reduced and modified for coarse grid
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   crs_dom_wri    : create and write mesh and mask file(s)
15   !!----------------------------------------------------------------------
16   USE timing          ! Timing
17   USE dom_oce         ! ocean space and time domain
18   USE in_out_manager  ! I/O manager
19   USE par_kind, ONLY: wp
20   USE lib_mpp         ! MPP library
21   USE iom_def
22   USE iom
23   USE crs         ! coarse grid domain
24   USE crsdom         ! coarse grid domain
25   USE crslbclnk       ! crs mediator to lbclnk
26   USE wrk_nemo        ! Working array
27
28
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC crs_dom_wri        ! routine called by crsini.F90
34
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE crs_dom_wri
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE crs_dom_wri  ***
45      !!
46      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
47      !!      ocean domain informations (mesh and mask arrays). This (these)
48      !!      file(s) is (are) used for visualisation (SAXO software) and
49      !!      diagnostic computation.
50      !!
51      !! ** Method  :   Write in a file all the arrays generated in routines
52      !!      crsini for meshes and mask. In three separate files:
53      !!      domain size, horizontal grid-point position,
54      !!      masks, depth and vertical scale factors
55      !!     
56      !! ** Output files :   mesh_hgr_crs.nc, mesh_zgr_crs.nc, mesh_mask.nc
57      !!----------------------------------------------------------------------
58      !!
59      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
60      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
61      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
62      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
63      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
64      INTEGER           ::   iif, iil, ijf, ijl
65      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
66      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
67      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
68      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
69      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
70      INTEGER           ::   ji, jj, jk   ! dummy loop indices
71      !                                   !  workspaces
72      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
73      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
74      REAL(wp), POINTER, DIMENSION(:,:  ) :: ze3tp, ze3wp
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )  CALL timing_start('crs_dom_wri')
78      !
79      CALL wrk_alloc( jpi_crs, jpj_crs,      zprt , zprw  )
80      CALL wrk_alloc( jpi_crs, jpj_crs,      ze3tp, ze3wp )
81      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zdepu, zdepv )
82
83      ze3tp(:,:) = 0.0
84      ze3wp(:,:) = 0.0
85
86      !
87      IF(lwp) WRITE(numout,*)
88      IF(lwp) WRITE(numout,*) 'crs_dom_wri : create NetCDF mesh and mask information file(s)'
89      IF(lwp) WRITE(numout,*) '~~~~~~~'
90     
91      clnam0 = 'mesh_mask_crs'  ! filename (mesh and mask informations)
92      clnam1 = 'mesh_crs'       ! filename (mesh informations)
93      clnam2 = 'mask_crs'       ! filename (mask informations)
94      clnam3 = 'mesh_hgr_crs'   ! filename (horizontal mesh informations)
95      clnam4 = 'mesh_zgr_crs'   ! filename (vertical   mesh informations)
96     
97
98      SELECT CASE ( MOD(nn_msh_crs, 3) )
99         !                                  ! ============================
100      CASE ( 1 )                            !  create 'mesh_mask.nc' file
101         !                                  ! ============================
102         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
103         inum2 = inum0                                            ! put all the informations
104         inum3 = inum0                                            ! in unit inum0
105         inum4 = inum0
106         
107         !                                  ! ============================
108      CASE ( 2 )                            !  create 'mesh.nc' and
109         !                                  !         'mask.nc' files
110         !                                  ! ============================
111         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
112         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
113         inum3 = inum1                                            ! put mesh informations
114         inum4 = inum1                                            ! in unit inum1
115         !                                  ! ============================
116      CASE ( 0 )                            !  create 'mesh_hgr.nc'
117         !                                  !         'mesh_zgr.nc' and
118         !                                  !         'mask.nc'     files
119         !                                  ! ============================
120         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
121         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
122         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
123         !
124      END SELECT
125 
126      !========================================================
127      !                                                         ! masks (inum2)
128      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask_crs, ktype = jp_i1 )     !    ! land-sea mask
129      CALL iom_rstput( 0, 0, inum2, 'umask', umask_crs, ktype = jp_i1 )
130      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask_crs, ktype = jp_i1 )
131      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask_crs, ktype = jp_i1 )
132     
133     
134      tmask_i_crs(:,:) = tmask_crs(:,:,1)
135      iif = nn_hls
136      iil = nlci_crs - nn_hls + 1
137      ijf = nn_hls
138      ijl = nlcj_crs - nn_hls + 1
139     
140      tmask_i_crs( 1:iif ,    :  ) = 0._wp
141      tmask_i_crs(iil:jpi_crs,    :  ) = 0._wp
142      tmask_i_crs(   :   , 1:ijf ) = 0._wp
143      tmask_i_crs(   :   ,ijl:jpj_crs) = 0._wp
144     
145     
146      tpol_crs(1:jpiglo_crs,:) = 1._wp
147      fpol_crs(1:jpiglo_crs,:) = 1._wp
148      IF( jperio == 3 .OR. jperio == 4 ) THEN
149         tpol_crs(jpiglo_crs/2+1:jpiglo_crs,:) = 0._wp
150         fpol_crs(       1      :jpiglo_crs,:) = 0._wp
151         IF( mjg_crs(nlej_crs) == jpiglo_crs ) THEN
152            DO ji = iif+1, iil-1
153               tmask_i_crs(ji,nlej_crs-1) = tmask_i_crs(ji,nlej_crs-1) &
154               & * tpol_crs(mig_crs(ji),1)
155            ENDDO
156         ENDIF
157      ENDIF
158      IF( jperio == 5 .OR. jperio == 6 ) THEN
159         tpol_crs(      1       :jpiglo_crs,:)=0._wp
160         fpol_crs(jpiglo_crs/2+1:jpiglo_crs,:)=0._wp
161      ENDIF
162     
163      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', tmask_i_crs, ktype = jp_i1 )
164                                   !    ! unique point mask
165      CALL dom_uniq_crs( zprw, 'U' )
166      zprt = umask_crs(:,:,1) * zprw
167      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
168      CALL dom_uniq_crs( zprw, 'V' )
169      zprt = vmask_crs(:,:,1) * zprw
170      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
171      CALL dom_uniq_crs( zprw, 'F' )
172      zprt = fmask_crs(:,:,1) * zprw
173      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
174      !========================================================
175      !                                                         ! horizontal mesh (inum3)
176      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt_crs, ktype = jp_r4 )     !    ! latitude
177      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu_crs, ktype = jp_r4 )
178      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv_crs, ktype = jp_r4 )
179      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf_crs, ktype = jp_r4 )
180     
181      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit_crs, ktype = jp_r4 )     !    ! longitude
182      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu_crs, ktype = jp_r4 )
183      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv_crs, ktype = jp_r4 )
184      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif_crs, ktype = jp_r4 )
185     
186      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t_crs, ktype = jp_r8 )         !    ! e1 scale factors
187      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u_crs, ktype = jp_r8 )
188      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v_crs, ktype = jp_r8 )
189      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f_crs, ktype = jp_r8 )
190     
191      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t_crs, ktype = jp_r8 )         !    ! e2 scale factors
192      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u_crs, ktype = jp_r8 )
193      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v_crs, ktype = jp_r8 )
194      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f_crs, ktype = jp_r8 )
195     
196      CALL iom_rstput( 0, 0, inum3, 'ff', ff_crs, ktype = jp_r8 )           !    ! coriolis factor
197
198      !========================================================
199      !                                                         ! vertical mesh (inum4)
200!     ! note that mbkt is set to 1 over land ==> use surface tmask_crs
201      zprt(:,:) = tmask_crs(:,:,1) * REAL( mbkt_crs(:,:) , wp )
202      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
203
204      IF( ln_zps ) THEN                       ! z-coordinate - partial steps
205
206           
207         IF ( nn_msh_crs <= 6 ) THEN
208            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t_crs )     
209            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w_crs )     
210            CALL iom_rstput( 0, 0, inum4, 'e3u', e3u_crs )     
211            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v_crs )     
212         ELSE
213            DO jj = 1,jpj_crs   
214               DO ji = 1,jpi_crs
215                  ze3tp(ji,jj) = e3t_crs(ji,jj,mbkt_crs(ji,jj)) * tmask_crs(ji,jj,1)
216                  ze3wp(ji,jj) = e3w_crs(ji,jj,mbkt_crs(ji,jj)) * tmask_crs(ji,jj,1)
217               END DO
218            END DO
219
220            CALL crs_lbc_lnk( ze3tp,'T', 1.0 )
221            CALL crs_lbc_lnk( ze3wp,'W', 1.0 )
222 
223            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', ze3tp )     
224            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', ze3wp )
225         ENDIF
226
227         IF ( nn_msh_crs <= 3 ) THEN
228            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept_crs, ktype = jp_r4 ) 
229            DO jk = 1,jpk   
230               DO jj = 1, jpj_crsm1   
231                  DO ji = 1, jpi_crsm1  ! jes what to do for fs_jpim1??vector opt.
232                     zdepu(ji,jj,jk) = MIN( gdept_crs(ji,jj,jk) , gdept_crs(ji+1,jj  ,jk) ) * umask_crs(ji,jj,jk)
233                     zdepv(ji,jj,jk) = MIN( gdept_crs(ji,jj,jk) , gdept_crs(ji  ,jj+1,jk) ) * vmask_crs(ji,jj,jk)
234                  END DO   
235               END DO   
236            END DO
237
238            CALL crs_lbc_lnk( zdepu,'U', 1. )   ;   CALL crs_lbc_lnk( zdepv,'V', 1. ) 
239            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
240            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
241            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw_crs, ktype = jp_r4 )
242         ELSE
243            DO jj = 1,jpj_crs   
244               DO ji = 1,jpi_crs
245                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
246                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
247               END DO
248            END DO
249            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
250            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
251         ENDIF
252
253         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )     !    ! reference z-coord.
254         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
255         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
256         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
257
258         CALL iom_rstput(  0, 0, inum4, 'ocean_volume_t', ocean_volume_crs_t ) 
259         CALL iom_rstput(  0, 0, inum4, 'facvol_t' , facvol_t  ) 
260         CALL iom_rstput(  0, 0, inum4, 'facvol_w' , facvol_w  ) 
261         CALL iom_rstput(  0, 0, inum4, 'facsurfu' , facsurfu  ) 
262         CALL iom_rstput(  0, 0, inum4, 'facsurfv' , facsurfv  ) 
263         CALL iom_rstput(  0, 0, inum4, 'e1e2w_msk', e1e2w_msk ) 
264         CALL iom_rstput(  0, 0, inum4, 'e2e3u_msk', e2e3u_msk ) 
265         CALL iom_rstput(  0, 0, inum4, 'e1e3v_msk', e1e3v_msk )
266         CALL iom_rstput(  0, 0, inum4, 'e1e2w'    , e1e2w_crs ) 
267         CALL iom_rstput(  0, 0, inum4, 'e2e3u'    , e2e3u_crs ) 
268         CALL iom_rstput(  0, 0, inum4, 'e1e3v'    , e1e3v_crs )
269         CALL iom_rstput(  0, 0, inum4, 'bt'       , bt_crs    )
270         CALL iom_rstput(  0, 0, inum4, 'r1_bt'    , r1_bt_crs )
271
272         CALL iom_rstput(  0, 0, inum4, 'crs_surfu_wgt', crs_surfu_wgt ) 
273         CALL iom_rstput(  0, 0, inum4, 'crs_surfv_wgt', crs_surfv_wgt ) 
274         CALL iom_rstput(  0, 0, inum4, 'crs_volt_wgt' , crs_volt_wgt  ) 
275
276      ENDIF
277     
278     IF( ln_zco ) THEN
279         !                                                      ! z-coordinate - full steps
280        CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )     !    ! depth
281        CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
282        CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )     !    ! scale factors
283        CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
284     ENDIF
285      !                                     ! ============================
286      !                                     !        close the files
287      !                                     ! ============================
288      SELECT CASE ( MOD(nn_msh_crs, 3) )
289      CASE ( 1 )               
290         CALL iom_close( inum0 )
291      CASE ( 2 )
292         CALL iom_close( inum1 )
293         CALL iom_close( inum2 )
294      CASE ( 0 )
295         CALL iom_close( inum2 )
296         CALL iom_close( inum3 )
297         CALL iom_close( inum4 )
298      END SELECT
299      !
300      CALL wrk_dealloc( jpi_crs, jpj_crs,      zprt , zprw  )
301      CALL wrk_dealloc( jpi_crs, jpj_crs,      ze3tp, ze3wp )
302      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zdepu, zdepv )
303      !
304      IF( nn_timing == 1 )  CALL timing_stop('crs_dom_wri')
305      !
306     
307   END SUBROUTINE crs_dom_wri
308
309
310   SUBROUTINE dom_uniq_crs( puniq, cdgrd )
311      !!----------------------------------------------------------------------
312      !!                  ***  ROUTINE crs_dom_uniq_crs  ***
313      !!                   
314      !! ** Purpose :   identify unique point of a grid (TUVF)
315      !!
316      !! ** Method  :   1) apply crs_lbc_lnk on an array with different values for each element
317      !!                2) check which elements have been changed
318      !!----------------------------------------------------------------------
319      !
320      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
321      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
322      !
323      REAL(wp) ::  zshift   ! shift value link to the process number
324      INTEGER  ::  ji       ! dummy loop indices
325      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
326      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
327      !!----------------------------------------------------------------------
328      !
329      IF( nn_timing == 1 )  CALL timing_start('crs_dom_uniq_crs')
330      !
331      CALL wrk_alloc( jpi_crs, jpj_crs, ztstref )
332      !
333      ! build an array with different values for each element
334      ! in mpp: make sure that these values are different even between process
335      ! -> apply a shift value according to the process number
336      zshift = jpi_crs * jpj_crs * ( narea - 1 )
337      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi_crs*jpj_crs) /), (/ jpi_crs, jpj_crs /) )
338      !
339      puniq(:,:) = ztstref(:,:)                   ! default definition
340      CALL crs_lbc_lnk( puniq,cdgrd, 1. )            ! apply boundary conditions
341      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
342      !
343      puniq(:,:) = 1.                             ! default definition
344      ! fill only the inner part of the cpu with llbl converted into real
345      puniq(nldi_crs:nlei_crs,nldj_crs:nlej_crs) = REAL( COUNT( lldbl(nldi_crs:nlei_crs,nldj_crs:nlej_crs,:), dim = 3 ) , wp )
346      !
347      CALL wrk_dealloc( jpi_crs, jpj_crs, ztstref )
348      !
349      IF( nn_timing == 1 )  CALL timing_stop('crs_dom_uniq_crs')
350      !
351     
352   END SUBROUTINE dom_uniq_crs
353
354   !!======================================================================
355
356END MODULE crsdomwri
357
358
Note: See TracBrowser for help on using the repository browser.