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

source: trunk/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 4528

Last change on this file since 4528 was 4294, checked in by cetlod, 10 years ago

dev_MERGE_2013 : changes to improve compilation

File size: 17.0 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         ELSE
208            DO jj = 1,jpj_crs   
209               DO ji = 1,jpi_crs
210                  ze3tp(ji,jj) = e3t_crs(ji,jj,mbkt_crs(ji,jj)) * tmask_crs(ji,jj,1)
211                  ze3wp(ji,jj) = e3w_crs(ji,jj,mbkt_crs(ji,jj)) * tmask_crs(ji,jj,1)
212               END DO
213            END DO
214
215            CALL crs_lbc_lnk( ze3tp,'T', 1.0 )
216            CALL crs_lbc_lnk( ze3wp,'W', 1.0 )
217 
218            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', ze3tp )     
219            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', ze3wp )
220         ENDIF
221
222         IF ( nn_msh_crs <= 3 ) THEN
223            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept_crs, ktype = jp_r4 ) 
224            DO jk = 1,jpk   
225               DO jj = 1, jpj_crsm1   
226                  DO ji = 1, jpi_crsm1  ! jes what to do for fs_jpim1??vector opt.
227                     zdepu(ji,jj,jk) = MIN( gdept_crs(ji,jj,jk) , gdept_crs(ji+1,jj  ,jk) ) * umask_crs(ji,jj,jk)
228                     zdepv(ji,jj,jk) = MIN( gdept_crs(ji,jj,jk) , gdept_crs(ji  ,jj+1,jk) ) * vmask_crs(ji,jj,jk)
229                  END DO   
230               END DO   
231            END DO
232
233            CALL crs_lbc_lnk( zdepu,'U', 1. )   ;   CALL crs_lbc_lnk( zdepv,'V', 1. ) 
234            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
235            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
236            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw_crs, ktype = jp_r4 )
237         ELSE
238            DO jj = 1,jpj_crs   
239               DO ji = 1,jpi_crs
240                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
241                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
242               END DO
243            END DO
244            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
245            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
246         ENDIF
247
248         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )     !    ! reference z-coord.
249         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
250         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
251         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
252
253         CALL iom_rstput(  0, 0, inum4, 'ocean_volume_t', ocean_volume_crs_t ) 
254         CALL iom_rstput(  0, 0, inum4, 'facvol_t' , facvol_t  ) 
255         CALL iom_rstput(  0, 0, inum4, 'facvol_w' , facvol_w  ) 
256         CALL iom_rstput(  0, 0, inum4, 'facsurfu' , facsurfu  ) 
257         CALL iom_rstput(  0, 0, inum4, 'facsurfv' , facsurfv  ) 
258         CALL iom_rstput(  0, 0, inum4, 'e1e2w_msk', e1e2w_msk ) 
259         CALL iom_rstput(  0, 0, inum4, 'e2e3u_msk', e2e3u_msk ) 
260         CALL iom_rstput(  0, 0, inum4, 'e1e3v_msk', e1e3v_msk )
261         CALL iom_rstput(  0, 0, inum4, 'e1e2w'    , e1e2w_crs ) 
262         CALL iom_rstput(  0, 0, inum4, 'e2e3u'    , e2e3u_crs ) 
263         CALL iom_rstput(  0, 0, inum4, 'e1e3v'    , e1e3v_crs )
264         CALL iom_rstput(  0, 0, inum4, 'bt'       , bt_crs    )
265         CALL iom_rstput(  0, 0, inum4, 'r1_bt'    , r1_bt_crs )
266
267         CALL iom_rstput(  0, 0, inum4, 'crs_surfu_wgt', crs_surfu_wgt ) 
268         CALL iom_rstput(  0, 0, inum4, 'crs_surfv_wgt', crs_surfv_wgt ) 
269         CALL iom_rstput(  0, 0, inum4, 'crs_volt_wgt' , crs_volt_wgt  ) 
270
271      ENDIF
272     
273     IF( ln_zco ) THEN
274         !                                                      ! z-coordinate - full steps
275        CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )     !    ! depth
276        CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
277        CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )     !    ! scale factors
278        CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
279     ENDIF
280      !                                     ! ============================
281      !                                     !        close the files
282      !                                     ! ============================
283      SELECT CASE ( MOD(nn_msh_crs, 3) )
284      CASE ( 1 )               
285         CALL iom_close( inum0 )
286      CASE ( 2 )
287         CALL iom_close( inum1 )
288         CALL iom_close( inum2 )
289      CASE ( 0 )
290         CALL iom_close( inum2 )
291         CALL iom_close( inum3 )
292         CALL iom_close( inum4 )
293      END SELECT
294      !
295      CALL wrk_dealloc( jpi_crs, jpj_crs,      zprt , zprw  )
296      CALL wrk_dealloc( jpi_crs, jpj_crs,      ze3tp, ze3wp )
297      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zdepu, zdepv )
298      !
299      IF( nn_timing == 1 )  CALL timing_stop('crs_dom_wri')
300      !
301     
302   END SUBROUTINE crs_dom_wri
303
304
305   SUBROUTINE dom_uniq_crs( puniq, cdgrd )
306      !!----------------------------------------------------------------------
307      !!                  ***  ROUTINE crs_dom_uniq_crs  ***
308      !!                   
309      !! ** Purpose :   identify unique point of a grid (TUVF)
310      !!
311      !! ** Method  :   1) apply crs_lbc_lnk on an array with different values for each element
312      !!                2) check which elements have been changed
313      !!----------------------------------------------------------------------
314      !
315      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
316      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
317      !
318      REAL(wp) ::  zshift   ! shift value link to the process number
319      INTEGER  ::  ji       ! dummy loop indices
320      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
321      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
322      !!----------------------------------------------------------------------
323      !
324      IF( nn_timing == 1 )  CALL timing_start('crs_dom_uniq_crs')
325      !
326      CALL wrk_alloc( jpi_crs, jpj_crs, ztstref )
327      !
328      ! build an array with different values for each element
329      ! in mpp: make sure that these values are different even between process
330      ! -> apply a shift value according to the process number
331      zshift = jpi_crs * jpj_crs * ( narea - 1 )
332      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi_crs*jpj_crs) /), (/ jpi_crs, jpj_crs /) )
333      !
334      puniq(:,:) = ztstref(:,:)                   ! default definition
335      CALL crs_lbc_lnk( puniq,cdgrd, 1. )            ! apply boundary conditions
336      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
337      !
338      puniq(:,:) = 1.                             ! default definition
339      ! fill only the inner part of the cpu with llbl converted into real
340      puniq(nldi_crs:nlei_crs,nldj_crs:nlej_crs) = REAL( COUNT( lldbl(nldi_crs:nlei_crs,nldj_crs:nlej_crs,:), dim = 3 ) , wp )
341      !
342      CALL wrk_dealloc( jpi_crs, jpj_crs, ztstref )
343      !
344      IF( nn_timing == 1 )  CALL timing_stop('crs_dom_uniq_crs')
345      !
346     
347   END SUBROUTINE dom_uniq_crs
348
349   !!======================================================================
350
351END MODULE crsdomwri
352
353
Note: See TracBrowser for help on using the repository browser.