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/2013/dev_r3940_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2013/dev_r3940_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 4015

Last change on this file since 4015 was 4015, checked in by cetlod, 11 years ago

2013/dev_r3940_CNRS4_IOCRS: 1st step, add new routines for outputs coarsening

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