Context Navigation

domwri.F90 @ 9012

Last change on this file since 9012 was 8882, checked in by flavoni, 6 years ago
dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864
Property svn:keywords set to `Id`
File size: 16.6 KB

Line
1	MODULE 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 ! 2016-01 (G. Madec) simplified mesh_mask.nc file
11	!!----------------------------------------------------------------------
12
13	!!----------------------------------------------------------------------
14	!! dom_wri : create and write mesh and mask file(s)
15	!! dom_uniq : identify unique point of a grid (TUVF)
16	!! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
17	!!----------------------------------------------------------------------
18	USE dom_oce ! ocean space and time domain
19	USE phycst , ONLY : rsmall
20	USE wet_dry, ONLY : ln_wd, ht_wd
21	!
22	USE in_out_manager ! I/O manager
23	USE iom ! I/O library
24	USE lbclnk ! lateral boundary conditions - mpp exchanges
25	USE lib_mpp ! MPP library
26	USE timing ! Timing
27
28	IMPLICIT NONE
29	PRIVATE
30
31	PUBLIC dom_wri ! routine called by inidom.F90
32	PUBLIC dom_stiff ! routine called by inidom.F90
33
34	!! * Substitutions
35	# include "vectopt_loop_substitute.h90"
36	!!----------------------------------------------------------------------
37	!! NEMO/OPA 4.0 , NEMO Consortium (2016)
38	!! $Id$
39	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
40	!!----------------------------------------------------------------------
41	CONTAINS
42
43	SUBROUTINE dom_wri
44	!!----------------------------------------------------------------------
45	!! * ROUTINE dom_wri *
46	!!
47	!! ** Purpose : Create the NetCDF file(s) which contain(s) all the
48	!! ocean domain informations (mesh and mask arrays). This (these)
49	!! file(s) is (are) used for visualisation (SAXO software) and
50	!! diagnostic computation.
51	!!
52	!! ** Method : Write in a file all the arrays generated in routines
53	!! domhgr, domzgr, and dommsk. Note: the file contain depends on
54	!! the vertical coord. used (z-coord, partial steps, s-coord)
55	!! MOD(nn_msh, 3) = 1 : 'mesh_mask.nc' file
56	!! = 2 : 'mesh.nc' and mask.nc' files
57	!! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and
58	!! 'mask.nc' files
59	!! For huge size domain, use option 2 or 3 depending on your
60	!! vertical coordinate.
61	!!
62	!! if nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
63	!! if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
64	!! corresponding to the depth of the bottom t- and w-points
65	!! if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the
66	!! thickness (e3[tw]_ps) of the bottom points
67	!!
68	!! ** output file : meshmask.nc : domain size, horizontal grid-point position,
69	!! masks, depth and vertical scale factors
70	!!----------------------------------------------------------------------
71	INTEGER :: inum ! temprary units for 'mesh_mask.nc' file
72	CHARACTER(len=21) :: clnam ! filename (mesh and mask informations)
73	INTEGER :: ji, jj, jk ! dummy loop indices
74	INTEGER :: izco, izps, isco, icav
75	!
76	REAL(wp), DIMENSION(jpi,jpj) :: zprt, zprw ! 2D workspace
77	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepu, zdepv ! 3D workspace
78	!!----------------------------------------------------------------------
79	!
80	IF( ln_timing ) CALL timing_start('dom_wri')
81	!
82	IF(lwp) WRITE(numout,*)
83	IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
84	IF(lwp) WRITE(numout,*) '~~~~~~~'
85
86	clnam = 'mesh_mask' ! filename (mesh and mask informations)
87
88	! ! ============================
89	! ! create 'mesh_mask.nc' file
90	! ! ============================
91	CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
92	!
93	! ! global domain size
94	CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
95	CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
96	CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
97
98	! ! domain characteristics
99	CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
100	! ! type of vertical coordinate
101	IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
102	IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
103	IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
104	CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
105	CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
106	CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
107	! ! ocean cavities under iceshelves
108	IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF
109	CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
110
111	! ! masks
112	CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask
113	CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
114	CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
115	CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
116
117	CALL dom_uniq( zprw, 'T' )
118	DO jj = 1, jpj
119	DO ji = 1, jpi
120	zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
121	END DO
122	END DO ! ! unique point mask
123	CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )
124	CALL dom_uniq( zprw, 'U' )
125	DO jj = 1, jpj
126	DO ji = 1, jpi
127	zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj) ! ! unique point mask
128	END DO
129	END DO
130	CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )
131	CALL dom_uniq( zprw, 'V' )
132	DO jj = 1, jpj
133	DO ji = 1, jpi
134	zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
135	END DO
136	END DO
137	CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )
138	!!gm ssfmask has been removed ==>> find another solution to defined fmaskutil
139	!! Here we just remove the output of fmaskutil.
140	! CALL dom_uniq( zprw, 'F' )
141	! DO jj = 1, jpj
142	! DO ji = 1, jpi
143	! zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
144	! END DO
145	! END DO
146	! CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )
147	!!gm
148
149	! ! horizontal mesh (inum3)
150	CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude
151	CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
152	CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
153	CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
154
155	CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude
156	CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
157	CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
158	CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
159
160	CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors
161	CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
162	CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
163	CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
164
165	CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors
166	CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
167	CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
168	CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
169
170	CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 ) ! ! coriolis factor
171	CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
172
173	! note that mbkt is set to 1 over land ==> use surface tmask
174	zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
175	CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
176	zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
177	CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
178	zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
179	CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 ) ! ! nb of ocean T-points
180	! ! vertical mesh
181	CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8 ) ! ! scale factors
182	CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8 )
183	CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8 )
184	CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8 )
185	!
186	CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 ) ! stretched system
187	CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
188	CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )
189	CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
190	!
191	IF( ln_sco ) THEN ! s-coordinate stiffness
192	CALL dom_stiff( zprt )
193	CALL iom_rstput( 0, 0, inum, 'stiffness', zprt ) ! Max. grid stiffness ratio
194	ENDIF
195	!
196	IF( ln_wd ) THEN ! wetting and drying domain
197	CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 )
198	CALL iom_rstput( 0, 0, inum, 'ht_wd' , ht_wd , ktype = jp_r8 )
199	ENDIF
200	! ! ============================
201	CALL iom_close( inum ) ! close the files
202	! ! ============================
203	!
204	IF( ln_timing ) CALL timing_stop('dom_wri')
205	!
206	END SUBROUTINE dom_wri
207
208
209	SUBROUTINE dom_uniq( puniq, cdgrd )
210	!!----------------------------------------------------------------------
211	!! * ROUTINE dom_uniq *
212	!!
213	!! ** Purpose : identify unique point of a grid (TUVF)
214	!!
215	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
216	!! 2) check which elements have been changed
217	!!----------------------------------------------------------------------
218	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
219	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
220	!
221	REAL(wp) :: zshift ! shift value link to the process number
222	INTEGER :: ji ! dummy loop indices
223	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
224	REAL(wp), DIMENSION(jpi,jpj) :: ztstref
225	!!----------------------------------------------------------------------
226	!
227	IF( ln_timing ) CALL timing_start('dom_uniq')
228	!
229	! build an array with different values for each element
230	! in mpp: make sure that these values are different even between process
231	! -> apply a shift value according to the process number
232	zshift = jpi * jpj * ( narea - 1 )
233	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
234	!
235	puniq(:,:) = ztstref(:,:) ! default definition
236	CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions
237	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
238	!
239	puniq(:,:) = 1. ! default definition
240	! fill only the inner part of the cpu with llbl converted into real
241	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
242	!
243	IF( ln_timing ) CALL timing_stop('dom_uniq')
244	!
245	END SUBROUTINE dom_uniq
246
247
248	SUBROUTINE dom_stiff( px1 )
249	!!----------------------------------------------------------------------
250	!! * ROUTINE dom_stiff *
251	!!
252	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
253	!!
254	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
255	!! Save the maximum in the vertical direction
256	!! (this number is only relevant in s-coordinates)
257	!!
258	!! Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
259	!!----------------------------------------------------------------------
260	REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness
261	!
262	INTEGER :: ji, jj, jk
263	REAL(wp) :: zrxmax
264	REAL(wp), DIMENSION(4) :: zr1
265	REAL(wp), DIMENSION(jpi,jpj) :: zx1
266	!!----------------------------------------------------------------------
267	zx1(:,:) = 0._wp
268	zrxmax = 0._wp
269	zr1(:) = 0._wp
270	!
271	DO ji = 2, jpim1
272	DO jj = 2, jpjm1
273	DO jk = 1, jpkm1
274	!!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation....
275	!! especially since it is gde3w which is used to compute the pressure gradient
276	!! furthermore, I think gdept_0 should be used below instead of w point in the numerator
277	!! so that the ratio is computed at the same point (i.e. uw and vw) ....
278	zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &
279	& +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &
280	& / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &
281	& -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)
282	zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &
283	& +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &
284	& / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &
285	& -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)
286	zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &
287	& +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &
288	& / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &
289	& -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)
290	zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &
291	& +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &
292	& / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &
293	& -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)
294	zrxmax = MAXVAL( zr1(1:4) )
295	zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
296	END DO
297	END DO
298	END DO
299	CALL lbc_lnk( zx1, 'T', 1. )
300	!
301	IF( PRESENT( px1 ) ) px1 = zx1
302	!
303	zrxmax = MAXVAL( zx1 )
304	!
305	IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain
306	!
307	IF(lwp) THEN
308	WRITE(numout,*)
309	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
310	WRITE(numout,*) '~~~~~~~~~'
311	ENDIF
312	!
313	END SUBROUTINE dom_stiff
314
315	!!======================================================================
316	END MODULE domwri

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 9012

Download in other formats: