Context Navigation

domwri.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 8 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

Property svn:keywords set to Id

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

Download in other formats: