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.

create_meshmask.f90 in utils/tools/SIREN/src – NEMO

Context Navigation

source: utils/tools/SIREN/src/create_meshmask.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 4 years ago
update nemo trunk
File size: 70.1 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! DESCRIPTION:
6	!> @file
7	!> @brief
8	!> This program creates the NetCDF file(s) which contain(s) all the
9	!> ocean domain informations.
10	!> it also permits to create the domain_cfg.nc file (needed to run NEMO v4.0
11	!> and upper), or the mesh_mask file(s).
12	!>
13	!> @details
14	!> @section sec1 method
15	!> bathymetry (and optionally ice shelf draft) is read on input file.<br/>
16	!> horizontal grid-point position, scale factors, and the coriolis factor
17	!> are read in coordinates file or computed.<br/>
18	!> vertical coordinate is defined, and the bathymetry recomputed to fit the
19	!> vertical grid.<br/>
20	!> finally the masks from the bathymetry are computed.
21	!>
22	!> all the variables read and or computed, are writen in one to three file(s) depending on
23	!> output option.
24	!> @note
25	!> the file contain depends on
26	!> the vertical coordinate used (z-coord, partial steps, s-coord)
27	!>
28	!> @section sec2 how to
29	!> USAGE: create_meshmask create_meshmask.nam [-v] [-h]<br/>
30	!> - positional arguments:<br/>
31	!> - create_meshmask.nam<br/>
32	!> namelist of create_meshmask
33	!> @note
34	!> a template of the namelist could be created running (in templates directory):
35	!> @code{.sh}
36	!> python create_templates.py create_meshmask
37	!> @endcode
38	!>
39	!> - optional arguments:<br/>
40	!> - -h, --help<br/>
41	!> show this help message (and exit)<br/>
42	!> - -v, --version<br/>
43	!> show Siren's version (and exit)
44	!>
45	!> @section sec_meshmask create_meshmask.nam
46	!> create_meshmask.nam contains 13 sub-namelists:<br/>
47	!> - namlog to set logger parameters
48	!> - namcfg to set configuration file parameters
49	!> - namsrc to set source files parameters
50	!> - namhgr to set horizontal grid parameters
51	!> - namzgr to set vertical grid parameters
52	!> - namdmin to set minimum depth parameters
53	!> - namzco to set vertical coordinate parameters
54	!> - namzps to set partial step parameters
55	!> - namsco to set sigma or hybrid parameters
56	!> - namlbc to set lateral boundary condition parameters
57	!> - namwd to set wetting and drying parameters
58	!> - namgrd to set grid parameters
59	!> - namout to set output parameters
60	!>
61	!> here after, each sub-namelist parameters is detailed.
62	!> @note
63	!> default values are specified between brackets
64	!>
65	!> @subsection sublog namlog
66	!> the logger sub-namelist parameters are :
67	!>
68	!> - cn_logfile [@a create_meshmask.log]<br/>
69	!> logger filename
70	!>
71	!> - cn_verbosity [@a warning]<br/>
72	!> verbosity level, choose between :
73	!> - trace
74	!> - debug
75	!> - info
76	!> - warning
77	!> - error
78	!> - fatal
79	!> - none
80	!>
81	!> - in_maxerror [@a 5]<br/>
82	!> maximum number of error allowed
83	!>
84	!> @subsection subcfg namcfg
85	!> the configuration sub-namelist parameters are :
86	!>
87	!> - cn_varcfg [@a ./cfg/variable.cfg]<br/>
88	!> path to the variable configuration file.<br/>
89	!> the variable configuration file defines standard name,
90	!> default interpolation method, axis,...
91	!> to be used for some known variables.<br/>
92	!>
93	!> - cn_dimcfg [@a ./cfg/dimension.cfg]<br/>
94	!> path to the dimension configuration file.<br/>
95	!> the dimension configuration file defines dimensions allowed.<br/>
96	!>
97	!> - cn_dumcfg [@a ./cfg/dummy.cfg]<br/>
98	!> path to the useless (dummy) configuration file.<br/>
99	!> the dummy configuration file defines useless
100	!> dimension or variable. these dimension(s) or variable(s) will not be
101	!> processed.<br/>
102	!>
103	!> @subsection subsrc namsrc
104	!> the source grid sub-namelist parameters are :
105	!>
106	!> - cn_bathy [@a ]<br/>
107	!> path to the bathymetry file
108	!> - cn_varbathy [@a ]<br/>
109	!> bathymetry variable name
110	!> - cn_coord [@a ]<br/>
111	!> path to the coordinate file (in_mshhgr=0)
112	!> - cn_isfdep [@a ]<br/>
113	!> iceshelf draft (ln_isfcav=true, see namzgr)
114	!> - cn_varisfdep [@a isfdraft]<br/>
115	!> iceshelf draft variable name (ln_isfcav=true, see namzgr)
116	!> - in_perio [@a ]<br/>
117	!> NEMO periodicity
118	!> - ln_closea [@a .TRUE.]<br/>
119	!> logical to fill closed sea or not
120	!>
121	!> @subsection subhgr namhgr
122	!> the grid sub-namelist parameters are :
123	!>
124	!> - in_mshhgr [@a 0]<br/>
125	!> type of horizontal mesh
126	!> - 0: curvilinear coordinate on the sphere read in coordinate.nc
127	!> - 1: geographical mesh on the sphere with regular grid-spacing
128	!> - 2: f-plane with regular grid-spacing
129	!> - 3: beta-plane with regular grid-spacing
130	!> - 4: Mercator grid with T/U point at the equator
131	!> - 5: beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
132	!> - dn_ppglam0 [@a ]<br/>
133	!> longitude of first raw and column T-point (in_mshhgr = 1 or 4)
134	!> - dn_ppgphi0 [@a ]<br/>
135	!> latitude of first raw and column T-point (in_mshhgr = 1 or 4)
136	!> - dn_ppe1_deg [@a ]<br/>
137	!> zonal grid-spacing (degrees) (in_mshhgr = 1,2,3 or 4)
138	!> - dn_ppe2_deg [@a ]<br/>
139	!> meridional grid-spacing (degrees) (in_mshhgr = 1,2,3 or 4)
140	!>
141	!>
142	!> @subsection subzgr namzgr
143	!> the vertical grid sub-namelist parameters are :
144	!>
145	!> - ln_zco [@a .FALSE.]<br/>
146	!> z-coordinate - full steps
147	!> - ln_zps [@a .FALSE.]<br/>
148	!> z-coordinate - partial steps
149	!> - ln_sco [@a .FALSE.]<br/>
150	!> s- or hybrid z-s-coordinate
151	!> - ln_isfcav [@a .FALSE.]<br/>
152	!> ice shelf cavities
153	!> - ln_iscpl [@a .FALSE.]<br/>
154	!> coupling with ice sheet
155	!> - ln_wd [@a .FALSE.]<br/>
156	!> Wetting/drying activation
157	!> - in_nlevel [@a 75]<br/>
158	!> number of vertical level
159	!>
160	!> @subsection subdmin namdmin
161	!> the minimum depth sub-namelist parameters are :
162	!>
163	!> - dn_hmin [@a ]<br/>
164	!> minimum ocean depth (>0) or minimum number of ocean levels (<0)
165	!> - dn_isfhmin [@a ]<br/>
166	!> threshold to discriminate grounded ice to floating ice
167	!>
168	!> @subsection subzco namzco
169	!> the vertical coordinate sub-namelist parameters are :
170	!>
171	!> - dn_pp_to_be_computed [@a 0]<br/>
172	!>
173	!> - dn_ppsur [@a -3958.951371276829]<br/>
174	!> coefficient to compute vertical grid
175	!>
176	!> - dn_ppa0 [@a 103.953009600000]<br/>
177	!> coefficient to compute vertical grid
178	!>
179	!> - dn_ppa1 [@a 2.415951269000]<br/>
180	!> coefficient to compute vertical grid
181	!>
182	!> - dn_ppkth [@a 15.351013700000]<br/>
183	!> coefficient to compute vertical grid
184	!>
185	!> - dn_ppacr [@a 7.000000000000]<br/>
186	!> coefficient to compute vertical grid
187	!>
188	!> - dn_ppdzmin [@a 6.]<br/>
189	!> minimum vertical spacing
190	!>
191	!> - dn_pphmax [@a 5750.]<br/>
192	!> maximum depth
193	!>
194	!> - @b ln_dbletanh [@a .TRUE.]<br/>
195	!> use double tanh to compute vartical grid
196	!>
197	!> - dn_ppa2 [@a 100.760928500000]<br/>
198	!> double tanh function parameter
199	!>
200	!> - dn_ppkth2 [@a 48.029893720000]<br/>
201	!> double tanh function parameter
202	!>
203	!> - dn_ppacr2 [@a 13.000000000000]<br/>
204	!> double tanh function parameter
205	!>
206	!> @note
207	!> If dn_ppa1, dn_ppa0 and dn_ppsur are undefined,
208	!> NEMO will compute them from dn_ppdzmin, dn_pphmax, dn_ppkth, dn_ppacr
209	!>
210	!> @warning
211	!> this namelist is also needed to define partial steps, sigma or hybrid coordinate.
212	!>
213	!> @subsection subzps namzps
214	!> the partial step sub-namelist parameters are :
215	!>
216	!> - dn_e3zps_min [@a 25.]<br/>
217	!> minimum thickness of partial step level (meters)
218	!> - dn_e3zps_rat [@a 0.2]<br/>
219	!> minimum thickness ratio of partial step level
220	!>
221	!>
222	!> @subsection subsco namsco
223	!> the sigma or hybrid sub-namelist parameters are :
224	!>
225	!> - ln_s_sh94 [@a .FALSE.]<br/>
226	!> use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1
227	!> - ln_s_sf12 [@a .FALSE.]<br/>
228	!> use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma
229	!> - dn_sbot_min [@a ]<br/>
230	!> minimum depth of s-bottom surface (>0) (m)
231	!> - dn_sbot_max [@a ]<br/>
232	!> maximum depth of s-bottom surface (= ocean depth) (>0) (m)
233	!> - dn_hc [@a ]<br/>
234	!> Critical depth for transition from sigma to stretched coordinates
235	!> <br/> <br/>
236	!> Song and Haidvogel 1994 stretching additional parameters
237	!> - dn_rmax [@a ]<br/>
238	!> maximum cut-off r-value allowed (0<dn_rmax<1)
239	!> - dn_theta [@a ]<br/>
240	!> surface control parameter (0<=dn_theta<=20)
241	!> - dn_thetb [@a ]<br/>
242	!> bottom control parameter (0<=dn_thetb<= 1)
243	!> - dn_bb [@a ]<br/>
244	!> stretching parameter ( dn_bb=0; top only, dn_bb =1; top and bottom)
245	!> <br/> <br/>
246	!> Siddorn and Furner stretching additional parameters
247	!> - ln_sigcrit [@a .FALSE.]<br/>
248	!> switching to sigma (T) or Z (F) at H<Hc
249	!> - dn_alpha [@a ]<br/>
250	!> stretchin parameter ( >1 surface; <1 bottom)
251	!> - dn_efold [@a ]<br/>
252	!> e-fold length scale for transition region
253	!> - dn_zs [@a ]<br/>
254	!> Surface cell depth (Zs) (m)
255	!> <br/>
256	!> Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
257	!> - dn_zb_a [@a ]<br/>
258	!> Bathymetry multiplier for Zb
259	!> - dn_zb_b [@a ]<br/>
260	!> Offset for Zb
261	!>
262	!> @subsection sublbc namlbc
263	!> the lateral boundary condition sub-namelist parameters are :
264	!>
265	!> - rn_shlat [@a 2.]<br/>
266	!> lateral boundary conditions at the coast (modify fmask)
267	!> - shlat = 0 : free slip
268	!> - 0 < shlat < 2 : partial slip
269	!> - shlat = 2 : no slip
270	!> - shlat > 2 : strong slip
271	!>
272	!> for more information see Boundary Condition at the Coast
273	!> in [NEMO documentation](https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/manual/pdf/NEMO_manual.pdf)
274	!>
275	!> @subsection subwd namwd
276	!> the wetting and drying sub-namelist parameters are :
277	!>
278	!> - dn_wdmin1 [@a ]<br/>
279	!> minimum water depth on dried cells
280	!> - dn_wdmin2 [@a ]<br/>
281	!> tolerrance of minimum water depth on dried cells
282	!> - dn_wdld [@a ]<br/>
283	!> land elevation below which wetting/drying
284	!>
285	!> @subsection subgrd namgrd
286	!> the grid sub-namelist parameters are :
287	!>
288	!> - in_cfg [@a 0]<br/>
289	!> inverse resolution of the configuration (1/4° => 4)
290	!> - ln_bench [@a .FALSE.]<br/>
291	!> GYRE (in_mshhgr = 5 ) used as Benchmark.<br/>
292	!> => forced the resolution to be about 100 km
293	!> - ln_c1d [@a .FALSE.]<br/>
294	!> use configuration 1D
295	!> - ln_e3_dep [@a .FALSE.]<br/>
296	!> vertical scale factors =T: e3.=dk[depth] =F: old definition
297	!>
298	!> @subsection subout namout
299	!> the output sub-namelist parameters are :
300	!>
301	!> - cn_domcfg [@a domain_cfg.nc]<br/>
302	!> output file name
303	!> - in_msh [@a 0]<br/>
304	!> number of output file and contain (0-9)
305	!> - in_nproc [@a 1]<br/>
306	!> number of processor to be used
307	!> - in_niproc [@a 1]<br/>
308	!> i-direction number of processor
309	!> - in_njproc [@a 1]<br/>
310	!> j-direction numebr of processor
311	!>
312	!> - if niproc, and njproc are provided : the program only look for land
313	!> processor to be removed
314	!> - if nproc is provided : the program compute each possible domain layout,
315	!> and save the one with the most land processor to be removed
316	!> - with no information about number of processors, the program
317	!> assume to use only one processor
318	!>
319	!> @note
320	!> - if in_msh = 0 : write 'domain_cfg.nc' file.
321	!> - if MOD(in_msh, 3) = 1 : write '<b>mesh_mask.nc</b>' file.
322	!> - if MOD(in_msh, 3) = 2 : write '<b>mesh.nc</b>' and '<b>mask.nc</b>' files.
323	!> - if MOD(in_msh, 3) = 0 : write '<b>mesh_hgr.nc</b>', '<b>mesh_zgr.nc</b>' and '<b>mask.nc</b>' files.<br/>
324	!> For huge size domain, use option 2 or 3 depending on your vertical coordinate.<br/>
325	!>
326	!> @note
327	!> - if 0 < in_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
328	!> - if 3 < in_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
329	!> corresponding to the depth of the bottom t- and w-points
330	!> - if 6 < in_msh <= 9: write 2D arrays corresponding to the depth and the
331	!> thickness (e3[tw]_ps) of the bottom points
332	!>
333	!> <hr>
334	!> @author J.Paul
335	!>
336	!> @date September, 2015 - Initial Version (based on domhgr.F90, domzgr.F90, domwri.F90)
337	!> @date October, 2016
338	!> - update from trunk (revision 6961): add wetting and drying, ice sheet coupling..
339	!> @date October, 2016
340	!> - dimension to be used select from configuration file
341	!> - do not use anymore special case for ORCA grid
342	!> - allow to write domain_cfg file
343	!> @date November, 2016
344	!> - choose vertical scale factors (e3.=dk[depth] or old definition)
345	!> @date January, 2019
346	!> - add url path to global attributes of output file(s)
347	!> @date February, 2019
348	!> - rename sub namelist namin to namsrc
349	!> @date Ocober, 2019
350	!> - add help and version optional arguments
351	!>
352	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
353	!----------------------------------------------------------------------
354	PROGRAM create_meshmask
355
356	USE global ! global variable
357	USE kind ! F90 kind parameter
358	USE phycst ! physical constant
359	USE logger ! log file manager
360	USE date ! date manager
361	USE fct ! basic useful function
362	USE att ! attribute manager
363	USE dim ! dimension manager
364	USE var ! variable manager
365	USE mpp ! MPP manager
366	USE iom_mpp ! I/O MPP manager
367	USE lbc ! lateral boundary conditions
368	USE grid
369	USE grid_hgr
370	USE grid_zgr
371
372	IMPLICIT NONE
373
374	! parameters
375	CHARACTER(LEN=lc), PARAMETER :: cp_myname = "create_meshmask"
376
377	! local variable
378	CHARACTER(LEN=lc) :: cl_arg
379	CHARACTER(LEN=lc) :: cl_namelist
380	CHARACTER(LEN=lc) :: cl_date
381	CHARACTER(LEN=lc) :: cl_url
382	CHARACTER(LEN=lc) :: cl_errormsg
383
384	INTEGER(i1), DIMENSION(:) , ALLOCATABLE :: bl_tmp
385
386	INTEGER(i4) :: il_narg
387	INTEGER(i4) :: il_status
388	INTEGER(i4) :: il_fileid
389	INTEGER(i4) :: il_attid
390	INTEGER(i4) :: il_ew
391	INTEGER(i4) :: jpi
392	INTEGER(i4) :: jpj
393	INTEGER(i4) :: jpk
394	INTEGER(i4), DIMENSION(:) , ALLOCATABLE :: il_tmp
395	INTEGER(i4), DIMENSION(:,:) , ALLOCATABLE :: il_mask
396
397	LOGICAL :: ll_exist
398	LOGICAL :: ll_domcfg
399
400	REAL(dp) , DIMENSION(:,:) , ALLOCATABLE :: dl_tmp2D
401	REAL(dp) , DIMENSION(:,:,:) , ALLOCATABLE :: dl_tmp3D
402
403	TYPE(TATT) :: tl_att
404	TYPE(TATT) , DIMENSION(:) , ALLOCATABLE :: tl_gatt
405
406	TYPE(TDIM) :: tl_dim
407
408	TYPE(TVAR) :: tl_bathy
409	TYPE(TVAR) :: tl_risfdep
410	TYPE(TVAR) :: tl_misf
411	TYPE(TVAR) :: tl_gdepu
412	TYPE(TVAR) :: tl_gdepv
413	TYPE(TVAR) :: tl_hdept
414	TYPE(TVAR) :: tl_hdepw
415	TYPE(TVAR) :: tl_scalar
416
417	TYPE(TNAMH) :: tl_namh
418	TYPE(TNAMZ) :: tl_namz
419
420	TYPE(TMPP) :: tl_mpp
421	TYPE(TMPP) , TARGET :: tl_mppout0
422	TYPE(TMPP) , TARGET :: tl_mppout1
423	TYPE(TMPP) , TARGET :: tl_mppout2
424	TYPE(TMPP) , POINTER :: tl_mppmsk
425	TYPE(TMPP) , POINTER :: tl_mpphgr
426	TYPE(TMPP) , POINTER :: tl_mppzgr
427
428	! parameter
429	INTEGER(i4) , PARAMETER :: ip_maxatt = 40
430
431	! loop indices
432	INTEGER(i4) :: ji
433	INTEGER(i4) :: jj
434	INTEGER(i4) :: jk
435
436	INTEGER(i4) :: ik
437
438	! namelist variable
439	! namlog
440	CHARACTER(LEN=lc) :: cn_logfile = 'create_meshmask.log'
441	CHARACTER(LEN=lc) :: cn_verbosity= 'warning'
442	INTEGER(i4) :: in_maxerror = 5
443
444	! namcfg
445	CHARACTER(LEN=lc) :: cn_varcfg = './cfg/variable.cfg'
446	CHARACTER(LEN=lc) :: cn_dimcfg = './cfg/dimension.cfg'
447	CHARACTER(LEN=lc) :: cn_dumcfg = './cfg/dummy.cfg'
448
449	! namsrc
450	CHARACTER(LEN=lc) :: cn_bathy = ''
451	CHARACTER(LEN=lc) :: cn_varbathy = ''
452	CHARACTER(LEN=lc) :: cn_coord = ''
453	CHARACTER(LEN=lc) :: cn_isfdep = ''
454	CHARACTER(LEN=lc) :: cn_varisfdep= 'isfdraft'
455	INTEGER(i4) :: in_perio = -1
456	LOGICAL :: ln_closea = .TRUE.
457
458	! namzgr
459	LOGICAL :: ln_zco = .FALSE.
460	LOGICAL :: ln_zps = .FALSE.
461	LOGICAL :: ln_sco = .FALSE.
462	LOGICAL :: ln_isfcav = .FALSE.
463	LOGICAL :: ln_iscpl = .FALSE.
464	LOGICAL :: ln_wd = .FALSE.
465	INTEGER(i4) :: in_nlevel = 75
466
467	! namlbc
468	REAL(sp) :: rn_shlat = 2.
469
470	! namout
471	CHARACTER(LEN=lc) :: cn_domcfg = 'domain_cfg.nc'
472	INTEGER(i4) :: in_msh = 0
473	CHARACTER(LEN=lc) :: cn_type = 'cdf'
474	INTEGER(i4) :: in_nproc = 0
475	INTEGER(i4) :: in_niproc = 0
476	INTEGER(i4) :: in_njproc = 0
477	!-------------------------------------------------------------------
478	NAMELIST /namlog/ & !< logger namelist
479	& cn_logfile, & !< log file
480	& cn_verbosity, & !< log verbosity
481	& in_maxerror !< logger maximum error
482
483	NAMELIST /namcfg/ & !< configuration namelist
484	& cn_varcfg, & !< variable configuration file
485	& cn_dimcfg, & !< dimension configuration file
486	& cn_dumcfg !< dummy configuration file
487
488	NAMELIST /namsrc/ & !< source namelist
489	& cn_bathy, & !< Bathymetry file
490	& cn_varbathy, & !< Bathymetry variable name
491	& cn_coord, & !< Coordinate file (in_mshhgr=0)
492	& cn_isfdep, & !< Iceshelf draft (ln_isfcav=true)
493	& cn_varisfdep, & !< Iceshelf draft variable name (ln_isfcav=true)
494	& in_perio, & !< NEMO periodicity
495	& ln_closea
496
497	NAMELIST /namzgr/ &
498	& ln_zco, & !< z-coordinate
499	& ln_zps, & !< z-coordinate with partial steps
500	& ln_sco, & !< s-coordinate
501	& ln_isfcav, & !< presence of ISF
502	& ln_iscpl, & !< coupling with ice sheet
503	& ln_wd, & !< Wetting/drying activation
504	& in_nlevel
505
506	NAMELIST /namlbc/ &
507	& rn_shlat !< lateral momentum boundary condition
508
509	NAMELIST /namout/ & !< output namlist
510	& cn_domcfg, & !< output file name
511	& in_msh, & !< number of output file (1,2,3)
512	& in_nproc, & !< number of processor to be used
513	& in_niproc, & !< i-direction number of processor
514	& in_njproc !< j-direction numebr of processor
515	!-------------------------------------------------------------------
516
517	!
518	! Initialisation
519	! --------------
520	!
521	il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec
522
523	! Traitement des arguments fournis
524	! --------------------------------
525	IF( il_narg /= 1 )THEN
526	WRITE(cl_errormsg,*) ' ERROR : one argument is needed '
527	CALL fct_help(cp_myname,cl_errormsg)
528	CALL EXIT(1)
529	ELSE
530
531	CALL GET_COMMAND_ARGUMENT(1,cl_arg) !f03 intrinsec
532	SELECT CASE (cl_arg)
533	CASE ('-v', '--version')
534
535	CALL fct_version(cp_myname)
536	CALL EXIT(0)
537
538	CASE ('-h', '--help')
539
540	CALL fct_help(cp_myname)
541	CALL EXIT(0)
542
543	CASE DEFAULT
544
545	cl_namelist=cl_arg
546
547	! read namelist
548	INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist)
549	IF( ll_exist )THEN
550
551	il_fileid=fct_getunit()
552
553	OPEN( il_fileid, FILE=TRIM(cl_namelist), &
554	& FORM='FORMATTED', &
555	& ACCESS='SEQUENTIAL', &
556	& STATUS='OLD', &
557	& ACTION='READ', &
558	& IOSTAT=il_status)
559	CALL fct_err(il_status)
560	IF( il_status /= 0 )THEN
561	WRITE(cl_errormsg,*) " ERROR : error opening "//TRIM(cl_namelist)
562	CALL fct_help(cp_myname,cl_errormsg)
563	CALL EXIT(1)
564	ENDIF
565
566	READ( il_fileid, NML = namlog )
567
568	! define logger file
569	CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror)
570	CALL logger_header()
571
572	READ( il_fileid, NML = namcfg )
573	! get variable extra information
574	CALL var_def_extra(TRIM(cn_varcfg))
575
576	! get dimension allowed
577	CALL dim_def_extra(TRIM(cn_dimcfg))
578
579	! get dummy variable
580	CALL var_get_dummy(TRIM(cn_dumcfg))
581	! get dummy dimension
582	CALL dim_get_dummy(TRIM(cn_dumcfg))
583	! get dummy attribute
584	CALL att_get_dummy(TRIM(cn_dumcfg))
585
586	READ( il_fileid, NML = namsrc )
587	READ( il_fileid, NML = namzgr )
588	READ( il_fileid, NML = namlbc )
589
590	READ( il_fileid, NML = namout )
591
592	CLOSE( il_fileid, IOSTAT=il_status )
593	CALL fct_err(il_status)
594	IF( il_status /= 0 )THEN
595	CALL logger_error("CREATE MASK: closing "//TRIM(cl_namelist))
596	ENDIF
597
598	ELSE
599
600	WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cl_namelist)
601	CALL fct_help(cp_myname,cl_errormsg)
602	CALL EXIT(1)
603
604	ENDIF
605
606	END SELECT
607	ENDIF
608
609	ll_domcfg=.FALSE.
610	IF( in_msh == 0 )THEN
611	ll_domcfg=.TRUE.
612	ENDIF
613
614	! open file
615	IF( cn_bathy /= '' )THEN
616	tl_mpp=mpp_init( file_init(TRIM(cn_bathy)), id_perio=in_perio)
617	CALL grid_get_info(tl_mpp)
618	ELSE
619	CALL logger_fatal("CREATE MESH MASK: no input bathymetry file found. "//&
620	& "check namelist")
621	ENDIF
622
623	! read bathymetry
624	WRITE(,) 'FILE TO BE USED:',TRIM(cn_bathy)
625	CALL iom_mpp_open(tl_mpp)
626
627	! get dimension
628	jpi=tl_mpp%t_dim(jp_I)%i_len
629	jpj=tl_mpp%t_dim(jp_J)%i_len
630	jpk=in_nlevel
631
632	WRITE(,) 'DIMENSION TO BE USED :',jpi,jpj,jpk
633
634	! read variable
635	IF( TRIM(cn_varbathy) == '' )THEN
636	IF( ln_zco )THEN
637	cn_varbathy='Bathy_level'
638	ELSEIF( ln_zps .OR. ln_sco )THEN
639	IF( ln_isfcav )THEN
640	cn_varbathy='Bathymetry_isf'
641	ELSE
642	cn_varbathy='Bathymetry'
643	ENDIF
644	ENDIF
645	ENDIF
646	WRITE(,) 'VARIABLE READ : '//TRIM(cn_varbathy)
647	tl_bathy=iom_mpp_read_var(tl_mpp, TRIM(cn_varbathy))
648	CALL iom_mpp_close(tl_mpp)
649
650	WHERE( tl_bathy%d_value(:,:,1,1) == tl_bathy%d_fill .OR. &
651	& tl_bathy%d_value(:,:,1,1) < 0._dp )
652	tl_bathy%d_value(:,:,1,1) = 0._dp
653	END WHERE
654
655	IF ( ln_isfcav ) THEN
656	WRITE(,) 'ICESHELF DRAFT FILE TO BE USED:',TRIM(cn_isfdep)
657	WRITE(,) 'ICESHELF VARIABLE READ : '//TRIM(cn_varisfdep)
658	! open Iceshelf draft
659	IF( cn_isfdep /= '' )THEN
660	tl_mpp=mpp_init( file_init(TRIM(cn_isfdep)), id_perio=in_perio)
661	CALL grid_get_info(tl_mpp)
662	ELSE
663	CALL logger_fatal("CREATE MESH MASK: no input Iceshelf draft '//&
664	& 'file found. check namelist")
665	ENDIF
666
667	! read Iceshelf draft
668	CALL iom_mpp_open(tl_mpp)
669	IF( ln_zps .OR. ln_sco ) THEN
670	tl_risfdep=iom_mpp_read_var(tl_mpp, cn_varisfdep)
671	ENDIF
672	CALL iom_mpp_close(tl_mpp)
673	ELSE
674	ALLOCATE(dl_tmp2D(jpi,jpj))
675	dl_tmp2D(:,:)=0._dp
676
677	tl_risfdep=var_init(cn_varisfdep, dl_tmp2D(:,:), id_type=NF90_DOUBLE)
678
679	DEALLOCATE(dl_tmp2D)
680	ENDIF
681
682	! fill closed sea
683	IF( ln_closea )THEN
684	WRITE(,) "CLOSED SEA"
685	ALLOCATE( il_mask(tl_bathy%t_dim(1)%i_len, &
686	& tl_bathy%t_dim(2)%i_len) )
687
688	! split domain in N sea subdomain
689	il_mask(:,:)=grid_split_domain(tl_bathy)
690
691	! fill smallest domain
692	CALL grid_fill_small_dom( tl_bathy, il_mask(:,:) )
693
694	DEALLOCATE( il_mask )
695	ENDIF
696
697	in_perio = tl_mpp%i_perio
698	il_ew = tl_mpp%i_ew
699
700	! clean
701	CALL mpp_clean(tl_mpp)
702
703	! Horizontal mesh (dom_hgr) -------------------------------------------------
704	tl_namh=grid_hgr_nam( cn_coord, in_perio, cl_namelist )
705
706	! init Horizontal grid global variable
707	CALL grid_hgr_init(jpi,jpj,jpk,ll_domcfg)
708
709	! compute horizontal mesh
710	WRITE(,) "COMPUTE HORIZONTAL MESH"
711	CALL grid_hgr_fill(tl_namh,jpi,jpj,ll_domcfg)
712
713	! Vertyical mesh (dom_zgr) -------------------------------------------------
714	tl_namz=grid_zgr_nam( cn_coord, in_perio, cl_namelist )
715
716	! init Vertical grid global variable
717	CALL grid_zgr_init(jpi,jpj,jpk,ln_sco)
718	IF( ln_zps ) CALL grid_zgr_zps_init(jpi,jpj)
719	IF( ln_sco ) CALL grid_zgr_sco_init(jpi,jpj)
720
721	! compute vertical mesh
722	WRITE(,) "COMPUTE VERTICAL MESH"
723	CALL grid_zgr_fill( tl_namz,jpi,jpj,jpk, tl_bathy, tl_risfdep )
724
725	! compute masks
726	WRITE(,) "COMPUTE MASK"
727	CALL create_meshmask__mask(tl_namh,jpi,jpj,jpk,ll_domcfg)
728
729	! Maximum stiffness ratio/hydrostatic consistency
730	IF( ln_sco ) CALL grid_zgr_sco_stiff(tl_namz,jpi,jpj,jpk)
731
732	! clean
733	CALL var_clean(tl_bathy)
734
735	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
736	! create ouptut structure
737	IF( in_niproc == 0 .AND. &
738	& in_njproc == 0 .AND. &
739	& in_nproc == 0 )THEN
740	in_niproc = 1
741	in_njproc = 1
742	in_nproc = 1
743	ENDIF
744
745	WRITE(,) "WRITE FILE(S)"
746	IF( ll_domcfg )THEN
747	! ! ============================
748	! ! create 'domain_cfg.nc' file
749	! ! ============================
750	tl_mppout0=mpp_init( cn_domcfg, tg_tmask, &
751	& in_niproc, in_njproc, in_nproc, &
752	& cd_type=cn_type )
753
754	tl_mppmsk=>tl_mppout0
755	tl_mpphgr=>tl_mppout0
756	tl_mppzgr=>tl_mppout0
757
758	ELSE
759	SELECT CASE ( MOD(in_msh, 3) )
760	! ! ============================
761	CASE ( 1 ) ! create 'mesh_mask.nc' file
762	! ! ============================
763	tl_mppout0=mpp_init( 'mesh_mask.nc', tg_tmask, &
764	& in_niproc, in_njproc, in_nproc, &
765	& cd_type=cn_type )
766
767	tl_mppmsk=>tl_mppout0
768	tl_mpphgr=>tl_mppout0
769	tl_mppzgr=>tl_mppout0
770
771	! ! ============================
772	CASE ( 2 ) ! create 'mesh.nc' and
773	! ! 'mask.nc' files
774	! ! ============================
775	tl_mppout0=mpp_init( 'mask.nc', tg_tmask, &
776	& in_niproc, in_njproc, in_nproc, &
777	& cd_type=cn_type )
778	tl_mppout1=mpp_init( 'mesh.nc', tg_tmask, &
779	& in_niproc, in_njproc, in_nproc, &
780	& cd_type=cn_type )
781
782	tl_mppmsk=>tl_mppout0
783	tl_mpphgr=>tl_mppout1
784	tl_mppzgr=>tl_mppout1
785
786	! ! ============================
787	CASE ( 0 ) ! create 'mesh_hgr.nc'
788	! ! 'mesh_zgr.nc' and
789	! ! 'mask.nc' files
790	! ! ============================
791	tl_mppout0=mpp_init( 'mask.nc', tg_tmask, &
792	& in_niproc, in_njproc, in_nproc, &
793	& cd_type=cn_type )
794	tl_mppout1=mpp_init( 'mesh_hgr.nc', tg_tmask, &
795	& in_niproc, in_njproc, in_nproc, &
796	& cd_type=cn_type )
797	tl_mppout2=mpp_init( 'mesh_zgr.nc', tg_tmask, &
798	& in_niproc, in_njproc, in_nproc, &
799	& cd_type=cn_type )
800	!
801
802	tl_mppmsk=>tl_mppout0
803	tl_mpphgr=>tl_mppout1
804	tl_mppzgr=>tl_mppout2
805
806	END SELECT
807	ENDIF
808
809	! add variables
810	IF( ll_domcfg )THEN
811	ALLOCATE(il_tmp(1))
812	tl_dim%l_use=.FALSE.
813
814	il_tmp(:)=jpi
815	tl_scalar=var_init('jpiglo', il_tmp(:), id_type=NF90_INT, td_dim=tl_dim)
816	CALL mpp_add_var(tl_mppmsk, tl_scalar)
817
818	il_tmp(:)=jpj
819	tl_scalar=var_init('jpjglo', il_tmp(:), id_type=NF90_INT, td_dim=tl_dim)
820	CALL mpp_add_var(tl_mppmsk, tl_scalar)
821
822	il_tmp(:)=jpk
823	tl_scalar=var_init('jpkglo', il_tmp(:), id_type=NF90_INT, td_dim=tl_dim)
824	CALL mpp_add_var(tl_mppmsk, tl_scalar)
825
826	il_tmp(:)=tl_mppout0%i_perio
827	tl_scalar=var_init('jperio', il_tmp(:), id_type=NF90_INT, td_dim=tl_dim)
828	CALL mpp_add_var(tl_mppmsk, tl_scalar)
829
830	DEALLOCATE(il_tmp)
831	ALLOCATE(bl_tmp(1))
832
833	bl_tmp(:)=0
834	IF( ln_zco ) bl_tmp(:)=1
835	tl_scalar=var_init('ln_zco',bl_tmp(:), id_type=NF90_BYTE, td_dim=tl_dim)
836	CALL mpp_add_var(tl_mppmsk, tl_scalar)
837
838	bl_tmp(:)=0
839	IF( ln_zps ) bl_tmp(:)=1
840	tl_scalar=var_init('ln_zps',bl_tmp(:), id_type=NF90_BYTE, td_dim=tl_dim)
841	CALL mpp_add_var(tl_mppmsk, tl_scalar)
842
843	bl_tmp(:)=0
844	IF( ln_sco ) bl_tmp(:)=1
845	tl_scalar=var_init('ln_sco',bl_tmp(:), id_type=NF90_BYTE, td_dim=tl_dim)
846	CALL mpp_add_var(tl_mppmsk, tl_scalar)
847
848	bl_tmp(:)=0
849	IF( ln_isfcav ) bl_tmp(:)=1
850	tl_scalar=var_init('ln_isfcav',bl_tmp(:), id_type=NF90_BYTE, td_dim=tl_dim)
851	CALL mpp_add_var(tl_mppmsk, tl_scalar)
852
853	DEALLOCATE(bl_tmp)
854	CALL var_clean(tl_scalar)
855	ENDIF
856
857	!!! mask (msk)
858	!!!----------------------
859	IF( .NOT. ll_domcfg )THEN
860	! tmask
861	CALL mpp_add_var(tl_mppmsk, tg_tmask)
862	CALL var_clean(tg_tmask)
863	! umask
864	CALL mpp_add_var(tl_mppmsk, tg_umask)
865	CALL var_clean(tg_umask)
866	! vmask
867	CALL mpp_add_var(tl_mppmsk, tg_vmask)
868	CALL var_clean(tg_vmask)
869	! fmask
870	CALL mpp_add_var(tl_mppmsk, tg_fmask)
871	CALL var_clean(tg_fmask)
872	ENDIF
873
874	!!! horizontal mesh (hgr)
875	!!!----------------------
876
877	! latitude
878	! glamt
879	CALL mpp_add_var(tl_mpphgr, tg_glamt)
880	CALL var_clean(tg_glamt)
881	! glamu
882	CALL mpp_add_var(tl_mpphgr, tg_glamu)
883	CALL var_clean(tg_glamu)
884	! glamV
885	CALL mpp_add_var(tl_mpphgr, tg_glamv)
886	CALL var_clean(tg_glamv)
887	! glamf
888	CALL mpp_add_var(tl_mpphgr, tg_glamf)
889	CALL var_clean(tg_glamf)
890
891	! longitude
892	! gphit
893	CALL mpp_add_var(tl_mpphgr, tg_gphit)
894	CALL var_clean(tg_gphit)
895	! gphiu
896	CALL mpp_add_var(tl_mpphgr, tg_gphiu)
897	CALL var_clean(tg_gphiu)
898	! gphiv
899	CALL mpp_add_var(tl_mpphgr, tg_gphiv)
900	CALL var_clean(tg_gphiv)
901	! gphif
902	CALL mpp_add_var(tl_mpphgr, tg_gphif)
903	CALL var_clean(tg_gphif)
904
905	! e1 scale factors
906	! e1t
907	CALL mpp_add_var(tl_mpphgr, tg_e1t)
908	CALL var_clean(tg_e1t)
909	! e1u
910	CALL mpp_add_var(tl_mpphgr, tg_e1u)
911	CALL var_clean(tg_e1u)
912	! e1v
913	CALL mpp_add_var(tl_mpphgr, tg_e1v)
914	CALL var_clean(tg_e1v)
915	! e1f
916	CALL mpp_add_var(tl_mpphgr, tg_e1f)
917	CALL var_clean(tg_e1f)
918
919	! e2 scale factors
920	! e2t
921	CALL mpp_add_var(tl_mpphgr, tg_e2t)
922	CALL var_clean(tg_e2t)
923	! e2u
924	CALL mpp_add_var(tl_mpphgr, tg_e2u)
925	CALL var_clean(tg_e2u)
926	! e2v
927	CALL mpp_add_var(tl_mpphgr, tg_e2v)
928	CALL var_clean(tg_e2v)
929	! e2f
930	CALL mpp_add_var(tl_mpphgr, tg_e2f)
931	CALL var_clean(tg_e2f)
932
933	! coriolis factor
934	! ff_t
935	CALL mpp_add_var(tl_mpphgr, tg_ff_t)
936	CALL var_clean(tg_ff_t)
937	! ff_f
938	CALL mpp_add_var(tl_mpphgr, tg_ff_f)
939	CALL var_clean(tg_ff_f)
940
941	! angles
942	IF( .NOT. ll_domcfg )THEN
943	! cost
944	CALL mpp_add_var(tl_mpphgr, tg_gcost)
945	CALL var_clean(tg_gcost)
946	! cosu
947	CALL mpp_add_var(tl_mpphgr, tg_gcosu)
948	CALL var_clean(tg_gcosu)
949	! cosv
950	CALL mpp_add_var(tl_mpphgr, tg_gcosv)
951	CALL var_clean(tg_gcosv)
952	! cosf
953	CALL mpp_add_var(tl_mpphgr, tg_gcosf)
954	CALL var_clean(tg_gcosf)
955
956	! sint
957	CALL mpp_add_var(tl_mpphgr, tg_gsint)
958	CALL var_clean(tg_gsint)
959	! sinu
960	CALL mpp_add_var(tl_mpphgr, tg_gsinu)
961	CALL var_clean(tg_gsinu)
962	! sinv
963	CALL mpp_add_var(tl_mpphgr, tg_gsinv)
964	CALL var_clean(tg_gsinv)
965	! sinf
966	CALL mpp_add_var(tl_mpphgr, tg_gsinf)
967	CALL var_clean(tg_gsinf)
968	ENDIF
969
970	!!! vertical mesh (zgr)
971	!!!----------------------
972	! note that mbkt is set to 1 over land ==> use surface tmask
973	!
974	! mbathy
975	tg_mbathy%d_value(:,:,:,:) = tg_ssmask%d_value(:,:,:,:) * &
976	& tg_mbkt%d_value(:,:,:,:)
977	!
978	IF( ll_domcfg )THEN
979	tg_mbathy%c_name='bottom_level'
980	ENDIF
981	CALL mpp_add_var(tl_mppzgr, tg_mbathy)
982	CALL var_clean(tg_mbathy)
983
984	! misf
985	ALLOCATE(dl_tmp2D(jpi,jpj))
986	dl_tmp2D(:,:)=dp_fill
987
988	tl_misf =var_init('misf ',dl_tmp2D(:,:), id_type=NF90_INT)
989
990	DEALLOCATE(dl_tmp2D)
991	tl_misf%d_value(:,:,1,1) = tg_ssmask%d_value(:,:,1,1) * &
992	& tg_mikt%d_value(:,:,1,1)
993	!
994	IF( ll_domcfg ) tl_misf%c_name='top_level'
995	CALL mpp_add_var(tl_mppzgr, tl_misf)
996	CALL var_clean(tl_misf)
997
998	IF( .NOT. ll_domcfg )THEN
999	! isfdraft
1000	tl_risfdep%d_value(:,:,:,:) = tl_risfdep%d_value(:,:,:,:) * &
1001	& tg_mikt%d_value(:,:,:,:)
1002
1003	CALL mpp_add_var(tl_mppzgr, tl_risfdep)
1004	CALL var_clean(tl_risfdep)
1005	ENDIF
1006
1007	IF( ln_sco ) THEN ! s-coordinate
1008
1009	IF( .NOT. ll_domcfg )THEN
1010	! hbatt
1011	CALL mpp_add_var(tl_mppzgr, tg_hbatt)
1012	CALL var_clean(tg_hbatt)
1013	! hbatu
1014	CALL mpp_add_var(tl_mppzgr, tg_hbatu)
1015	CALL var_clean(tg_hbatu)
1016	! hbatv
1017	CALL mpp_add_var(tl_mppzgr, tg_hbatv)
1018	CALL var_clean(tg_hbatv)
1019	! hbatf
1020	CALL mpp_add_var(tl_mppzgr, tg_hbatf)
1021	CALL var_clean(tg_hbatf)
1022
1023	! scaling coef.
1024	IF( .NOT. (tl_namz%l_s_sh94 .OR. tl_namz%l_s_sf12) )THEN
1025	! gsigt
1026	CALL mpp_add_var(tl_mppzgr, tg_gsigt)
1027	CALL var_clean(tg_gsigt)
1028	! gsigw
1029	CALL mpp_add_var(tl_mppzgr, tg_gsigw)
1030	CALL var_clean(tg_gsigw)
1031	! gsi3w
1032	CALL mpp_add_var(tl_mppzgr, tg_gsi3w)
1033	CALL var_clean(tg_gsi3w)
1034	! esigt
1035	CALL mpp_add_var(tl_mppzgr, tg_esigt)
1036	CALL var_clean(tg_esigt)
1037	! esigw
1038	CALL mpp_add_var(tl_mppzgr, tg_esigw)
1039	CALL var_clean(tg_esigw)
1040	ENDIF
1041	ENDIF
1042
1043	! scale factors
1044	! e3t_0
1045	CALL mpp_add_var(tl_mppzgr, tg_e3t_0)
1046	CALL var_clean(tg_e3t_0)
1047	! e3u_0
1048	CALL mpp_add_var(tl_mppzgr, tg_e3u_0)
1049	CALL var_clean(tg_e3u_0)
1050	! e3v_0
1051	CALL mpp_add_var(tl_mppzgr, tg_e3v_0)
1052	CALL var_clean(tg_e3v_0)
1053	! e3f_0
1054	CALL mpp_add_var(tl_mppzgr, tg_e3f_0)
1055	CALL var_clean(tg_e3f_0)
1056	! e3w_0
1057	CALL mpp_add_var(tl_mppzgr, tg_e3w_0)
1058	CALL var_clean(tg_e3w_0)
1059	! e3uw_0
1060	CALL mpp_add_var(tl_mppzgr, tg_e3uw_0)
1061	CALL var_clean(tg_e3uw_0)
1062	! e3vw_0
1063	CALL mpp_add_var(tl_mppzgr, tg_e3vw_0)
1064	CALL var_clean(tg_e3vw_0)
1065
1066	! Max. grid stiffness ratio
1067	! rx1
1068	IF( ll_domcfg ) tg_rx1%c_name='stiffness'
1069	CALL mpp_add_var(tl_mppzgr, tg_rx1)
1070	CALL var_clean(tg_rx1)
1071
1072	! stretched system
1073	IF( .NOT. tl_namz%l_e3_dep )THEN
1074	! gdept_1d
1075	CALL mpp_add_var(tl_mppzgr, tg_gdept_1d)
1076	CALL var_clean(tg_gdept_1d)
1077	! gdepw_1d
1078	CALL mpp_add_var(tl_mppzgr, tg_gdepw_1d)
1079	CALL var_clean(tg_gdepw_1d)
1080
1081	! gdept_0
1082	CALL mpp_add_var(tl_mppzgr, tg_gdept_0)
1083	CALL var_clean(tg_gdept_0)
1084	! gdepw_0
1085	CALL mpp_add_var(tl_mppzgr, tg_gdepw_0)
1086	CALL var_clean(tg_gdepw_0)
1087	ENDIF
1088
1089	ENDIF
1090
1091	IF( ln_zps ) THEN ! z-coordinate - partial steps
1092
1093	IF( ll_domcfg .OR. in_msh <= 6 ) THEN ! 3D vertical scale factors
1094
1095	! e3t_0
1096	CALL mpp_add_var(tl_mppzgr, tg_e3t_0)
1097	CALL var_clean(tg_e3t_0)
1098	! e3u_0
1099	CALL mpp_add_var(tl_mppzgr, tg_e3u_0)
1100	CALL var_clean(tg_e3u_0)
1101	! e3v_0
1102	CALL mpp_add_var(tl_mppzgr, tg_e3v_0)
1103	CALL var_clean(tg_e3v_0)
1104	! e3w_0
1105	CALL mpp_add_var(tl_mppzgr, tg_e3w_0)
1106	CALL var_clean(tg_e3w_0)
1107
1108	ELSE
1109
1110	DO jj = 1,jpj
1111	DO ji = 1,jpi
1112	ik=tg_mbkt%d_value(ji,jj,1,1)
1113	tg_e3tp%d_value(ji,jj,1,1) = tg_e3t_0%d_value(ji,jj,ik,1) * &
1114	& tg_ssmask%d_value(ji,jj,1,1)
1115	tg_e3wp%d_value(ji,jj,1,1) = tg_e3w_0%d_value(ji,jj,ik,1) * &
1116	& tg_ssmask%d_value(ji,jj,1,1)
1117	END DO
1118	END DO
1119	! e3t_ps
1120	CALL mpp_add_var(tl_mppzgr, tg_e3tp)
1121	CALL var_clean(tg_e3tp)
1122	! e3w_ps
1123	CALL mpp_add_var(tl_mppzgr, tg_e3wp)
1124	CALL var_clean(tg_e3wp)
1125
1126	ENDIF ! 3D vertical scale factors
1127
1128	IF( ll_domcfg .OR. in_msh <= 3 ) THEN ! 3D depth
1129
1130	IF( .NOT. tl_namz%l_e3_dep )THEN
1131
1132	! gdepu, gdepv
1133	IF( .NOT. ll_domcfg )THEN
1134	ALLOCATE(dl_tmp3D(jpi,jpj,jpk))
1135	dl_tmp3D(:,:,:)=dp_fill
1136
1137	tl_gdepu=var_init('gdepu',dl_tmp3D(:,:,:), id_type=NF90_FLOAT)
1138	tl_gdepv=var_init('gdepv',dl_tmp3D(:,:,:), id_type=NF90_FLOAT)
1139
1140	DEALLOCATE(dl_tmp3D)
1141	DO jk = 1,jpk
1142	DO jj = 1, jpj-1
1143	DO ji = 1, jpi-1 ! vector opt.
1144	tl_gdepu%d_value(ji,jj,jk,1) = MIN( tg_gdept_0%d_value(ji ,jj ,jk,1) , &
1145	& tg_gdept_0%d_value(ji+1,jj ,jk,1) )
1146
1147	tl_gdepv%d_value(ji,jj,jk,1) = MIN( tg_gdept_0%d_value(ji ,jj ,jk,1) , &
1148	& tg_gdept_0%d_value(ji ,jj+1,jk,1) )
1149	END DO
1150	END DO
1151	END DO
1152	CALL lbc_lnk( tl_gdepu%d_value(:,:,:,1), 'U', in_perio, 1._dp )
1153	CALL lbc_lnk( tl_gdepv%d_value(:,:,:,1), 'V', in_perio, 1._dp )
1154
1155	! gdepu
1156	CALL mpp_add_var(tl_mppzgr, tl_gdepu)
1157	CALL var_clean(tl_gdepu)
1158	! gdepv
1159	CALL mpp_add_var(tl_mppzgr, tl_gdepv)
1160	CALL var_clean(tl_gdepv)
1161	ENDIF
1162
1163	! gdept_0
1164	CALL mpp_add_var(tl_mppzgr, tg_gdept_0)
1165	CALL var_clean(tg_gdept_0)
1166
1167	! gdepw_0
1168	CALL mpp_add_var(tl_mppzgr, tg_gdepw_0)
1169	CALL var_clean(tg_gdepw_0)
1170	ENDIF
1171
1172	ELSE ! 2D bottom depth
1173	ALLOCATE(dl_tmp2D(jpi,jpj))
1174	dl_tmp2D(:,:)=dp_fill
1175
1176	tl_hdept=var_init('hdept',dl_tmp2D(:,:), id_type=NF90_INT)
1177	tl_hdepw=var_init('hdepw',dl_tmp2D(:,:), id_type=NF90_INT)
1178
1179	DEALLOCATE(dl_tmp2D)
1180	DO jj = 1,jpj
1181	DO ji = 1,jpi
1182	ik=tg_mbkt%d_value(ji,jj,1,1)
1183	tl_hdept%d_value(ji,jj,1,1) = tg_gdept_0%d_value(ji,jj,ik ,1) * tg_ssmask%d_value(ji,jj,1,1)
1184	tl_hdepw%d_value(ji,jj,1,1) = tg_gdepw_0%d_value(ji,jj,ik+1,1) * tg_ssmask%d_value(ji,jj,1,1)
1185	END DO
1186	END DO
1187	! hdept
1188	CALL mpp_add_var(tl_mppzgr, tl_hdept)
1189	CALL var_clean(tl_hdept)
1190	! hdepw
1191	CALL mpp_add_var(tl_mppzgr, tl_hdepw)
1192	CALL var_clean(tl_hdepw)
1193
1194	! clean
1195	CALL var_clean(tg_gdept_0)
1196
1197	ENDIF ! 3D depth
1198
1199	ENDIF
1200
1201	! scale factors
1202	IF( ll_domcfg )THEN
1203	! e3t_1d
1204	CALL mpp_add_var(tl_mppzgr, tg_e3t_1d)
1205	CALL var_clean(tg_e3t_1d)
1206	! e3w_1d
1207	CALL mpp_add_var(tl_mppzgr, tg_e3w_1d)
1208	CALL var_clean(tg_e3w_1d)
1209	ENDIF
1210
1211	IF( ln_zps .OR. ln_zco )THEN ! z-coordinate
1212	IF( .NOT. tl_namz%l_e3_dep )THEN
1213	! depth
1214	! gdept_1d
1215	CALL mpp_add_var(tl_mppzgr, tg_gdept_1d)
1216	CALL var_clean(tg_gdept_1d)
1217	! gdepw_1d
1218	CALL mpp_add_var(tl_mppzgr, tg_gdepw_1d)
1219	CALL var_clean(tg_gdepw_1d)
1220	ENDIF
1221	ENDIF
1222
1223	! define global attributes
1224	ALLOCATE(tl_gatt(ip_maxatt))
1225
1226	tl_gatt(:) = create_meshmask__gloatt(cn_bathy,cn_coord,cn_isfdep,tl_namh,tl_namz)
1227
1228
1229	IF( in_msh == 0 ) in_msh=1
1230	SELECT CASE ( MOD(in_msh, 3) )
1231	CASE ( 1 )
1232	! add some attribute
1233	tl_att=att_init("Created_by","SIREN create_meshmask")
1234	CALL mpp_add_att(tl_mppmsk, tl_att)
1235
1236	!add source url
1237	cl_url=fct_split(fct_split(cp_url,2,'$'),2,'URL:')
1238	tl_att=att_init("SIREN_url",cl_url)
1239	CALL mpp_add_att(tl_mppmsk, tl_att)
1240
1241	! add date of creation
1242	cl_date=date_print(date_now())
1243	tl_att=att_init("Creation_date",TRIM(cl_date))
1244	CALL mpp_add_att(tl_mppmsk, tl_att)
1245
1246	! add attribute periodicity
1247	il_attid=0
1248	IF( ASSOCIATED(tl_mppmsk%t_proc(1)%t_att) )THEN
1249	il_attid=att_get_id(tl_mppmsk%t_proc(1)%t_att(:),'periodicity')
1250	ENDIF
1251	IF( in_perio >= 0 .AND. il_attid == 0 )THEN
1252	tl_att=att_init('periodicity',in_perio)
1253	CALL mpp_add_att(tl_mppmsk,tl_att)
1254	ENDIF
1255
1256	il_attid=0
1257	IF( ASSOCIATED(tl_mppmsk%t_proc(1)%t_att) )THEN
1258	il_attid=att_get_id(tl_mppmsk%t_proc(1)%t_att(:),'ew_overlap')
1259	ENDIF
1260	IF( il_ew >= 0 .AND. il_attid == 0 )THEN
1261	tl_att=att_init('ew_overlap',il_ew)
1262	CALL mpp_add_att(tl_mppmsk,tl_att)
1263	ENDIF
1264
1265	ji=1
1266	DO WHILE( tl_gatt(ji)%c_name /= '' )
1267	CALL mpp_add_att(tl_mppmsk,tl_gatt(ji))
1268	ji=ji+1
1269	ENDDO
1270	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1271	! create file
1272	CALL iom_mpp_create(tl_mppmsk)
1273
1274	! write file
1275	CALL iom_mpp_write_file(tl_mppmsk)
1276	! close file
1277	CALL iom_mpp_close(tl_mppmsk)
1278
1279	! clean
1280	CALL mpp_clean(tl_mppmsk)
1281
1282	CASE ( 2 )
1283	! add some attribute
1284	tl_att=att_init("Created_by","SIREN create_meshmask")
1285	CALL mpp_add_att(tl_mppmsk, tl_att)
1286	CALL mpp_add_att(tl_mpphgr, tl_att)
1287
1288	!add source url
1289	cl_url=fct_split(fct_split(cp_url,2,'$'),2,'URL:')
1290	tl_att=att_init("SIREN_url",cl_url)
1291	CALL mpp_add_att(tl_mppmsk, tl_att)
1292	CALL mpp_add_att(tl_mpphgr, tl_att)
1293
1294	! add date of creation
1295	cl_date=date_print(date_now())
1296	tl_att=att_init("Creation_date",TRIM(cl_date))
1297	CALL mpp_add_att(tl_mppmsk, tl_att)
1298	CALL mpp_add_att(tl_mpphgr, tl_att)
1299
1300	! add attribute periodicity
1301	il_attid=0
1302	IF( ASSOCIATED(tl_mppmsk%t_proc(1)%t_att) )THEN
1303	il_attid=att_get_id(tl_mppmsk%t_proc(1)%t_att(:),'periodicity')
1304	ENDIF
1305	IF( in_perio >= 0 .AND. il_attid == 0 )THEN
1306	tl_att=att_init('periodicity',in_perio)
1307	CALL mpp_add_att(tl_mppmsk,tl_att)
1308	CALL mpp_add_att(tl_mpphgr,tl_att)
1309	ENDIF
1310
1311	il_attid=0
1312	IF( ASSOCIATED(tl_mppmsk%t_proc(1)%t_att) )THEN
1313	il_attid=att_get_id(tl_mppmsk%t_proc(1)%t_att(:),'ew_overlap')
1314	ENDIF
1315	IF( il_ew >= 0 .AND. il_attid == 0 )THEN
1316	tl_att=att_init('ew_overlap',il_ew)
1317	CALL mpp_add_att(tl_mppmsk,tl_att)
1318	CALL mpp_add_att(tl_mpphgr,tl_att)
1319	ENDIF
1320
1321	ji=1
1322	DO WHILE( tl_gatt(ji)%c_name /= '' )
1323	CALL mpp_add_att(tl_mppmsk,tl_gatt(ji))
1324	CALL mpp_add_att(tl_mpphgr,tl_gatt(ji))
1325	ji=ji+1
1326	ENDDO
1327
1328	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1329	! create mask file
1330	!-----------------
1331	CALL iom_mpp_create(tl_mppmsk)
1332
1333	! write file
1334	CALL iom_mpp_write_file(tl_mppmsk)
1335	! close file
1336	CALL iom_mpp_close(tl_mppmsk)
1337
1338	! clean
1339	CALL mpp_clean(tl_mppmsk)
1340
1341	! create mesh file
1342	!-----------------
1343	CALL iom_mpp_create(tl_mpphgr)
1344
1345	! write file
1346	CALL iom_mpp_write_file(tl_mpphgr)
1347	! close file
1348	CALL iom_mpp_close(tl_mpphgr)
1349
1350	! clean
1351	CALL mpp_clean(tl_mpphgr)
1352
1353	CASE(0)
1354	! add some attribute
1355	tl_att=att_init("Created_by","SIREN create_meshmask")
1356	CALL mpp_add_att(tl_mppmsk, tl_att)
1357	CALL mpp_add_att(tl_mpphgr, tl_att)
1358	CALL mpp_add_att(tl_mppzgr, tl_att)
1359
1360	!add source url
1361	cl_url=fct_split(fct_split(cp_url,2,'$'),2,'URL:')
1362	tl_att=att_init("SIREN_url",cl_url)
1363	CALL mpp_add_att(tl_mppmsk, tl_att)
1364	CALL mpp_add_att(tl_mpphgr, tl_att)
1365	CALL mpp_add_att(tl_mppzgr, tl_att)
1366
1367	! add date of creation
1368	cl_date=date_print(date_now())
1369	tl_att=att_init("Creation_date",TRIM(cl_date))
1370	CALL mpp_add_att(tl_mppmsk, tl_att)
1371	CALL mpp_add_att(tl_mpphgr, tl_att)
1372	CALL mpp_add_att(tl_mppzgr, tl_att)
1373
1374	! add attribute periodicity
1375	il_attid=0
1376	IF( ASSOCIATED(tl_mppmsk%t_proc(1)%t_att) )THEN
1377	il_attid=att_get_id(tl_mppmsk%t_proc(1)%t_att(:),'periodicity')
1378	ENDIF
1379	IF( in_perio >= 0 .AND. il_attid == 0 )THEN
1380	tl_att=att_init('periodicity',in_perio)
1381	CALL mpp_add_att(tl_mppmsk,tl_att)
1382	CALL mpp_add_att(tl_mpphgr,tl_att)
1383	CALL mpp_add_att(tl_mppzgr,tl_att)
1384	ENDIF
1385
1386	il_attid=0
1387	IF( ASSOCIATED(tl_mppmsk%t_proc(1)%t_att) )THEN
1388	il_attid=att_get_id(tl_mppmsk%t_proc(1)%t_att(:),'ew_overlap')
1389	ENDIF
1390	IF( il_ew >= 0 .AND. il_attid == 0 )THEN
1391	tl_att=att_init('ew_overlap',il_ew)
1392	CALL mpp_add_att(tl_mppmsk,tl_att)
1393	CALL mpp_add_att(tl_mpphgr,tl_att)
1394	CALL mpp_add_att(tl_mppzgr,tl_att)
1395	ENDIF
1396
1397	ji=1
1398	DO WHILE( tl_gatt(ji)%c_name /= '' )
1399	CALL mpp_add_att(tl_mppmsk,tl_gatt(ji))
1400	CALL mpp_add_att(tl_mpphgr,tl_gatt(ji))
1401	CALL mpp_add_att(tl_mppzgr,tl_gatt(ji))
1402	ji=ji+1
1403	ENDDO
1404
1405	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1406	! create mask file
1407	!-----------------
1408	CALL iom_mpp_create(tl_mppmsk)
1409
1410	! write file
1411	CALL iom_mpp_write_file(tl_mppmsk)
1412	! close file
1413	CALL iom_mpp_close(tl_mppmsk)
1414
1415	! clean
1416	WRITE(,) "CLEAN MSK"
1417	CALL mpp_clean(tl_mppmsk)
1418
1419	! create mesh_hgr file
1420	!-----------------
1421	CALL iom_mpp_create(tl_mpphgr)
1422
1423	! write file
1424	CALL iom_mpp_write_file(tl_mpphgr)
1425	! close file
1426	CALL iom_mpp_close(tl_mpphgr)
1427
1428	! clean
1429	WRITE(,) "CLEAN HGR"
1430	CALL mpp_clean(tl_mpphgr)
1431
1432	! create mesh_zgr file
1433	!-----------------
1434	WRITE(,) "CREATE ZGR"
1435	CALL iom_mpp_create(tl_mppzgr)
1436
1437	! write file
1438	WRITE(,) "WRITE ZGR"
1439	CALL iom_mpp_write_file(tl_mppzgr)
1440	! close file
1441	WRITE(,) "CLOSE ZGR"
1442	CALL iom_mpp_close(tl_mppzgr)
1443
1444	! clean
1445	WRITE(,) "CLEAN ZGR"
1446	CALL mpp_clean(tl_mppzgr)
1447
1448	CASE DEFAULT
1449	CALL logger_fatal("CREATE MESHMASK : something wrong with in_msh("//&
1450	& TRIM(fct_str(in_msh))//"), check your namelist.")
1451	END SELECT
1452
1453	! clean
1454	WRITE(,) "CLEAN"
1455	CALL att_clean(tl_att)
1456	CALL att_clean(tl_gatt)
1457
1458	DEALLOCATE(tl_gatt)
1459
1460	CALL grid_hgr_clean(ll_domcfg)
1461	CALL grid_zgr_clean(ln_sco)
1462	IF( ln_zps ) CALL grid_zgr_zps_clean()
1463	IF( ln_sco ) CALL grid_zgr_sco_clean()
1464	CALL var_clean_extra()
1465
1466	! close log file
1467	CALL logger_footer()
1468	CALL logger_close()
1469	CONTAINS
1470	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1471	SUBROUTINE create_meshmask__mask(td_nam,jpi,jpj,jpk,ld_domcfg)
1472	!-------------------------------------------------------------------
1473	!> @brief This subroutine compute land/ocean mask arrays at tracer points,
1474	!> horizontal velocity points (u & v), vorticity points (f) and
1475	!> barotropic stream function points (b).
1476	!>
1477	!> @details
1478	!>
1479	!> ** Method : The ocean/land mask is computed from the basin bathymetry in level (mbathy)
1480	!> which is defined or read in dommba.
1481	!> mbathy equals 0 over continental T-point and the number of ocean level over the ocean.
1482	!>
1483	!> At a given position (ji,jj,jk) the ocean/land mask is given by:
1484	!> - t-point :
1485	!> - 0. IF mbathy( ji ,jj) =< 0
1486	!> - 1. IF mbathy( ji ,jj) >= jk
1487	!> - u-point :
1488	!> - 0. IF mbathy( ji ,jj) or mbathy(ji+1, jj ) =< 0
1489	!> - 1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
1490	!> - v-point :
1491	!> - 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) =< 0
1492	!> - 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
1493	!> - f-point :
1494	!> - 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) or mbathy(ji+1,jj) or mbathy(ji+1,jj+1) =< 0
1495	!> - 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
1496	!> - b-point : the same definition as for f-point of the first ocean
1497	!> level (surface level) but with 0 along coastlines.
1498	!> - tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
1499	!> rows/lines due to cyclic or North Fold boundaries as well
1500	!> as MPP halos.
1501	!>
1502	!> @warning do not set the lateral friction through the value of fmask along
1503	!> the coast and topography.
1504	!>
1505	!> @note If nperio not equal to 0, the land/ocean mask arrays
1506	!> are defined with the proper value at lateral domain boundaries,
1507	!> but bmask. indeed, bmask defined the domain over which the
1508	!> barotropic stream function is computed. this domain cannot
1509	!> contain identical columns because the matrix associated with
1510	!> the barotropic stream function equation is then no more inverti-
1511	!> ble. therefore bmask is set to 0 along lateral domain boundaries
1512	!> even IF nperio is not zero.
1513	!>
1514	!> In case of open boundaries (lk_bdy=T):
1515	!> - tmask is set to 1 on the points to be computed bay the open
1516	!> boundaries routines.
1517	!> - bmask is set to 0 on the open boundaries.
1518	!>
1519	!> ** Action :
1520	!> - tmask : land/ocean mask at t-point (=0. or 1.)
1521	!> - umask : land/ocean mask at u-point (=0. or 1.)
1522	!> - vmask : land/ocean mask at v-point (=0. or 1.)
1523	!> - fmask : land/ocean mask at f-point (=0. or 1.)
1524	!> - bmask : land/ocean mask at barotropic stream
1525	!> function point (=0. or 1.) and set to 0 along lateral boundaries
1526	!> - tmask_i : interior ocean mask
1527	!>
1528	!> @author J.Paul
1529	!> @date September, 2015 - rewrite from dom_msk
1530	!> @date October, 2016
1531	!> - do not use anymore special case for ORCA grid
1532	!>
1533	!> @param[in] td_nam
1534	!> @param[in] jpi
1535	!> @param[in] jpj
1536	!> @param[in] jpk
1537	!-------------------------------------------------------------------
1538
1539	IMPLICIT NONE
1540
1541	! Argument
1542	TYPE(TNAMH), INTENT(IN) :: td_nam
1543	INTEGER(i4), INTENT(IN) :: jpi
1544	INTEGER(i4), INTENT(IN) :: jpj
1545	INTEGER(i4), INTENT(IN) :: jpk
1546	LOGICAL , INTENT(IN) :: ld_domcfg
1547
1548	! local variable
1549	! INTEGER(i4) :: ii0, ii1 ! local integers
1550	! INTEGER(i4) :: ij0, ij1
1551	! INTEGER(i4) :: isrow
1552
1553	! INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: imsk
1554	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: zwf
1555
1556	! REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_tpol
1557	! REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_fpol
1558
1559	! loop indices
1560	INTEGER(i4) :: ji
1561	INTEGER(i4) :: jj
1562	INTEGER(i4) :: jk
1563	!----------------------------------------------------------------
1564
1565	! ALLOCATE( imsk(jpi,jpj) )
1566
1567	! ALLOCATE( dl_tpol(jpi) )
1568	! ALLOCATE( dl_fpol(jpi) )
1569
1570	CALL logger_info('dommsk : ocean mask ')
1571	CALL logger_info('~~~~~~')
1572	IF ( rn_shlat == 0. )THEN ; CALL logger_info(' ocean lateral free-slip ')
1573	ELSEIF( rn_shlat == 2. )THEN ; CALL logger_info(' ocean lateral no-slip ')
1574	ELSEIF( 0. < rn_shlat .AND. rn_shlat < 2. )THEN ; CALL logger_info(' ocean lateral partial-slip ')
1575	ELSEIF( 2. < rn_shlat )THEN ; CALL logger_info(' ocean lateral strong-slip ')
1576	ELSE ; CALL logger_info(' rn_shlat is negative = '//TRIM(fct_str(rn_shlat)))
1577	ENDIF
1578
1579	! 1. Ocean/land mask at t-point (computed from mbathy)
1580	! -----------------------------
1581	! N.B. tmask has already the right boundary conditions since mbathy is ok
1582	!
1583	tg_tmask%d_value(:,:,:,1) = 0._dp
1584	DO jk = 1, jpk
1585	DO jj = 1, jpj
1586	DO ji = 1, jpi
1587	IF( tg_mbathy%d_value(ji,jj,1,1) - REAL(jk,dp) + 0.1_dp >= 0._dp )THEN
1588	tg_tmask%d_value(ji,jj,jk,1) = 1._dp
1589	ENDIF
1590	ENDDO
1591	ENDDO
1592	ENDDO
1593
1594	! (ISF) define barotropic mask and mask the ice shelf point
1595	tg_ssmask%d_value(:,:,1,1)=tg_tmask%d_value(:,:,1,1) ! at this stage ice shelf is not masked
1596
1597	DO jk = 1, jpk
1598	DO jj = 1, jpj
1599	DO ji = 1, jpi
1600	IF( tg_misfdep%d_value(ji,jj,1,1) - REAL(jk,dp) - 0.1_dp >= 0._dp ) THEN
1601	tg_tmask%d_value(ji,jj,jk,1) = 0._dp
1602	END IF
1603	ENDDO
1604	ENDDO
1605	ENDDO
1606
1607	! ! Interior domain mask (used for global sum)
1608	! ! --------------------
1609	! tg_tmask_i%d_value(:,:,1,1) = tg_ssmask%d_value(:,:,1,1) ! (ISH) tmask_i = 1 even on the ice shelf
1610	!
1611	! tg_tmask_h%d_value(:,:,1,1) = 1._dp ! 0 on the halo and 1 elsewhere
1612	!
1613	! tg_tmask_h%d_value( 1 , : ,1,1) = 0._dp ! first columns
1614	! tg_tmask_h%d_value(jpi, : ,1,1) = 0._dp ! last columns
1615	! tg_tmask_h%d_value( : , 1 ,1,1) = 0._dp ! first rows
1616	! tg_tmask_h%d_value( : ,jpj,1,1) = 0._dp ! last rows
1617	!
1618	! ! north fold mask
1619	! ! ---------------
1620	! dl_tpol(1:jpi) = 1._dp
1621	! dl_fpol(1:jpi) = 1._dp
1622	! IF( td_nam%i_perio == 3 .OR. td_nam%i_perio == 4 )THEN ! T-point pivot
1623	! dl_tpol(jpi/2+1:jpi) = 0._dp
1624	! dl_fpol( 1 :jpi) = 0._dp
1625	! IF( mjg(nlej) == jpj ) THEN ! only half of the nlcj-1 row
1626	! DO ji = iif+1, iil-1
1627	! tg_tmask_h%d_value(ji,nlej-1,1,1) = tg_tmask_h%d_value(ji,nlej-1,1,1) * dl_tpol(mig(ji))
1628	! END DO
1629	! ENDIF
1630	! ENDIF
1631	!
1632	! tg_tmask_i%d_value(:,:,1,1) = tg_tmask_i%d_value(:,:,1,1) * tg_tmask_h%d_value(:,:,1,1)
1633	!
1634	! IF( td_nam%i_perio == 5 .OR. td_nam%i_perio == 6 )THEN ! F-point pivot
1635	! dl_tpol( 1 :jpi) = 0._dp
1636	! dl_fpol(jpi/2+1:jpi) = 0._dp
1637	! ENDIF
1638
1639	! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask)
1640	! -------------------------------------------
1641	DO jk = 1, jpk
1642	DO jj = 1, jpj-1
1643	DO ji = 1, jpi-1 ! vector loop
1644	tg_umask%d_value(ji,jj,jk,1) = tg_tmask%d_value(ji ,jj ,jk,1) * &
1645	& tg_tmask%d_value(ji+1,jj ,jk,1)
1646
1647	tg_vmask%d_value(ji,jj,jk,1) = tg_tmask%d_value(ji ,jj ,jk,1) * &
1648	& tg_tmask%d_value(ji ,jj+1,jk,1)
1649	!END DO
1650	!DO ji = 1, jpi-1 ! NO vector opt.
1651	IF( .NOT. ld_domcfg )THEN
1652	tg_fmask%d_value(ji,jj,jk,1) = tg_tmask%d_value(ji ,jj ,jk,1) * &
1653	& tg_tmask%d_value(ji+1,jj ,jk,1) * &
1654	& tg_tmask%d_value(ji ,jj+1,jk,1) * &
1655	& tg_tmask%d_value(ji+1,jj+1,jk,1)
1656	ENDIF
1657	ENDDO
1658	ENDDO
1659	ENDDO
1660
1661	! ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
1662	! DO jj = 1, jpjm1
1663	! DO ji = 1, fs_jpim1 ! vector loop
1664	! tg_ssumask%d_value(ji,jj,1,1) = tg_ssmask%d_value(ji ,jj ,1,1) * &
1665	! & tg_ssmask%d_value(ji+1,jj ,1,1) * &
1666	! & MIN( 1._wp, SUM(tg_umask%d_value(ji,jj,:,1)) )
1667	! tg_ssvmask%d_value(ji,jj,1,1) = tg_ssmask%d_value(ji ,jj ,1,1) * &
1668	! & tg_ssmask%d_value(ji ,jj+1,1,1) * &
1669	! & MIN( 1._wp, SUM(tg_vmask%d_value(ji,jj,:,1)) )
1670	! END DO
1671	! DO ji = 1, jpim1 ! NO vector opt.
1672	! tg_ssfmask%d_value(ji,jj,1,1) = tg_ssmask%d_value(ji ,jj ,1,1) * &
1673	! & tg_ssmask%d_value(ji+1,jj ,1,1) * &
1674	! & tg_ssmask%d_value(ji ,jj+1,1,1) * &
1675	! & tg_ssmask%d_value(ji+1,jj+1,1,1) * &
1676	! & MIN( 1._wp, SUM(tg_fmask%d_value(ji,jj,:,1)) )
1677	! END DO
1678	! END DO
1679	CALL lbc_lnk( tg_umask%d_value (:,:,:,1), 'U', td_nam%i_perio, 1._dp ) ! Lateral boundary conditions
1680	CALL lbc_lnk( tg_vmask%d_value (:,:,:,1), 'V', td_nam%i_perio, 1._dp )
1681	IF( .NOT. ld_domcfg )THEN
1682	CALL lbc_lnk( tg_fmask%d_value (:,:,:,1), 'F', td_nam%i_perio, 1._dp )
1683	ENDIF
1684	! CALL lbc_lnk( tg_ssumask%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp ) ! Lateral boundary conditions
1685	! CALL lbc_lnk( tg_ssvmask%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
1686	! CALL lbc_lnk( tg_ssfmask%d_value(:,:,:,1), 'F', td_nam%i_perio, 1._dp )
1687
1688	! 3. Ocean/land mask at wu-, wv- and w points
1689	!----------------------------------------------
1690	! tg_wmask%d_value (:,:,1,1) = tg_tmask%d_value(:,:,1,1) ! surface
1691	! tg_wumask%d_value(:,:,1,1) = tg_umask%d_value(:,:,1,1)
1692	! tg_wvmask%d_value(:,:,1,1) = tg_vmask%d_value(:,:,1,1)
1693	!
1694	! DO jk=2,jpk ! interior values
1695	! tg_wmask%d_value (:,:,jk,1) = tg_tmask%d_value(:,:,jk ,1) * &
1696	! & tg_tmask%d_value(:,:,jk-1,1)
1697	! tg_wumask%d_value(:,:,jk,1) = tg_umask%d_value(:,:,jk ,1) * &
1698	! & tg_umask%d_value(:,:,jk-1,1)
1699	! tg_wvmask%d_value(:,:,jk,1) = tg_vmask%d_value(:,:,jk ,1) * &
1700	! & tg_vmask%d_value(:,:,jk-1,1)
1701	! ENDDO
1702
1703	! Lateral boundary conditions on velocity (modify fmask)
1704	! ---------------------------------------
1705	IF( .NOT. ld_domcfg )THEN
1706	ALLOCATE( zwf(jpi,jpj) )
1707	DO jk = 1, jpk
1708	zwf(:,:) = tg_fmask%d_value(:,:,jk,1)
1709	DO jj = 2, jpj-1
1710	DO ji = 2, jpi-1 ! vector opt.
1711	IF( tg_fmask%d_value(ji,jj,jk,1) == 0._dp )THEN
1712	tg_fmask%d_value(ji,jj,jk,1) = rn_shlat * &
1713	& MIN(1._dp , MAX(zwf(ji+1,jj), zwf(ji,jj+1), &
1714	& zwf(ji-1,jj), zwf(ji,jj-1)) )
1715	ENDIF
1716	END DO
1717	END DO
1718	DO jj = 2, jpj-1
1719	IF( tg_fmask%d_value(1,jj,jk,1) == 0._dp )THEN
1720	tg_fmask%d_value(1 ,jj,jk,1) = rn_shlat * &
1721	& MIN(1._dp, MAX(zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1)))
1722	ENDIF
1723	IF( tg_fmask%d_value(jpi,jj,jk,1) == 0._dp )THEN
1724	tg_fmask%d_value(jpi,jj,jk,1) = rn_shlat * &
1725	& MIN(1._wp, MAX(zwf(jpi,jj+1), zwf(jpi-1,jj), zwf(jpi,jj-1)))
1726	ENDIF
1727	END DO
1728	DO ji = 2, jpi-1
1729	IF( tg_fmask%d_value(ji,1,jk,1) == 0._dp )THEN
1730	tg_fmask%d_value(ji, 1 ,jk,1) = rn_shlat * &
1731	& MIN(1._dp, MAX(zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1)))
1732	ENDIF
1733	IF( tg_fmask%d_value(ji,jpj,jk,1) == 0._dp )THEN
1734	tg_fmask%d_value(ji,jpj,jk,1) = rn_shlat * &
1735	& MIN(1._dp, MAX(zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpj-1)))
1736	ENDIF
1737	END DO
1738	END DO
1739	DEALLOCATE( zwf )
1740
1741	! IF( td_nam%c_cfg == "orca" .AND. td_nam%i_cfg == 2 )THEN ! ORCA_R2 configuration
1742	! ! ! Increased lateral friction near of some straits
1743	! IF( td_nam%i_cla == 0 ) THEN
1744	! ! Gibraltar strait : partial slip (fmask=0.5)
1745	! ij0 = 101 ; ij1 = 101
1746	! ii0 = 139 ; ii1 = 140
1747	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1) = 0.5_dp
1748	!
1749	! ij0 = 102 ; ij1 = 102
1750	! ii0 = 139 ; ii1 = 140
1751	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1) = 0.5_dp
1752	! !
1753	! !Bab el Mandeb : partial slip (fmask=1)
1754	! ij0 = 87 ; ij1 = 88
1755	! ii0 = 160 ; ii1 = 160
1756	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1) = 1._dp
1757	!
1758	! ij0 = 88 ; ij1 = 88
1759	! ii0 = 159 ; ii1 = 159
1760	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1) = 1._dp
1761	! !
1762	! ENDIF
1763	! ENDIF
1764	! !
1765	! IF( td_nam%c_cfg == "orca" .AND. td_nam%i_cfg == 1 )THEN ! ORCA R1 configuration
1766	! ! Increased lateral friction near of some straits
1767	! ! This dirty section will be suppressed by simplification process:
1768	! ! all this will come back in input files
1769	! ! Currently these hard-wired indices relate to configuration with
1770	! ! extend grid (jpjglo=332)
1771	! !
1772	! isrow = 332 - jpj
1773	! ! Gibraltar Strait
1774	! ii0 = 282 ; ii1 = 283
1775	! ij0 = 201 + isrow ; ij1 = 241 - isrow
1776	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1) = 2._dp
1777	!
1778	! ! Bhosporus Strait
1779	! ii0 = 314 ; ii1 = 315
1780	! ij0 = 208 + isrow ; ij1 = 248 - isrow
1781	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 2._dp
1782	!
1783	! ! Makassar Strait (Top)
1784	! ii0 = 48 ; ii1 = 48
1785	! ij0 = 149 + isrow ; ij1 = 190 - isrow
1786	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 3._dp
1787	!
1788	! ! Lombok Strait
1789	! ii0 = 44 ; ii1 = 44
1790	! ij0 = 124 + isrow ; ij1 = 165 - isrow
1791	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 2._dp
1792	!
1793	! ! Ombai Strait
1794	! ii0 = 53 ; ii1 = 53
1795	! ij0 = 124 + isrow ; ij1 = 165 - isrow
1796	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 2._dp
1797	!
1798	! ! Timor Passage
1799	! ii0 = 56 ; ii1 = 56
1800	! ij0 = 124 + isrow ; ij1 = 165 - isrow
1801	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 2._dp
1802	!
1803	! ! West Halmahera Strait
1804	! ii0 = 58 ; ii1 = 58
1805	! ij0 = 141 + isrow ; ij1 = 182 - isrow
1806	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 3._dp
1807	!
1808	! ! East Halmahera Strait
1809	! ii0 = 55 ; ii1 = 55
1810	! ij0 = 141 + isrow ; ij1 = 182 - isrow
1811	! tg_fmask%d_value(ii0:ii1,ij0:ij1,1:jpk,1 ) = 3._dp
1812	! !
1813	! ENDIF
1814	!
1815	CALL lbc_lnk( tg_fmask%d_value(:,:,:,1), 'F', td_nam%i_perio, 1._dp ) ! Lateral boundary conditions on fmask
1816	ENDIF
1817
1818	! DEALLOCATE( imsk )
1819
1820	! DEALLOCATE( dl_tpol )
1821	! DEALLOCATE( dl_fpol )
1822
1823	END SUBROUTINE create_meshmask__mask
1824	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1825	FUNCTION create_meshmask__gloatt(cd_bathy,cd_coord,cd_isfdep,td_namh,td_namz) &
1826	& RESULT(tf_att)
1827	!-------------------------------------------------------------------
1828	!> @brief
1829	!> this function create array of global attributes.
1830	!>
1831	!> @author J.Paul
1832	!> @date October, 2016 - initial release
1833	!>
1834	!> @param[in] cd_bathy
1835	!> @param[in] cd_coord
1836	!> @param[in] cd_isfdep
1837	!> @param[in] td_namh
1838	!> @param[in] td_namz
1839	!-------------------------------------------------------------------
1840
1841	IMPLICIT NONE
1842
1843	! Argument
1844	CHARACTER(LEN=*), INTENT(IN ) :: cd_bathy
1845	CHARACTER(LEN=*), INTENT(IN ) :: cd_coord
1846	CHARACTER(LEN=*), INTENT(IN ) :: cd_isfdep
1847	TYPE(TNAMH) , INTENT(IN ) :: td_namh
1848	TYPE(TNAMZ) , INTENT(IN ) :: td_namz
1849
1850	! function
1851	TYPE(TATT), DIMENSION(ip_maxatt) :: tf_att
1852
1853	! loop indices
1854	INTEGER(i4) :: ji
1855	!----------------------------------------------------------------
1856
1857	ji=1 ; tf_att(ji)=att_init("src_bathy",TRIM(cd_bathy))
1858	! horizontal grid
1859	ji=ji+1 ; tf_att(ji)=att_init("in_mshhgr",td_namh%i_mshhgr)
1860	SELECT CASE(td_namh%i_mshhgr)
1861	CASE(0)
1862	ji=ji+1 ; tf_att(ji)=att_init("src_coord",TRIM(cd_coord))
1863	CASE(1,4)
1864	ji=ji+1 ; tf_att(ji)=att_init("ppglam0",td_namh%d_ppglam0)
1865	ji=ji+1 ; tf_att(ji)=att_init("ppgphi0",td_namh%d_ppgphi0)
1866	CASE(2,3)
1867	ji=ji+1 ; tf_att(ji)=att_init("ppglam0",td_namh%d_ppglam0)
1868	ji=ji+1 ; tf_att(ji)=att_init("ppgphi0",td_namh%d_ppgphi0)
1869	ji=ji+1 ; tf_att(ji)=att_init("ppe1_deg",td_namh%d_ppe1_deg)
1870	ji=ji+1 ; tf_att(ji)=att_init("ppe2_deg",td_namh%d_ppe2_deg)
1871	END SELECT
1872
1873	IF( td_namz%l_isfcav )THEN
1874	ji=ji+1 ; tf_att(ji)=att_init("ice_shelf_cavities","activated")
1875	ji=ji+1 ; tf_att(ji)=att_init("src_isfdep",TRIM(cd_isfdep))
1876	ENDIF
1877	IF( td_namz%l_iscpl )THEN
1878	ji=ji+1 ; tf_att(ji)=att_init("ice_sheet_coupling","activated")
1879	ENDIF
1880
1881	! vertical grid
1882	IF( td_namz%l_zco )THEN
1883	ji=ji+1 ; tf_att(ji)=att_init("z_coord","full steps")
1884	ENDIF
1885	IF( td_namz%l_zps )THEN
1886	ji=ji+1 ; tf_att(ji)=att_init("z_coord","partial steps")
1887	ENDIF
1888	IF( td_namz%l_sco )THEN
1889	IF( td_namz%l_s_sh94 )THEN
1890	ji=ji+1 ; tf_att(ji)=att_init("z_coord","hybrid Song and Haidvogel 1994")
1891	ELSEIF( td_namz%l_s_sf12 )THEN
1892	ji=ji+1 ; tf_att(ji)=att_init("z_coord","hybrid Siddorn and Furner 2012")
1893	ELSE
1894	ji=ji+1 ; tf_att(ji)=att_init("z_coord","sigma")
1895	ENDIF
1896	ENDIF
1897	ji=ji+1 ; tf_att(ji)=att_init("hmin",td_namz%d_hmin)
1898	IF( td_namz%l_isfcav ) ji=ji+1 ; tf_att(ji)=att_init("isfhmin",td_namz%d_isfhmin)
1899
1900	! zco
1901	IF( td_namz%d_ppsur /= NF90_FILL_DOUBLE )THEN
1902	ji=ji+1 ; tf_att(ji)=att_init("ppsur",td_namz%d_ppsur)
1903	ELSE
1904	ji=ji+1 ; tf_att(ji)=att_init("ppsur","to_be_computed")
1905	ENDIF
1906	IF( td_namz%d_ppa0 /= NF90_FILL_DOUBLE )THEN
1907	ji=ji+1 ; tf_att(ji)=att_init("ppa0",td_namz%d_ppa0)
1908	ELSE
1909	ji=ji+1 ; tf_att(ji)=att_init("ppa0","to_be_computed")
1910	ENDIF
1911	IF( td_namz%d_ppa1 /= NF90_FILL_DOUBLE )THEN
1912	ji=ji+1 ; tf_att(ji)=att_init("ppa1",td_namz%d_ppa1)
1913	ELSE
1914	ji=ji+1 ; tf_att(ji)=att_init("ppa1","to_be_computed")
1915	ENDIF
1916
1917	ji=ji+1 ; tf_att(ji)=att_init("ppkth",td_namz%d_ppkth)
1918	ji=ji+1 ; tf_att(ji)=att_init("ppacr",td_namz%d_ppacr)
1919	ji=ji+1 ; tf_att(ji)=att_init("ppdzmin",td_namz%d_ppdzmin)
1920	ji=ji+1 ; tf_att(ji)=att_init("pphmax",td_namz%d_pphmax)
1921
1922	IF( td_namz%l_dbletanh )THEN
1923	ji=ji+1 ; tf_att(ji)=att_init("ppa2",td_namz%d_ppa2)
1924	ji=ji+1 ; tf_att(ji)=att_init("ppkth2",td_namz%d_ppkth2)
1925	ji=ji+1 ; tf_att(ji)=att_init("ppacr2",td_namz%d_ppacr2)
1926	ENDIF
1927
1928	IF( td_namz%l_zps )THEN
1929	ji=ji+1 ; tf_att(ji)=att_init("e3zps_min",td_namz%d_e3zps_min)
1930	ji=ji+1 ; tf_att(ji)=att_init("e3zps_rat",td_namz%d_e3zps_rat)
1931	ENDIF
1932
1933	IF( td_namz%l_sco )THEN
1934	ji=ji+1 ; tf_att(ji)=att_init("sbot_min",td_namz%d_sbot_min)
1935	ji=ji+1 ; tf_att(ji)=att_init("sbot_max",td_namz%d_sbot_max)
1936	ji=ji+1 ; tf_att(ji)=att_init("hc",td_namz%d_hc)
1937	IF( td_namz%l_s_sh94 )THEN
1938	ji=ji+1 ; tf_att(ji)=att_init("rmax",td_namz%d_rmax)
1939	ji=ji+1 ; tf_att(ji)=att_init("theta",td_namz%d_theta)
1940	ji=ji+1 ; tf_att(ji)=att_init("thetb",td_namz%d_thetb)
1941	ji=ji+1 ; tf_att(ji)=att_init("bb",td_namz%d_bb)
1942	ELSEIF( td_namz%l_s_sf12 )THEN
1943	IF( td_namz%l_sigcrit ) ji=ji+1 ; tf_att(ji)=att_init("sigma_below_critical_depth","activated")
1944	ji=ji+1 ; tf_att(ji)=att_init("alpha",td_namz%d_alpha)
1945	ji=ji+1 ; tf_att(ji)=att_init("efold",td_namz%d_efold)
1946	ji=ji+1 ; tf_att(ji)=att_init("zs",td_namz%d_zs)
1947	ji=ji+1 ; tf_att(ji)=att_init("zb_a",td_namz%d_zb_a)
1948	ji=ji+1 ; tf_att(ji)=att_init("zb_b",td_namz%d_zb_b)
1949	ENDIF
1950	ENDIF
1951
1952	IF( td_namz%l_wd )THEN
1953	ji=ji+1 ; tf_att(ji)=att_init("wetting_drying","activated")
1954	ENDIF
1955
1956	IF( td_namz%l_wd )THEN
1957	ji=ji+1 ; tf_att(ji)=att_init("wdmin1",td_namz%d_wdmin1)
1958	ji=ji+1 ; tf_att(ji)=att_init("wdmin2",td_namz%d_wdmin2)
1959	ji=ji+1 ; tf_att(ji)=att_init("wdld",td_namz%d_wdld)
1960	ENDIF
1961
1962	END FUNCTION create_meshmask__gloatt
1963	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1964	END PROGRAM

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

Download in other formats: