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 @ 5601

Last change on this file since 5601 was 5601, checked in by cbricaud, 9 years ago

commit changes/bugfix/... for crs ; ok with time-splitting/fixed volume

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