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

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months 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!----------------------------------------------------------------------
354PROGRAM 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()
1469CONTAINS
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   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1964END PROGRAM
Note: See TracBrowser for help on using the repository browser.