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.
Changeset 5870 for branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Ignore:
Timestamp:
2015-11-09T18:33:54+01:00 (9 years ago)
Author:
acc
Message:

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.

Location:
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r5506 r5870  
    7171      !!                                   =2 put at location runoff 
    7272      !!---------------------------------------------------------------------- 
    73       INTEGER ::   jc            ! dummy loop indices 
    74       INTEGER :: isrow           ! local index 
    75       !!---------------------------------------------------------------------- 
    76        
     73      INTEGER ::   jc      ! dummy loop indices 
     74      INTEGER ::   isrow   ! local index 
     75      !!---------------------------------------------------------------------- 
     76      ! 
    7777      IF(lwp) WRITE(numout,*) 
    7878      IF(lwp) WRITE(numout,*)'dom_clo : closed seas ' 
    7979      IF(lwp) WRITE(numout,*)'~~~~~~~' 
    80  
     80      ! 
    8181      ! initial values 
    8282      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1 
    8383      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1 
    84  
     84      ! 
    8585      ! set the closed seas (in data domain indices) 
    8686      ! ------------------- 
    87  
     87      ! 
    8888      IF( cp_cfg == "orca" ) THEN 
    8989         ! 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5842 r5870  
    77   !! History :  1.0  ! 2005-10  (A. Beckmann, G. Madec)  reactivate s-coordinate  
    88   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
    9    !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
     9   !!            3.4  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
    1010   !!            3.5  ! 2012     (S. Mocavero, I. Epicoco) Add arrays associated 
    1111   !!                             to the optimization of BDY communications 
    12    !!            3.6.?! 2014     (H. Liu) Add arrays associated with Wetting and Drying 
     12   !!            3.7  ! 2015-11  (G. Madec) introduce surface and scale factor ratio 
     13   !!                 !          (H. Liu) Add arrays associated with Wetting and Drying 
    1314   !!---------------------------------------------------------------------- 
    1415 
    … …  
    2122 
    2223   IMPLICIT NONE 
    23    PUBLIC             ! allows the acces to par_oce when dom_oce is used 
    24    !                  ! exception to coding rules... to be suppressed ??? 
     24   PUBLIC             ! allows the acces to par_oce when dom_oce is used (exception to coding rules) 
    2525 
    2626   PUBLIC dom_oce_alloc  ! Called from nemogcm.F90 
    … …  
    108108   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dtra  !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 
    109109 
    110    !                                         !!* Namelist namcla : cross land advection 
    111    INTEGER, PUBLIC ::   nn_cla               !: =1 cross land advection for exchanges through some straits (ORCA2) 
    112  
    113110   !!---------------------------------------------------------------------- 
    114111   !! space domain parameters 
    … …  
    159156   !! horizontal curvilinear coordinate and scale factors 
    160157   !! --------------------------------------------------------------------- 
    161    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  glamt, glamu   !: longitude of t-, u-, v- and f-points (degre) 
    162    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  glamv, glamf   !: 
    163    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphit, gphiu   !: latitude  of t-, u-, v- and f-points (degre) 
    164    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphiv, gphif   !: 
    165    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1t, e2t, r1_e1t, r1_e2t   !: horizontal scale factors and inverse at t-point (m) 
    166    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1u, e2u, r1_e1u, r1_e2u   !: horizontal scale factors and inverse at u-point (m) 
    167    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1v, e2v, r1_e1v, r1_e2v   !: horizontal scale factors and inverse at v-point (m) 
    168    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1f, e2f, r1_e1f, r1_e2f   !: horizontal scale factors and inverse at f-point (m) 
    169    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  e1e2t          !: surface at t-point (m2) 
    170    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 
     158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   glamt , glamu, glamv , glamf    !: longitude at t, u, v, f-points [degree] 
     159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gphit , gphiu, gphiv , gphif    !: latitude  at t, u, v, f-points [degree] 
     160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::   e1t   , e2t  , r1_e1t, r1_e2t   !: t-point horizontal scale factors    [m] 
     161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::   e1u   , e2u  , r1_e1u, r1_e2u   !: horizontal scale factors at u-point [m] 
     162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::   e1v   , e2v  , r1_e1v, r1_e2v   !: horizontal scale factors at v-point [m] 
     163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::   e1f   , e2f  , r1_e1f, r1_e2f   !: horizontal scale factors at f-point [m] 
     164   ! 
     165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2t , r1_e1e2t                !: associated metrics at t-point 
     166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2u , r1_e1e2u , e2_e1u       !: associated metrics at u-point 
     167   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2v , r1_e1e2v , e1_e2v       !: associated metrics at v-point 
     168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point 
     169   ! 
     170   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff                              !: coriolis factor                   [1/s] 
    171171 
    172172   !!---------------------------------------------------------------------- 
    … …  
    217217   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0           !: reference depth at t-       points (meters) 
    218218   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters) 
    219    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re2u_e1u       !: scale factor coeffs at u points (e2u/e1u) 
    220    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re1v_e2v       !: scale factor coeffs at v points (e1v/e2v) 
    221    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12t , r1_e12t !: horizontal cell surface and inverse at t points 
    222    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12u , r1_e12u !: horizontal cell surface and inverse at u points 
    223    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12v , r1_e12v !: horizontal cell surface and inverse at v points 
    224    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12f , r1_e12f !: horizontal cell surface and inverse at f points 
    225219 
    226220   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
    … …  
    266260 
    267261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    268  
    269 #if defined key_noslip_accurate 
    270    INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  )  :: npcoa              !: ??? 
    271    INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ??? 
    272 #endif 
    273  
    274    !!---------------------------------------------------------------------- 
    275    !! critical depths,filters, limiters,and masks for  Wetting and Drying 
    276    !! --------------------------------------------------------------------- 
    277  
    278    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wduflt, wdvflt     !: u- and v- filter 
    279    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter  
    280  
    281    LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off) 
    282    REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells 
    283    REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
    284    REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
    285    INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter 
    286  
    287    !LOGICAL,  PUBLIC, SAVE ::   ln_wd     =  .FALSE.  !: turn on wetting/drying (T: on, F: off) 
    288    !REAL(wp), PUBLIC, SAVE ::   rn_wdmin1 = 0.10_wp   !: minimum water depth on dried cells 
    289    !REAL(wp), PUBLIC, SAVE ::   rn_wdmin2 = 0.01_wp   !: tolerrance of minimum water depth on dried cells 
    290    !REAL(wp), PUBLIC, SAVE ::   rn_wdld   = 20.0_wp   !: land elevation below which wetting/drying will be considered 
    291    !INTEGER , PUBLIC, SAVE ::   nn_wdit   =   10      !: maximum number of iteration for W/D limiter 
    292262 
    293263   !!---------------------------------------------------------------------- 
    … …  
    353323   INTEGER FUNCTION dom_oce_alloc() 
    354324      !!---------------------------------------------------------------------- 
    355       INTEGER, DIMENSION(12) :: ierr 
     325      INTEGER, DIMENSION(13) :: ierr 
    356326      !!---------------------------------------------------------------------- 
    357327      ierr(:) = 0 
    … …  
    366336         &      tpol(jpiglo)  , fpol(jpiglo)                               , STAT=ierr(2) ) 
    367337         ! 
    368       ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) ,   &  
    369          &      glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) ,   &   
    370          &      glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) ,   &   
    371          &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) ,   & 
    372          &      e1e2t(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )      
     338      ALLOCATE( glamt(jpi,jpj) ,    glamu(jpi,jpj) ,  glamv(jpi,jpj) ,  glamf(jpi,jpj) ,     & 
     339         &      gphit(jpi,jpj) ,    gphiu(jpi,jpj) ,  gphiv(jpi,jpj) ,  gphif(jpi,jpj) ,     & 
     340         &       e1t (jpi,jpj) ,     e2t (jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) ,     & 
     341         &       e1u (jpi,jpj) ,     e2u (jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) ,     & 
     342         &       e1v (jpi,jpj) ,     e2v (jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) ,     & 
     343         &       e1f (jpi,jpj) ,     e2f (jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) ,     & 
     344         &      e1e2t(jpi,jpj) , r1_e1e2t(jpi,jpj)                                     ,     & 
     345         &      e1e2u(jpi,jpj) , r1_e1e2u(jpi,jpj) , e2_e1u(jpi,jpj)                   ,     & 
     346         &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     & 
     347         &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     & 
     348         &        ff (jpi,jpj)                                                         , STAT=ierr(3) ) 
    373349         ! 
    374350      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         & 
    … …  
    384360         &      gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) ,                           & 
    385361         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) ,                           & 
    386          &      ehu_a    (jpi,jpj)    , ehv_a (jpi,jpj),                                                     & 
    387          &      ehur_a   (jpi,jpj)    , ehvr_a (jpi,jpj),                                                     & 
    388          &      ehu_b    (jpi,jpj)    , ehv_b (jpi,jpj),                                                     & 
    389          &      ehur_b   (jpi,jpj)    , ehvr_b (jpi,jpj),                                  STAT=ierr(5) )                           
     362         &      ehu_a   (jpi,jpj)     , ehv_a (jpi,jpj),                                                     & 
     363         &      ehur_a  (jpi,jpj)     , ehvr_a(jpi,jpj),                                                     & 
     364         &      ehu_b   (jpi,jpj)     , ehv_b (jpi,jpj),                                                     & 
     365         &      ehur_b  (jpi,jpj)     , ehvr_b(jpi,jpj),                                  STAT=ierr(5) )                           
    390366#endif 
    391367         ! 
    392       ALLOCATE( hu      (jpi,jpj) , hur     (jpi,jpj) , hu_0(jpi,jpj) , ht_0  (jpi,jpj) ,     & 
    393          &      hv      (jpi,jpj) , hvr     (jpi,jpj) , hv_0(jpi,jpj) , ht    (jpi,jpj) ,     & 
    394          &      re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) ,                                       & 
    395          &      e12t    (jpi,jpj) , r1_e12t (jpi,jpj) ,                                       & 
    396          &      e12u    (jpi,jpj) , r1_e12u (jpi,jpj) ,                                       & 
    397          &      e12v    (jpi,jpj) , r1_e12v (jpi,jpj) ,                                       & 
    398          &      e12f    (jpi,jpj) , r1_e12f (jpi,jpj) ,                                   STAT=ierr(6)  ) 
     368      ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , ht_0(jpi,jpj) ,     & 
     369         &      hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , ht  (jpi,jpj) , STAT=ierr(6)  ) 
    399370         ! 
    400371      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) ,                                     & 
    … …  
    407378         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     & 
    408379         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     & 
    409          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
     380         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 
    410381 
    411382      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      & 
    412383         &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 
    413          &     bmask(jpi,jpj)  ,                                                       & 
     384         &     bmask  (jpi,jpj) ,                                                       & 
    414385         &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
    415386 
    416387! (ISF) Allocation of basic array    
    417       ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),     & 
    418          &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
    419          &     mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
     388      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),                   & 
     389         &      mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
     390         &      mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
    420391 
    421392      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  
    … …  
    424395      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
    425396 
    426 #if defined key_noslip_accurate 
    427       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
    428 #endif 
    429  
    430       IF(ln_wd) & 
    431       ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr(12) ) 
    432397      ! 
    433398      dom_oce_alloc = MAXVAL(ierr) 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5363 r5870  
    1919   !!   dom_nam        : read and contral domain namelists 
    2020   !!   dom_ctl        : control print for the ocean domain 
     21   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 
    2122   !!---------------------------------------------------------------------- 
    2223   USE oce             ! ocean variables 
    … …  
    2526   USE phycst          ! physical constants 
    2627   USE closea          ! closed seas 
    27    USE in_out_manager  ! I/O manager 
    28    USE lib_mpp         ! distributed memory computing library 
    29  
    3028   USE domhgr          ! domain: set the horizontal mesh 
    3129   USE domzgr          ! domain: set the vertical mesh 
    … …  
    3634   USE c1d             ! 1D vertical configuration 
    3735   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine) 
     36   ! 
     37   USE in_out_manager  ! I/O manager 
     38   USE lib_mpp         ! distributed memory computing library 
     39   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    3840   USE timing          ! Timing 
    39    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    4041 
    4142   IMPLICIT NONE 
    … …  
    8182      ENDIF 
    8283      ! 
    83                              CALL dom_nam      ! read namelist ( namrun, namdom, namcla ) 
     84                             CALL dom_nam      ! read namelist ( namrun, namdom ) 
    8485                             CALL dom_clo      ! Closed seas and lake 
    8586                             CALL dom_hgr      ! Horizontal mesh 
    … …  
    8889      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency 
    8990      ! 
    90       ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points 
    91       hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points 
    92       hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points 
     91      ht_0(:,:) = 0._wp                        ! Reference ocean depth at T-points 
     92      hu_0(:,:) = 0._wp                        ! Reference ocean depth at U-points 
     93      hv_0(:,:) = 0._wp                        ! Reference ocean depth at V-points 
    9394      DO jk = 1, jpk 
    9495         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
    … …  
    9798      END DO 
    9899      ! 
    99       IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh 
     100      IF( lk_vvl         )   CALL dom_vvl_init ! Vertical variable mesh 
    100101      ! 
    101102      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point 
    … …  
    131132      !! ** input   : - namrun namelist 
    132133      !!              - namdom namelist 
    133       !!              - namcla namelist 
    134134      !!              - namnc4 namelist   ! "key_netcdf4" only 
    135135      !!---------------------------------------------------------------------- 
    … …  
    146146         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 
    147147         &             ppa2, ppkth2, ppacr2 
    148       NAMELIST/namcla/ nn_cla 
    149148#if defined key_netcdf4 
    150149      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
    … …  
    155154      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run 
    156155      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901) 
    157 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp ) 
     156901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp ) 
    158157 
    159158      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run 
    160159      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 ) 
    161 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 
     160902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 
    162161      IF(lwm) WRITE ( numond, namrun ) 
    163162      ! 
    … …  
    251250904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    252251      IF(lwm) WRITE ( numond, namdom ) 
    253  
     252      ! 
    254253      IF(lwp) THEN 
    255254         WRITE(numout,*) 
    … …  
    293292         WRITE(numout,*) '                                      ppacr2            = ', ppacr2 
    294293      ENDIF 
    295  
     294      ! 
    296295      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon) 
    297296      e3zps_min = rn_e3zps_min 
    … …  
    304303      rdtmax    = rn_rdtmin 
    305304      rdth      = rn_rdth 
    306  
    307       REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection 
    308       READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905) 
    309 905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp ) 
    310  
    311       REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection 
    312       READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 ) 
    313 906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp ) 
    314       IF(lwm) WRITE( numond, namcla ) 
    315  
    316       IF(lwp) THEN 
    317          WRITE(numout,*) 
    318          WRITE(numout,*) '   Namelist namcla' 
    319          WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla 
    320       ENDIF 
    321       IF ( nn_cla .EQ. 1 ) THEN 
    322          IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2  
    323             CONTINUE 
    324          ELSE 
    325             CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' ) 
    326          ENDIF 
    327       ENDIF 
    328305 
    329306#if defined key_netcdf4 
    … …  
    409386   END SUBROUTINE dom_ctl 
    410387 
     388 
    411389   SUBROUTINE dom_stiff 
    412390      !!---------------------------------------------------------------------- 
    … …  
    427405      REAL(wp), DIMENSION(4) :: zr1 
    428406      !!---------------------------------------------------------------------- 
    429       rx1(:,:) = 0.e0 
    430       zrxmax   = 0.e0 
    431       zr1(:)   = 0.e0 
    432        
     407      rx1(:,:) = 0._wp 
     408      zrxmax   = 0._wp 
     409      zr1(:)   = 0._wp 
     410      ! 
    433411      DO ji = 2, jpim1 
    434412         DO jj = 2, jpjm1 
    … …  
    455433         END DO 
    456434      END DO 
    457  
    458435      CALL lbc_lnk( rx1, 'T', 1. ) 
    459  
    460       zrxmax = MAXVAL(rx1) 
    461  
     436      ! 
     437      zrxmax = MAXVAL( rx1 ) 
     438      ! 
    462439      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain 
    463  
     440      ! 
    464441      IF(lwp) THEN 
    465442         WRITE(numout,*) 
    … …  
    467444         WRITE(numout,*) '~~~~~~~~~' 
    468445      ENDIF 
    469  
     446      ! 
    470447   END SUBROUTINE dom_stiff 
    471  
    472  
    473448 
    474449   !!====================================================================== 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5656 r5870  
    1414   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist 
    1515   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration  
    16    !!            4.0  ! 2011-02  (G. Madec) add cell surface (e1e2t) 
     16   !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse 
     17   !!                                       add optional read of e1e2u & e1e2v 
    1718   !!---------------------------------------------------------------------- 
    1819 
    … …  
    2324   USE dom_oce        ! ocean space and time domain 
    2425   USE phycst         ! physical constants 
     26   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files 
     27   ! 
    2528   USE in_out_manager ! I/O manager 
    2629   USE lib_mpp        ! MPP library 
    … …  
    3538 
    3639   !!---------------------------------------------------------------------- 
    37    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     40   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    3841   !! $Id$  
    3942   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    106109      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
    107110      INTEGER  ::   isrow                ! index for ORCA1 starting row 
    108  
     111      INTEGER  ::   ie1e2u_v             ! fag for u- & v-surface read in coordinate file or not 
    109112      !!---------------------------------------------------------------------- 
    110113      ! 
    … …  
    122125         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m   
    123126      ENDIF 
    124  
    125  
    126       SELECT CASE( jphgr_msh )   ! type of horizontal mesh 
    127  
    128       CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file 
    129  
     127      ! 
     128      ! 
     129      SELECT CASE( jphgr_msh )   !  type of horizontal mesh   
     130      ! 
     131      CASE ( 0 )                     !==  read in coordinate.nc file  ==! 
     132         ! 
    130133         IF(lwp) WRITE(numout,*) 
    131134         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file' 
    132  
    133          CALL hgr_read           ! Defaultl option  :   NetCDF file 
    134  
    135          !                                                ! ===================== 
    136          IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    137             !                                             ! ===================== 
    138             IF( nn_cla == 0 ) THEN 
    139                ! 
    140                ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km) 
    141                ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    142                IF(lwp) WRITE(numout,*) 
    143                IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km' 
    144                ! 
    145                ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km) 
    146                ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3 
    147                                                e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3 
    148                IF(lwp) WRITE(numout,*) 
    149                IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km' 
    150                IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km' 
    151             ENDIF 
    152  
    153             ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km) 
    154             ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    155             IF(lwp) WRITE(numout,*) 
    156             IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km' 
    157             ! 
     135         ! 
     136         ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v 
     137         ! 
     138         CALL hgr_read( ie1e2u_v )     ! read the coordinate.nc file 
     139         ! 
     140         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them 
     141            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction 
     142            e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
     143            e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
    158144         ENDIF 
    159  
    160             !                                             ! ===================== 
    161          IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    162             !                                             ! ===================== 
    163             ! This dirty section will be suppressed by simplification process: all this will come back in input files 
    164             ! Currently these hard-wired indices relate to configuration with 
    165             ! extend grid (jpjglo=332) 
    166             ! which had a grid-size of 362x292. 
    167             !  
    168             isrow = 332 - jpjglo 
    169             ! 
    170             ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km) 
    171             ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    172             IF(lwp) WRITE(numout,*) 
    173             IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km' 
    174  
    175             ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km) 
    176             ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    177             IF(lwp) WRITE(numout,*) 
    178             IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km' 
    179  
    180             ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km) 
    181             ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3 
    182             IF(lwp) WRITE(numout,*) 
    183             IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km' 
    184  
    185             ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 
    186             ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3 
    187             IF(lwp) WRITE(numout,*) 
    188             IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km' 
    189  
    190             ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km) 
    191             ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 
    192             IF(lwp) WRITE(numout,*) 
    193             IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km' 
    194  
    195             ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km) 
    196             ij0 = 164 - isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 
    197             IF(lwp) WRITE(numout,*) 
    198             IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km' 
    199  
    200             ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km) 
    201             ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 
    202             IF(lwp) WRITE(numout,*) 
    203             IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km' 
    204  
    205             ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km) 
    206             ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 
    207             IF(lwp) WRITE(numout,*) 
    208             IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km' 
    209             ! 
    210             ! 
    211          ENDIF 
    212  
    213          !                                                ! ====================== 
    214          IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration 
    215             !                                             ! ====================== 
    216             ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km) 
    217             ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    218             IF(lwp) WRITE(numout,*) 
    219             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait' 
    220             ! 
    221             ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km) 
    222             ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    223             IF(lwp) WRITE(numout,*) 
    224             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait' 
    225             ! 
    226             ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km) 
    227             ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3 
    228             IF(lwp) WRITE(numout,*) 
    229             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait' 
    230             ! 
    231             ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km) 
    232             ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3 
    233             IF(lwp) WRITE(numout,*) 
    234             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait' 
    235             ! 
    236             ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km) 
    237             ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    238             IF(lwp) WRITE(numout,*) 
    239             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait' 
    240             ! 
    241             ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km) 
    242             ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    243             IF(lwp) WRITE(numout,*) 
    244             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait' 
    245             ! 
    246             ! 
    247             ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km) 
    248             ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3 
    249             IF(lwp) WRITE(numout,*) 
    250             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb' 
    251             ! 
    252          ENDIF 
    253  
    254  
    255          ! N.B. :  General case, lat and long function of both i and j indices: 
    256          !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   & 
    257          !                                  + (                           fsdiph( zti, ztj ) )**2  ) 
    258          !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   & 
    259          !                                  + (                           fsdiph( zui, zuj ) )**2  ) 
    260          !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   & 
    261          !                                  + (                           fsdiph( zvi, zvj ) )**2  ) 
    262          !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   & 
    263          !                                  + (                           fsdiph( zfi, zfj ) )**2  ) 
    264          ! 
    265          !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   & 
    266          !                                  + (                           fsdjph( zti, ztj ) )**2  ) 
    267          !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   & 
    268          !                                  + (                           fsdjph( zui, zuj ) )**2  ) 
    269          !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   & 
    270          !                                  + (                           fsdjph( zvi, zvj ) )**2  ) 
    271          !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   & 
    272          !                                  + (                           fsdjph( zfi, zfj ) )**2  ) 
    273  
    274  
    275       CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing 
    276  
     145         ! 
     146      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==! 
     147         ! 
    277148         IF(lwp) WRITE(numout,*) 
    278149         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing' 
    279150         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg'  
    280  
     151         ! 
    281152         DO jj = 1, jpj 
    282153            DO ji = 1, jpi 
    283                zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    284                zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 ) 
    285                zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5 
    286                zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5 
     154               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 ) 
     155               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 ) 
     156               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
     157               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
    287158         ! Longitude 
    288159               glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
    … …  
    307178            END DO 
    308179         END DO 
    309  
    310  
    311       CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing 
    312  
     180         ! 
     181      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==! 
     182         ! 
    313183         IF(lwp) WRITE(numout,*) 
    314184         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing' 
    315185         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'  
    316  
     186         ! 
    317187         ! Position coordinates (in kilometers) 
    318188         !                          ========== 
    319          glam0 = 0.e0 
     189         glam0 = 0._wp 
    320190         gphi0 = - ppe2_m * 1.e-3 
    321           
     191         ! 
    322192#if defined key_agrif  
    323193         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
    … …  
    332202         DO jj = 1, jpj 
    333203            DO ji = 1, jpi 
    334                glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       ) 
    335                glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ) 
     204               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       ) 
     205               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 ) 
    336206               glamv(ji,jj) = glamt(ji,jj) 
    337207               glamf(ji,jj) = glamu(ji,jj) 
    338     
    339                gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       ) 
     208               ! 
     209               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       ) 
    340210               gphiu(ji,jj) = gphit(ji,jj) 
    341                gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 ) 
     211               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 
    342212               gphif(ji,jj) = gphiv(ji,jj) 
    343213            END DO 
    344214         END DO 
    345  
     215         ! 
    346216         ! Horizontal scale factors (in meters) 
    347217         !                              ====== 
    … …  
    350220         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m 
    351221         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m 
    352  
    353       CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type 
    354  
     222         ! 
     223      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==! 
     224         ! 
    355225         IF(lwp) WRITE(numout,*) 
    356226         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type' 
    357227         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg' 
    358228         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 
    359  
     229         ! 
    360230         !  Find index corresponding to the equator, given the grid spacing e1_deg 
    361231         !  and the (approximate) southern latitude ppgphi0. 
    … …  
    365235         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 
    366236         IF(  ppgphi0 > 0 )  ijeq = -ijeq 
    367  
     237         ! 
    368238         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq 
    369  
     239         ! 
    370240         DO jj = 1, jpj 
    371241            DO ji = 1, jpi 
    372                zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 ) 
    373                zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 ) 
    374                zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5 
    375                zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5 
     242               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 ) 
     243               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 ) 
     244               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
     245               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
    376246         ! Longitude 
    377247               glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
    … …  
    396266            END DO 
    397267         END DO 
    398  
    399       CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration) 
    400  
     268         ! 
     269      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 
     270         ! 
    401271         IF(lwp) WRITE(numout,*) 
    402272         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 
    403273         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
    404  
     274         ! 
    405275         ! Position coordinates (in kilometers) 
    406276         !                          ========== 
    407  
     277         ! 
    408278         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 
    409          zlam1 = -85 
    410          zphi1 = 29 
     279         zlam1 = -85._wp 
     280         zphi1 =  29._wp 
    411281         ! resolution in meters 
    412          ze1 = 106000. / FLOAT(jp_cfg)             
     282         ze1 = 106000. / REAL( jp_cfg , wp )             
    413283         ! benchmark: forced the resolution to be about 100 km 
    414          IF( nbench /= 0 )   ze1 = 106000.e0      
    415          zsin_alpha = - SQRT( 2. ) / 2. 
    416          zcos_alpha =   SQRT( 2. ) / 2. 
     284         IF( nbench /= 0 )   ze1 = 106000._wp      
     285         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 
     286         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp 
    417287         ze1deg = ze1 / (ra * rad) 
    418          IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon 
    419          !                                                          ! at the right jp_cfg resolution 
    420          glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 ) 
    421          gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 ) 
    422  
     288         IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon 
     289         !                                                           ! at the right jp_cfg resolution 
     290         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     291         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     292         ! 
    423293         IF( nprint==1 .AND. lwp )   THEN 
    424294            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
    425295            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 
    426296         ENDIF 
    427  
     297         ! 
    428298         DO jj = 1, jpj 
    429            DO ji = 1, jpi 
    430              zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5 
    431              zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5 
    432  
    433              glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
    434              gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
    435  
    436              glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
    437              gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
    438  
    439              glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
    440              gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
    441  
    442              glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
    443              gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
    444            END DO 
    445           END DO 
    446  
     299            DO ji = 1, jpi 
     300               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5 
     301               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5 
     302               ! 
     303               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     304               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     305               ! 
     306               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     307               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     308               ! 
     309               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     310               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     311               ! 
     312               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     313               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     314            END DO 
     315         END DO 
     316         ! 
    447317         ! Horizontal scale factors (in meters) 
    448318         !                              ====== 
    … …  
    451321         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1 
    452322         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1 
    453  
     323         ! 
    454324      CASE DEFAULT 
    455325         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh 
    456326         CALL ctl_stop( ctmp1 ) 
    457  
     327         ! 
    458328      END SELECT 
    459329       
    460       ! T-cell surface 
    461       ! -------------- 
    462       e1e2t(:,:) = e1t(:,:) * e2t(:,:) 
    463      
    464       ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning) 
    465       ! --------------------------------------------------------------------------- 
    466       e12t    (:,:) = e1t(:,:) * e2t(:,:) 
    467       e12u    (:,:) = e1u(:,:) * e2u(:,:) 
    468       e12v    (:,:) = e1v(:,:) * e2v(:,:) 
    469       e12f    (:,:) = e1f(:,:) * e2f(:,:) 
    470       r1_e12t (:,:) = 1._wp    / e12t(:,:) 
    471       r1_e12u (:,:) = 1._wp    / e12u(:,:) 
    472       r1_e12v (:,:) = 1._wp    / e12v(:,:) 
    473       r1_e12f (:,:) = 1._wp    / e12f(:,:) 
    474       re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 
    475       re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 
    476       r1_e1t  (:,:) = 1._wp    / e1t(:,:) 
    477       r1_e1u  (:,:) = 1._wp    / e1u(:,:) 
    478       r1_e1v  (:,:) = 1._wp    / e1v(:,:) 
    479       r1_e1f  (:,:) = 1._wp    / e1f(:,:) 
    480       r1_e2t  (:,:) = 1._wp    / e2t(:,:) 
    481       r1_e2u  (:,:) = 1._wp    / e2u(:,:) 
    482       r1_e2v  (:,:) = 1._wp    / e2v(:,:) 
    483       r1_e2f  (:,:) = 1._wp    / e2f(:,:) 
    484  
    485       ! Control printing : Grid informations (if not restart) 
    486       ! ---------------- 
    487  
    488       IF( lwp .AND. .NOT.ln_rstart ) THEN 
     330      ! associated horizontal metrics 
     331      ! ----------------------------- 
     332      ! 
     333      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:) 
     334      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:) 
     335      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:) 
     336      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:) 
     337      ! 
     338      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 
     339      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 
     340      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them 
     341         e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
     342         e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
     343      ENDIF 
     344      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases 
     345      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) 
     346      !    
     347      e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 
     348      e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 
     349 
     350      IF( lwp .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart) 
    489351         WRITE(numout,*) 
    490352         WRITE(numout,*) '          longitude and e1 scale factors' 
    … …  
    4963589300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
    497359            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 
    498           
     360            ! 
    499361         WRITE(numout,*) 
    500362         WRITE(numout,*) '          latitude and e2 scale factors' 
    … …  
    506368      ENDIF 
    507369 
    508        
    509       IF( nprint == 1 .AND. lwp ) THEN 
    510          WRITE(numout,*) '          e1u e2u ' 
    511          CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    512          CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    513          WRITE(numout,*) '          e1v e2v  ' 
    514          CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    515          CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    516          WRITE(numout,*) '          e1f e2f  ' 
    517          CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    518          CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    519       ENDIF 
    520  
    521370 
    522371      ! ================= ! 
    … …  
    525374 
    526375      SELECT CASE( jphgr_msh )   ! type of horizontal mesh 
    527  
     376      ! 
    528377      CASE ( 0, 1, 4 )               ! mesh on the sphere 
    529  
     378         ! 
    530379         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
    531  
     380         ! 
    532381      CASE ( 2 )                     ! f-plane at ppgphi0  
    533  
     382         ! 
    534383         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
    535  
     384         ! 
    536385         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1) 
    537  
     386         ! 
    538387      CASE ( 3 )                     ! beta-plane 
    539  
     388         ! 
    540389         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0 
    541          zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points 
    542           
     390         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points 
     391         ! 
    543392#if defined key_agrif 
    544393         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
    545394            IF( .NOT. Agrif_Root() ) THEN 
    546               zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m)   &  
    547                     &           / (ra * rad) 
     395              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 
    548396            ENDIF 
    549397         ENDIF 
    550398#endif          
    551399         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south 
    552  
     400         ! 
    553401         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
    554           
     402         ! 
    555403         IF(lwp) THEN 
    556404            WRITE(numout,*)  
    … …  
    565413            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    566414         END IF 
    567  
     415         ! 
    568416      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration) 
    569  
     417         ! 
    570418         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0 
    571          zphi0 = 15.e0                                                      ! latitude of the first row F-points 
     419         zphi0 = 15._wp                                                     ! latitude of the first row F-points 
    572420         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south 
    573  
     421         ! 
    574422         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
    575  
     423         ! 
    576424         IF(lwp) THEN 
    577425            WRITE(numout,*)  
    … …  
    579427            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 
    580428         ENDIF 
    581  
     429         ! 
    582430         IF( lk_mpp ) THEN  
    583431            zminff=ff(nldi,nldj) 
    … …  
    587435            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    588436         END IF 
    589  
     437         ! 
    590438      END SELECT 
    591439 
    … …  
    596444 
    597445      IF( nperio == 2 ) THEN 
    598          znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi ) 
     446         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
    599447         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 
    600448      ENDIF 
    … …  
    605453 
    606454 
    607    SUBROUTINE hgr_read 
     455   SUBROUTINE hgr_read( ke1e2u_v ) 
    608456      !!--------------------------------------------------------------------- 
    609457      !!              ***  ROUTINE hgr_read  *** 
    610458      !! 
    611       !! ** Purpose :   Read a coordinate file in NetCDF format  
    612       !! 
    613       !! ** Method  :   The mesh file has been defined trough a analytical  
    614       !!      or semi-analytical method. It is read in a NetCDF file.  
    615       !!      
     459      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM 
     460      !! 
    616461      !!---------------------------------------------------------------------- 
    617462      USE iom 
    618  
     463      !! 
     464      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 
     465      ! 
    619466      INTEGER ::   inum   ! temporary logical unit 
    620467      !!---------------------------------------------------------------------- 
    621  
     468      ! 
    622469      IF(lwp) THEN 
    623470         WRITE(numout,*) 
    … …  
    625472         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 
    626473      ENDIF 
    627        
     474      ! 
    628475      CALL iom_open( 'coordinates', inum ) 
    629        
     476      ! 
    630477      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 
    631478      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 
    632479      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 
    633480      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 
    634        
     481      ! 
    635482      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 
    636483      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 
    637484      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 
    638485      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 
    639        
    640       CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 
    641       CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 
    642       CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 
    643       CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 
    644        
    645       CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 
    646       CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 
    647       CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 
    648       CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 
    649        
     486      ! 
     487      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr ) 
     488      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr ) 
     489      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr ) 
     490      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr ) 
     491      ! 
     492      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr ) 
     493      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr ) 
     494      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr ) 
     495      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr ) 
     496      ! 
     497      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN 
     498         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file' 
     499         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr ) 
     500         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr ) 
     501         ke1e2u_v = 1 
     502      ELSE 
     503         ke1e2u_v = 0 
     504      ENDIF 
     505      ! 
    650506      CALL iom_close( inum ) 
    651507       
     508!!gm   THIS is TO BE REMOVED !!!!!!! 
     509 
    652510! need to be define for the extended grid south of -80S 
    653511! some point are undefined but you need to have e1 and e2 .NE. 0 
    … …  
    676534         e2f=1.0e2 
    677535      END WHERE 
     536!!gm end 
    678537        
    679538    END SUBROUTINE hgr_read 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5552 r5870  
    1717   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization 
    1818   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
     19   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask 
    1920   !!---------------------------------------------------------------------- 
    2021 
    2122   !!---------------------------------------------------------------------- 
    2223   !!   dom_msk        : compute land/ocean mask 
    23    !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used. 
    2424   !!---------------------------------------------------------------------- 
    2525   USE oce             ! ocean dynamics and tracers 
    … …  
    3636 
    3737   PUBLIC   dom_msk         ! routine called by inidom.F90 
    38    PUBLIC   dom_msk_alloc   ! routine called by nemogcm.F90 
    3938 
    4039   !                            !!* Namelist namlbc : lateral boundary condition * 
    … …  
    4241   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition  
    4342   !                                            with analytical eqs. 
    44  
    45  
    46    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  icoord ! Workspace for dom_msk_nsa() 
    4743 
    4844   !! * Substitutions 
    … …  
    5450   !!---------------------------------------------------------------------- 
    5551CONTAINS 
    56     
    57    INTEGER FUNCTION dom_msk_alloc() 
    58       !!--------------------------------------------------------------------- 
    59       !!                 ***  FUNCTION dom_msk_alloc  *** 
    60       !!--------------------------------------------------------------------- 
    61       dom_msk_alloc = 0 
    62 #if defined key_noslip_accurate 
    63       ALLOCATE(icoord(jpi*jpj*jpk,3), STAT=dom_msk_alloc) 
    64 #endif 
    65       IF( dom_msk_alloc /= 0 )   CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array') 
    66       ! 
    67    END FUNCTION dom_msk_alloc 
    68  
    6952 
    7053   SUBROUTINE dom_msk 
    … …  
    129112      !!               tmask_i  : interior ocean mask 
    130113      !!---------------------------------------------------------------------- 
    131       ! 
    132       INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     114      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    133115      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers 
    134116      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
    … …  
    199181      END DO   
    200182 
    201 !!gm  ???? 
    202 #if defined key_zdfkpp 
    203       IF( cp_cfg == 'orca' ) THEN 
    204          IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section 
    205             ij0 =  87   ;   ij1 =  88 
    206             ii0 = 160   ;   ii1 = 161 
    207             tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp 
    208          ELSE 
    209             IF(lwp) WRITE(numout,*) 
    210             IF(lwp) WRITE(numout,cform_war) 
    211             IF(lwp) WRITE(numout,*) 
    212             IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait' 
    213             IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations' 
    214             IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved' 
    215             IF(lwp) WRITE(numout,*) 
    216          ENDIF 
    217       ENDIF 
    218 #endif 
    219 !!gm end 
    220  
    221183      ! Interior domain mask (used for global sum) 
    222184      ! -------------------- 
    … …  
    284246      ! 3. Ocean/land mask at wu-, wv- and w points  
    285247      !---------------------------------------------- 
    286       wmask (:,:,1) = tmask(:,:,1) ! ???????? 
    287       wumask(:,:,1) = umask(:,:,1) ! ???????? 
    288       wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
    289       DO jk=2,jpk 
    290          wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
    291          wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
    292          wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
     248      wmask (:,:,1) = tmask(:,:,1)     ! surface 
     249      wumask(:,:,1) = umask(:,:,1) 
     250      wvmask(:,:,1) = vmask(:,:,1) 
     251      DO jk = 2, jpk                   ! interior values 
     252         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 
     253         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)    
     254         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 
    293255      END DO 
    294256 
    … …  
    339301      ENDIF 
    340302 
    341  
    342       ! mask for second order calculation of vorticity 
    343       ! ---------------------------------------------- 
    344       CALL dom_msk_nsa 
    345  
    346        
    347303      ! Lateral boundary conditions on velocity (modify fmask) 
    348304      ! ---------------------------------------      
    … …  
    377333      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration 
    378334         !                                                 ! Increased lateral friction near of some straits 
    379          IF( nn_cla == 0 ) THEN 
    380             !                                ! Gibraltar strait  : partial slip (fmask=0.5) 
    381             ij0 = 101   ;   ij1 = 101 
    382             ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    383             ij0 = 102   ;   ij1 = 102 
    384             ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    385             ! 
    386             !                                ! Bab el Mandeb : partial slip (fmask=1) 
    387             ij0 =  87   ;   ij1 =  88 
    388             ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    389             ij0 =  88   ;   ij1 =  88 
    390             ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    391             ! 
    392          ENDIF 
     335         !                                ! Gibraltar strait  : partial slip (fmask=0.5) 
     336         ij0 = 101   ;   ij1 = 101 
     337         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
     338         ij0 = 102   ;   ij1 = 102 
     339         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
     340         ! 
     341         !                                ! Bab el Mandeb : partial slip (fmask=1) 
     342         ij0 =  87   ;   ij1 =  88 
     343         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
     344         ij0 =  88   ;   ij1 =  88 
     345         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
     346         ! 
    393347         !                                ! Danish straits  : strong slip (fmask > 2) 
    394348! We keep this as an example but it is instable in this case  
    … …  
    500454      ! 
    501455   END SUBROUTINE dom_msk 
    502  
    503 #if defined key_noslip_accurate 
    504    !!---------------------------------------------------------------------- 
    505    !!   'key_noslip_accurate' :         accurate no-slip boundary condition 
    506    !!---------------------------------------------------------------------- 
    507     
    508    SUBROUTINE dom_msk_nsa 
    509       !!--------------------------------------------------------------------- 
    510       !!                 ***  ROUTINE dom_msk_nsa  *** 
    511       !!  
    512       !! ** Purpose : 
    513       !! 
    514       !! ** Method  : 
    515       !! 
    516       !! ** Action : 
    517       !!---------------------------------------------------------------------- 
    518       INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices 
    519       INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd 
    520       REAL(wp) ::   zaa 
    521       !!--------------------------------------------------------------------- 
    522       ! 
    523       IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa') 
    524       ! 
    525       IF(lwp) WRITE(numout,*) 
    526       IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition' 
    527       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme' 
    528       IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' ) 
    529  
    530       ! mask for second order calculation of vorticity 
    531       ! ---------------------------------------------- 
    532       ! noslip boundary condition: fmask=1  at convex corner, store 
    533       ! index of straight coast meshes ( 'west', refering to a coast, 
    534       ! means west of the ocean, aso) 
    535        
    536       DO jk = 1, jpk 
    537          DO jl = 1, 4 
    538             npcoa(jl,jk) = 0 
    539             DO ji = 1, 2*(jpi+jpj) 
    540                nicoa(ji,jl,jk) = 0 
    541                njcoa(ji,jl,jk) = 0 
    542             END DO 
    543          END DO 
    544       END DO 
    545        
    546       IF( jperio == 2 ) THEN 
    547          WRITE(numout,*) ' ' 
    548          WRITE(numout,*) ' symetric boundary conditions need special' 
    549          WRITE(numout,*) ' treatment not implemented. we stop.' 
    550          STOP 
    551       ENDIF 
    552        
    553       ! convex corners 
    554        
    555       DO jk = 1, jpkm1 
    556          DO jj = 1, jpjm1 
    557             DO ji = 1, jpim1 
    558                zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   & 
    559                   &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) 
    560                IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp 
    561             END DO 
    562          END DO 
    563       END DO 
    564  
    565       ! north-south straight coast 
    566  
    567       DO jk = 1, jpkm1 
    568          inw = 0 
    569          ine = 0 
    570          DO jj = 2, jpjm1 
    571             DO ji = 2, jpim1 
    572                zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) 
    573                IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    574                   inw = inw + 1 
    575                   nicoa(inw,1,jk) = ji 
    576                   njcoa(inw,1,jk) = jj 
    577                   IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj 
    578                ENDIF 
    579                zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk) 
    580                IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    581                   ine = ine + 1 
    582                   nicoa(ine,2,jk) = ji 
    583                   njcoa(ine,2,jk) = jj 
    584                   IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj 
    585                ENDIF 
    586             END DO 
    587          END DO 
    588          npcoa(1,jk) = inw 
    589          npcoa(2,jk) = ine 
    590       END DO 
    591  
    592       ! west-east straight coast 
    593  
    594       DO jk = 1, jpkm1 
    595          ins = 0 
    596          inn = 0 
    597          DO jj = 2, jpjm1 
    598             DO ji =2, jpim1 
    599                zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) 
    600                IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    601                   ins = ins + 1 
    602                   nicoa(ins,3,jk) = ji 
    603                   njcoa(ins,3,jk) = jj 
    604                   IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj 
    605                ENDIF 
    606                zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk) 
    607                IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    608                   inn = inn + 1 
    609                   nicoa(inn,4,jk) = ji 
    610                   njcoa(inn,4,jk) = jj 
    611                   IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj 
    612                ENDIF 
    613             END DO 
    614          END DO 
    615          npcoa(3,jk) = ins 
    616          npcoa(4,jk) = inn 
    617       END DO 
    618  
    619       itest = 2 * ( jpi + jpj ) 
    620       DO jk = 1, jpk 
    621          IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   & 
    622              npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN 
    623              
    624             WRITE(ctmp1,*) ' level jk = ',jk 
    625             WRITE(ctmp2,*) ' straight coast index arraies are too small.:' 
    626             WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   & 
    627                 &                                     npcoa(3,jk), npcoa(4,jk) 
    628             WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.' 
    629             CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 ) 
    630         ENDIF 
    631       END DO 
    632  
    633       ierror = 0 
    634       iind = 0 
    635       ijnd = 0 
    636       IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2 
    637       IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2 
    638       DO jk = 1, jpk 
    639          DO jl = 1, npcoa(1,jk) 
    640             IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN 
    641                ierror = ierror+1 
    642                icoord(ierror,1) = nicoa(jl,1,jk) 
    643                icoord(ierror,2) = njcoa(jl,1,jk) 
    644                icoord(ierror,3) = jk 
    645             ENDIF 
    646          END DO 
    647          DO jl = 1, npcoa(2,jk) 
    648             IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN 
    649                ierror = ierror + 1 
    650                icoord(ierror,1) = nicoa(jl,2,jk) 
    651                icoord(ierror,2) = njcoa(jl,2,jk) 
    652                icoord(ierror,3) = jk 
    653             ENDIF 
    654          END DO 
    655          DO jl = 1, npcoa(3,jk) 
    656             IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN 
    657                ierror = ierror + 1 
    658                icoord(ierror,1) = nicoa(jl,3,jk) 
    659                icoord(ierror,2) = njcoa(jl,3,jk) 
    660                icoord(ierror,3) = jk 
    661             ENDIF 
    662          END DO 
    663          DO jl = 1, npcoa(4,jk) 
    664             IF( njcoa(jl,4,jk)-2 < 1) THEN 
    665                ierror=ierror + 1 
    666                icoord(ierror,1) = nicoa(jl,4,jk) 
    667                icoord(ierror,2) = njcoa(jl,4,jk) 
    668                icoord(ierror,3) = jk 
    669             ENDIF 
    670          END DO 
    671       END DO 
    672        
    673       IF( ierror > 0 ) THEN 
    674          IF(lwp) WRITE(numout,*) 
    675          IF(lwp) WRITE(numout,*) '              Problem on lateral conditions' 
    676          IF(lwp) WRITE(numout,*) '                 Bad marking off at points:' 
    677          DO jl = 1, ierror 
    678             IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   & 
    679                &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')' 
    680          END DO 
    681          CALL ctl_stop( 'We stop...' ) 
    682       ENDIF 
    683       ! 
    684       IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa') 
    685       ! 
    686    END SUBROUTINE dom_msk_nsa 
    687  
    688 #else 
    689    !!---------------------------------------------------------------------- 
    690    !!   Default option :                                      Empty routine 
    691    !!---------------------------------------------------------------------- 
    692    SUBROUTINE dom_msk_nsa        
    693    END SUBROUTINE dom_msk_nsa 
    694 #endif 
    695456    
    696457   !!====================================================================== 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

    r3294 r5870  
    3737      !!                -> not good if located at too high latitude... 
    3838      !!---------------------------------------------------------------------- 
    39       ! 
    4039      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point 
    4140      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point 
    … …  
    4948      IF( nn_timing == 1 )  CALL timing_start('dom_ngb') 
    5049      ! 
    51       CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist ) 
     50      CALL wrk_alloc( jpi,jpj,  zglam, zgphi, zmask, zdist ) 
    5251      ! 
    5352      zmask(:,:) = 0._wp 
    … …  
    7675      ENDIF 
    7776      ! 
    78       CALL wrk_dealloc( jpi, jpj, zglam, zgphi, zmask, zdist ) 
     77      CALL wrk_dealloc( jpi,jpj,  zglam, zgphi, zmask, zdist ) 
    7978      ! 
    8079      IF( nn_timing == 1 )  CALL timing_stop('dom_ngb') 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5842 r5870  
    1010   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_vvl'                              variable volume 
    13    !!---------------------------------------------------------------------- 
     12 
    1413   !!---------------------------------------------------------------------- 
    1514   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
    … …  
    1918   !!   dom_vvl_rst      : read/write restart file 
    2019   !!   dom_vvl_ctl      : Check the vvl options 
    21    !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors  
    22    !!                    : to account for manual changes to e[1,2][u,v] in some Straits  
    2320   !!---------------------------------------------------------------------- 
    24    !! * Modules used 
    2521   USE oce             ! ocean dynamics and tracers 
    2622   USE dom_oce         ! ocean space and time domain 
    2723   USE sbc_oce         ! ocean surface boundary condition 
     24   USE wet_dry         ! wetting and drying 
    2825   USE in_out_manager  ! I/O manager 
    2926   USE iom             ! I/O manager library 
    … …  
    3734   PRIVATE 
    3835 
    39    !! * Routine accessibility 
    4036   PUBLIC  dom_vvl_init       ! called by domain.F90 
    4137   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    4238   PUBLIC  dom_vvl_sf_swp     ! called by step.F90 
    4339   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    44    PRIVATE dom_vvl_orca_fix   ! called by dom_vvl_interpol 
    45  
    46    !!* Namelist nam_vvl 
    47    LOGICAL , PUBLIC                                      :: ln_vvl_zstar = .FALSE.              ! zstar  vertical coordinate 
    48    LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.             ! ztilde vertical coordinate 
    49    LOGICAL , PUBLIC                                      :: ln_vvl_layer = .FALSE.              ! level  vertical coordinate 
    50    LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate 
    51    LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor = .FALSE.     ! ztilde vertical coordinate 
    52    LOGICAL , PUBLIC                                      :: ln_vvl_kepe = .FALSE.               ! kinetic/potential energy transfer 
    53    !                                                                                            ! conservation: not used yet 
    54    REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient 
    55    REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days] 
    56    REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days] 
    57    REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation 
    58    LOGICAL , PUBLIC                                      :: ln_vvl_dbg = .FALSE.      ! debug control prints 
    59  
    60    !! * Module variables 
    61    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport 
    62    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence 
    63    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors 
    64    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors 
    65    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors 
    66    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence 
     40 
     41   !                                                      !!* Namelist nam_vvl 
     42   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate 
     43   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate 
     44   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate 
     45   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate 
     46   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate 
     47   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer 
     48   !                                                       ! conservation: not used yet 
     49   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient 
     50   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days] 
     51   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days] 
     52   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation 
     53   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints 
     54 
     55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport 
     56   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence 
     57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors 
     58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors 
     59   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors 
     60   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6761 
    6862   !! * Substitutions 
    … …  
    146140      CALL dom_vvl_rst( nit000, 'READ' ) 
    147141      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
    148  
    149       !IF(ln_wd) THEN 
    150       !  DO jj = 1, jpj 
    151       !    DO ji = 1, jpi 
    152       !      IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
    153       !       fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
    154       !       fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
    155       !       fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
    156       !      END IF 
    157       !    ENDDO 
    158       !  ENDDO 
    159       !END IF 
    160142 
    161143      ! Reconstruction of all vertical scale factors at now and before time steps 
    … …  
    384366            DO jj = 1, jpjm1 
    385367               DO ji = 1, fs_jpim1   ! vector opt. 
    386                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
    387                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    388                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
    389                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     368                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)          & 
     369                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     370                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)          &  
     371                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    390372                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    391373                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    … …  
    406388                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    407389                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    408                      &                                            ) * r1_e12t(ji,jj) 
     390                     &                                            ) * r1_e1e2t(ji,jj) 
    409391               END DO 
    410392            END DO 
    … …  
    707689      !!                - vertical interpolation: simple averaging 
    708690      !!---------------------------------------------------------------------- 
    709       !! * Arguments 
    710       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated 
    711       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3 
    712       CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors 
    713       !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    714       !! * Local declarations 
    715       REAL(wp) ::  zwad                                                ! = 1.0 when ln_wd = .true. 
    716                                                                        ! = 0.0 when ln_wd = .false. 
    717                                                                        !  
    718       INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    719       LOGICAL ::   l_is_orca                                           ! local logical 
    720       !!---------------------------------------------------------------------- 
     691      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated 
     692      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3 
     693      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors 
     694      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
     695      ! 
     696      INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
     697      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd = T/F 
     698      !!---------------------------------------------------------------------- 
     699      ! 
    721700      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol') 
    722          ! 
     701      ! 
    723702      IF(ln_wd) THEN 
    724         zwad = 1.0_wp 
     703        zlnwd = 1.0_wp 
    725704      ELSE 
    726         zwad = 0.0_wp 
     705        zlnwd = 0.0_wp 
    727706      END IF 
    728  
    729       l_is_orca = .FALSE. 
    730       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations 
    731  
    732       SELECT CASE ( pout ) 
    733          !               ! ------------------------------------- ! 
    734       CASE( 'U' )        ! interpolation from T-point to U-point ! 
    735          !               ! ------------------------------------- ! 
    736          ! horizontal surface weighted interpolation 
     707      ! 
     708      SELECT CASE ( pout )    !==  type of interpolation  ==! 
     709         ! 
     710      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    737711         DO jk = 1, jpk 
    738712            DO jj = 1, jpjm1 
    739713               DO ji = 1, fs_jpim1   ! vector opt. 
    740                   !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
    741                   pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12u(ji,jj)        & 
    742                      &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    743                      &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     714                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     715                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     716                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    744717               END DO 
    745718            END DO 
    746719         END DO 
    747          ! 
    748          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    749          ! boundary conditions 
    750720         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 
    751721         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    752          !               ! ------------------------------------- ! 
    753       CASE( 'V' )        ! interpolation from T-point to V-point ! 
    754          !               ! ------------------------------------- ! 
    755          ! horizontal surface weighted interpolation 
     722         ! 
     723      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    756724         DO jk = 1, jpk 
    757725            DO jj = 1, jpjm1 
    758726               DO ji = 1, fs_jpim1   ! vector opt. 
    759                   !pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
    760                   pe3_out(ji,jj,jk) = 0.5_wp * (vmask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12v(ji,jj)        & 
    761                      &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    762                      &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     727                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     728                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     729                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    763730               END DO 
    764731            END DO 
    765732         END DO 
    766          ! 
    767          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    768          ! boundary conditions 
    769733         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 
    770734         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    771          !               ! ------------------------------------- ! 
    772       CASE( 'F' )        ! interpolation from U-point to F-point ! 
    773          !               ! ------------------------------------- ! 
    774          ! horizontal surface weighted interpolation 
     735         ! 
     736      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    775737         DO jk = 1, jpk 
    776738            DO jj = 1, jpjm1 
    777739               DO ji = 1, fs_jpim1   ! vector opt. 
    778                   !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
    779                   pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zwad) + zwad)     & 
    780                      &                       * r1_e12f(ji,jj)                                                     & 
    781                      &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    782                      &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     740                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     741                     &                       *    r1_e1e2f(ji,jj)                                                  & 
     742                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     743                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    783744               END DO 
    784745            END DO 
    785746         END DO 
    786          ! 
    787          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    788          ! boundary conditions 
    789747         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 
    790748         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
    791          !               ! ------------------------------------- ! 
    792       CASE( 'W' )        ! interpolation from T-point to W-point ! 
    793          !               ! ------------------------------------- ! 
    794          ! vertical simple interpolation 
     749         ! 
     750      CASE( 'W' )                   !* from T- to W-point : vertical simple mean 
     751         ! 
    795752         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 
    796          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
     753         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 
     754!!gm BUG? use here wmask in case of ISF ?  to be checked 
    797755         DO jk = 2, jpk 
    798             !pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   & 
    799             !   &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
    800             pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 
    801                &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                         & 
    802                &                            +            0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad)   & 
     756            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   & 
     757               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               & 
     758               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     & 
    803759               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
    804760         END DO 
    805          !               ! -------------------------------------- ! 
    806       CASE( 'UW' )       ! interpolation from U-point to UW-point ! 
    807          !               ! -------------------------------------- ! 
    808          ! vertical simple interpolation 
     761         ! 
     762      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean 
     763         ! 
    809764         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 
    810765         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
     766!!gm BUG? use here wumask in case of ISF ?  to be checked 
    811767         DO jk = 2, jpk 
    812             !pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   & 
    813             !   &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
    814             pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 
    815                &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                         & 
    816                &                             +            0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad)   & 
     768            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
     769               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              & 
     770               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
    817771               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
    818772         END DO 
    819          !               ! -------------------------------------- ! 
    820       CASE( 'VW' )       ! interpolation from V-point to VW-point ! 
    821          !               ! -------------------------------------- ! 
    822          ! vertical simple interpolation 
     773         ! 
     774      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean 
     775         ! 
    823776         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 
    824777         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
     778!!gm BUG? use here wvmask in case of ISF ?  to be checked 
    825779         DO jk = 2, jpk 
    826             !pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   & 
    827             !   &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
    828             pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 
    829                &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                         & 
    830                &                             +            0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad)   & 
     780            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
     781               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              & 
     782               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
    831783               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
    832784         END DO 
    833785      END SELECT 
    834786      ! 
    835  
    836787      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol') 
    837  
     788      ! 
    838789   END SUBROUTINE dom_vvl_interpol 
     790 
    839791 
    840792   SUBROUTINE dom_vvl_rst( kt, cdrw ) 
    … …  
    850802      !!                they are set to 0. 
    851803      !!---------------------------------------------------------------------- 
    852       !! * Arguments 
    853804      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    854805      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    855       !! * Local declarations 
     806      ! 
    856807      INTEGER ::   ji, jj, jk 
    857808      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    … …  
    943894            sshn(:,:) = 0.0_wp 
    944895 
    945             IF(ln_wd) THEN 
     896            IF( ln_wd ) THEN 
    946897              DO jj = 1, jpj 
    947898                DO ji = 1, jpi 
    948                   !IF(e3t_0(ji,jj,1) < 0._wp) THEN 
    949                   !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
    950                   IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
    951                     fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1  
    952                     fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1  
    953                     fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
    954                     sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    955                     sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    956                     ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     899                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 
     900                     fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1  
     901                     fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1  
     902                     fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
     903                     sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     904                     sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     905                     ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    957906                  ENDIF 
    958907                ENDDO 
    … …  
    966915            END IF 
    967916         ENDIF 
    968  
     917         ! 
    969918      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    970919         !                                   ! =================== 
    … …  
    986935            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 
    987936         ENDIF 
    988  
    989       ENDIF 
     937         ! 
     938      ENDIF 
     939      ! 
    990940      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst') 
    991  
     941      ! 
    992942   END SUBROUTINE dom_vvl_rst 
    993943 
    … …  
    1000950      !!                for vertical coordinate 
    1001951      !!---------------------------------------------------------------------- 
    1002       INTEGER ::   ioptio 
    1003       INTEGER ::   ios 
    1004  
     952      INTEGER ::   ioptio, ios 
     953      !! 
    1005954      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 
    1006                       & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
    1007                       & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
     955         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
     956         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
    1008957      !!----------------------------------------------------------------------  
    1009  
     958      ! 
    1010959      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :  
    1011960      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 
    1012961901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 
    1013  
     962      ! 
    1014963      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run 
    1015964      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 
    1016965902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp ) 
    1017966      IF(lwm) WRITE ( numond, nam_vvl ) 
    1018  
     967      ! 
    1019968      IF(lwp) THEN                    ! Namelist print 
    1020969         WRITE(numout,*) 
    … …  
    1049998         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg 
    1050999      ENDIF 
    1051  
     1000      ! 
    10521001      ioptio = 0                      ! Parameter control 
    1053       IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. 
    1054       IF( ln_vvl_zstar           )        ioptio = ioptio + 1 
    1055       IF( ln_vvl_ztilde          )        ioptio = ioptio + 1 
    1056       IF( ln_vvl_layer           )        ioptio = ioptio + 1 
    1057  
     1002      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true. 
     1003      IF( ln_vvl_zstar           )   ioptio = ioptio + 1 
     1004      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1 
     1005      IF( ln_vvl_layer           )   ioptio = ioptio + 1 
     1006      ! 
    10581007      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    10591008      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    1060  
     1009      ! 
    10611010      IF(lwp) THEN                   ! Print the choice 
    10621011         WRITE(numout,*) 
    … …  
    10691018         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used' 
    10701019      ENDIF 
    1071  
     1020      ! 
    10721021#if defined key_agrif 
    10731022      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' ) 
    10741023#endif 
    1075  
     1024      ! 
    10761025   END SUBROUTINE dom_vvl_ctl 
    1077  
    1078    SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    1079       !!--------------------------------------------------------------------- 
    1080       !!                   ***  ROUTINE dom_vvl_orca_fix  *** 
    1081       !!                      
    1082       !! ** Purpose :   Correct surface weighted, horizontally interpolated,  
    1083       !!                scale factors at locations that have been individually 
    1084       !!                modified in domhgr. Such modifications break the 
    1085       !!                relationship between e12t and e1u*e2u etc. 
    1086       !!                Recompute some scale factors ignoring the modified metric. 
    1087       !!---------------------------------------------------------------------- 
    1088       !! * Arguments 
    1089       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated 
    1090       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3 
    1091       CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors 
    1092       !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    1093       !! * Local declarations 
    1094       INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    1095       INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices 
    1096       INTEGER ::   isrow                                               ! index for ORCA1 starting row 
    1097       !! acc 
    1098       !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for 
    1099       !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations 
    1100       !!  
    1101       !                                                ! ===================== 
    1102       IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration 
    1103          !                                             ! ===================== 
    1104       !! acc 
    1105          IF( nn_cla == 0 ) THEN 
    1106             ! 
    1107             ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified) 
    1108             ij0 = 102   ;   ij1 = 102 
    1109             DO jk = 1, jpkm1 
    1110                DO jj = mj0(ij0), mj1(ij1) 
    1111                   DO ji = mi0(ii0), mi1(ii1) 
    1112                      SELECT CASE ( pout ) 
    1113                      CASE( 'U' ) 
    1114                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1115                        &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1116                        &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1117                        &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1118                      CASE( 'F' ) 
    1119                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1120                        &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1121                        &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1122                        &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1123                      END SELECT 
    1124                   END DO 
    1125                END DO 
    1126             END DO 
    1127             ! 
    1128             ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified) 
    1129             ij0 =  88   ;   ij1 =  88 
    1130             DO jk = 1, jpkm1 
    1131                DO jj = mj0(ij0), mj1(ij1) 
    1132                   DO ji = mi0(ii0), mi1(ii1) 
    1133                      SELECT CASE ( pout ) 
    1134                      CASE( 'U' ) 
    1135                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1136                        &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1137                        &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1138                        &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1139                      CASE( 'V' ) 
    1140                         pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1141                        &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1142                        &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1143                        &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1144                      CASE( 'F' ) 
    1145                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1146                        &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1147                        &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1148                        &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1149                      END SELECT 
    1150                   END DO 
    1151                END DO 
    1152             END DO 
    1153          ENDIF 
    1154  
    1155          ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified) 
    1156          ij0 = 116   ;   ij1 = 116 
    1157          DO jk = 1, jpkm1 
    1158             DO jj = mj0(ij0), mj1(ij1) 
    1159                DO ji = mi0(ii0), mi1(ii1) 
    1160                   SELECT CASE ( pout ) 
    1161                   CASE( 'U' ) 
    1162                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1163                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1164                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1165                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1166                   CASE( 'F' ) 
    1167                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1168                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1169                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1170                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1171                   END SELECT 
    1172                END DO 
    1173             END DO 
    1174          END DO 
    1175       ENDIF 
    1176       ! 
    1177          !                                             ! ===================== 
    1178       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    1179          !                                             ! ===================== 
    1180          ! This dirty section will be suppressed by simplification process: 
    1181          ! all this will come back in input files 
    1182          ! Currently these hard-wired indices relate to configuration with 
    1183          ! extend grid (jpjglo=332) 
    1184          ! which had a grid-size of 362x292. 
    1185          isrow = 332 - jpjglo 
    1186          ! 
    1187          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified) 
    1188          ij0 = 241 - isrow   ;   ij1 = 241 - isrow 
    1189          DO jk = 1, jpkm1 
    1190             DO jj = mj0(ij0), mj1(ij1) 
    1191                DO ji = mi0(ii0), mi1(ii1) 
    1192                   SELECT CASE ( pout ) 
    1193                   CASE( 'U' ) 
    1194                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1195                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1196                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1197                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1198                   CASE( 'F' ) 
    1199                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1200                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1201                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1202                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1203                   END SELECT 
    1204                END DO 
    1205             END DO 
    1206          END DO 
    1207          ! 
    1208          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
    1209          ij0 = 248 - isrow   ;   ij1 = 248 - isrow 
    1210          DO jk = 1, jpkm1 
    1211             DO jj = mj0(ij0), mj1(ij1) 
    1212                DO ji = mi0(ii0), mi1(ii1) 
    1213                   SELECT CASE ( pout ) 
    1214                   CASE( 'U' ) 
    1215                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &   
    1216                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1217                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1218                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1219                   CASE( 'F' ) 
    1220                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &   
    1221                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1222                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1223                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1224                   END SELECT 
    1225                END DO 
    1226             END DO 
    1227          END DO 
    1228          ! 
    1229          ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
    1230          ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    1231          DO jk = 1, jpkm1 
    1232             DO jj = mj0(ij0), mj1(ij1) 
    1233                DO ji = mi0(ii0), mi1(ii1) 
    1234                   SELECT CASE ( pout ) 
    1235                   CASE( 'V' ) 
    1236                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1237                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1238                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1239                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1240                   END SELECT 
    1241                END DO 
    1242             END DO 
    1243          END DO 
    1244          ! 
    1245          ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
    1246          ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    1247          DO jk = 1, jpkm1 
    1248             DO jj = mj0(ij0), mj1(ij1) 
    1249                DO ji = mi0(ii0), mi1(ii1) 
    1250                   SELECT CASE ( pout ) 
    1251                   CASE( 'V' ) 
    1252                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1253                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1254                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1255                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1256                   END SELECT 
    1257                END DO 
    1258             END DO 
    1259          END DO 
    1260          ! 
    1261          ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
    1262          ij0 = 164 - isrow  ;   ij1 = 165  - isrow   
    1263          DO jk = 1, jpkm1 
    1264             DO jj = mj0(ij0), mj1(ij1) 
    1265                DO ji = mi0(ii0), mi1(ii1) 
    1266                   SELECT CASE ( pout ) 
    1267                   CASE( 'V' ) 
    1268                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1269                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1270                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1271                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1272                   END SELECT 
    1273                END DO 
    1274             END DO 
    1275          END DO 
    1276          ! 
    1277          ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified) 
    1278          ij0 = 164 - isrow    ;   ij1 = 165  - isrow   
    1279          DO jk = 1, jpkm1 
    1280             DO jj = mj0(ij0), mj1(ij1) 
    1281                DO ji = mi0(ii0), mi1(ii1) 
    1282                   SELECT CASE ( pout ) 
    1283                   CASE( 'V' ) 
    1284                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1285                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1286                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1287                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1288                   END SELECT 
    1289                END DO 
    1290             END DO 
    1291          END DO 
    1292          ! 
    1293          ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
    1294          ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    1295          DO jk = 1, jpkm1 
    1296             DO jj = mj0(ij0), mj1(ij1) 
    1297                DO ji = mi0(ii0), mi1(ii1) 
    1298                   SELECT CASE ( pout ) 
    1299                   CASE( 'V' ) 
    1300                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1301                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1302                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1303                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1304                   END SELECT 
    1305                END DO 
    1306             END DO 
    1307          END DO 
    1308          ! 
    1309          ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
    1310          ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    1311          DO jk = 1, jpkm1 
    1312             DO jj = mj0(ij0), mj1(ij1) 
    1313                DO ji = mi0(ii0), mi1(ii1) 
    1314                   SELECT CASE ( pout ) 
    1315                   CASE( 'V' ) 
    1316                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1317                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1318                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1319                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1320                   END SELECT 
    1321                END DO 
    1322             END DO 
    1323          END DO 
    1324       ENDIF 
    1325          !                                             ! ===================== 
    1326       IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration 
    1327          !                                             ! ===================== 
    1328          ! 
    1329          ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified) 
    1330          ij0 = 327   ;   ij1 = 327 
    1331          DO jk = 1, jpkm1 
    1332             DO jj = mj0(ij0), mj1(ij1) 
    1333                DO ji = mi0(ii0), mi1(ii1) 
    1334                   SELECT CASE ( pout ) 
    1335                   CASE( 'U' ) 
    1336                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1337                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1338                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1339                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1340                   CASE( 'F' ) 
    1341                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1342                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1343                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1344                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1345                   END SELECT 
    1346                END DO 
    1347             END DO 
    1348          END DO 
    1349          ! 
    1350          ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified) 
    1351          ij0 = 343   ;   ij1 = 343 
    1352          DO jk = 1, jpkm1 
    1353             DO jj = mj0(ij0), mj1(ij1) 
    1354                DO ji = mi0(ii0), mi1(ii1) 
    1355                   SELECT CASE ( pout ) 
    1356                   CASE( 'U' ) 
    1357                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &   
    1358                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1359                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1360                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1361                   CASE( 'F' ) 
    1362                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &   
    1363                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1364                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1365                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1366                   END SELECT 
    1367                END DO 
    1368             END DO 
    1369          END DO 
    1370          ! 
    1371          ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified) 
    1372          ij0 = 232   ;   ij1 = 232 
    1373          DO jk = 1, jpkm1 
    1374             DO jj = mj0(ij0), mj1(ij1) 
    1375                DO ji = mi0(ii0), mi1(ii1) 
    1376                   SELECT CASE ( pout ) 
    1377                   CASE( 'U' ) 
    1378                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1379                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1380                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1381                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1382                   CASE( 'F' ) 
    1383                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1384                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1385                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1386                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1387                   END SELECT 
    1388                END DO 
    1389             END DO 
    1390          END DO 
    1391          ! 
    1392          ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified) 
    1393          ij0 = 232   ;   ij1 = 232 
    1394          DO jk = 1, jpkm1 
    1395             DO jj = mj0(ij0), mj1(ij1) 
    1396                DO ji = mi0(ii0), mi1(ii1) 
    1397                   SELECT CASE ( pout ) 
    1398                   CASE( 'U' ) 
    1399                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1400                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1401                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1402                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1403                   CASE( 'F' ) 
    1404                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1405                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1406                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1407                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1408                   END SELECT 
    1409                END DO 
    1410             END DO 
    1411          END DO 
    1412          ! 
    1413          ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified) 
    1414          ij0 = 270   ;   ij1 = 270 
    1415          DO jk = 1, jpkm1 
    1416             DO jj = mj0(ij0), mj1(ij1) 
    1417                DO ji = mi0(ii0), mi1(ii1) 
    1418                   SELECT CASE ( pout ) 
    1419                   CASE( 'U' ) 
    1420                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1421                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1422                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1423                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1424                   CASE( 'F' ) 
    1425                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1426                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1427                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1428                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1429                   END SELECT 
    1430                END DO 
    1431             END DO 
    1432          END DO 
    1433          ! 
    1434          ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified) 
    1435          ij0 = 232   ;   ij1 = 233 
    1436          DO jk = 1, jpkm1 
    1437             DO jj = mj0(ij0), mj1(ij1) 
    1438                DO ji = mi0(ii0), mi1(ii1) 
    1439                   SELECT CASE ( pout ) 
    1440                   CASE( 'V' ) 
    1441                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1442                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1443                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1444                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1445                   END SELECT 
    1446                END DO 
    1447             END DO 
    1448          END DO 
    1449          ! 
    1450          ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified) 
    1451          ij0 = 276   ;   ij1 = 276 
    1452          DO jk = 1, jpkm1 
    1453             DO jj = mj0(ij0), mj1(ij1) 
    1454                DO ji = mi0(ii0), mi1(ii1) 
    1455                   SELECT CASE ( pout ) 
    1456                   CASE( 'V' ) 
    1457                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1458                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1459                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1460                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1461                   END SELECT 
    1462                END DO 
    1463             END DO 
    1464          END DO 
    1465       ENDIF 
    1466    END SUBROUTINE dom_vvl_orca_fix 
    14671026 
    14681027   !!====================================================================== 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r5603 r5870  
    77   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL 
    88   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file 
    9    !!            3.0  ! 2008-01  (S. Masson) add dom_uniq  
    10    !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
     9   !!            3.0  ! 2008-01  (S. Masson)  add dom_uniq  
    1110   !!---------------------------------------------------------------------- 
    1211 
    1312   !!---------------------------------------------------------------------- 
    1413   !!   dom_wri        : create and write mesh and mask file(s) 
    15    !!   dom_uniq       : 
     14   !!   dom_uniq       : identify unique point of a grid (TUVF) 
    1615   !!---------------------------------------------------------------------- 
    1716   USE dom_oce         ! ocean space and time domain 
    … …  
    2625   PRIVATE 
    2726 
    28    PUBLIC dom_wri        ! routine called by inidom.F90 
    29  
     27   PUBLIC   dom_wri              ! routine called by inidom.F90 
     28   PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90 
    3029   !! * Substitutions 
    3130#  include "vectopt_loop_substitute.h90" 
    … …  
    3635   !!---------------------------------------------------------------------- 
    3736CONTAINS 
     37 
     38   SUBROUTINE dom_wri_coordinate 
     39      !!---------------------------------------------------------------------- 
     40      !!                  ***  ROUTINE dom_wri_coordinate  *** 
     41      !!                    
     42      !! ** Purpose :   Create the NetCDF file which contains all the 
     43      !!              standard coordinate information plus the surface, 
     44      !!              e1e2u and e1e2v. By doing so, those surface will 
     45      !!              not be changed by the reduction of e1u or e2v scale  
     46      !!              factors in some straits.  
     47      !!                 NB: call just after the read of standard coordinate 
     48      !!              and the reduction of scale factors in some straits 
     49      !! 
     50      !! ** output file :   coordinate_e1e2u_v.nc 
     51      !!---------------------------------------------------------------------- 
     52      INTEGER           ::   inum0    ! temprary units for 'coordinate_e1e2u_v.nc' file 
     53      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations) 
     54      !                                   !  workspaces 
     55      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw  
     56      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv 
     57      !!---------------------------------------------------------------------- 
     58      ! 
     59      IF( nn_timing == 1 )  CALL timing_start('dom_wri_coordinate') 
     60      ! 
     61      IF(lwp) WRITE(numout,*) 
     62      IF(lwp) WRITE(numout,*) 'dom_wri_coordinate : create NetCDF coordinate file' 
     63      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 
     64       
     65      clnam0 = 'coordinate_e1e2u_v'  ! filename (mesh and mask informations) 
     66       
     67      !  create 'coordinate_e1e2u_v.nc' file 
     68      ! ============================ 
     69      ! 
     70      CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
     71      ! 
     72      !                                                         ! horizontal mesh (inum3) 
     73      CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude 
     74      CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r4 ) 
     75      CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r4 ) 
     76      CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r4 ) 
     77       
     78      CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude 
     79      CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r4 ) 
     80      CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r4 ) 
     81      CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r4 ) 
     82       
     83      CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors 
     84      CALL iom_rstput( 0, 0, inum0, 'e1u', e1u, ktype = jp_r8 ) 
     85      CALL iom_rstput( 0, 0, inum0, 'e1v', e1v, ktype = jp_r8 ) 
     86      CALL iom_rstput( 0, 0, inum0, 'e1f', e1f, ktype = jp_r8 ) 
     87       
     88      CALL iom_rstput( 0, 0, inum0, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors 
     89      CALL iom_rstput( 0, 0, inum0, 'e2u', e2u, ktype = jp_r8 ) 
     90      CALL iom_rstput( 0, 0, inum0, 'e2v', e2v, ktype = jp_r8 ) 
     91      CALL iom_rstput( 0, 0, inum0, 'e2f', e2f, ktype = jp_r8 ) 
     92       
     93      CALL iom_rstput( 0, 0, inum0, 'e1e2u', e1e2u, ktype = jp_r8 ) 
     94      CALL iom_rstput( 0, 0, inum0, 'e1e2v', e1e2v, ktype = jp_r8 ) 
     95 
     96      CALL iom_close( inum0 ) 
     97      ! 
     98      IF( nn_timing == 1 )  CALL timing_stop('dom_wri_coordinate') 
     99      ! 
     100   END SUBROUTINE dom_wri_coordinate 
     101 
    38102 
    39103   SUBROUTINE dom_wri 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5842 r5870  
    3737   USE oce               ! ocean variables 
    3838   USE dom_oce           ! ocean domain 
     39   USE wet_dry           ! wetting and drying 
    3940   USE closea            ! closed seas 
    4041   USE c1d               ! 1D vertical configuration 
    … …  
    502503            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    503504               !                                             ! ===================== 
    504                IF( nn_cla == 0 ) THEN 
    505                   ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
    506                   ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995) 
    507                   DO ji = mi0(ii0), mi1(ii1) 
    508                      DO jj = mj0(ij0), mj1(ij1) 
    509                         mbathy(ji,jj) = 15 
    510                      END DO 
     505               ! 
     506               ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
     507               ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995) 
     508               DO ji = mi0(ii0), mi1(ii1) 
     509                  DO jj = mj0(ij0), mj1(ij1) 
     510                     mbathy(ji,jj) = 15 
    511511                  END DO 
    512                   IF(lwp) WRITE(numout,*) 
    513                   IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
    514                   ! 
    515                   ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open 
    516                   ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995) 
    517                   DO ji = mi0(ii0), mi1(ii1) 
    518                      DO jj = mj0(ij0), mj1(ij1) 
    519                         mbathy(ji,jj) = 12 
    520                      END DO 
     512               END DO 
     513               IF(lwp) WRITE(numout,*) 
     514               IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
     515               ! 
     516               ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open 
     517               ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995) 
     518               DO ji = mi0(ii0), mi1(ii1) 
     519                  DO jj = mj0(ij0), mj1(ij1) 
     520                     mbathy(ji,jj) = 12 
    521521                  END DO 
    522                   IF(lwp) WRITE(numout,*) 
    523                   IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    524                ENDIF 
     522               END DO 
     523               IF(lwp) WRITE(numout,*) 
     524               IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    525525               ! 
    526526            ENDIF 
    … …  
    546546            !        
    547547            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    548                ! 
    549               IF( nn_cla == 0 ) THEN 
    550                  ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open  
    551                  ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995) 
    552                  DO ji = mi0(ii0), mi1(ii1) 
    553                     DO jj = mj0(ij0), mj1(ij1) 
    554                        bathy(ji,jj) = 284._wp 
    555                     END DO 
     548            ! 
     549              ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open  
     550              ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995) 
     551              DO ji = mi0(ii0), mi1(ii1) 
     552                 DO jj = mj0(ij0), mj1(ij1) 
     553                    bathy(ji,jj) = 284._wp 
    556554                 END DO 
    557                  IF(lwp) WRITE(numout,*)      
    558                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
    559                  ! 
    560                  ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open 
    561                  ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995) 
    562                  DO ji = mi0(ii0), mi1(ii1) 
    563                     DO jj = mj0(ij0), mj1(ij1) 
    564                        bathy(ji,jj) = 137._wp 
    565                     END DO 
     555               END DO 
     556              IF(lwp) WRITE(numout,*)      
     557              IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
     558              ! 
     559              ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open 
     560              ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995) 
     561               DO ji = mi0(ii0), mi1(ii1) 
     562                 DO jj = mj0(ij0), mj1(ij1) 
     563                    bathy(ji,jj) = 137._wp 
    566564                 END DO 
    567                  IF(lwp) WRITE(numout,*) 
    568                  IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    569               ENDIF 
     565              END DO 
     566              IF(lwp) WRITE(numout,*) 
     567               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    570568              ! 
    571569           ENDIF 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r5215 r5870  
    174174            END DO 
    175175         END DO 
    176          IF( nn_cla == 1 ) THEN                          ! Cross Land advection 
    177             il0 = 138   ;   il1 = 138                          ! set T & S profile at Gibraltar Strait 
    178             ij0 = 101   ;   ij1 = 102 
    179             ii0 = 139   ;   ii1 = 139 
    180             DO jl = mi0(il0), mi1(il1) 
    181                DO jj = mj0(ij0), mj1(ij1) 
    182                   DO ji = mi0(ii0), mi1(ii1) 
    183                      sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 
    184                      sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 
    185                   END DO 
    186                END DO 
    187             END DO 
    188             il0 = 164   ;   il1 = 164                          ! set T & S profile at Bab el Mandeb Strait 
    189             ij0 =  87   ;   ij1 =  88 
    190             ii0 = 161   ;   ii1 = 163 
    191             DO jl = mi0(il0), mi1(il1) 
    192                DO jj = mj0(ij0), mj1(ij1) 
    193                   DO ji = mi0(ii0), mi1(ii1) 
    194                      sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 
    195                      sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 
    196                   END DO 
    197                END DO 
    198             END DO 
    199          ELSE                                            ! No Cross Land advection 
    200             ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea 
    201             ii0 = 148   ;   ii1 = 160 
    202             sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp 
    203             sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 
    204             sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 
    205          ENDIF 
     176         ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea 
     177         ii0 = 148   ;   ii1 = 160 
     178         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp 
     179         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 
     180         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 
    206181      ENDIF 
    207182      ! 

  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5332 r5870  
    2929   USE daymod          ! calendar 
    3030   USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine) 
    31    USE ldftra_oce      ! ocean active tracers: lateral physics 
     31   USE ldftra          ! lateral physics: ocean active tracers 
    3232   USE zdf_oce         ! ocean vertical physics 
    3333   USE phycst          ! physical constants 
    3434   USE dtatsd          ! data temperature and salinity   (dta_tsd routine) 
    3535   USE dtauvd          ! data: U & V current             (dta_uvd routine) 
    36    USE zpshde          ! partial step: hor. derivative (zps_hde routine) 
    37    USE eosbn2          ! equation of state            (eos bn2 routine) 
    3836   USE domvvl          ! varying vertical mesh 
    3937   USE dynspg_oce      ! pressure gradient schemes 
    … …  
    7674      ! 
    7775 
    78       IF(lwp) WRITE(numout,*) ' ' 
     76      IF(lwp) WRITE(numout,*) 
    7977      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 
    8078      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    8179 
    82       CALL dta_tsd_init                       ! Initialisation of T & S input data 
    83       IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data 
     80                     CALL dta_tsd_init        ! Initialisation of T & S input data 
     81      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data 
    8482 
    8583      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    … …  
    103101         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp 
    104102         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp   
    105          rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp 
    106          hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp 
     103                                    hdivn(:,:,:) = 0._wp 
    107104         ! 
    108105         IF( cp_cfg == 'eel' ) THEN 
    … …  
    119116            ENDIF 
    120117            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 
    121                CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd ) 
     118               CALL wrk_alloc( jpi,jpj,jpk,2,  zuvd ) 
    122119               CALL dta_uvd( nit000, zuvd ) 
    123120               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:) 
    124121               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:) 
    125                CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd ) 
     122               CALL wrk_dealloc( jpi,jpj,jpk,2,  zuvd ) 
    126123            ENDIF 
    127124         ENDIF 
    128          ! 
    129          CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    130 #if ! defined key_c1d 
    131          IF( ln_zps .AND. .NOT. ln_isfcav)                                 & 
    132             &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    133             &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    134          IF( ln_zps .AND.       ln_isfcav)                                 & 
    135             &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    136             &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    137             &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    138 #endif 
    139125         !    
    140126         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 
    … …  
    142128            DO jk = 1, jpk 
    143129               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
    144             ENDDO 
     130            END DO 
    145131         ENDIF 
    146132         !  
    … …  
    163149      ! 
    164150      ! 
    165       un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 
    166       ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 
     151      un_b(:,:) = 0._wp   ;  vn_b(:,:) = 0._wp 
     152      ub_b(:,:) = 0._wp   ;  vb_b(:,:) = 0._wp 
    167153      ! 
    168154      DO jk = 1, jpkm1 
    … …  
    202188      !! References :  Philander ??? 
    203189      !!---------------------------------------------------------------------- 
    204       INTEGER  :: ji, jj, jk 
    205       REAL(wp) ::   zsal = 35.50 
     190      INTEGER  ::   ji, jj, jk 
     191      REAL(wp) ::   zsal = 35.50_wp 
    206192      !!---------------------------------------------------------------------- 
    207193      ! 
    … …  
    233219      !!                and relative vorticity fields 
    234220      !!---------------------------------------------------------------------- 
    235       USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine) 
     221      USE divhor     ! hor. divergence      (div_hor routine) 
    236222      USE iom 
    237223      ! 
    … …  
    282268            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 
    283269            ! 
    284             ! set the dynamics: U,V, hdiv, rot (and ssh if necessary) 
     270            ! set the dynamics: U,V, hdiv (and ssh if necessary) 
    285271            ! ---------------- 
    286272            ! Start EEL5 configuration with barotropic geostrophic velocities  
    … …  
    318304            ENDIF 
    319305            ! 
    320             CALL div_cur( nit000 )                  ! horizontal divergence and relative vorticity (curl) 
     306!!gm  Check  here call to div_hor should not be necessary 
     307!!gm         div_hor call runoffs  not sure they are defined at that level 
     308            CALL div_hor( nit000 )                  ! horizontal divergence and relative vorticity (curl) 
    321309            ! N.B. the vertical velocity will be computed from the horizontal divergence field 
    322310            ! in istate by a call to wzv routine 
    … …  
    371359      !! 
    372360      !! ** Method  : - set temprature field 
    373       !!              - set salinity field 
     361      !!              - set salinity   field 
    374362      !!---------------------------------------------------------------------- 
    375363      INTEGER :: ji, jj, jk  ! dummy loop indices 
    … …  
    458446      !!---------------------------------------------------------------------- 
    459447      USE dynspg          ! surface pressure gradient             (dyn_spg routine) 
    460       USE divcur          ! hor. divergence & rel. vorticity      (div_cur routine) 
     448      USE divhor          ! hor. divergence                       (div_hor routine) 
    461449      USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    462450      ! 
    … …  
    467455      !!---------------------------------------------------------------------- 
    468456      ! 
    469       CALL wrk_alloc( jpi, jpj, jpk, zprn) 
     457      CALL wrk_alloc( jpi,jpj,jpk,  zprn) 
    470458      ! 
    471459      IF(lwp) WRITE(numout,*)  
    … …  
    544532      indic = 0 
    545533      CALL dyn_spg( nit000, indic )       ! surface pressure gradient 
    546  
     534      ! 
    547535      ! the new velocity is ua*rdt 
    548  
     536      ! 
    549537      CALL lbc_lnk( ua, 'U', -1. ) 
    550538      CALL lbc_lnk( va, 'V', -1. ) 
    … …  
    556544      un(:,:,:) = ub(:,:,:) 
    557545      vn(:,:,:) = vb(:,:,:) 
    558         
    559       ! Compute the divergence and curl 
    560  
    561       CALL div_cur( nit000 )            ! now horizontal divergence and curl 
    562  
    563       hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value 
    564       rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value 
    565       ! 
    566       CALL wrk_dealloc( jpi, jpj, jpk, zprn) 
     546      ! 
     547!!gm  Check  here call to div_hor should not be necessary 
     548!!gm         div_hor call runoffs  not sure they are defined at that level 
     549      CALL div_hor( nit000 )            ! now horizontal divergence 
     550      ! 
     551      CALL wrk_dealloc( jpi,jpj,jpk,   zprn) 
    567552      ! 
    568553   END SUBROUTINE istate_uvg 
Note: See TracChangeset for help on using the changeset viewer.