Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domwri.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.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( nn_timing == 1 ) CALL timing_start('dom_wri')
81	!
82	!
83	IF(lwp) WRITE(numout,*)
84	IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
85	IF(lwp) WRITE(numout,*) '~~~~~~~'
86
87	clnam = 'mesh_mask' ! filename (mesh and mask informations)
88
89	! ! ============================
90	! ! create 'mesh_mask.nc' file
91	! ! ============================
92	CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
93	!
94	! ! global domain size
95	CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
96	CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
97	CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
98
99	! ! domain characteristics
100	CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
101	! ! type of vertical coordinate
102	IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
103	IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
104	IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
105	CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
106	CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
107	CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
108	! ! ocean cavities under iceshelves
109	IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF
110	CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
111
112	! ! masks
113	CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask
114	CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
115	CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
116	CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
117
118	CALL dom_uniq( zprw, 'T' )
119	DO jj = 1, jpj
120	DO ji = 1, jpi
121	zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
122	END DO
123	END DO ! ! unique point mask
124	CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )
125	CALL dom_uniq( zprw, 'U' )
126	DO jj = 1, jpj
127	DO ji = 1, jpi
128	zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj) ! ! unique point mask
129	END DO
130	END DO
131	CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )
132	CALL dom_uniq( zprw, 'V' )
133	DO jj = 1, jpj
134	DO ji = 1, jpi
135	zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
136	END DO
137	END DO
138	CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )
139	!!gm ssfmask has been removed ==>> find another solution to defined fmaskutil
140	!! Here we just remove the output of fmaskutil.
141	! CALL dom_uniq( zprw, 'F' )
142	! DO jj = 1, jpj
143	! DO ji = 1, jpi
144	! zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
145	! END DO
146	! END DO
147	! CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )
148	!!gm
149
150	! ! horizontal mesh (inum3)
151	CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude
152	CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
153	CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
154	CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
155
156	CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude
157	CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
158	CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
159	CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
160
161	CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors
162	CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
163	CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
164	CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
165
166	CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors
167	CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
168	CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
169	CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
170
171	CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 ) ! ! coriolis factor
172	CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
173
174	! note that mbkt is set to 1 over land ==> use surface tmask
175	zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
176	CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
177	zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
178	CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
179	zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
180	CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 ) ! ! nb of ocean T-points
181	! ! vertical mesh
182	CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8 ) ! ! scale factors
183	CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8 )
184	CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8 )
185	CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8 )
186	!
187	CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 ) ! stretched system
188	CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
189	CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )
190	CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
191	!
192	IF( ln_sco ) THEN ! s-coordinate stiffness
193	CALL dom_stiff( zprt )
194	CALL iom_rstput( 0, 0, inum, 'stiffness', zprt ) ! Max. grid stiffness ratio
195	ENDIF
196	!
197	IF( ln_wd ) THEN ! wetting and drying domain
198	CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 )
199	CALL iom_rstput( 0, 0, inum, 'ht_wd' , ht_wd , ktype = jp_r8 )
200	ENDIF
201	! ! ============================
202	CALL iom_close( inum ) ! close the files
203	! ! ============================
204	!
205	!
206	IF( nn_timing == 1 ) CALL timing_stop('dom_wri')
207	!
208	END SUBROUTINE dom_wri
209
210
211	SUBROUTINE dom_uniq( puniq, cdgrd )
212	!!----------------------------------------------------------------------
213	!! * ROUTINE dom_uniq *
214	!!
215	!! ** Purpose : identify unique point of a grid (TUVF)
216	!!
217	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
218	!! 2) check which elements have been changed
219	!!----------------------------------------------------------------------
220	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
221	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
222	!
223	REAL(wp) :: zshift ! shift value link to the process number
224	INTEGER :: ji ! dummy loop indices
225	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
226	REAL(wp), DIMENSION(jpi,jpj) :: ztstref
227	!!----------------------------------------------------------------------
228	!
229	IF( nn_timing == 1 ) CALL timing_start('dom_uniq')
230	!
231	!
232	! build an array with different values for each element
233	! in mpp: make sure that these values are different even between process
234	! -> apply a shift value according to the process number
235	zshift = jpi * jpj * ( narea - 1 )
236	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
237	!
238	puniq(:,:) = ztstref(:,:) ! default definition
239	CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions
240	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
241	!
242	puniq(:,:) = 1. ! default definition
243	! fill only the inner part of the cpu with llbl converted into real
244	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
245	!
246	!
247	IF( nn_timing == 1 ) CALL timing_stop('dom_uniq')
248	!
249	END SUBROUTINE dom_uniq
250
251
252	SUBROUTINE dom_stiff( px1 )
253	!!----------------------------------------------------------------------
254	!! * ROUTINE dom_stiff *
255	!!
256	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
257	!!
258	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
259	!! Save the maximum in the vertical direction
260	!! (this number is only relevant in s-coordinates)
261	!!
262	!! Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
263	!!----------------------------------------------------------------------
264	REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness
265	!
266	INTEGER :: ji, jj, jk
267	REAL(wp) :: zrxmax
268	REAL(wp), DIMENSION(4) :: zr1
269	REAL(wp), DIMENSION(jpi,jpj) :: zx1
270	!!----------------------------------------------------------------------
271	zx1(:,:) = 0._wp
272	zrxmax = 0._wp
273	zr1(:) = 0._wp
274	!
275	DO ji = 2, jpim1
276	DO jj = 2, jpjm1
277	DO jk = 1, jpkm1
278	!!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation....
279	!! especially since it is gde3w which is used to compute the pressure gradient
280	!! furthermore, I think gdept_0 should be used below instead of w point in the numerator
281	!! so that the ratio is computed at the same point (i.e. uw and vw) ....
282	zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &
283	& +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &
284	& / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &
285	& -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)
286	zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &
287	& +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &
288	& / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &
289	& -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)
290	zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &
291	& +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &
292	& / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &
293	& -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)
294	zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &
295	& +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &
296	& / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &
297	& -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)
298	zrxmax = MAXVAL( zr1(1:4) )
299	zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
300	END DO
301	END DO
302	END DO
303	CALL lbc_lnk( zx1, 'T', 1. )
304	!
305	IF( PRESENT( px1 ) ) px1 = zx1
306	!
307	zrxmax = MAXVAL( zx1 )
308	!
309	IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain
310	!
311	IF(lwp) THEN
312	WRITE(numout,*)
313	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
314	WRITE(numout,*) '~~~~~~~~~'
315	ENDIF
316	!
317	END SUBROUTINE dom_stiff
318
319	!!======================================================================
320	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_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 7910

Download in other formats: