Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domwri.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago
GMED 462 add flush
File size: 15.3 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 : ll_wd ! Wetting and drying
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
27	IMPLICIT NONE
28	PRIVATE
29
30	PUBLIC dom_wri ! routine called by inidom.F90
31	PUBLIC dom_stiff ! routine called by inidom.F90
32
33	!! * Substitutions
34	# include "vectopt_loop_substitute.h90"
35	!!----------------------------------------------------------------------
36	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
37	!! $Id$
38	!! Software governed by the CeCILL license (see ./LICENSE)
39	!!----------------------------------------------------------------------
40	CONTAINS
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 : create a file with all domain related arrays
52	!!
53	!! ** output file : meshmask.nc : domain size, horizontal grid-point position,
54	!! masks, depth and vertical scale factors
55	!!----------------------------------------------------------------------
56	INTEGER :: inum ! temprary units for 'mesh_mask.nc' file
57	CHARACTER(len=21) :: clnam ! filename (mesh and mask informations)
58	INTEGER :: ji, jj, jk ! dummy loop indices
59	INTEGER :: izco, izps, isco, icav
60	!
61	REAL(wp), DIMENSION(jpi,jpj) :: zprt, zprw ! 2D workspace
62	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepu, zdepv ! 3D workspace
63	!!----------------------------------------------------------------------
64	!
65	IF(lwp) THEN
66	WRITE(numout,*)
67	WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
68	WRITE(numout,*) '~~~~~~~'
69	IF(lflush) CALL FLUSH(numout)
70	ENDIF
71
72	clnam = 'mesh_mask' ! filename (mesh and mask informations)
73
74	! ! ============================
75	! ! create 'mesh_mask.nc' file
76	! ! ============================
77	CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE. )
78	!
79	! ! global domain size
80	CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
81	CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
82	CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
83
84	! ! domain characteristics
85	CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
86	! ! type of vertical coordinate
87	IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
88	IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
89	IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
90	CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
91	CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
92	CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
93	! ! ocean cavities under iceshelves
94	IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF
95	CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
96
97	! ! masks
98	CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask
99	CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
100	CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
101	CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
102
103	CALL dom_uniq( zprw, 'T' )
104	DO jj = 1, jpj
105	DO ji = 1, jpi
106	zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
107	END DO
108	END DO ! ! unique point mask
109	CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )
110	CALL dom_uniq( zprw, 'U' )
111	DO jj = 1, jpj
112	DO ji = 1, jpi
113	zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj) ! ! unique point mask
114	END DO
115	END DO
116	CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )
117	CALL dom_uniq( zprw, 'V' )
118	DO jj = 1, jpj
119	DO ji = 1, jpi
120	zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
121	END DO
122	END DO
123	CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )
124	!!gm ssfmask has been removed ==>> find another solution to defined fmaskutil
125	!! Here we just remove the output of fmaskutil.
126	! CALL dom_uniq( zprw, 'F' )
127	! DO jj = 1, jpj
128	! DO ji = 1, jpi
129	! zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
130	! END DO
131	! END DO
132	! CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )
133	!!gm
134
135	! ! horizontal mesh (inum3)
136	CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude
137	CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
138	CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
139	CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
140
141	CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude
142	CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
143	CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
144	CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
145
146	CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors
147	CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
148	CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
149	CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
150
151	CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors
152	CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
153	CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
154	CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
155
156	CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 ) ! ! coriolis factor
157	CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
158
159	! note that mbkt is set to 1 over land ==> use surface tmask
160	zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
161	CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
162	zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
163	CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
164	zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
165	CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 ) ! ! nb of ocean T-points
166	! ! vertical mesh
167	CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8 ) ! ! scale factors
168	CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8 )
169	CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8 )
170	CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8 )
171	!
172	CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 ) ! stretched system
173	CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
174	CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )
175	CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
176	!
177	IF( ln_sco ) THEN ! s-coordinate stiffness
178	CALL dom_stiff( zprt )
179	CALL iom_rstput( 0, 0, inum, 'stiffness', zprt ) ! Max. grid stiffness ratio
180	ENDIF
181	!
182	IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 )
183
184	! ! ============================
185	CALL iom_close( inum ) ! close the files
186	! ! ============================
187	END SUBROUTINE dom_wri
188
189
190	SUBROUTINE dom_uniq( puniq, cdgrd )
191	!!----------------------------------------------------------------------
192	!! * ROUTINE dom_uniq *
193	!!
194	!! ** Purpose : identify unique point of a grid (TUVF)
195	!!
196	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
197	!! 2) check which elements have been changed
198	!!----------------------------------------------------------------------
199	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
200	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
201	!
202	REAL(wp) :: zshift ! shift value link to the process number
203	INTEGER :: ji ! dummy loop indices
204	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
205	REAL(wp), DIMENSION(jpi,jpj) :: ztstref
206	!!----------------------------------------------------------------------
207	!
208	! build an array with different values for each element
209	! in mpp: make sure that these values are different even between process
210	! -> apply a shift value according to the process number
211	zshift = jpi * jpj * ( narea - 1 )
212	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
213	!
214	puniq(:,:) = ztstref(:,:) ! default definition
215	CALL lbc_lnk( 'domwri', puniq, cdgrd, 1. ) ! apply boundary conditions
216	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
217	!
218	puniq(:,:) = 1. ! default definition
219	! fill only the inner part of the cpu with llbl converted into real
220	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
221	!
222	END SUBROUTINE dom_uniq
223
224
225	SUBROUTINE dom_stiff( px1 )
226	!!----------------------------------------------------------------------
227	!! * ROUTINE dom_stiff *
228	!!
229	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
230	!!
231	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
232	!! Save the maximum in the vertical direction
233	!! (this number is only relevant in s-coordinates)
234	!!
235	!! Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
236	!!----------------------------------------------------------------------
237	REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness
238	!
239	INTEGER :: ji, jj, jk
240	REAL(wp) :: zrxmax
241	REAL(wp), DIMENSION(4) :: zr1
242	REAL(wp), DIMENSION(jpi,jpj) :: zx1
243	!!----------------------------------------------------------------------
244	zx1(:,:) = 0._wp
245	zrxmax = 0._wp
246	zr1(:) = 0._wp
247	!
248	DO ji = 2, jpim1
249	DO jj = 2, jpjm1
250	DO jk = 1, jpkm1
251	!!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation....
252	!! especially since it is gde3w which is used to compute the pressure gradient
253	!! furthermore, I think gdept_0 should be used below instead of w point in the numerator
254	!! so that the ratio is computed at the same point (i.e. uw and vw) ....
255	zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &
256	& +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &
257	& / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &
258	& -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)
259	zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &
260	& +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &
261	& / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &
262	& -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)
263	zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &
264	& +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &
265	& / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &
266	& -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)
267	zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &
268	& +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &
269	& / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &
270	& -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)
271	zrxmax = MAXVAL( zr1(1:4) )
272	zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
273	END DO
274	END DO
275	END DO
276	CALL lbc_lnk( 'domwri', zx1, 'T', 1. )
277	!
278	IF( PRESENT( px1 ) ) px1 = zx1
279	!
280	zrxmax = MAXVAL( zx1 )
281	!
282	CALL mpp_max( 'domwri', zrxmax ) ! max over the global domain
283	!
284	IF(lwp) THEN
285	WRITE(numout,*)
286	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
287	WRITE(numout,*) '~~~~~~~~~'
288	IF(lflush) CALL FLUSH(numout)
289	ENDIF
290	!
291	END SUBROUTINE dom_stiff
292
293	!!======================================================================
294	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: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DOM/domwri.F90 @ 10986

Download in other formats: