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/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 6772

Last change on this file since 6772 was 6772, checked in by cbricaud, 8 years ago

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