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 5737 for branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Ignore:
Timestamp:
2015-09-13T09:42:41+02:00 (9 years ago)
Author:
gm
Message:

#1593: LDF-ADV, step I: Phasing of horizontal scale factors correct 2

Location:
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5123 r5737  
    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.7  ! 2015-11  (G. Madec) introduce surface and scale factor ratio 
    1213   !!---------------------------------------------------------------------- 
    1314 
     
    158159   !! horizontal curvilinear coordinate and scale factors 
    159160   !! --------------------------------------------------------------------- 
    160    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  glamt, glamu   !: longitude of t-, u-, v- and f-points (degre) 
    161    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  glamv, glamf   !: 
    162    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphit, gphiu   !: latitude  of t-, u-, v- and f-points (degre) 
    163    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphiv, gphif   !: 
    164    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1t, e2t, r1_e1t, r1_e2t   !: horizontal scale factors and inverse at t-point (m) 
    165    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1u, e2u, r1_e1u, r1_e2u   !: horizontal scale factors and inverse at u-point (m) 
    166    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1v, e2v, r1_e1v, r1_e2v   !: horizontal scale factors and inverse at v-point (m) 
    167    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1f, e2f, r1_e1f, r1_e2f   !: horizontal scale factors and inverse at f-point (m) 
    168    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  e1e2t          !: surface at t-point (m2) 
    169    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 
     161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   glamt , glamu, glamv , glamf    !: longitude at t, u, v, f-points [degree] 
     162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gphit , gphiu, gphiv , gphif    !: latitude  at t, u, v, f-points [degree] 
     163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1t   , e2t  , r1_e1t, r1_e2t   !: t-point horizontal scale factors    [m] 
     164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1u   , e2u  , r1_e1u, r1_e2u   !: horizontal scale factors at u-point [m] 
     165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1v   , e2v  , r1_e1v, r1_e2v   !: horizontal scale factors at v-point [m] 
     166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1f   , e2f  , r1_e1f, r1_e2f   !: horizontal scale factors at f-point [m] 
     167   ! 
     168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2t , r1_e1e2t                !: associated metrics at t-point 
     169   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2u , r1_e1e2u , e2_e1u       !: associated metrics at u-point 
     170   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2v , r1_e1e2v , e1_e2v       !: associated metrics at v-point 
     171   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point 
     172   ! 
     173   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff                              !: coriolis factor                   [1/s] 
    170174 
    171175   !!---------------------------------------------------------------------- 
     
    216220   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0           !: reference depth at t-       points (meters) 
    217221   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters) 
    218    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re2u_e1u       !: scale factor coeffs at u points (e2u/e1u) 
    219    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re1v_e2v       !: scale factor coeffs at v points (e1v/e2v) 
    220    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12t , r1_e12t !: horizontal cell surface and inverse at t points 
    221    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12u , r1_e12u !: horizontal cell surface and inverse at u points 
    222    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12v , r1_e12v !: horizontal cell surface and inverse at v points 
    223    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12f , r1_e12f !: horizontal cell surface and inverse at f points 
    224222 
    225223   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
     
    333331   INTEGER FUNCTION dom_oce_alloc() 
    334332      !!---------------------------------------------------------------------- 
    335       INTEGER, DIMENSION(12) :: ierr 
     333      INTEGER, DIMENSION(13) :: ierr 
    336334      !!---------------------------------------------------------------------- 
    337335      ierr(:) = 0 
     
    346344         &      tpol(jpiglo)  , fpol(jpiglo)                               , STAT=ierr(2) ) 
    347345         ! 
    348       ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) ,   &  
    349          &      glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) ,   &   
    350          &      glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) ,   &   
    351          &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) ,   & 
    352          &      e1e2t(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )      
     346      ALLOCATE( glamt(jpi,jpj) ,    glamu(jpi,jpj) ,  glamv(jpi,jpj) ,  glamf(jpi,jpj) ,     & 
     347         &      gphit(jpi,jpj) ,    gphiu(jpi,jpj) ,  gphiv(jpi,jpj) ,  gphif(jpi,jpj) ,     & 
     348         &       e1t (jpi,jpj) ,     e2t (jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) ,     & 
     349         &       e1u (jpi,jpj) ,     e2u (jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) ,     & 
     350         &       e1v (jpi,jpj) ,     e2v (jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) ,     & 
     351         &       e1f (jpi,jpj) ,     e2f (jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) ,     & 
     352         &      e1e2t(jpi,jpj) , r1_e1e2t(jpi,jpj)                                     ,     & 
     353         &      e1e2u(jpi,jpj) , r1_e1e2u(jpi,jpj) , e2_e1u(jpi,jpj)                   ,     & 
     354         &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     & 
     355         &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     & 
     356         &        ff (jpi,jpj)                                                         , STAT=ierr(3) ) 
    353357         ! 
    354358      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         & 
     
    364368         &      gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) ,                           & 
    365369         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) ,                           & 
    366          &      ehu_a    (jpi,jpj)    , ehv_a  (jpi,jpj),                                                     & 
    367          &      ehur_a   (jpi,jpj)    , ehvr_a (jpi,jpj),                                                     & 
    368          &      ehu_b    (jpi,jpj)    , ehv_b  (jpi,jpj),                                                     & 
    369          &      ehur_b   (jpi,jpj)    , ehvr_b (jpi,jpj),                                  STAT=ierr(5) )                           
    370 #endif 
    371          ! 
    372       ALLOCATE( hu      (jpi,jpj) , hur     (jpi,jpj) , hu_0(jpi,jpj) , ht_0  (jpi,jpj) ,     & 
    373          &      hv      (jpi,jpj) , hvr     (jpi,jpj) , hv_0(jpi,jpj) , ht    (jpi,jpj) ,     & 
    374          &      re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) ,                                       & 
    375          &      e12t    (jpi,jpj) , r1_e12t (jpi,jpj) ,                                       & 
    376          &      e12u    (jpi,jpj) , r1_e12u (jpi,jpj) ,                                       & 
    377          &      e12v    (jpi,jpj) , r1_e12v (jpi,jpj) ,                                       & 
    378          &      e12f    (jpi,jpj) , r1_e12f (jpi,jpj) ,                                   STAT=ierr(6)  ) 
     370         &      ehu_a   (jpi,jpj)     , ehv_a (jpi,jpj),                                                     & 
     371         &      ehur_a  (jpi,jpj)     , ehvr_a(jpi,jpj),                                                     & 
     372         &      ehu_b   (jpi,jpj)     , ehv_b (jpi,jpj),                                                     & 
     373         &      ehur_b  (jpi,jpj)     , ehvr_b(jpi,jpj),                                  STAT=ierr(5) )                           
     374#endif 
     375         ! 
     376      ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , ht_0(jpi,jpj) ,     & 
     377         &      hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , ht  (jpi,jpj) , STAT=ierr(6)  ) 
    379378         ! 
    380379      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) ,                                     & 
     
    387386         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     & 
    388387         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     & 
    389          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
     388         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 
    390389 
    391390      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      & 
    392391         &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 
    393          &     bmask(jpi,jpj)  ,                                                       & 
     392         &     bmask  (jpi,jpj) ,                                                       & 
    394393         &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
    395394 
    396395! (ISF) Allocation of basic array    
    397       ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),     & 
    398          &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
    399          &     mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
     396      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),                   & 
     397         &      mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
     398         &      mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
    400399 
    401400      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  
     
    405404 
    406405#if defined key_noslip_accurate 
    407       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
     406      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(13) ) 
    408407#endif 
    409408      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5656 r5737  
    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) 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      ie1e2u_v = 0               !  set to unread e1e2u and e1e2v 
     129      ! 
     130      SELECT CASE( jphgr_msh )   !  type of horizontal mesh   
     131      ! 
     132      CASE ( 0 )                     !==  read in coordinate.nc file  ==! 
     133         ! 
    130134         IF(lwp) WRITE(numout,*) 
    131135         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file' 
    132  
    133          CALL hgr_read           ! Defaultl option  :   NetCDF file 
    134  
     136         ! 
     137         CALL hgr_read( ie1e2u_v ) 
     138         ! 
     139         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them 
     140            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction 
     141            e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
     142            e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
     143 
     144         ! 
    135145         !                                                ! ===================== 
    136146         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
     
    157167            ! 
    158168         ENDIF 
    159  
    160             !                                             ! ===================== 
     169         ! 
     170         !                                                ! ===================== 
    161171         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    162172            !                                             ! ===================== 
    163173            ! This dirty section will be suppressed by simplification process: all this will come back in input files 
    164174            ! Currently these hard-wired indices relate to configuration with 
    165             ! extend grid (jpjglo=332) 
    166             ! which had a grid-size of 362x292. 
     175            ! extend grid (jpjglo=332)  which had a grid-size of 362x292. 
    167176            !  
    168177            isrow = 332 - jpjglo 
     
    208217            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km' 
    209218            ! 
    210             ! 
    211          ENDIF 
    212  
     219         ENDIF 
     220         ! 
    213221         !                                                ! ====================== 
    214222         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration 
     
    251259            ! 
    252260         ENDIF 
    253  
    254  
     261          
     262            !                       ! create 'coordinate_e1e2u_v.nc' file  that contains 
     263            !                       ! reduced scale factor in some strait but full e1e2u and e1e2v surfaces          
     264            IF( ie1e2u_v == 0 )   CALL dom_wri_coordinate 
     265            ! 
     266            ! 
     267         ENDIF 
     268 
     269 
     270         ! 
    255271         ! N.B. :  General case, lat and long function of both i and j indices: 
    256272         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   & 
     
    271287         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   & 
    272288         !                                  + (                           fsdjph( zfi, zfj ) )**2  ) 
    273  
    274  
    275       CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing 
    276  
     289         ! 
     290         ! 
     291      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==! 
     292         ! 
    277293         IF(lwp) WRITE(numout,*) 
    278294         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing' 
    279295         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg'  
    280  
     296         ! 
    281297         DO jj = 1, jpj 
    282298            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 
     299               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 ) 
     300               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 ) 
     301               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
     302               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
    287303         ! Longitude 
    288304               glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
     
    307323            END DO 
    308324         END DO 
    309  
    310  
    311       CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing 
    312  
     325         ! 
     326      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==! 
     327         ! 
    313328         IF(lwp) WRITE(numout,*) 
    314329         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing' 
    315330         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'  
    316  
     331         ! 
    317332         ! Position coordinates (in kilometers) 
    318333         !                          ========== 
    319334         glam0 = 0.e0 
    320335         gphi0 = - ppe2_m * 1.e-3 
    321           
     336         ! 
    322337#if defined key_agrif  
    323338         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
     
    332347         DO jj = 1, jpj 
    333348            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 ) 
     349               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       ) 
     350               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 ) 
    336351               glamv(ji,jj) = glamt(ji,jj) 
    337352               glamf(ji,jj) = glamu(ji,jj) 
    338     
    339                gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       ) 
     353               ! 
     354               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       ) 
    340355               gphiu(ji,jj) = gphit(ji,jj) 
    341                gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 ) 
     356               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 
    342357               gphif(ji,jj) = gphiv(ji,jj) 
    343358            END DO 
    344359         END DO 
    345  
     360         ! 
    346361         ! Horizontal scale factors (in meters) 
    347362         !                              ====== 
     
    350365         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m 
    351366         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m 
    352  
    353       CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type 
    354  
     367         ! 
     368      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==! 
     369         ! 
    355370         IF(lwp) WRITE(numout,*) 
    356371         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type' 
    357372         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg' 
    358373         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 
    359  
     374         ! 
    360375         !  Find index corresponding to the equator, given the grid spacing e1_deg 
    361376         !  and the (approximate) southern latitude ppgphi0. 
     
    365380         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 
    366381         IF(  ppgphi0 > 0 )  ijeq = -ijeq 
    367  
     382         ! 
    368383         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq 
    369  
     384         ! 
    370385         DO jj = 1, jpj 
    371386            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 
     387               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 ) 
     388               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 ) 
     389               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
     390               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
    376391         ! Longitude 
    377392               glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
     
    396411            END DO 
    397412         END DO 
    398  
    399       CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration) 
    400  
     413         ! 
     414      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 
     415         ! 
    401416         IF(lwp) WRITE(numout,*) 
    402417         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 
    403418         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
    404  
     419         ! 
    405420         ! Position coordinates (in kilometers) 
    406421         !                          ========== 
    407  
     422         ! 
    408423         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 
    409          zlam1 = -85 
    410          zphi1 = 29 
     424         zlam1 = -85._wp 
     425         zphi1 =  29._wp 
    411426         ! resolution in meters 
    412          ze1 = 106000. / FLOAT(jp_cfg)             
     427         ze1 = 106000. / REAL( jp_cfg , wp )             
    413428         ! benchmark: forced the resolution to be about 100 km 
    414429         IF( nbench /= 0 )   ze1 = 106000.e0      
    415          zsin_alpha = - SQRT( 2. ) / 2. 
    416          zcos_alpha =   SQRT( 2. ) / 2. 
     430         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 
     431         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp 
    417432         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  
     433         IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon 
     434         !                                                           ! at the right jp_cfg resolution 
     435         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     436         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     437         ! 
    423438         IF( nprint==1 .AND. lwp )   THEN 
    424439            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
    425440            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 
    426441         ENDIF 
    427  
     442         ! 
    428443         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  
     444            DO ji = 1, jpi 
     445               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5 
     446               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5 
     447               ! 
     448               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     449               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     450               ! 
     451               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     452               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     453               ! 
     454               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     455               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     456               ! 
     457               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     458               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     459            END DO 
     460         END DO 
     461         ! 
    447462         ! Horizontal scale factors (in meters) 
    448463         !                              ====== 
     
    451466         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1 
    452467         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1 
    453  
     468         ! 
    454469      CASE DEFAULT 
    455470         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh 
    456471         CALL ctl_stop( ctmp1 ) 
    457  
     472         ! 
    458473      END SELECT 
    459474       
    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 
     475      ! associated horizontal metrics 
     476      ! ----------------------------- 
     477      ! 
     478      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:) 
     479      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:) 
     480      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:) 
     481      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:) 
     482      ! 
     483      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 
     484      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 
     485      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them 
     486         e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
     487         e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
     488      ENDIF 
     489      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases 
     490      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) 
     491      !    
     492      e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 
     493      e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 
     494 
     495      IF( lwp .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart) 
    489496         WRITE(numout,*) 
    490497         WRITE(numout,*) '          longitude and e1 scale factors' 
     
    4965039300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
    497504            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 
    498           
     505            ! 
    499506         WRITE(numout,*) 
    500507         WRITE(numout,*) '          latitude and e2 scale factors' 
     
    506513      ENDIF 
    507514 
    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  
    521515 
    522516      ! ================= ! 
     
    525519 
    526520      SELECT CASE( jphgr_msh )   ! type of horizontal mesh 
    527  
     521      ! 
    528522      CASE ( 0, 1, 4 )               ! mesh on the sphere 
    529  
     523         ! 
    530524         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
    531  
     525         ! 
    532526      CASE ( 2 )                     ! f-plane at ppgphi0  
    533  
     527         ! 
    534528         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
    535  
     529         ! 
    536530         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1) 
    537  
     531         ! 
    538532      CASE ( 3 )                     ! beta-plane 
    539  
     533         ! 
    540534         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           
     535         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points 
     536         ! 
    543537#if defined key_agrif 
    544538         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
    545539            IF( .NOT. Agrif_Root() ) THEN 
    546               zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m)   &  
    547                     &           / (ra * rad) 
     540              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 
    548541            ENDIF 
    549542         ENDIF 
    550543#endif          
    551544         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south 
    552  
     545         ! 
    553546         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
    554           
     547         ! 
    555548         IF(lwp) THEN 
    556549            WRITE(numout,*)  
     
    565558            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    566559         END IF 
    567  
     560         ! 
    568561      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration) 
    569  
     562         ! 
    570563         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0 
    571564         zphi0 = 15.e0                                                      ! latitude of the first row F-points 
    572565         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south 
    573  
     566         ! 
    574567         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
    575  
     568         ! 
    576569         IF(lwp) THEN 
    577570            WRITE(numout,*)  
     
    579572            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 
    580573         ENDIF 
    581  
     574         ! 
    582575         IF( lk_mpp ) THEN  
    583576            zminff=ff(nldi,nldj) 
     
    587580            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    588581         END IF 
    589  
     582         ! 
    590583      END SELECT 
    591584 
     
    596589 
    597590      IF( nperio == 2 ) THEN 
    598          znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi ) 
     591         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
    599592         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 
    600593      ENDIF 
     
    605598 
    606599 
    607    SUBROUTINE hgr_read 
     600   SUBROUTINE hgr_read( ke1e2u_v ) 
    608601      !!--------------------------------------------------------------------- 
    609602      !!              ***  ROUTINE hgr_read  *** 
    610603      !! 
    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       !!      
     604      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM 
     605      !! 
    616606      !!---------------------------------------------------------------------- 
    617607      USE iom 
    618  
     608      !! 
     609      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 
     610      ! 
    619611      INTEGER ::   inum   ! temporary logical unit 
    620612      !!---------------------------------------------------------------------- 
    621  
     613      ! 
    622614      IF(lwp) THEN 
    623615         WRITE(numout,*) 
     
    625617         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 
    626618      ENDIF 
    627        
     619      ! 
    628620      CALL iom_open( 'coordinates', inum ) 
    629        
     621      ! 
    630622      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 
    631623      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 
    632624      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 
    633625      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 
    634        
     626      ! 
    635627      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 
    636628      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 
    637629      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 
    638630      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        
     631      ! 
     632      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr ) 
     633      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr ) 
     634      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr ) 
     635      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr ) 
     636      ! 
     637      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr ) 
     638      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr ) 
     639      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr ) 
     640      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr ) 
     641      ! 
     642      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN 
     643         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file' 
     644         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr ) 
     645         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr ) 
     646         ke1e2u_v = 1 
     647      ELSE 
     648         ke1e2u_v = 0 
     649      ENDIF 
     650      ! 
    650651      CALL iom_close( inum ) 
    651652       
     653!!gm   THIS is TO BE REMOVED !!!!!!! 
     654 
    652655! need to be define for the extended grid south of -80S 
    653656! some point are undefined but you need to have e1 and e2 .NE. 0 
     
    676679         e2f=1.0e2 
    677680      END WHERE 
     681!!gm end 
    678682        
    679683    END SUBROUTINE hgr_read 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5506 r5737  
    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 
     
    3733   PRIVATE 
    3834 
    39    !! * Routine accessibility 
    4035   PUBLIC  dom_vvl_init       ! called by domain.F90 
    4136   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    4237   PUBLIC  dom_vvl_sf_swp     ! called by step.F90 
    4338   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 
     39 
     40   !                                                      !!* Namelist nam_vvl 
     41   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate 
     42   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate 
     43   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate 
     44   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate 
     45   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate 
     46   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer 
     47   !                                                       ! conservation: not used yet 
     48   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient 
     49   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days] 
     50   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days] 
     51   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation 
     52   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints 
     53 
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport 
     55   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence 
     56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors 
     57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors 
     58   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors 
     59   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6760 
    6861   !! * Substitutions 
     
    372365            DO jj = 1, jpjm1 
    373366               DO ji = 1, fs_jpim1   ! vector opt. 
    374                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
    375                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    376                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
    377                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     367                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)          & 
     368                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     369                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)          &  
     370                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    378371                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    379372                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     
    394387                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    395388                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    396                      &                                            ) * r1_e12t(ji,jj) 
     389                     &                                            ) * r1_e1e2t(ji,jj) 
    397390               END DO 
    398391            END DO 
     
    695688      !!                - vertical interpolation: simple averaging 
    696689      !!---------------------------------------------------------------------- 
    697       !! * Arguments 
    698690      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated 
    699691      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3 
    700692      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors 
    701693      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    702       !! * Local declarations 
     694      ! 
    703695      INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    704696      LOGICAL ::   l_is_orca                                           ! local logical 
     
    717709            DO jj = 1, jpjm1 
    718710               DO ji = 1, fs_jpim1   ! vector opt. 
    719                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
    720                      &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    721                      &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     711                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj)                                   & 
     712                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     713                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    722714               END DO 
    723715            END DO 
    724716         END DO 
    725717         ! 
    726          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    727718         ! boundary conditions 
    728719         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 
     
    735726            DO jj = 1, jpjm1 
    736727               DO ji = 1, fs_jpim1   ! vector opt. 
    737                   pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
    738                      &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    739                      &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     728                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj)                                   & 
     729                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     730                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    740731               END DO 
    741732            END DO 
    742733         END DO 
    743734         ! 
    744          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    745735         ! boundary conditions 
    746736         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 
     
    753743            DO jj = 1, jpjm1 
    754744               DO ji = 1, fs_jpim1   ! vector opt. 
    755                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
    756                      &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    757                      &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     745                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e1e2f(ji,jj)               & 
     746                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     747                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    758748               END DO 
    759749            END DO 
    760750         END DO 
    761751         ! 
    762          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    763752         ! boundary conditions 
    764753         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 
     
    10211010   END SUBROUTINE dom_vvl_ctl 
    10221011 
    1023    SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
    1024       !!--------------------------------------------------------------------- 
    1025       !!                   ***  ROUTINE dom_vvl_orca_fix  *** 
    1026       !!                      
    1027       !! ** Purpose :   Correct surface weighted, horizontally interpolated,  
    1028       !!                scale factors at locations that have been individually 
    1029       !!                modified in domhgr. Such modifications break the 
    1030       !!                relationship between e12t and e1u*e2u etc. 
    1031       !!                Recompute some scale factors ignoring the modified metric. 
    1032       !!---------------------------------------------------------------------- 
    1033       !! * Arguments 
    1034       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated 
    1035       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3 
    1036       CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors 
    1037       !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    1038       !! * Local declarations 
    1039       INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    1040       INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices 
    1041       INTEGER ::   isrow                                               ! index for ORCA1 starting row 
    1042       !! acc 
    1043       !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for 
    1044       !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations 
    1045       !!  
    1046       !                                                ! ===================== 
    1047       IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration 
    1048          !                                             ! ===================== 
    1049       !! acc 
    1050          IF( nn_cla == 0 ) THEN 
    1051             ! 
    1052             ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified) 
    1053             ij0 = 102   ;   ij1 = 102 
    1054             DO jk = 1, jpkm1 
    1055                DO jj = mj0(ij0), mj1(ij1) 
    1056                   DO ji = mi0(ii0), mi1(ii1) 
    1057                      SELECT CASE ( pout ) 
    1058                      CASE( 'U' ) 
    1059                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1060                        &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1061                        &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1062                        &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1063                      CASE( 'F' ) 
    1064                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1065                        &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1066                        &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1067                        &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1068                      END SELECT 
    1069                   END DO 
    1070                END DO 
    1071             END DO 
    1072             ! 
    1073             ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified) 
    1074             ij0 =  88   ;   ij1 =  88 
    1075             DO jk = 1, jpkm1 
    1076                DO jj = mj0(ij0), mj1(ij1) 
    1077                   DO ji = mi0(ii0), mi1(ii1) 
    1078                      SELECT CASE ( pout ) 
    1079                      CASE( 'U' ) 
    1080                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1081                        &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1082                        &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1083                        &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1084                      CASE( 'V' ) 
    1085                         pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1086                        &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1087                        &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1088                        &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1089                      CASE( 'F' ) 
    1090                         pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1091                        &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1092                        &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1093                        &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1094                      END SELECT 
    1095                   END DO 
    1096                END DO 
    1097             END DO 
    1098          ENDIF 
    1099  
    1100          ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified) 
    1101          ij0 = 116   ;   ij1 = 116 
    1102          DO jk = 1, jpkm1 
    1103             DO jj = mj0(ij0), mj1(ij1) 
    1104                DO ji = mi0(ii0), mi1(ii1) 
    1105                   SELECT CASE ( pout ) 
    1106                   CASE( 'U' ) 
    1107                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1108                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1109                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1110                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1111                   CASE( 'F' ) 
    1112                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1113                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1114                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1115                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1116                   END SELECT 
    1117                END DO 
    1118             END DO 
    1119          END DO 
    1120       ENDIF 
    1121       ! 
    1122          !                                             ! ===================== 
    1123       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    1124          !                                             ! ===================== 
    1125          ! This dirty section will be suppressed by simplification process: 
    1126          ! all this will come back in input files 
    1127          ! Currently these hard-wired indices relate to configuration with 
    1128          ! extend grid (jpjglo=332) 
    1129          ! which had a grid-size of 362x292. 
    1130          isrow = 332 - jpjglo 
    1131          ! 
    1132          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified) 
    1133          ij0 = 241 - isrow   ;   ij1 = 241 - isrow 
    1134          DO jk = 1, jpkm1 
    1135             DO jj = mj0(ij0), mj1(ij1) 
    1136                DO ji = mi0(ii0), mi1(ii1) 
    1137                   SELECT CASE ( pout ) 
    1138                   CASE( 'U' ) 
    1139                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1140                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1141                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1142                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1143                   CASE( 'F' ) 
    1144                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1145                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1146                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1147                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1148                   END SELECT 
    1149                END DO 
    1150             END DO 
    1151          END DO 
    1152          ! 
    1153          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
    1154          ij0 = 248 - isrow   ;   ij1 = 248 - isrow 
    1155          DO jk = 1, jpkm1 
    1156             DO jj = mj0(ij0), mj1(ij1) 
    1157                DO ji = mi0(ii0), mi1(ii1) 
    1158                   SELECT CASE ( pout ) 
    1159                   CASE( 'U' ) 
    1160                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &   
    1161                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1162                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1163                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1164                   CASE( 'F' ) 
    1165                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &   
    1166                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1167                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1168                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1169                   END SELECT 
    1170                END DO 
    1171             END DO 
    1172          END DO 
    1173          ! 
    1174          ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
    1175          ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    1176          DO jk = 1, jpkm1 
    1177             DO jj = mj0(ij0), mj1(ij1) 
    1178                DO ji = mi0(ii0), mi1(ii1) 
    1179                   SELECT CASE ( pout ) 
    1180                   CASE( 'V' ) 
    1181                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1182                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1183                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1184                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1185                   END SELECT 
    1186                END DO 
    1187             END DO 
    1188          END DO 
    1189          ! 
    1190          ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
    1191          ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    1192          DO jk = 1, jpkm1 
    1193             DO jj = mj0(ij0), mj1(ij1) 
    1194                DO ji = mi0(ii0), mi1(ii1) 
    1195                   SELECT CASE ( pout ) 
    1196                   CASE( 'V' ) 
    1197                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1198                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1199                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1200                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1201                   END SELECT 
    1202                END DO 
    1203             END DO 
    1204          END DO 
    1205          ! 
    1206          ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
    1207          ij0 = 164 - isrow  ;   ij1 = 165  - isrow   
    1208          DO jk = 1, jpkm1 
    1209             DO jj = mj0(ij0), mj1(ij1) 
    1210                DO ji = mi0(ii0), mi1(ii1) 
    1211                   SELECT CASE ( pout ) 
    1212                   CASE( 'V' ) 
    1213                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1214                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1215                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1216                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1217                   END SELECT 
    1218                END DO 
    1219             END DO 
    1220          END DO 
    1221          ! 
    1222          ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified) 
    1223          ij0 = 164 - isrow    ;   ij1 = 165  - isrow   
    1224          DO jk = 1, jpkm1 
    1225             DO jj = mj0(ij0), mj1(ij1) 
    1226                DO ji = mi0(ii0), mi1(ii1) 
    1227                   SELECT CASE ( pout ) 
    1228                   CASE( 'V' ) 
    1229                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1230                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1231                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1232                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1233                   END SELECT 
    1234                END DO 
    1235             END DO 
    1236          END DO 
    1237          ! 
    1238          ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
    1239          ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    1240          DO jk = 1, jpkm1 
    1241             DO jj = mj0(ij0), mj1(ij1) 
    1242                DO ji = mi0(ii0), mi1(ii1) 
    1243                   SELECT CASE ( pout ) 
    1244                   CASE( 'V' ) 
    1245                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1246                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1247                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1248                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1249                   END SELECT 
    1250                END DO 
    1251             END DO 
    1252          END DO 
    1253          ! 
    1254          ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
    1255          ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    1256          DO jk = 1, jpkm1 
    1257             DO jj = mj0(ij0), mj1(ij1) 
    1258                DO ji = mi0(ii0), mi1(ii1) 
    1259                   SELECT CASE ( pout ) 
    1260                   CASE( 'V' ) 
    1261                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1262                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1263                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1264                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1265                   END SELECT 
    1266                END DO 
    1267             END DO 
    1268          END DO 
    1269       ENDIF 
    1270          !                                             ! ===================== 
    1271       IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration 
    1272          !                                             ! ===================== 
    1273          ! 
    1274          ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified) 
    1275          ij0 = 327   ;   ij1 = 327 
    1276          DO jk = 1, jpkm1 
    1277             DO jj = mj0(ij0), mj1(ij1) 
    1278                DO ji = mi0(ii0), mi1(ii1) 
    1279                   SELECT CASE ( pout ) 
    1280                   CASE( 'U' ) 
    1281                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1282                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1283                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1284                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1285                   CASE( 'F' ) 
    1286                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1287                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1288                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1289                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1290                   END SELECT 
    1291                END DO 
    1292             END DO 
    1293          END DO 
    1294          ! 
    1295          ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified) 
    1296          ij0 = 343   ;   ij1 = 343 
    1297          DO jk = 1, jpkm1 
    1298             DO jj = mj0(ij0), mj1(ij1) 
    1299                DO ji = mi0(ii0), mi1(ii1) 
    1300                   SELECT CASE ( pout ) 
    1301                   CASE( 'U' ) 
    1302                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &   
    1303                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1304                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1305                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1306                   CASE( 'F' ) 
    1307                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &   
    1308                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1309                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1310                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1311                   END SELECT 
    1312                END DO 
    1313             END DO 
    1314          END DO 
    1315          ! 
    1316          ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified) 
    1317          ij0 = 232   ;   ij1 = 232 
    1318          DO jk = 1, jpkm1 
    1319             DO jj = mj0(ij0), mj1(ij1) 
    1320                DO ji = mi0(ii0), mi1(ii1) 
    1321                   SELECT CASE ( pout ) 
    1322                   CASE( 'U' ) 
    1323                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1324                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1325                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1326                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1327                   CASE( 'F' ) 
    1328                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1329                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1330                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1331                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1332                   END SELECT 
    1333                END DO 
    1334             END DO 
    1335          END DO 
    1336          ! 
    1337          ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified) 
    1338          ij0 = 232   ;   ij1 = 232 
    1339          DO jk = 1, jpkm1 
    1340             DO jj = mj0(ij0), mj1(ij1) 
    1341                DO ji = mi0(ii0), mi1(ii1) 
    1342                   SELECT CASE ( pout ) 
    1343                   CASE( 'U' ) 
    1344                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1345                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1346                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1347                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1348                   CASE( 'F' ) 
    1349                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1350                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1351                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1352                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1353                   END SELECT 
    1354                END DO 
    1355             END DO 
    1356          END DO 
    1357          ! 
    1358          ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified) 
    1359          ij0 = 270   ;   ij1 = 270 
    1360          DO jk = 1, jpkm1 
    1361             DO jj = mj0(ij0), mj1(ij1) 
    1362                DO ji = mi0(ii0), mi1(ii1) 
    1363                   SELECT CASE ( pout ) 
    1364                   CASE( 'U' ) 
    1365                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
    1366                     &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) & 
    1367                     &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 
    1368                     &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk) 
    1369                   CASE( 'F' ) 
    1370                      pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
    1371                     &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) & 
    1372                     &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 
    1373                     &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk) 
    1374                   END SELECT 
    1375                END DO 
    1376             END DO 
    1377          END DO 
    1378          ! 
    1379          ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified) 
    1380          ij0 = 232   ;   ij1 = 233 
    1381          DO jk = 1, jpkm1 
    1382             DO jj = mj0(ij0), mj1(ij1) 
    1383                DO ji = mi0(ii0), mi1(ii1) 
    1384                   SELECT CASE ( pout ) 
    1385                   CASE( 'V' ) 
    1386                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1387                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1388                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1389                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1390                   END SELECT 
    1391                END DO 
    1392             END DO 
    1393          END DO 
    1394          ! 
    1395          ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified) 
    1396          ij0 = 276   ;   ij1 = 276 
    1397          DO jk = 1, jpkm1 
    1398             DO jj = mj0(ij0), mj1(ij1) 
    1399                DO ji = mi0(ii0), mi1(ii1) 
    1400                   SELECT CASE ( pout ) 
    1401                   CASE( 'V' ) 
    1402                      pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        & 
    1403                     &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) & 
    1404                     &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 
    1405                     &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk) 
    1406                   END SELECT 
    1407                END DO 
    1408             END DO 
    1409          END DO 
    1410       ENDIF 
    1411    END SUBROUTINE dom_vvl_orca_fix 
    1412  
    14131012   !!====================================================================== 
    14141013END MODULE domvvl 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r5603 r5737  
    2626   PRIVATE 
    2727 
    28    PUBLIC dom_wri        ! routine called by inidom.F90 
    29  
     28   PUBLIC   dom_wri              ! routine called by inidom.F90 
     29   PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90 
    3030   !! * Substitutions 
    3131#  include "vectopt_loop_substitute.h90" 
     
    3636   !!---------------------------------------------------------------------- 
    3737CONTAINS 
     38 
     39 
     40 
     41   SUBROUTINE dom_wri_coordinate 
     42      !!---------------------------------------------------------------------- 
     43      !!                  ***  ROUTINE dom_wri_coordinate  *** 
     44      !!                    
     45      !! ** Purpose :   Create the NetCDF file which contains all the 
     46      !!              standard coordinate information plus the surface, 
     47      !!              e1e2u and e1e2v. By doing so, those surface will 
     48      !!              not be changed by the reduction of e1u or e2v scale  
     49      !!              factors in some straits.  
     50      !!                 NB: call just after the read of standard coordinate 
     51      !!              and the reduction of scale factors in some straits 
     52      !! 
     53      !! ** output file :   coordinate_e1e2u_v.nc 
     54      !!---------------------------------------------------------------------- 
     55      INTEGER           ::   inum0    ! temprary units for 'coordinate_e1e2u_v.nc' file 
     56      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations) 
     57      !                                   !  workspaces 
     58      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw  
     59      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv 
     60      !!---------------------------------------------------------------------- 
     61      ! 
     62      IF( nn_timing == 1 )  CALL timing_start('dom_wri_coordinate') 
     63      ! 
     64      IF(lwp) WRITE(numout,*) 
     65      IF(lwp) WRITE(numout,*) 'dom_wri_coordinate : create NetCDF coordinate file' 
     66      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 
     67       
     68      clnam0 = 'coordinate_e1e2u_v'  ! filename (mesh and mask informations) 
     69       
     70      !  create 'coordinate_e1e2u_v.nc' file 
     71      ! ============================ 
     72      ! 
     73      CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
     74      ! 
     75      !                                                         ! horizontal mesh (inum3) 
     76      CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude 
     77      CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r4 ) 
     78      CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r4 ) 
     79      CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r4 ) 
     80       
     81      CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude 
     82      CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r4 ) 
     83      CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r4 ) 
     84      CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r4 ) 
     85       
     86      CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors 
     87      CALL iom_rstput( 0, 0, inum0, 'e1u', e1u, ktype = jp_r8 ) 
     88      CALL iom_rstput( 0, 0, inum0, 'e1v', e1v, ktype = jp_r8 ) 
     89      CALL iom_rstput( 0, 0, inum0, 'e1f', e1f, ktype = jp_r8 ) 
     90       
     91      CALL iom_rstput( 0, 0, inum0, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors 
     92      CALL iom_rstput( 0, 0, inum0, 'e2u', e2u, ktype = jp_r8 ) 
     93      CALL iom_rstput( 0, 0, inum0, 'e2v', e2v, ktype = jp_r8 ) 
     94      CALL iom_rstput( 0, 0, inum0, 'e2f', e2f, ktype = jp_r8 ) 
     95       
     96      CALL iom_rstput( 0, 0, inum0, 'e1e2u', e1e2u, ktype = jp_r8 ) 
     97      CALL iom_rstput( 0, 0, inum0, 'e1e2v', e1e2v, ktype = jp_r8 ) 
     98 
     99      CALL iom_close( inum0 ) 
     100      ! 
     101      IF( nn_timing == 1 )  CALL timing_stop('dom_wri_coordinate') 
     102      ! 
     103   END SUBROUTINE dom_wri_coordinate 
     104 
    38105 
    39106   SUBROUTINE dom_wri 
Note: See TracChangeset for help on using the changeset viewer.