source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

File size: 19.6 KB
Line 
1MODULE domwri
2   !!======================================================================
3   !!                       ***  MODULE domwri  ***
4   !! Ocean initialization : write the ocean domain mesh file(s)
5   !!======================================================================
6   !! History :  OPA  ! 1997-02  (G. Madec)  Original code
7   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL
8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file
9   !!            3.0  ! 2008-01  (S. Masson) add dom_uniq
10   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dom_wri        : create and write mesh and mask file(s)
15   !!   dom_uniq       :
16   !!----------------------------------------------------------------------
17   USE dom_oce         ! ocean space and time domain
18   USE in_out_manager  ! I/O manager
19   USE iom             ! I/O library
20   USE lbclnk          ! lateral boundary conditions - mpp exchanges
21   USE lib_mpp         ! MPP library
22   USE wrk_nemo        ! Memory allocation
23   USE timing          ! Timing
24
25   USE yomhook, ONLY: lhook, dr_hook
26   USE parkind1, ONLY: jprb, jpim
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC dom_wri, dom_uniq  ! routines called by inidom.F90 and iom.F90
32
33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dom_wri
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dom_wri  ***
45      !!                   
46      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
47      !!      ocean domain informations (mesh and mask arrays). This (these)
48      !!      file(s) is (are) used for visualisation (SAXO software) and
49      !!      diagnostic computation.
50      !!
51      !! ** Method  :   Write in a file all the arrays generated in routines
52      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
53      !!      the vertical coord. used (z-coord, partial steps, s-coord)
54      !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file
55      !!                         = 2  :   'mesh.nc' and mask.nc' files
56      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
57      !!                                  'mask.nc' files
58      !!      For huge size domain, use option 2 or 3 depending on your
59      !!      vertical coordinate.
60      !!
61      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
62      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
63      !!                        corresponding to the depth of the bottom t- and w-points
64      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
65      !!                        thickness (e3[tw]_ps) of the bottom points
66      !!
67      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
68      !!                                   masks, depth and vertical scale factors
69      !!----------------------------------------------------------------------
70      !!
71      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
72      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
73      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
74      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
75      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
76      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
77      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
78      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
79      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
80      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
81      INTEGER           ::   ji, jj, jk   ! dummy loop indices
82      !                                   !  workspaces
83      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
84      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
85      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
86      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
87      REAL(KIND=jprb)               :: zhook_handle
88
89      CHARACTER(LEN=*), PARAMETER :: RoutineName='DOM_WRI'
90
91      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
92
93      !!----------------------------------------------------------------------
94      !
95      IF( nn_timing == 1 )  CALL timing_start('dom_wri')
96      !
97      CALL wrk_alloc( jpi, jpj, zprt, zprw )
98      CALL wrk_alloc( jpi, jpj, jpk, zdepu, zdepv )
99      !
100      IF(lwp) WRITE(numout,*)
101      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
102      IF(lwp) WRITE(numout,*) '~~~~~~~'
103     
104      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
105      clnam1 = 'mesh'       ! filename (mesh informations)
106      clnam2 = 'mask'       ! filename (mask informations)
107      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
108      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
109     
110      SELECT CASE ( MOD(nmsh, 3) )
111         !                                  ! ============================
112      CASE ( 1 )                            !  create 'mesh_mask.nc' file
113         !                                  ! ============================
114         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
115         inum2 = inum0                                            ! put all the informations
116         inum3 = inum0                                            ! in unit inum0
117         inum4 = inum0
118         
119         !                                  ! ============================
120      CASE ( 2 )                            !  create 'mesh.nc' and
121         !                                  !         'mask.nc' files
122         !                                  ! ============================
123         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
124         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
125         inum3 = inum1                                            ! put mesh informations
126         inum4 = inum1                                            ! in unit inum1
127         !                                  ! ============================
128      CASE ( 0 )                            !  create 'mesh_hgr.nc'
129         !                                  !         'mesh_zgr.nc' and
130         !                                  !         'mask.nc'     files
131         !                                  ! ============================
132         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
133         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
134         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
135         !
136      END SELECT
137     
138      !                                                         ! masks (inum2)
139      IF(nn_timing == 2)  CALL timing_start('iom_rstput')
140      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
141      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
142      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
143      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
144      IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
145     
146      CALL dom_uniq( zprw, 'T' )
147      DO jj = 1, jpj
148         DO ji = 1, jpi
149            jk=mikt(ji,jj) 
150            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
151         END DO
152      END DO                             !    ! unique point mask
153      IF(nn_timing == 2)  CALL timing_start('iom_rstput')
154      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
155      IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
156      CALL dom_uniq( zprw, 'U' )
157      DO jj = 1, jpj
158         DO ji = 1, jpi
159            jk=miku(ji,jj) 
160            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
161         END DO
162      END DO
163      IF(nn_timing == 2)  CALL timing_start('iom_rstput')
164      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
165      IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
166      CALL dom_uniq( zprw, 'V' )
167      DO jj = 1, jpj
168         DO ji = 1, jpi
169            jk=mikv(ji,jj) 
170            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
171         END DO
172      END DO
173      IF(nn_timing == 2)  CALL timing_start('iom_rstput')
174      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
175      IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
176      CALL dom_uniq( zprw, 'F' )
177      DO jj = 1, jpj
178         DO ji = 1, jpi
179            jk=mikf(ji,jj) 
180            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
181         END DO
182      END DO
183      IF(nn_timing == 2)  CALL timing_start('iom_rstput')
184      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
185
186      !                                                         ! horizontal mesh (inum3)
187      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
188      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
189      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
190      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
191     
192      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
193      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
194      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
195      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
196     
197      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
198      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
199      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
200      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
201     
202      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
203      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
204      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
205      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
206     
207      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
208     
209      ! note that mbkt is set to 1 over land ==> use surface tmask
210      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
211      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
212      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
213      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points
214      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
215      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points
216           
217      IF( ln_sco ) THEN                                         ! s-coordinate
218         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )
219         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
220         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
221         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
222         !
223         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
224         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
225         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
226         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
227         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
228         !
229         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
230         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
231         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
232         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
233         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio
234         !
235         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
236         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
237         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
238         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )     
239      ENDIF
240      IF(nn_timing == 2)  CALL timing_stop('iom_rstput') 
241      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
242         !
243         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
244            IF(nn_timing == 2)  CALL timing_start('iom_rstput')
245            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
246            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
247            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
248            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
249            IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
250         ELSE                                                   !    ! 2D masked bottom ocean scale factors
251            DO jj = 1,jpj   
252               DO ji = 1,jpi
253                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
254                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
255               END DO
256            END DO
257            IF(nn_timing == 2)  CALL timing_start('iom_rstput')
258            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
259            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
260            IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
261         END IF
262         !
263         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
264            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
265            DO jk = 1,jpk   
266               DO jj = 1, jpjm1   
267                  DO ji = 1, fs_jpim1   ! vector opt.
268                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
269                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
270                  END DO   
271               END DO   
272            END DO
273            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
274            IF(nn_timing == 2)  CALL timing_start('iom_rstput')
275            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
276            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
277            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
278            IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
279         ELSE                                                   !    ! 2D bottom depth
280            DO jj = 1,jpj   
281               DO ji = 1,jpi
282                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj)
283                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj)
284               END DO
285            END DO
286            IF(nn_timing == 2)  CALL timing_start('iom_rstput')
287            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
288            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
289            IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
290         ENDIF
291         !
292         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
293         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord.
294         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
295         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
296         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
297         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
298      ENDIF
299     
300      IF( ln_zco ) THEN
301         !                                                      ! z-coordinate - full steps
302         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
303         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
304         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
305         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
306         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
307         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
308      ENDIF
309      !                                     ! ============================
310      !                                     !        close the files
311      !                                     ! ============================
312      SELECT CASE ( MOD(nmsh, 3) )
313      CASE ( 1 )               
314         CALL iom_close( inum0 )
315      CASE ( 2 )
316         CALL iom_close( inum1 )
317         CALL iom_close( inum2 )
318      CASE ( 0 )
319         CALL iom_close( inum2 )
320         CALL iom_close( inum3 )
321         CALL iom_close( inum4 )
322      END SELECT
323      !
324      CALL wrk_dealloc( jpi, jpj, zprt, zprw )
325      CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
326      !
327      IF( nn_timing == 1 )  CALL timing_stop('dom_wri')
328      !
329      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
330   END SUBROUTINE dom_wri
331
332
333   SUBROUTINE dom_uniq( puniq, cdgrd )
334      !!----------------------------------------------------------------------
335      !!                  ***  ROUTINE dom_uniq  ***
336      !!                   
337      !! ** Purpose :   identify unique point of a grid (TUVF)
338      !!
339      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
340      !!                2) check which elements have been changed
341      !!----------------------------------------------------------------------
342      !
343      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
344      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
345      !
346      REAL(wp) ::  zshift   ! shift value link to the process number
347      INTEGER  ::  ji       ! dummy loop indices
348      LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) ::  lldbl  ! store whether each point is unique or not
349      REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
350      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
351      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
352      REAL(KIND=jprb)               :: zhook_handle
353
354      CHARACTER(LEN=*), PARAMETER :: RoutineName='DOM_UNIQ'
355
356      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
357
358      !!----------------------------------------------------------------------
359      !
360      IF( nn_timing == 1 )  CALL timing_start('dom_uniq')
361      !
362      CALL wrk_alloc( jpi, jpj, ztstref )
363      !
364      ! build an array with different values for each element
365      ! in mpp: make sure that these values are different even between process
366      ! -> apply a shift value according to the process number
367      zshift = jpi * jpj * ( narea - 1 )
368      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
369      !
370      puniq(:,:) = ztstref(:,:)                   ! default definition
371      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
372      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
373      !
374      puniq(:,:) = 1.                             ! default definition
375      ! fill only the inner part of the cpu with llbl converted into real
376      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
377      !
378      CALL wrk_dealloc( jpi, jpj, ztstref )
379      !
380      IF( nn_timing == 1 )  CALL timing_stop('dom_uniq')
381      !
382      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
383   END SUBROUTINE dom_uniq
384
385   !!======================================================================
386END MODULE domwri
Note: See TracBrowser for help on using the repository browser.