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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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