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.
crsdomwri.F90 in branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 8895

Last change on this file since 8895 was 8895, checked in by andmirek, 6 years ago

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