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.

domain.F90 in branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 6624

Last change on this file since 6624 was 6624, checked in by gm, 9 years ago
#1692 - branch SIMPLIF_2_usrdef: add domain_cfg.nc file which includes jperio, and vert. coord. logicals. +code cleaning
Property svn:keywords set to `Id`
File size: 34.2 KB

Line
1	MODULE domain
2	!!==============================================================================
3	!! * MODULE domain *
4	!! Ocean initialization : domain initialization
5	!!==============================================================================
6	!! History : OPA ! 1990-10 (C. Levy - G. Madec) Original code
7	!! ! 1992-01 (M. Imbard) insert time step initialization
8	!! ! 1996-06 (G. Madec) generalized vertical coordinate
9	!! ! 1997-02 (G. Madec) creation of domwri.F
10	!! ! 2001-05 (E.Durand - G. Madec) insert closed sea
11	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
12	!! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization
13	!! 3.3 ! 2010-11 (G. Madec) initialisation in C1D configuration
14	!! 3.6 ! 2013 ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
15	!! 3.7 ! 2015-11 (G. Madec, A. Coward) time varying zgr by default
16	!!----------------------------------------------------------------------
17
18	!!----------------------------------------------------------------------
19	!! dom_init : initialize the space and time domain
20	!! dom_nam : read and contral domain namelists
21	!! dom_ctl : control print for the ocean domain
22	!! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
23	!! cfg_wri : create the "domain_cfg.nc" file containing all required configuration information
24	!!----------------------------------------------------------------------
25	USE oce ! ocean variables
26	USE dom_oce ! domain: ocean
27	USE sbc_oce ! surface boundary condition: ocean
28	USE phycst ! physical constants
29	USE closea ! closed seas
30	USE domhgr ! domain: set the horizontal mesh
31	USE domzgr ! domain: set the vertical mesh
32	USE dommsk ! domain: set the mask system
33	USE domwri ! domain: write the meshmask file
34	USE domvvl ! variable volume
35	USE c1d ! 1D vertical configuration
36	USE dyncor_c1d ! Coriolis term (c1d case) (cor_c1d routine)
37	!
38	USE in_out_manager ! I/O manager
39	USE iom ! I/O library
40	USE lbclnk ! ocean lateral boundary condition (or mpp link)
41	USE lib_mpp ! distributed memory computing library
42	USE wrk_nemo ! Memory Allocation
43	USE timing ! Timing
44
45	IMPLICIT NONE
46	PRIVATE
47
48	PUBLIC dom_init ! called by opa.F90
49
50	!!-------------------------------------------------------------------------
51	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
52	!! $Id$
53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54	!!-------------------------------------------------------------------------
55	CONTAINS
56
57	SUBROUTINE dom_init
58	!!----------------------------------------------------------------------
59	!! * ROUTINE dom_init *
60	!!
61	!! ** Purpose : Domain initialization. Call the routines that are
62	!! required to create the arrays which define the space
63	!! and time domain of the ocean model.
64	!!
65	!! ** Method : - dom_msk: compute the masks from the bathymetry file
66	!! - dom_hgr: compute or read the horizontal grid-point position
67	!! and scale factors, and the coriolis factor
68	!! - dom_zgr: define the vertical coordinate and the bathymetry
69	!! - dom_wri: create the meshmask file if nmsh=1
70	!! - 1D configuration, move Coriolis, u and v at T-point
71	!!----------------------------------------------------------------------
72	INTEGER :: jk ! dummy loop indices
73	INTEGER :: iconf = 0 ! local integers
74	CHARACTER (len=64) :: cform = "(A12, 3(A13, I7) )"
75	REAL(wp), POINTER, DIMENSION(:,:) :: z1_hu_0, z1_hv_0
76	!!----------------------------------------------------------------------
77	!
78	IF( nn_timing == 1 ) CALL timing_start('dom_init')
79	!
80	IF(lwp) THEN
81	WRITE(numout,*)
82	WRITE(numout,*) 'dom_init : domain initialization'
83	WRITE(numout,*) '~~~~~~~~'
84	ENDIF
85	!
86	! Ocean domain Parameters (control print)
87	! -----------------------
88	IF(lwp) THEN
89	WRITE(numout,*) ' Domain info'
90	WRITE(numout,*) ' dimension of model'
91	WRITE(numout,*) ' Local domain Global domain Data domain '
92	WRITE(numout,cform) ' ',' jpi : ', jpi, ' jpiglo : ', jpiglo, ' jpidta : ', jpidta
93	WRITE(numout,cform) ' ',' jpj : ', jpj, ' jpjglo : ', jpjglo, ' jpjdta : ', jpjdta
94	WRITE(numout,cform) ' ',' jpk : ', jpk, ' jpkglo : ', jpkglo, ' jpkdta : ', jpkdta
95	WRITE(numout,cform) ' ' ,' jpij : ', jpij
96	WRITE(numout,*) ' mpp local domain info (mpp)'
97	WRITE(numout,*) ' jpni : ', jpni, ' jpreci : ', jpreci
98	WRITE(numout,*) ' jpnj : ', jpnj, ' jprecj : ', jprecj
99	WRITE(numout,*) ' jpnij : ', jpnij
100	WRITE(numout,*) ' lateral boundary of the Global domain : jperio = ', jperio
101	ENDIF
102
103	! !== Reference coordinate system ==!
104	!
105	CALL dom_nam ! read namelist ( namrun, namdom )
106	CALL dom_clo ! Closed seas and lake
107	CALL dom_hgr ! Horizontal mesh
108	CALL dom_zgr ! Vertical mesh and bathymetry
109	CALL dom_msk ! Masks
110	!
111	IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency
112	!
113	ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1) ! Reference ocean thickness
114	hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
115	hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
116	DO jk = 2, jpk
117	ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
118	hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
119	hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
120	END DO
121	!
122	! !== time varying part of coordinate system ==!
123	!
124	IF( ln_linssh ) THEN ! Fix in time : set to the reference one for all
125	! before ! now ! after !
126	; gdept_b = gdept_0 ; gdept_n = gdept_0 ! --- ! depth of grid-points
127	; gdepw_b = gdepw_0 ; gdepw_n = gdepw_0 ! --- !
128	; ; gde3w_n = gde3w_0 ! --- !
129	!
130	; e3t_b = e3t_0 ; e3t_n = e3t_0 ; e3t_a = e3t_0 ! scale factors
131	; e3u_b = e3u_0 ; e3u_n = e3u_0 ; e3u_a = e3u_0 !
132	; e3v_b = e3v_0 ; e3v_n = e3v_0 ; e3v_a = e3v_0 !
133	; ; e3f_n = e3f_0 ! --- !
134	; e3w_b = e3w_0 ; e3w_n = e3w_0 ! --- !
135	; e3uw_b = e3uw_0 ; e3uw_n = e3uw_0 ! --- !
136	; e3vw_b = e3vw_0 ; e3vw_n = e3vw_0 ! --- !
137	!
138	CALL wrk_alloc( jpi,jpj, z1_hu_0, z1_hv_0 )
139	!
140	z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF
141	z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
142	!
143	! before ! now ! after !
144	; ; ht_n = ht_0 ! ! water column thickness
145	; hu_b = hu_0 ; hu_n = hu_0 ; hu_a = hu_0 !
146	; hv_b = hv_0 ; hv_n = hv_0 ; hv_a = hv_0 !
147	; r1_hu_b = z1_hu_0 ; r1_hu_n = z1_hu_0 ; r1_hu_a = z1_hu_0 ! inverse of water column thickness
148	; r1_hv_b = z1_hv_0 ; r1_hv_n = z1_hv_0 ; r1_hv_a = z1_hv_0 !
149	!
150	CALL wrk_dealloc( jpi,jpj, z1_hu_0, z1_hv_0 )
151	!
152	ELSE ! time varying : initialize before/now/after variables
153	!
154	CALL dom_vvl_init
155	!
156	ENDIF
157	!
158	IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point
159	!
160	IF( nmsh > 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file
161	IF( nmsh > 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file
162	IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control
163	!
164
165	IF(lwp) THEN
166	WRITE(numout,*)
167	WRITE(numout,*) 'dom_init : end of domain initialization nmsh=', nmsh
168	WRITE(numout,*)
169	ENDIF
170	!
171	IF( nmsh == -1 ) CALL cfg_wri ! create the configuration file
172	!
173	IF( nn_timing == 1 ) CALL timing_stop('dom_init')
174	!
175	END SUBROUTINE dom_init
176
177
178	SUBROUTINE dom_nam
179	!!----------------------------------------------------------------------
180	!! * ROUTINE dom_nam *
181	!!
182	!! ** Purpose : read domaine namelists and print the variables.
183	!!
184	!! ** input : - namrun namelist
185	!! - namdom namelist
186	!! - namnc4 namelist ! "key_netcdf4" only
187	!!----------------------------------------------------------------------
188	USE ioipsl
189	NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, &
190	nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , &
191	& nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , &
192	& nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , &
193	& ln_cfmeta, ln_iscpl
194	NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
195	& rn_atfp , rn_rdt , nn_closea , ln_crs , &
196	& ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
197	& ppa2, ppkth2, ppacr2
198	#if defined key_netcdf4
199	NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
200	#endif
201	INTEGER :: ios ! Local integer output status for namelist read
202	!!----------------------------------------------------------------------
203
204	REWIND( numnam_ref ) ! Namelist namrun in reference namelist : Parameters of the run
205	READ ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
206	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
207
208	REWIND( numnam_cfg ) ! Namelist namrun in configuration namelist : Parameters of the run
209	READ ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
210	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
211	IF(lwm) WRITE ( numond, namrun )
212	!
213	IF(lwp) THEN ! control print
214	WRITE(numout,*)
215	WRITE(numout,*) 'dom_nam : domain initialization through namelist read'
216	WRITE(numout,*) '~~~~~~~ '
217	WRITE(numout,*) ' Namelist namrun'
218	WRITE(numout,*) ' job number nn_no = ', nn_no
219	WRITE(numout,*) ' experiment name for output cn_exp = ', cn_exp
220	WRITE(numout,*) ' file prefix restart input cn_ocerst_in= ', cn_ocerst_in
221	WRITE(numout,*) ' restart input directory cn_ocerst_indir= ', cn_ocerst_indir
222	WRITE(numout,*) ' file prefix restart output cn_ocerst_out= ', cn_ocerst_out
223	WRITE(numout,*) ' restart output directory cn_ocerst_outdir= ', cn_ocerst_outdir
224	WRITE(numout,*) ' restart logical ln_rstart = ', ln_rstart
225	WRITE(numout,*) ' start with forward time step nn_euler = ', nn_euler
226	WRITE(numout,*) ' control of time step nn_rstctl = ', nn_rstctl
227	WRITE(numout,*) ' number of the first time step nn_it000 = ', nn_it000
228	WRITE(numout,*) ' number of the last time step nn_itend = ', nn_itend
229	WRITE(numout,*) ' initial calendar date aammjj nn_date0 = ', nn_date0
230	WRITE(numout,*) ' initial time of day in hhmm nn_time0 = ', nn_time0
231	WRITE(numout,*) ' leap year calendar (0/1) nn_leapy = ', nn_leapy
232	WRITE(numout,*) ' initial state output nn_istate = ', nn_istate
233	IF( ln_rst_list ) THEN
234	WRITE(numout,*) ' list of restart dump times nn_stocklist =', nn_stocklist
235	ELSE
236	WRITE(numout,*) ' frequency of restart file nn_stock = ', nn_stock
237	ENDIF
238	WRITE(numout,*) ' frequency of output file nn_write = ', nn_write
239	WRITE(numout,*) ' mask land points ln_mskland = ', ln_mskland
240	WRITE(numout,*) ' additional CF standard metadata ln_cfmeta = ', ln_cfmeta
241	WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber
242	WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz
243	WRITE(numout,*) ' IS coupling at the restart step ln_iscpl = ', ln_iscpl
244	ENDIF
245
246	no = nn_no ! conversion DOCTOR names into model names (this should disappear soon)
247	cexper = cn_exp
248	nrstdt = nn_rstctl
249	nit000 = nn_it000
250	nitend = nn_itend
251	ndate0 = nn_date0
252	nleapy = nn_leapy
253	ninist = nn_istate
254	nstock = nn_stock
255	nstocklist = nn_stocklist
256	nwrite = nn_write
257	neuler = nn_euler
258	IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
259	WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
260	CALL ctl_warn( ctmp1 )
261	neuler = 0
262	ENDIF
263
264	! ! control of output frequency
265	IF ( nstock == 0 .OR. nstock > nitend ) THEN
266	WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
267	CALL ctl_warn( ctmp1 )
268	nstock = nitend
269	ENDIF
270	IF ( nwrite == 0 ) THEN
271	WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
272	CALL ctl_warn( ctmp1 )
273	nwrite = nitend
274	ENDIF
275
276	#if defined key_agrif
277	IF( Agrif_Root() ) THEN
278	#endif
279	SELECT CASE ( nleapy ) ! Choose calendar for IOIPSL
280	CASE ( 1 )
281	CALL ioconf_calendar('gregorian')
282	IF(lwp) WRITE(numout,*) ' The IOIPSL calendar is "gregorian", i.e. leap year'
283	CASE ( 0 )
284	CALL ioconf_calendar('noleap')
285	IF(lwp) WRITE(numout,*) ' The IOIPSL calendar is "noleap", i.e. no leap year'
286	CASE ( 30 )
287	CALL ioconf_calendar('360d')
288	IF(lwp) WRITE(numout,*) ' The IOIPSL calendar is "360d", i.e. 360 days in a year'
289	END SELECT
290	#if defined key_agrif
291	ENDIF
292	#endif
293
294	REWIND( numnam_ref ) ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
295	READ ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
296	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
297
298	!
299	REWIND( numnam_cfg ) ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
300	READ ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
301	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
302	IF(lwm) WRITE ( numond, namdom )
303	!
304	IF(lwp) THEN
305	WRITE(numout,*)
306	WRITE(numout,*) ' Namelist namdom : space & time domain'
307	WRITE(numout,*) ' flag read/compute bathymetry nn_bathy = ', nn_bathy
308	WRITE(numout,*) ' Depth (if =0 bathy=jpkm1) rn_bathy = ', rn_bathy
309	WRITE(numout,*) ' min depth of the ocean (>0) or rn_hmin = ', rn_hmin
310	WRITE(numout,*) ' min number of ocean level (<0) '
311	WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)'
312	WRITE(numout,*) ' minimum thickness of partial rn_e3zps_min = ', rn_e3zps_min, ' (m)'
313	WRITE(numout,*) ' step level rn_e3zps_rat = ', rn_e3zps_rat
314	WRITE(numout,*) ' create mesh/mask file(s) nn_msh = ', nn_msh
315	WRITE(numout,*) ' = 0 no file created '
316	WRITE(numout,*) ' = 1 mesh_mask '
317	WRITE(numout,*) ' = 2 mesh and mask '
318	WRITE(numout,*) ' = 3 mesh_hgr, msh_zgr and mask'
319	WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt
320	WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp
321	WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea
322	WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs
323
324	WRITE(numout,*) ' ORCA r4, r2 and r05 coefficients ppsur = ', ppsur
325	WRITE(numout,*) ' ppa0 = ', ppa0
326	WRITE(numout,*) ' ppa1 = ', ppa1
327	WRITE(numout,*) ' ppkth = ', ppkth
328	WRITE(numout,*) ' ppacr = ', ppacr
329	WRITE(numout,*) ' Minimum vertical spacing ppdzmin = ', ppdzmin
330	WRITE(numout,*) ' Maximum depth pphmax = ', pphmax
331	WRITE(numout,*) ' Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
332	WRITE(numout,*) ' Double tanh function parameters ppa2 = ', ppa2
333	WRITE(numout,*) ' ppkth2 = ', ppkth2
334	WRITE(numout,*) ' ppacr2 = ', ppacr2
335	ENDIF
336
337	call flush( numout )
338	!
339	ntopo = nn_bathy ! conversion DOCTOR names into model names (this should disappear soon)
340	e3zps_min = rn_e3zps_min
341	e3zps_rat = rn_e3zps_rat
342	nmsh = nn_msh
343	atfp = rn_atfp
344	rdt = rn_rdt
345
346	#if defined key_netcdf4
347	! ! NetCDF 4 case ("key_netcdf4" defined)
348	REWIND( numnam_ref ) ! Namelist namnc4 in reference namelist : NETCDF
349	READ ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
350	907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
351
352	REWIND( numnam_cfg ) ! Namelist namnc4 in configuration namelist : NETCDF
353	READ ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
354	908 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
355	IF(lwm) WRITE( numond, namnc4 )
356
357	IF(lwp) THEN ! control print
358	WRITE(numout,*)
359	WRITE(numout,*) ' Namelist namnc4 - Netcdf4 chunking parameters'
360	WRITE(numout,*) ' number of chunks in i-dimension nn_nchunks_i = ', nn_nchunks_i
361	WRITE(numout,*) ' number of chunks in j-dimension nn_nchunks_j = ', nn_nchunks_j
362	WRITE(numout,*) ' number of chunks in k-dimension nn_nchunks_k = ', nn_nchunks_k
363	WRITE(numout,*) ' apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
364	ENDIF
365
366	! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
367	! Note the chunk size in the unlimited (time) dimension will be fixed at 1
368	snc4set%ni = nn_nchunks_i
369	snc4set%nj = nn_nchunks_j
370	snc4set%nk = nn_nchunks_k
371	snc4set%luse = ln_nc4zip
372	#else
373	snc4set%luse = .FALSE. ! No NetCDF 4 case
374	#endif
375
376	call flush( numout )
377
378	!
379	END SUBROUTINE dom_nam
380
381
382	SUBROUTINE dom_ctl
383	!!----------------------------------------------------------------------
384	!! * ROUTINE dom_ctl *
385	!!
386	!! ** Purpose : Domain control.
387	!!
388	!! ** Method : compute and print extrema of masked scale factors
389	!!----------------------------------------------------------------------
390	INTEGER :: iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
391	INTEGER, DIMENSION(2) :: iloc !
392	REAL(wp) :: ze1min, ze1max, ze2min, ze2max
393	!!----------------------------------------------------------------------
394	!
395	IF(lk_mpp) THEN
396	CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
397	CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
398	CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
399	CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
400	ELSE
401	ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
402	ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
403	ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
404	ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
405
406	iloc = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
407	iimi1 = iloc(1) + nimpp - 1
408	ijmi1 = iloc(2) + njmpp - 1
409	iloc = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
410	iimi2 = iloc(1) + nimpp - 1
411	ijmi2 = iloc(2) + njmpp - 1
412	iloc = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
413	iima1 = iloc(1) + nimpp - 1
414	ijma1 = iloc(2) + njmpp - 1
415	iloc = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
416	iima2 = iloc(1) + nimpp - 1
417	ijma2 = iloc(2) + njmpp - 1
418	ENDIF
419	IF(lwp) THEN
420	WRITE(numout,*)
421	WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
422	WRITE(numout,*) '~~~~~~~'
423	WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
424	WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
425	WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
426	WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
427	ENDIF
428	!
429	END SUBROUTINE dom_ctl
430
431
432	SUBROUTINE dom_stiff
433	!!----------------------------------------------------------------------
434	!! * ROUTINE dom_stiff *
435	!!
436	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
437	!!
438	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
439	!! Save the maximum in the vertical direction
440	!! (this number is only relevant in s-coordinates)
441	!!
442	!! Haney, R. L., 1991: On the pressure gradient force
443	!! over steep topography in sigma coordinate ocean models.
444	!! J. Phys. Oceanogr., 21, 610-619.
445	!!----------------------------------------------------------------------
446	INTEGER :: ji, jj, jk
447	REAL(wp) :: zrxmax
448	REAL(wp), DIMENSION(4) :: zr1
449	!!----------------------------------------------------------------------
450	rx1(:,:) = 0._wp
451	zrxmax = 0._wp
452	zr1(:) = 0._wp
453	!
454	DO ji = 2, jpim1
455	DO jj = 2, jpjm1
456	DO jk = 1, jpkm1
457	zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &
458	& +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &
459	& / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &
460	& -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)
461	zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &
462	& +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &
463	& / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &
464	& -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)
465	zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &
466	& +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &
467	& / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &
468	& -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)
469	zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &
470	& +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &
471	& / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &
472	& -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)
473	zrxmax = MAXVAL( zr1(1:4) )
474	rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
475	END DO
476	END DO
477	END DO
478	CALL lbc_lnk( rx1, 'T', 1. )
479	!
480	zrxmax = MAXVAL( rx1 )
481	!
482	IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain
483	!
484	IF(lwp) THEN
485	WRITE(numout,*)
486	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
487	WRITE(numout,*) '~~~~~~~~~'
488	ENDIF
489	!
490	END SUBROUTINE dom_stiff
491
492
493	SUBROUTINE cfg_wri
494	!!----------------------------------------------------------------------
495	!! * ROUTINE cfg_wri *
496	!!
497	!! ** Purpose : Create the NetCDF file(s) which contain(s) all the
498	!! ocean domain informations (mesh and mask arrays). This (these)
499	!! file(s) is (are) used for visualisation (SAXO software) and
500	!! diagnostic computation.
501	!!
502	!! ** Method : Write in a file all the arrays generated in routines
503	!! domhgr, domzgr, and dommsk. Note: the file contain depends on
504	!! the vertical coord. used (z-coord, partial steps, s-coord)
505	!! MOD(nmsh, 3) = 1 : 'mesh_mask.nc' file
506	!! = 2 : 'mesh.nc' and mask.nc' files
507	!! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and
508	!! 'mask.nc' files
509	!! For huge size domain, use option 2 or 3 depending on your
510	!! vertical coordinate.
511	!!
512	!! if nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
513	!! if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
514	!! corresponding to the depth of the bottom t- and w-points
515	!! if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
516	!! thickness (e3[tw]_ps) of the bottom points
517	!!
518	!! ** output file : meshmask.nc : domain size, horizontal grid-point position,
519	!! masks, depth and vertical scale factors
520	!!----------------------------------------------------------------------
521	INTEGER :: ji, jj, jk ! dummy loop indices
522	INTEGER :: izco, izps, isco, icav
523	INTEGER :: inum ! temprary units for 'domain_cfg.nc' file
524	CHARACTER(len=21) :: clnam ! filename (mesh and mask informations)
525	!!----------------------------------------------------------------------
526	!
527	IF(lwp) WRITE(numout,*)
528	IF(lwp) WRITE(numout,*) 'cfg_wri : create the "domain_cfg.nc" file containing all required configuration information'
529	IF(lwp) WRITE(numout,*) '~~~~~~~'
530	!
531	! ! ============================= !
532	! ! create 'domain_cfg.nc' file !
533	! ! ============================= !
534	!
535	clnam = 'domain_cfg' ! filename (configuration information)
536	CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
537
538	! !== global domain size ==!
539	CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
540	CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
541	CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk , wp), ktype = jp_i4 )
542
543	! !== domain characteristics ==!
544	! ! lateral boundary of the global domain
545	CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
546	! ! type of vertical coordinate
547	IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
548	IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
549	IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
550	CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
551	CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
552	CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
553	! ! ocean cavities under iceshelves
554	IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF
555	CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
556
557	! !== horizontal mesh !
558	!
559	CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! latitude
560	CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
561	CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
562	CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
563	!
564	CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! longitude
565	CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
566	CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
567	CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
568	!
569	CALL iom_rstput( 0, 0, inum, 'e1t' , e1t , ktype = jp_r8 ) ! i-scale factors (e1.)
570	CALL iom_rstput( 0, 0, inum, 'e1u' , e1u , ktype = jp_r8 )
571	CALL iom_rstput( 0, 0, inum, 'e1v' , e1v , ktype = jp_r8 )
572	CALL iom_rstput( 0, 0, inum, 'e1f' , e1f , ktype = jp_r8 )
573	!
574	CALL iom_rstput( 0, 0, inum, 'e2t' , e2t , ktype = jp_r8 ) ! j-scale factors (e2.)
575	CALL iom_rstput( 0, 0, inum, 'e2u' , e2u , ktype = jp_r8 )
576	CALL iom_rstput( 0, 0, inum, 'e2v' , e2v , ktype = jp_r8 )
577	CALL iom_rstput( 0, 0, inum, 'e2f' , e2f , ktype = jp_r8 )
578	!
579	CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 ) ! coriolis factor
580	CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 )
581
582
583	! !== vertical mesh - 3D mask ==!
584	!
585	CALL iom_rstput( 0, 0, inum, 'gdept_1d', gdept_1d, ktype = jp_r8 ) ! reference 1D-coordinate
586	CALL iom_rstput( 0, 0, inum, 'gdepw_1d', gdepw_1d, ktype = jp_r8 )
587	CALL iom_rstput( 0, 0, inum, 'e3t_1d' , e3t_1d , ktype = jp_r8 )
588	CALL iom_rstput( 0, 0, inum, 'e3w_1d' , e3w_1d , ktype = jp_r8 )
589	!
590	CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 ) ! depth (t- & w-points)
591	CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
592	!
593	CALL iom_rstput( 0, 0, inum, 'e3t_0' , e3t_0 , ktype = jp_r8 ) ! vertical scale factors (e
594	CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8 )
595	CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8 )
596	CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8 )
597	!
598	CALL iom_rstput( 0, 0, inum, 'tmask' , tmask , ktype = jp_i1 ) ! masks (in bytes)
599	CALL iom_rstput( 0, 0, inum, 'umask' , umask , ktype = jp_i1 )
600	CALL iom_rstput( 0, 0, inum, 'vmask' , vmask , ktype = jp_i1 )
601	CALL iom_rstput( 0, 0, inum, 'fmask' , fmask , ktype = jp_i1 )
602
603	!!gm Probably not required fields :
604	CALL iom_rstput( 0, 0, inum, 'bathy' , bathy , ktype = jp_r8 ) ! depth of the ocean at T-points
605	CALL iom_rstput( 0, 0, inum, 'mbathy' , REAL( mbathy , wp ) , ktype = jp_i4 ) ! nb of ocean T-points
606	CALL iom_rstput( 0, 0, inum, 'mbkt' , REAL( mbkt , wp ) , ktype = jp_i4 ) ! nb of ocean T-points
607
608	!
609	CALL iom_rstput( 0, 0, inum, 'bathy' , risfdep , ktype = jp_r8 ) ! depth of the iceshelves at T-points
610	CALL iom_rstput( 0, 0, inum, 'misfdep', REAL( misfdep, wp ) , ktype = jp_i4 ) ! nb of ocean T-points (ISF)
611
612	!!gm end
613
614	!!gm ?
615	! CALL iom_rstput( 0, 0, inum, 'hbatt', hbatt )
616	! CALL iom_rstput( 0, 0, inum, 'hbatu', hbatu )
617	! CALL iom_rstput( 0, 0, inum, 'hbatv', hbatv )
618	! CALL iom_rstput( 0, 0, inum, 'hbatf', hbatf )
619	!!gm ?
620	!!gm ?
621	! CALL iom_rstput( 0, 0, inum, 'rx1' , rx1 ) ! Max. grid stiffness ratio
622	!!gm ?
623	!
624	! DO jk = 1,jpk
625	! DO jj = 1, jpjm1
626	! DO ji = 1, fs_jpim1 ! vector opt.
627	! zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk) )
628	! zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk) )
629	! END DO
630	! END DO
631	! END DO
632	! CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. )
633	! CALL iom_rstput( 0, 0, inum, 'gdepu' , zdepu , ktype = jp_r8 )
634	! CALL iom_rstput( 0, 0, inum, 'gdepv' , zdepv , ktype = jp_r8 )
635
636	! ! ============================
637	! ! close the files
638	! ! ============================
639	CALL iom_close( inum )
640	!
641	END SUBROUTINE cfg_wri
642
643	!!======================================================================
644	END MODULE domain

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

Download in other formats: