source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

  • Property svn:keywords set to Id
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
35   !! $Id$
36CONTAINS
37
38   SUBROUTINE crs_dom_wri
39      !!----------------------------------------------------------------------
40      !!                  ***  ROUTINE crs_dom_wri  ***
41      !!
42      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
43      !!      ocean domain informations (mesh and mask arrays). This (these)
44      !!      file(s) is (are) used for visualisation (SAXO software) and
45      !!      diagnostic computation.
46      !!
47      !! ** Method  :   Write in a file all the arrays generated in routines
48      !!      crsini for meshes and mask. In three separate files:
49      !!      domain size, horizontal grid-point position,
50      !!      masks, depth and vertical scale factors
51      !!     
52      !! ** Output files :   mesh_hgr_crs.nc, mesh_zgr_crs.nc, mesh_mask.nc
53      !!----------------------------------------------------------------------
54      !!
55      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
56      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
57      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
58      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
59      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
60      INTEGER           ::   iif, iil, ijf, ijl
61      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
62      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
63      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
64      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
65      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
66      INTEGER           ::   ji, jj, jk   ! dummy loop indices
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      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      CALL wrk_dealloc( jpi_crs, jpj_crs,      zprt , zprw  )
297      CALL wrk_dealloc( jpi_crs, jpj_crs,      ze3tp, ze3wp )
298      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zdepu, zdepv )
299      !
300      IF( nn_timing == 1 )  CALL timing_stop('crs_dom_wri')
301      !
302     
303   END SUBROUTINE crs_dom_wri
304
305
306   SUBROUTINE dom_uniq_crs( puniq, cdgrd )
307      !!----------------------------------------------------------------------
308      !!                  ***  ROUTINE crs_dom_uniq_crs  ***
309      !!                   
310      !! ** Purpose :   identify unique point of a grid (TUVF)
311      !!
312      !! ** Method  :   1) apply crs_lbc_lnk on an array with different values for each element
313      !!                2) check which elements have been changed
314      !!----------------------------------------------------------------------
315      !
316      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
317      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
318      !
319      REAL(wp) ::  zshift   ! shift value link to the process number
320      INTEGER  ::  ji       ! dummy loop indices
321      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
322      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
323      !!----------------------------------------------------------------------
324      !
325      IF( nn_timing == 1 )  CALL timing_start('crs_dom_uniq_crs')
326      !
327      CALL wrk_alloc( jpi_crs, jpj_crs, ztstref )
328      !
329      ! build an array with different values for each element
330      ! in mpp: make sure that these values are different even between process
331      ! -> apply a shift value according to the process number
332      zshift = jpi_crs * jpj_crs * ( narea - 1 )
333      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi_crs*jpj_crs) /), (/ jpi_crs, jpj_crs /) )
334      !
335      puniq(:,:) = ztstref(:,:)                   ! default definition
336      CALL crs_lbc_lnk( puniq,cdgrd, 1. )            ! apply boundary conditions
337      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
338      !
339      puniq(:,:) = 1.                             ! default definition
340      ! fill only the inner part of the cpu with llbl converted into real
341      puniq(nldi_crs:nlei_crs,nldj_crs:nlej_crs) = REAL( COUNT( lldbl(nldi_crs:nlei_crs,nldj_crs:nlej_crs,:), dim = 3 ) , wp )
342      !
343      CALL wrk_dealloc( jpi_crs, jpj_crs, ztstref )
344      !
345      IF( nn_timing == 1 )  CALL timing_stop('crs_dom_uniq_crs')
346      !
347     
348   END SUBROUTINE dom_uniq_crs
349
350   !!======================================================================
351
352END MODULE crsdomwri
353
354
Note: See TracBrowser for help on using the repository browser.