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 5620 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Ignore:
Timestamp:
2015-07-21T10:55:28+02:00 (9 years ago)
Author:
jamesharle
Message:

Merge with r5619 of trunk, update to unstructured BDY interpolation in
fldread.F90. Structured BDY interpolation incomplete.

Location:
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r4162 r5620  
    7272      !!---------------------------------------------------------------------- 
    7373      INTEGER ::   jc            ! dummy loop indices 
     74      INTEGER :: isrow           ! local index 
    7475      !!---------------------------------------------------------------------- 
    7576       
     
    9192         CASE ( 1 )                                  ! ORCA_R1 configuration 
    9293            !                                        ! ======================= 
     94            ! This dirty section will be suppressed by simplification process: 
     95            ! all this will come back in input files 
     96            ! Currently these hard-wired indices relate to configuration with 
     97            ! extend grid (jpjglo=332) 
     98            isrow = 332 - jpjglo 
     99            ! 
    93100            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea 
    94             ncsi1(1)   = 332  ; ncsj1(1)   = 203 
    95             ncsi2(1)   = 344  ; ncsj2(1)   = 235 
     101            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow 
     102            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow 
    96103            ncsir(1,1) = 1    ; ncsjr(1,1) = 1 
    97104            !                                         
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

    r5038 r5620  
    7373      !!---------------------------------------------------------------------- 
    7474      ! 
     75      ! max number of seconds between each restart 
     76      IF( REAL( nitend - nit000 + 1 ) * rdt > REAL( HUGE( nsec1jan000 ) ) ) THEN 
     77         CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ',   & 
     78            &           'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 
     79      ENDIF 
    7580      ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0 
    7681      IF( MOD( rday     , rdttra(1) ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' ) 
     
    238243               nday_year = 1 
    239244               nsec_year = ndt05 
    240                IF( nsec1jan000 >= 2 * (2**30 - nsecd * nyear_len(1) / 2 ) ) THEN   ! test integer 4 max value 
    241                   CALL ctl_stop( 'The number of seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h ',   & 
    242                      &           'of the current year is exceeding the INTEGER 4 max VALUE: 2^31-1 -> 68.09 years in seconds', & 
    243                      & 'You must do a restart at higher frequency (or remove this STOP and recompile everything in I8)' ) 
    244                ENDIF 
    245245               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1) 
    246246               IF( nleapy == 1 )   CALL day_mth 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5038 r5620  
    162162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphit, gphiu   !: latitude  of t-, u-, v- and f-points (degre) 
    163163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphiv, gphif   !: 
    164    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1t, e2t       !: horizontal scale factors at t-point (m) 
    165    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1u, e2u       !: horizontal scale factors at u-point (m) 
    166    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1v, e2v       !: horizontal scale factors at v-point (m) 
    167    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1f, e2f       !: horizontal scale factors at f-point (m) 
     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) 
    168168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  e1e2t          !: surface at t-point (m2) 
    169169   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 
     
    262262 
    263263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
     264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
    264265 
    265266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
     
    332333   INTEGER FUNCTION dom_oce_alloc() 
    333334      !!---------------------------------------------------------------------- 
    334       INTEGER, DIMENSION(11) :: ierr 
     335      INTEGER, DIMENSION(12) :: ierr 
    335336      !!---------------------------------------------------------------------- 
    336337      ierr(:) = 0 
     
    345346         &      tpol(jpiglo)  , fpol(jpiglo)                               , STAT=ierr(2) ) 
    346347         ! 
    347       ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) ,                      &  
    348          &      glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) ,                      &   
    349          &      glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , e1e2t(jpi,jpj) ,     &   
    350          &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )      
     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) )      
    351353         ! 
    352354      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         & 
     
    400402         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 
    401403 
     404      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
     405 
    402406#if defined key_noslip_accurate 
    403       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) ) 
     407      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
    404408#endif 
    405409      ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5038 r5620  
    135135      !!---------------------------------------------------------------------- 
    136136      USE ioipsl 
    137       NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     137      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
     138         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
    138139         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    139          &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz, nn_euler 
     140         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
    140141      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   & 
    141142         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  & 
     
    169170         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp 
    170171         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in 
     172         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir 
    171173         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out 
     174         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    172175         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
    173176         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
     
    178181         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy 
    179182         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate 
    180          WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock 
     183         IF( ln_rst_list ) THEN 
     184            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist 
     185         ELSE 
     186            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock 
     187         ENDIF 
    181188         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write 
    182189         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn 
    183190         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland 
     191         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta 
    184192         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
    185193         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
     
    195203      ninist = nn_istate 
    196204      nstock = nn_stock 
     205      nstocklist = nn_stocklist 
    197206      nwrite = nn_write 
    198207      neuler = nn_euler 
    199       IF ( neuler == 1 .AND. .NOT.ln_rstart ) THEN 
     208      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN 
    200209         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 ' 
    201210         CALL ctl_warn( ctmp1 ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5038 r5620  
    105105      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 
    106106      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
     107      INTEGER  ::   isrow                ! index for ORCA1 starting row 
     108 
    107109      !!---------------------------------------------------------------------- 
    108110      ! 
     
    159161         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    160162            !                                             ! ===================== 
    161  
    162             ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u = 20 km) 
    163             ij0 = 200   ;   ij1 = 200   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
     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 
    164172            IF(lwp) WRITE(numout,*) 
    165173            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km' 
    166174 
    167             ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km) 
    168             ij0 = 208   ;   ij1 = 208   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
     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 
    169177            IF(lwp) WRITE(numout,*) 
    170178            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km' 
    171179 
    172             ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v = 13 km) 
    173             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3 
     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 
    174182            IF(lwp) WRITE(numout,*) 
    175183            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km' 
    176184 
    177             ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 
    178             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3 
     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 
    179187            IF(lwp) WRITE(numout,*) 
    180188            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km' 
    181189 
    182             ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v = 13 km) 
    183             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 
     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 
    184192            IF(lwp) WRITE(numout,*) 
    185193            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km' 
    186194 
    187             ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v = 20 km) 
    188             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 
     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 
    189197            IF(lwp) WRITE(numout,*) 
    190198            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km' 
    191199 
    192             ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km) 
    193             ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 
     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 
    194202            IF(lwp) WRITE(numout,*) 
    195203            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km' 
    196204 
    197             ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km) 
    198             ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 
     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 
    199207            IF(lwp) WRITE(numout,*) 
    200208            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km' 
    201  
    202             ! 
    203  
    204             ! 
    205             ! 
    206209            ! 
    207210            ! 
     
    471474      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 
    472475      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(:,:) 
    473484 
    474485      ! Control printing : Grid informations (if not restart) 
     
    616627      CALL iom_open( 'coordinates', inum ) 
    617628       
    618       CALL iom_get( inum, jpdom_data, 'glamt', glamt ) 
    619       CALL iom_get( inum, jpdom_data, 'glamu', glamu ) 
    620       CALL iom_get( inum, jpdom_data, 'glamv', glamv ) 
    621       CALL iom_get( inum, jpdom_data, 'glamf', glamf ) 
    622        
    623       CALL iom_get( inum, jpdom_data, 'gphit', gphit ) 
    624       CALL iom_get( inum, jpdom_data, 'gphiu', gphiu ) 
    625       CALL iom_get( inum, jpdom_data, 'gphiv', gphiv ) 
    626       CALL iom_get( inum, jpdom_data, 'gphif', gphif ) 
    627        
    628       CALL iom_get( inum, jpdom_data, 'e1t', e1t ) 
    629       CALL iom_get( inum, jpdom_data, 'e1u', e1u ) 
    630       CALL iom_get( inum, jpdom_data, 'e1v', e1v ) 
    631       CALL iom_get( inum, jpdom_data, 'e1f', e1f ) 
    632        
    633       CALL iom_get( inum, jpdom_data, 'e2t', e2t ) 
    634       CALL iom_get( inum, jpdom_data, 'e2u', e2u ) 
    635       CALL iom_get( inum, jpdom_data, 'e2v', e2v ) 
    636       CALL iom_get( inum, jpdom_data, 'e2f', e2f ) 
     629      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 
     630      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 
     631      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 
     632      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 
     633       
     634      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 
     635      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 
     636      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 
     637      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 
     638       
     639      CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 
     640      CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 
     641      CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 
     642      CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 
     643       
     644      CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 
     645      CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 
     646      CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 
     647      CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 
    637648       
    638649      CALL iom_close( inum ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5038 r5620  
    134134      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
    135135      INTEGER  ::   ios 
     136      INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    136137      INTEGER , POINTER, DIMENSION(:,:) ::  imsk 
    137138      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf 
     
    281282      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    282283 
     284      ! 3. Ocean/land mask at wu-, wv- and w points  
     285      !---------------------------------------------- 
     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) 
     293      END DO 
    283294 
    284295      ! 4. ocean/land mask for the elliptic equation 
     
    391402      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration 
    392403         !                                                 ! Increased lateral friction near of some straits 
     404         ! This dirty section will be suppressed by simplification process: 
     405         ! all this will come back in input files 
     406         ! Currently these hard-wired indices relate to configuration with 
     407         ! extend grid (jpjglo=332) 
     408         ! 
     409         isrow = 332 - jpjglo 
     410         ! 
    393411         IF(lwp) WRITE(numout,*) 
    394412         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : ' 
    395413         IF(lwp) WRITE(numout,*) '      Gibraltar ' 
    396          ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait  
    397          ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     414         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait  
     415         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    398416 
    399417         IF(lwp) WRITE(numout,*) '      Bhosporus ' 
    400          ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait  
    401          ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     418         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait  
     419         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    402420 
    403421         IF(lwp) WRITE(numout,*) '      Makassar (Top) ' 
    404          ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top)  
    405          ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     422         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)  
     423         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    406424 
    407425         IF(lwp) WRITE(numout,*) '      Lombok ' 
    408          ii0 =  44   ;   ii1 =  44        ! Lombok Strait  
    409          ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     426         ii0 =  44           ;   ii1 =  44        ! Lombok Strait  
     427         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    410428 
    411429         IF(lwp) WRITE(numout,*) '      Ombai ' 
    412          ii0 =  53   ;   ii1 =  53        ! Ombai Strait  
    413          ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     430         ii0 =  53           ;   ii1 =  53        ! Ombai Strait  
     431         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    414432 
    415433         IF(lwp) WRITE(numout,*) '      Timor Passage ' 
    416          ii0 =  56   ;   ii1 =  56        ! Timor Passage  
    417          ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     434         ii0 =  56           ;   ii1 =  56        ! Timor Passage  
     435         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    418436 
    419437         IF(lwp) WRITE(numout,*) '      West Halmahera ' 
    420          ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait  
    421          ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     438         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait  
     439         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    422440 
    423441         IF(lwp) WRITE(numout,*) '      East Halmahera ' 
    424          ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait  
    425          ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     442         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait  
     443         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    426444         ! 
    427445      ENDIF 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5038 r5620  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: 
    99   !!                                          vvl option includes z_star and z_tilde coordinates 
     10   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   'key_vvl'                              variable volume 
     
    125126      INTEGER ::   ji,jj,jk 
    126127      INTEGER ::   ii0, ii1, ij0, ij1 
     128      REAL(wp)::   zcoef 
    127129      !!---------------------------------------------------------------------- 
    128130      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init') 
     
    164166      ! t- and w- points depth 
    165167      ! ---------------------- 
     168      ! set the isf depth as it is in the initial step 
    166169      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    167170      fsdepw_n(:,:,1) = 0.0_wp 
     
    169172      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 
    170173      fsdepw_b(:,:,1) = 0.0_wp 
    171       DO jj = 1,jpj 
    172          DO ji = 1,jpi 
    173             DO jk = 2,mikt(ji,jj)-1 
    174                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    175                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    176                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    177                fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 
    178                fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    179             END DO 
    180             IF (mikt(ji,jj) .GT. 1) THEN 
    181                jk = mikt(ji,jj) 
    182                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    183                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    184                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    185                fsdept_b(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_b(ji,jj,jk) 
    186                fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    187             END IF 
    188             DO jk = mikt(ji,jj)+1, jpk 
    189                fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     174 
     175      DO jk = 2, jpk 
     176         DO jj = 1,jpj 
     177            DO ji = 1,jpi 
     178              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     179                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     180                                                     ! 0.5 where jk = mikt   
     181               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    190182               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    191                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    192                fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 
     183               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     184                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
     185               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    193186               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 
     187               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  & 
     188                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk))  
    194189            END DO 
    195190         END DO 
     
    588583      INTEGER, INTENT( in )               :: kt       ! time step 
    589584      !! * Local declarations 
    590       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 
    591585      INTEGER                             :: ji,jj,jk       ! dummy loop indices 
     586      REAL(wp)                            :: zcoef 
    592587      !!---------------------------------------------------------------------- 
    593588 
    594589      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp') 
    595       ! 
    596       CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                ) 
    597590      ! 
    598591      IF( kt == nit000 )   THEN 
     
    638631      ! t- and w- points depth 
    639632      ! ---------------------- 
     633      ! set the isf depth as it is in the initial step 
    640634      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    641635      fsdepw_n(:,:,1) = 0.0_wp 
    642636      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
    643       DO jj = 1,jpj 
    644          DO ji = 1,jpi 
    645             DO jk = 2,mikt(ji,jj)-1 
    646                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    647                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    648                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    649             END DO 
    650             IF (mikt(ji,jj) .GT. 1) THEN 
    651                jk = mikt(ji,jj) 
    652                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    653                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    654                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    655             END IF 
    656             DO jk = mikt(ji,jj)+1, jpk 
    657                fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     637 
     638      DO jk = 2, jpk 
     639         DO jj = 1,jpj 
     640            DO ji = 1,jpi 
     641              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     642                                                                 ! 1 for jk = mikt 
     643               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    658644               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    659                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
     645               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     646                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
     647               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    660648            END DO 
    661649         END DO 
    662650      END DO 
     651 
    663652      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
    664653      ! ---------------------------------------------------------------------------------- 
     
    679668      ! Write outputs 
    680669      ! ============= 
    681       z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    682       CALL iom_put( "cellthc" , fse3t_n  (:,:,:) ) 
     670      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
     671      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
     672      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
     673      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
    683674      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
    684       CALL iom_put( "e3tdef"  , z_e3t_def(:,:,:) ) 
     675      IF( iom_use("e3tdef") )   & 
     676         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    685677 
    686678      ! write restart file 
    687679      ! ================== 
    688680      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 
    689       ! 
    690       CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 
    691681      ! 
    692682      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp') 
     
    10491039      INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    10501040      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices 
     1041      INTEGER ::   isrow                                               ! index for ORCA1 starting row 
    10511042      !! acc 
    10521043      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for 
     
    11321123      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    11331124         !                                             ! ===================== 
    1134          ! 
    1135          ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified) 
    1136          ij0 = 200   ;   ij1 = 200 
     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 
    11371134         DO jk = 1, jpkm1 
    11381135            DO jj = mj0(ij0), mj1(ij1) 
     
    11541151         END DO 
    11551152         ! 
    1156          ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
    1157          ij0 = 208   ;   ij1 = 208 
     1153         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
     1154         ij0 = 248 - isrow   ;   ij1 = 248 - isrow 
    11581155         DO jk = 1, jpkm1 
    11591156            DO jj = mj0(ij0), mj1(ij1) 
     
    11751172         END DO 
    11761173         ! 
    1177          ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
    1178          ij0 = 124   ;   ij1 = 125 
     1174         ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
     1175         ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    11791176         DO jk = 1, jpkm1 
    11801177            DO jj = mj0(ij0), mj1(ij1) 
     
    11911188         END DO 
    11921189         ! 
    1193          ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
    1194          ij0 = 124   ;   ij1 = 125 
     1190         ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
     1191         ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    11951192         DO jk = 1, jpkm1 
    11961193            DO jj = mj0(ij0), mj1(ij1) 
     
    12071204         END DO 
    12081205         ! 
    1209          ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
    1210          ij0 = 124   ;   ij1 = 125 
     1206         ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
     1207         ij0 = 164 - isrow  ;   ij1 = 165  - isrow   
    12111208         DO jk = 1, jpkm1 
    12121209            DO jj = mj0(ij0), mj1(ij1) 
     
    12231220         END DO 
    12241221         ! 
    1225          ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified) 
    1226          ij0 = 124   ;   ij1 = 125 
     1222         ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified) 
     1223         ij0 = 164 - isrow    ;   ij1 = 165  - isrow   
    12271224         DO jk = 1, jpkm1 
    12281225            DO jj = mj0(ij0), mj1(ij1) 
     
    12391236         END DO 
    12401237         ! 
    1241          ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
    1242          ij0 = 141   ;   ij1 = 142 
     1238         ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
     1239         ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    12431240         DO jk = 1, jpkm1 
    12441241            DO jj = mj0(ij0), mj1(ij1) 
     
    12551252         END DO 
    12561253         ! 
    1257          ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
    1258          ij0 = 141   ;   ij1 = 142 
     1254         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
     1255         ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    12591256         DO jk = 1, jpkm1 
    12601257            DO jj = mj0(ij0), mj1(ij1) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r5038 r5620  
    215215         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system 
    216216         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d ) 
     217         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )      
     218         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 ) 
    217219      ENDIF 
    218220       
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5038 r5620  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
     19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye   
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    3536   USE oce               ! ocean variables 
    3637   USE dom_oce           ! ocean domain 
    37    USE sbc_oce           ! surface variable (isf) 
    3838   USE closea            ! closed seas 
    3939   USE c1d               ! 1D vertical configuration 
     
    298298      ENDIF 
    299299 
     300      IF ( ln_isfcav ) THEN 
    300301! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    301302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
    302       DO jk = 1, jpkm1 
    303          e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
    304       END DO 
    305       e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
    306  
    307       DO jk = 2, jpk 
    308          e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
    309       END DO 
    310       e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     303         DO jk = 1, jpkm1 
     304            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
     305         END DO 
     306         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
     307 
     308         DO jk = 2, jpk 
     309            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
     310         END DO 
     311         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     312      END IF 
    311313 
    312314!!gm BUG in s-coordinate this does not work! 
     
    365367      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
    366368      INTEGER  ::   inum                      ! temporary logical unit 
     369      INTEGER  ::   ierror                    ! error flag 
    367370      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position 
    368371      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices 
    369372      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics  
    370373      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars 
    371       INTEGER , POINTER, DIMENSION(:,:) ::   idta   ! global domain integer data 
    372       REAL(wp), POINTER, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
     374      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data 
     375      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
    373376      !!---------------------------------------------------------------------- 
    374377      ! 
    375378      IF( nn_timing == 1 )  CALL timing_start('zgr_bat') 
    376       ! 
    377       CALL wrk_alloc( jpidta, jpjdta, idta ) 
    378       CALL wrk_alloc( jpidta, jpjdta, zdta ) 
    379379      ! 
    380380      IF(lwp) WRITE(numout,*) 
     
    385385         !                                            ! ================== ! 
    386386         !                                            ! global domain level and meter bathymetry (idta,zdta) 
     387         ! 
     388         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 
     389         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 
     390         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 
     391         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 
    387392         ! 
    388393         IF( ntopo == 0 ) THEN                        ! flat basin 
     
    468473         misfdep(:,:)=1 
    469474         ! 
    470          ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    471          IF( cp_cfg == "isomip" ) THEN  
    472            !  
    473            risfdep(:,:)=200.e0  
    474            misfdep(:,:)=1  
    475            ij0 = 1 ; ij1 = 40  
    476            DO jj = mj0(ij0), mj1(ij1)  
    477               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
    478                 END DO  
    479             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
    480            !  
    481          ELSEIF ( cp_cfg == "isomip2" ) THEN 
    482          !  
    483             risfdep(:,:)=0.e0 
    484             misfdep(:,:)=1 
    485             ij0 = 1 ; ij1 = 40 
    486             DO jj = mj0(ij0), mj1(ij1) 
    487                risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
    488             END DO 
    489             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    490          END IF 
     475         DEALLOCATE( idta, zdta ) 
    491476         ! 
    492477         !                                            ! ================ ! 
     
    529514         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry 
    530515            CALL iom_open ( 'bathy_meter.nc', inum )  
    531             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
     516            IF ( ln_isfcav ) THEN 
     517               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 
     518            ELSE 
     519               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  ) 
     520            END IF 
    532521            CALL iom_close( inum ) 
    533             !   
     522            !                                                 
    534523            risfdep(:,:)=0._wp          
    535524            misfdep(:,:)=1              
     
    579568      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    580569         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
    581          WHERE (bathy == risfdep) 
    582             bathy   = 0.0_wp ; risfdep = 0.0_wp 
    583          END WHERE 
     570         IF ( ln_isfcav ) THEN 
     571            WHERE (bathy == risfdep) 
     572               bathy   = 0.0_wp ; risfdep = 0.0_wp 
     573            END WHERE 
     574         END IF 
    584575         ! end patch 
    585576         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
     
    592583         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 
    593584      ENDIF 
    594       ! 
    595       CALL wrk_dealloc( jpidta, jpjdta, idta ) 
    596       CALL wrk_dealloc( jpidta, jpjdta, zdta ) 
    597585      ! 
    598586      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat') 
     
    959947      !!---------------------------------------------------------------------- 
    960948      !! 
     949      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     950      INTEGER  ::   ik, it, ikb, ikt ! temporary integers 
     951      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     952      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     953      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     954      REAL(wp) ::   zmax             ! Maximum depth 
     955      REAL(wp) ::   zdiff            ! temporary scalar 
     956      REAL(wp) ::   zrefdep          ! temporary scalar 
     957      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     958      !!--------------------------------------------------------------------- 
     959      ! 
     960      IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
     961      ! 
     962      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     963      ! 
     964      IF(lwp) WRITE(numout,*) 
     965      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
     966      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
     967      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     968 
     969      ll_print = .FALSE.                   ! Local variable for debugging 
     970       
     971      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
     972         WRITE(numout,*) 
     973         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     974         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     975      ENDIF 
     976 
     977 
     978      ! bathymetry in level (from bathy_meter) 
     979      ! =================== 
     980      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
     981      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     982      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     983      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
     984      END WHERE 
     985 
     986      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     987      ! find the number of ocean levels such that the last level thickness 
     988      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
     989      ! e3t_1d is the reference level thickness 
     990      DO jk = jpkm1, 1, -1 
     991         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     992         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     993      END DO 
     994 
     995      IF ( ln_isfcav ) CALL zgr_isf 
     996 
     997      ! Scale factors and depth at T- and W-points 
     998      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     999         gdept_0(:,:,jk) = gdept_1d(jk) 
     1000         gdepw_0(:,:,jk) = gdepw_1d(jk) 
     1001         e3t_0  (:,:,jk) = e3t_1d  (jk) 
     1002         e3w_0  (:,:,jk) = e3w_1d  (jk) 
     1003      END DO 
     1004      !  
     1005      DO jj = 1, jpj 
     1006         DO ji = 1, jpi 
     1007            ik = mbathy(ji,jj) 
     1008            IF( ik > 0 ) THEN               ! ocean point only 
     1009               ! max ocean level case 
     1010               IF( ik == jpkm1 ) THEN 
     1011                  zdepwp = bathy(ji,jj) 
     1012                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1013                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1014                  e3t_0(ji,jj,ik  ) = ze3tp 
     1015                  e3t_0(ji,jj,ik+1) = ze3tp 
     1016                  e3w_0(ji,jj,ik  ) = ze3wp 
     1017                  e3w_0(ji,jj,ik+1) = ze3tp 
     1018                  gdepw_0(ji,jj,ik+1) = zdepwp 
     1019                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1020                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1021                  ! 
     1022               ELSE                         ! standard case 
     1023                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1024                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1025                  ENDIF 
     1026!gm Bug?  check the gdepw_1d 
     1027                  !       ... on ik 
     1028                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1029                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1030                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1031                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1032                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1033                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1034                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1035                  !       ... on ik+1 
     1036                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1037                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1038                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     1039               ENDIF 
     1040            ENDIF 
     1041         END DO 
     1042      END DO 
     1043      ! 
     1044      it = 0 
     1045      DO jj = 1, jpj 
     1046         DO ji = 1, jpi 
     1047            ik = mbathy(ji,jj) 
     1048            IF( ik > 0 ) THEN               ! ocean point only 
     1049               e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1050               e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1051               ! test 
     1052               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1053               IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1054                  it = it + 1 
     1055                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1056                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1057                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1058                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1059               ENDIF 
     1060            ENDIF 
     1061         END DO 
     1062      END DO 
     1063      ! 
     1064      IF ( ln_isfcav ) THEN 
     1065      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1066         DO jj = 1, jpj  
     1067            DO ji = 1, jpi  
     1068               ik = misfdep(ji,jj)  
     1069               IF( ik > 1 ) THEN               ! ice shelf point only  
     1070                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1071                  gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1072!gm Bug?  check the gdepw_0  
     1073               !       ... on ik  
     1074                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1075                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1076                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1077                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1078                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1079 
     1080                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1081                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1082                  ENDIF  
     1083               !       ... on ik / ik-1  
     1084                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1085                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1086! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1087                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1088               ENDIF  
     1089            END DO  
     1090         END DO  
     1091      !  
     1092         it = 0  
     1093         DO jj = 1, jpj  
     1094            DO ji = 1, jpi  
     1095               ik = misfdep(ji,jj)  
     1096               IF( ik > 1 ) THEN               ! ice shelf point only  
     1097                  e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1098                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1099               ! test  
     1100                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1101                  IF( zdiff <= 0. .AND. lwp ) THEN   
     1102                     it = it + 1  
     1103                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1104                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1105                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1106                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1107                  ENDIF  
     1108               ENDIF  
     1109            END DO  
     1110         END DO  
     1111      END IF 
     1112      ! END (ISF) 
     1113 
     1114      ! Scale factors and depth at U-, V-, UW and VW-points 
     1115      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1116         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1117         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1118         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1119         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1120      END DO 
     1121      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1122         DO jj = 1, jpjm1 
     1123            DO ji = 1, fs_jpim1   ! vector opt. 
     1124               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1125               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1126               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1127               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1128            END DO 
     1129         END DO 
     1130      END DO 
     1131      IF ( ln_isfcav ) THEN 
     1132      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1133         DO jj = 2, jpjm1  
     1134            DO ji = 2, fs_jpim1   ! vector opt.  
     1135               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     1136               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     1137               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     1138                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     1139               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     1140               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     1141               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     1142                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     1143            END DO 
     1144         END DO 
     1145      END IF 
     1146 
     1147      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1148      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     1149      ! 
     1150      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1151         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1152         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1153         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1154         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1155      END DO 
     1156       
     1157      ! Scale factor at F-point 
     1158      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1159         e3f_0(:,:,jk) = e3t_1d(jk) 
     1160      END DO 
     1161      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1162         DO jj = 1, jpjm1 
     1163            DO ji = 1, fs_jpim1   ! vector opt. 
     1164               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1165            END DO 
     1166         END DO 
     1167      END DO 
     1168      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1169      ! 
     1170      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1171         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1172      END DO 
     1173!!gm  bug ? :  must be a do loop with mj0,mj1 
     1174      !  
     1175      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1176      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1177      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1178      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1179      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1180 
     1181      ! Control of the sign 
     1182      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1183      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1184      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1185      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1186      
     1187      ! Compute gdep3w_0 (vertical sum of e3w) 
     1188      IF ( ln_isfcav ) THEN ! if cavity 
     1189         WHERE (misfdep == 0) misfdep = 1 
     1190         DO jj = 1,jpj 
     1191            DO ji = 1,jpi 
     1192               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1193               DO jk = 2, misfdep(ji,jj) 
     1194                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1195               END DO 
     1196               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     1197               DO jk = misfdep(ji,jj) + 1, jpk 
     1198                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1199               END DO 
     1200            END DO 
     1201         END DO 
     1202      ELSE ! no cavity 
     1203         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1204         DO jk = 2, jpk 
     1205            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1206         END DO 
     1207      END IF 
     1208      !                                               ! ================= ! 
     1209      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     1210         !                                            ! ================= ! 
     1211         DO jj = 1,jpj 
     1212            DO ji = 1, jpi 
     1213               ik = MAX( mbathy(ji,jj), 1 ) 
     1214               zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
     1215               zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
     1216               zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
     1217               zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
     1218               zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
     1219               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
     1220            END DO 
     1221         END DO 
     1222         WRITE(numout,*) 
     1223         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1224         WRITE(numout,*) 
     1225         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1226         WRITE(numout,*) 
     1227         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1228         WRITE(numout,*) 
     1229         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1230         WRITE(numout,*) 
     1231         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1232         WRITE(numout,*) 
     1233         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1234      ENDIF   
     1235      ! 
     1236      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1237      ! 
     1238      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
     1239      ! 
     1240   END SUBROUTINE zgr_zps 
     1241 
     1242   SUBROUTINE zgr_isf 
     1243      !!---------------------------------------------------------------------- 
     1244      !!                    ***  ROUTINE zgr_isf  *** 
     1245      !!    
     1246      !! ** Purpose :   check the bathymetry in levels 
     1247      !!    
     1248      !! ** Method  :   THe water column have to contained at least 2 cells 
     1249      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1250      !!                this criterion. 
     1251      !!                  
     1252      !!    
     1253      !! ** Action  : - test compatibility between isfdraft and bathy  
     1254      !!              - bathy and isfdraft are modified 
     1255      !!---------------------------------------------------------------------- 
     1256      !!    
    9611257      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    9621258      INTEGER  ::   ik, it           ! temporary integers 
     
    9691265      REAL(wp) ::   zdiff            ! temporary scalar 
    9701266      REAL(wp) ::   zrefdep          ! temporary scalar 
    971       REAL(wp) ::   zbathydiff, zrisfdepdiff  
    972       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH) 
    973       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH) 
    974       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     1267      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1268      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1269      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    9751270      !!--------------------------------------------------------------------- 
    9761271      ! 
    977       IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
    978       ! 
    979       CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     1272      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1273      ! 
    9801274      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    981       CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 
    982       ! 
    983       IF(lwp) WRITE(numout,*) 
    984       IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
    985       IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
    986       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    987  
    988       ll_print = .FALSE.                   ! Local variable for debugging 
    989        
    990       IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
    991          WRITE(numout,*) 
    992          WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
    993          CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
    994       ENDIF 
    995  
    996       ! bathymetry in level (from bathy_meter) 
    997       ! =================== 
    998       zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
    999       bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    1000       WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
    1001       ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
    1002       END WHERE 
    1003  
    1004       ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    1005       ! find the number of ocean levels such that the last level thickness 
    1006       ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
    1007       ! e3t_1d is the reference level thickness 
    1008       DO jk = jpkm1, 1, -1 
    1009          zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1010          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    1011       END DO 
     1275      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1276 
     1277 
    10121278      ! (ISF) compute misfdep 
    10131279      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
     
    10531319            misfdep(jpi,:) = misfdep(  2  ,:)  
    10541320         ENDIF 
    1055   
     1321 
    10561322         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    10571323            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    10581324            mbathy(jpi,:) = mbathy(  2  ,:) 
    10591325         ENDIF 
    1060   
     1326 
    10611327         ! split last cell if possible (only where water column is 2 cell or less) 
    10621328         DO jk = jpkm1, 1, -1 
     
    10761342            END WHERE 
    10771343         END DO 
    1078   
     1344 
    10791345  
    10801346 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
     
    12521518  
    12531519 ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1254          DO jk = 1, jpk 
     1520         DO jk = 2, jpk 
    12551521            WHERE (misfdep==0) misfdep=jpk 
    12561522            zmask=0 
     
    13571623               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    13581624               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1359                IF( ibtest == 0 ) THEN 
     1625               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    13601626                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    13611627               END IF 
     
    14731739      ENDIF  
    14741740 
    1475       ! Scale factors and depth at T- and W-points 
    1476       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    1477          gdept_0(:,:,jk) = gdept_1d(jk) 
    1478          gdepw_0(:,:,jk) = gdepw_1d(jk) 
    1479          e3t_0  (:,:,jk) = e3t_1d  (jk) 
    1480          e3w_0  (:,:,jk) = e3w_1d  (jk) 
    1481       END DO 
    1482       !  
    1483       DO jj = 1, jpj 
    1484          DO ji = 1, jpi 
    1485             ik = mbathy(ji,jj) 
    1486             IF( ik > 0 ) THEN               ! ocean point only 
    1487                ! max ocean level case 
    1488                IF( ik == jpkm1 ) THEN 
    1489                   zdepwp = bathy(ji,jj) 
    1490                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1491                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1492                   e3t_0(ji,jj,ik  ) = ze3tp 
    1493                   e3t_0(ji,jj,ik+1) = ze3tp 
    1494                   e3w_0(ji,jj,ik  ) = ze3wp 
    1495                   e3w_0(ji,jj,ik+1) = ze3tp 
    1496                   gdepw_0(ji,jj,ik+1) = zdepwp 
    1497                   gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1498                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1499                   ! 
    1500                ELSE                         ! standard case 
    1501                   IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1502                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1503                   ENDIF 
    1504 !gm Bug?  check the gdepw_1d 
    1505                   !       ... on ik 
    1506                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1507                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1508                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1509                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1510                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1511                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1512                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1513                   !       ... on ik+1 
    1514                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1515                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1516                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1517                ENDIF 
    1518             ENDIF 
    1519          END DO 
    1520       END DO 
    1521       ! 
    1522       it = 0 
    1523       DO jj = 1, jpj 
    1524          DO ji = 1, jpi 
    1525             ik = mbathy(ji,jj) 
    1526             IF( ik > 0 ) THEN               ! ocean point only 
    1527                e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1528                e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1529                ! test 
    1530                zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1531                IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1532                   it = it + 1 
    1533                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1534                   WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1535                   WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1536                   WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1537                ENDIF 
    1538             ENDIF 
    1539          END DO 
    1540       END DO 
    1541       ! 
    1542       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1543       DO jj = 1, jpj  
    1544          DO ji = 1, jpi  
    1545             ik = misfdep(ji,jj)  
    1546             IF( ik > 1 ) THEN               ! ice shelf point only  
    1547                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1548                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1549 !gm Bug?  check the gdepw_0  
    1550                !       ... on ik  
    1551                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1552                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1553                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1554                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1555                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1556  
    1557                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1558                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1559                 ENDIF  
    1560                !       ... on ik / ik-1  
    1561                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1562                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1563 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1564                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    1565             ENDIF  
    1566          END DO  
    1567       END DO  
    1568       !  
    1569       it = 0  
    1570       DO jj = 1, jpj  
    1571          DO ji = 1, jpi  
    1572             ik = misfdep(ji,jj)  
    1573             IF( ik > 1 ) THEN               ! ice shelf point only  
    1574                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1575                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1576                ! test  
    1577                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1578                IF( zdiff <= 0. .AND. lwp ) THEN   
    1579                   it = it + 1  
    1580                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1581                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1582                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1583                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1584                ENDIF  
    1585             ENDIF  
    1586          END DO  
    1587       END DO  
    1588       ! END (ISF) 
    1589  
    1590       ! Scale factors and depth at U-, V-, UW and VW-points 
    1591       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1592          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1593          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1594          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1595          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1596       END DO 
    1597       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1598          DO jj = 1, jpjm1 
    1599             DO ji = 1, fs_jpim1   ! vector opt. 
    1600                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1601                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1602                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1603                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1604             END DO 
    1605          END DO 
    1606       END DO 
    1607       ! (ISF) define e3uw 
    1608       DO jk = 2,jpk                           
    1609          DO jj = 1, jpjm1  
    1610             DO ji = 1, fs_jpim1   ! vector opt.  
    1611                e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
    1612                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
    1613                e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
    1614                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
    1615             END DO  
    1616          END DO  
    1617       END DO 
    1618       !End (ISF) 
    1619        
    1620       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1621       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    1622       ! 
    1623       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1624          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1625          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1626          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1627          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1628       END DO 
    1629        
    1630       ! Scale factor at F-point 
    1631       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1632          e3f_0(:,:,jk) = e3t_1d(jk) 
    1633       END DO 
    1634       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1635          DO jj = 1, jpjm1 
    1636             DO ji = 1, fs_jpim1   ! vector opt. 
    1637                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1638             END DO 
    1639          END DO 
    1640       END DO 
    1641       CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1642       ! 
    1643       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1644          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1645       END DO 
    1646 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1647       !  
    1648       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1649       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1650       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1651       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1652       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1653  
    1654       ! Control of the sign 
    1655       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1656       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1657       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1658       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1659       
    1660       ! Compute gdep3w_0 (vertical sum of e3w) 
    1661       WHERE (misfdep == 0) misfdep = 1 
    1662       DO jj = 1,jpj 
    1663          DO ji = 1,jpi 
    1664             gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1665             DO jk = 2, misfdep(ji,jj) 
    1666                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1667             END DO 
    1668             IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
    1669             DO jk = misfdep(ji,jj) + 1, jpk 
    1670                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1671             END DO 
    1672         END DO 
    1673       END DO 
    1674       !                                               ! ================= ! 
    1675       IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    1676          !                                            ! ================= ! 
    1677          DO jj = 1,jpj 
    1678             DO ji = 1, jpi 
    1679                ik = MAX( mbathy(ji,jj), 1 ) 
    1680                zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
    1681                zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
    1682                zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
    1683                zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
    1684                zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
    1685                zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
    1686             END DO 
    1687          END DO 
    1688          WRITE(numout,*) 
    1689          WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1690          WRITE(numout,*) 
    1691          WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1692          WRITE(numout,*) 
    1693          WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1694          WRITE(numout,*) 
    1695          WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1696          WRITE(numout,*) 
    1697          WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1698          WRITE(numout,*) 
    1699          WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1700       ENDIF   
    1701       ! 
    1702       CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
    17031741      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    17041742      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1705       ! 
    1706       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1707       ! 
    1708    END SUBROUTINE zgr_zps 
     1743 
     1744      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
     1745       
     1746   END SUBROUTINE 
    17091747 
    17101748   SUBROUTINE zgr_sco 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    3939   !!---------------------------------------------------------------------- 
    4040   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    41    !! $Id: dtatem.F90 2392 2010-11-15 21:20:05Z gm $  
     41   !! $Id$  
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5038 r5620  
    6969      !! ** Purpose :   Initialization of the dynamics and tracer fields. 
    7070      !!---------------------------------------------------------------------- 
    71       ! - ML - needed for initialization of e3t_b 
    72       INTEGER  ::  ji,jj,jk     ! dummy loop indices 
    73       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::  zuvd    ! U & V data workspace 
     71      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     72      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace 
    7473      !!---------------------------------------------------------------------- 
    7574      ! 
     
    8483      IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data 
    8584 
    86       rhd  (:,:,:  ) = 0._wp 
    87       rhop (:,:,:  ) = 0._wp 
    88       rn2  (:,:,:  ) = 0._wp 
    89       tsa  (:,:,:,:) = 0._wp    
    90       rab_b(:,:,:,:) = 0._wp 
    91       rab_n(:,:,:,:) = 0._wp 
     85      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
     86      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
     87      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk 
     88      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    9289 
    9390      IF( ln_rstart ) THEN                    ! Restart from a file 
     
    113110         ELSEIF( cp_cfg == 'gyre' ) THEN          
    114111            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields 
    115         ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN 
    116             IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain'  
    117             tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields  
    118             tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:) 
    119             tsb(:,:,:,:)=tsn(:,:,:,:)   
    120112         ELSE                                    ! Initial T-S, U-V fields read in files 
    121113            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000 
     
    137129         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    138130#if ! defined key_c1d 
    139          IF( ln_zps )    CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    140             &                                      rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
    141             &                                      gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     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 
    142138#endif 
    143139         !    
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r5038 r5620  
    4141   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin] 
    4242#if defined key_lim3 
    43    REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp        !: melting point of snow         [Kelvin] 
    44    REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp        !: melting point of ice          [Kelvin] 
     43   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
     44   REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin] 
    4545#else 
    4646   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
     
    5151   REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat           [J/Kelvin] 
    5252   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J] 
     53   REAL(wp), PUBLIC ::   rau0_rcp                    !: = rau0 * rcp  
    5354   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp ) 
    5455 
     
    8283   REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3] 
    8384   REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3] 
     85#endif 
     86#if defined key_lim3 
     87   REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic 
     88   REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn 
    8489#endif 
    8590   !!---------------------------------------------------------------------- 
     
    166171      lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice 
    167172#endif 
    168  
     173#if defined key_lim3 
     174      r1_rhoic = 1._wp / rhoic 
     175      r1_rhosn = 1._wp / rhosn 
     176#endif 
    169177      IF(lwp) THEN 
    170178         WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.