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.
domwri.F90 in tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/DOM – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/DOM/domwri.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 15.7 KB
Line 
1MODULE domwri
2   !!======================================================================
3   !!                       ***  MODULE domwri  ***
4   !! Ocean initialization : write the ocean domain mesh ask file(s)
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_wri        : create and write mesh and mask file(s)
9   !!                    nmsh = 1  :   mesh_mask file
10   !!                         = 2  :   mesh and mask file
11   !!                         = 3  :   mesh_hgr, mesh_zgr and mask
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE dom_oce         ! ocean space and time domain
15   USE in_out_manager
16   USE iom
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18   USE lib_mpp
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC dom_wri        ! routine called by inidom.F90
24
25   !! * Substitutions
26#  include "vectopt_loop_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LOCEAN-IPSL (2005)
29   !! $Id: domwri.F90 1590 2009-08-06 10:18:30Z smasson $
30   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
31   !!----------------------------------------------------------------------
32
33CONTAINS
34
35   SUBROUTINE dom_wri
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE dom_wri  ***
38      !!                   
39      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
40      !!      ocean domain informations (mesh and mask arrays). This (these)
41      !!      file(s) is (are) used for visualisation (SAXO software) and
42      !!      diagnostic computation.
43      !!
44      !! ** Method  :   Write in a file all the arrays generated in routines
45      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
46      !!      the vertical coord. used (z-coord, partial steps, s-coord)
47      !!                    nmsh = 1  :   'mesh_mask.nc' file
48      !!                         = 2  :   'mesh.nc' and mask.nc' files
49      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
50      !!                                  'mask.nc' files
51      !!      For huge size domain, use option 2 or 3 depending on your
52      !!      vertical coordinate.
53      !!
54      !! ** output file :
55      !!      meshmask.nc  : domain size, horizontal grid-point position,
56      !!                     masks, depth and vertical scale factors
57      !!
58      !! History :
59      !!        !  97-02  (G. Madec)  Original code
60      !!        !  99-11  (M. Imbard)  NetCDF FORMAT with IOIPSL
61      !!   9.0  !  02-08  (G. Madec)  F90 and several file
62      !!----------------------------------------------------------------------
63      INTEGER                          ::   inum0    ! temprary units for 'mesh_mask.nc' file
64      INTEGER                          ::   inum1    ! temprary units for 'mesh.nc'      file
65      INTEGER                          ::   inum2    ! temprary units for 'mask.nc'      file
66      INTEGER                          ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
67      INTEGER                          ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
68      INTEGER                          ::   ji, jj, jk, ik
69      REAL(wp), DIMENSION(jpi,jpj)     ::   zprt     ! temporary array for bathymetry
70      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu    ! 3D depth of U point
71      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepv    ! 3D depth of V point
72      CHARACTER(len=21)                ::   clnam0   ! filename (mesh and mask informations)
73      CHARACTER(len=21)                ::   clnam1   ! filename (mesh informations)
74      CHARACTER(len=21)                ::   clnam2   ! filename (mask informations)
75      CHARACTER(len=21)                ::   clnam3   ! filename (horizontal mesh informations)
76      CHARACTER(len=21)                ::   clnam4   ! filename (vertical   mesh informations)
77      !!----------------------------------------------------------------------
78
79      IF(lwp) WRITE(numout,*)
80      IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
81      IF(lwp) WRITE(numout,*) '~~~~~~~'
82     
83      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
84      clnam1 = 'mesh'       ! filename (mesh informations)
85      clnam2 = 'mask'       ! filename (mask informations)
86      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
87      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
88     
89      SELECT CASE ( MOD(nmsh, 3) )
90         !                                  ! ============================
91      CASE ( 1 )                            !  create 'mesh_mask.nc' file
92         !                                  ! ============================
93         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
94         inum2 = inum0                                            ! put all the informations
95         inum3 = inum0                                            ! in unit inum0
96         inum4 = inum0
97         
98         !                                  ! ============================
99      CASE ( 2 )                            !  create 'mesh.nc' and
100         !                                  !         'mask.nc' files
101         !                                  ! ============================
102         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
103         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
104         inum3 = inum1                                            ! put mesh informations
105         inum4 = inum1                                            ! in unit inum1
106         !                                  ! ============================
107      CASE ( 0 )                            !  create 'mesh_hgr.nc'
108         !                                  !         'mesh_zgr.nc' and
109         !                                  !         'mask.nc'     files
110         !                                  ! ============================
111         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
112         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
113         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
114         
115      END SELECT
116     
117      !                                                         ! masks (inum2)
118      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
119      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
120      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
121      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
122     
123     
124      zprt = tmask(:,:,1) * dom_uniq('T')                               !    ! unique point mask
125      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 
126      zprt = umask(:,:,1) * dom_uniq('U')
127      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 
128      zprt = vmask(:,:,1) * dom_uniq('V')
129      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 
130      zprt = fmask(:,:,1) * dom_uniq('F')
131      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 
132
133      !                                                         ! horizontal mesh (inum3)
134      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
135      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
136      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
137      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
138     
139      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
140      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
141      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
142      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
143     
144      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
145      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
146      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
147      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
148     
149      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
150      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
151      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
152      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
153     
154      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
155     
156!       note that mbathy has been modified in dommsk or in solver.
157!       it is the number of non-zero "w" levels in the water, and the minimum
158!       value (on land) is 2. We define zprt as the number of "T" points in the ocean
159!       at any location, and zero on land.
160!
161      zprt = tmask(:,:,1)*(mbathy-1)
162      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )
163           
164#if ! defined key_zco
165      IF( ln_sco ) THEN                                         ! s-coordinate
166         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth
167         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
168         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
169         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
170         
171         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef.
172         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 
173         CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
174         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
175         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
176         
177         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors
178         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
179         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
180         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
181         
182         CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system
183         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 )
184      ENDIF
185     
186      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
187
188         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors
189            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )         
190            CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
191            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
192            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
193         ELSE                                                   !    ! 2D bottom scale factors
194            DO jj = 1,jpj   ;   DO ji = 1,jpi
195               ik = NINT( zprt(ji,jj) )   ! take care that mbathy is not what you think it is here !
196               IF ( ik /= 0 ) THEN   ;   e3tp(ji,jj) = e3t(ji,jj,ik)   ;   e3wp(ji,jj) = e3w(ji,jj,ik)
197               ELSE                  ;   e3tp(ji,jj) = 0.              ;   e3wp(ji,jj) = 0.
198               ENDIF
199            END DO   ;   END DO
200            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )     
201            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
202         END IF
203
204         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
205            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )     
206            DO jk = 1,jpk   ;   DO jj = 1, jpjm1   ;   DO ji = 1, fs_jpim1   ! vector opt.
207               zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji+1,jj  ,jk) )
208               zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji  ,jj+1,jk) )
209            END DO   ;   END DO   ;   END DO
210            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
211            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
212            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
213            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
214         ELSE                                                   !    ! 2D bottom depth
215            DO jj = 1,jpj   ;   DO ji = 1,jpi
216               ik = NINT( zprt(ji,jj) )   ! take care that mbathy is not what you think it is here !
217               IF ( ik /= 0 ) THEN   ;   hdept(ji,jj) = gdept(ji,jj,ik)   ;   hdepw(ji,jj) = gdepw(ji,jj,ik+1)
218               ELSE                  ;   hdept(ji,jj) = 0.                ;   hdepw(ji,jj) = 0.
219               ENDIF
220            END DO   ;   END DO
221            CALL iom_rstput( 0, 0, inum4, 'hdept' , hdept, ktype = jp_r4 )     
222            CALL iom_rstput( 0, 0, inum4, 'hdepw' , hdepw, ktype = jp_r4 ) 
223         ENDIF
224
225         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! reference z-coord.
226         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
227         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )
228         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
229      ENDIF
230     
231#endif
232     
233      IF( ln_zco ) THEN
234         !                                                      ! z-coordinate - full steps
235         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! depth
236         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
237         CALL iom_rstput( 0, 0, inum4, 'e3t_0'  , e3t_0   )     !    ! scale factors
238         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   )
239      ENDIF
240      !                                     ! ============================
241      !                                     !        close the files
242      !                                     ! ============================
243      SELECT CASE ( nmsh )
244      CASE ( 1 )               
245         CALL iom_close( inum0 )
246      CASE ( 2 )
247         CALL iom_close( inum1 )
248         CALL iom_close( inum2 )
249      CASE ( 3 )
250         CALL iom_close( inum2 )
251         CALL iom_close( inum3 )
252         CALL iom_close( inum4 )
253      END SELECT
254     
255   END SUBROUTINE dom_wri
256
257
258   FUNCTION dom_uniq( cdgrd )   RESULT( puniq )
259      !!----------------------------------------------------------------------
260      !!                  ***  ROUTINE dom_uniq  ***
261      !!                   
262      !! ** Purpose :   identify unique point of a grid (TUVF)
263      !!
264      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
265      !!                2) check which elements have been changed
266      !!
267      !!----------------------------------------------------------------------
268      CHARACTER(len=1)            , INTENT(in   ) ::  cdgrd   !
269      REAL(wp), DIMENSION(jpi,jpj)                ::  puniq   !
270
271      REAL(wp), DIMENSION(jpi,jpj  ) ::  ztstref   ! array with different values for each element
272      REAL(wp)                       ::  zshift    ! shift value link to the process number
273      LOGICAL , DIMENSION(jpi,jpj,1) ::  lldbl     ! is the point unique or not?
274      INTEGER                        ::  ji        ! dummy loop indices
275      !!----------------------------------------------------------------------
276
277      ! build an array with different values for each element
278      ! in mpp: make sure that these values are different even between process
279      ! -> apply a shift value according to the process number
280      zshift = jpi * jpj * ( narea - 1 )
281      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji, wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
282   
283      puniq(:,:) = ztstref(:,:)                   ! default definition
284      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions
285      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed
286     
287      puniq(:,:) = 1.                             ! default definition
288      ! fill only the inner part of the cpu with llbl converted into real
289      puniq(nldi:nlei,nldj:nlej) = REAL(COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ), wp)
290
291   END FUNCTION dom_uniq
292
293
294   !!======================================================================
295END MODULE domwri
Note: See TracBrowser for help on using the repository browser.