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 branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/create_meshmask.f90 @ 7233

Last change on this file since 7233 was 7233, checked in by jpaul, 7 years ago
see ticket #1781
File size: 64.9 KB

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

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

Download in other formats: