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

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

Last change on this file since 7235 was 7235, checked in by jpaul, 7 years ago

see ticket #1781

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