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/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 9124

Last change on this file since 9124 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

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