Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (9 months ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
3 deleted
12 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
  • NEMO/trunk/src/OCE/DOM/closea.F90

    r10425 r12377  
    1111   !!             4.0  !  2016-06  (G. Madec)  move to usrdef_closea, remove clo_ups 
    1212   !!             4.0  !  2017-12  (D. Storkey) new formulation based on masks read from file 
     13   !!             4.1  !  2019-07  (P. Mathiot) update to the new domcfg.nc input file 
    1314   !!---------------------------------------------------------------------- 
    1415 
    1516   !!---------------------------------------------------------------------- 
    1617   !!   dom_clo    : read in masks which define closed seas and runoff areas 
    17    !!   sbc_clo    : Special handling of freshwater fluxes over closed seas 
    1818   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf) 
    19    !!   clo_bat    : set to zero a field over closed sea (see domzgr) 
    20    !!---------------------------------------------------------------------- 
    21    USE oce             ! dynamics and tracers 
    22    USE dom_oce         ! ocean space and time domain 
    23    USE phycst          ! physical constants 
    24    USE sbc_oce         ! ocean surface boundary conditions 
    25    USE iom             ! I/O routines 
     19   !!   clo_msk    : set to zero a field over closed sea (see domzgr) 
     20   !!---------------------------------------------------------------------- 
     21   USE in_out_manager  ! I/O manager 
    2622   ! 
    27    USE in_out_manager  ! I/O manager 
    28    USE lib_fortran,    ONLY: glob_sum 
    29    USE lbclnk          ! lateral boundary condition - MPP exchanges 
    30    USE lib_mpp         ! MPP library 
    31    USE timing          ! Timing 
     23   USE diu_bulk    , ONLY: ln_diurnal_only            ! used for sanity check 
     24   USE iom         , ONLY: iom_open, iom_get, iom_close, jpdom_data ! I/O routines 
     25   USE lib_fortran , ONLY: glob_sum                   ! fortran library 
     26   USE lib_mpp     , ONLY: mpp_max, ctl_nam, ctl_stop ! MPP library 
    3227 
    3328   IMPLICIT NONE 
     29 
    3430   PRIVATE 
    3531 
    3632   PUBLIC dom_clo      ! called by domain module 
    37    PUBLIC sbc_clo      ! called by sbcmod module 
    3833   PUBLIC clo_rnf      ! called by sbcrnf module 
    39    PUBLIC clo_bat      ! called in domzgr module 
    40  
    41    LOGICAL, PUBLIC :: ln_closea  !:  T => keep closed seas (defined by closea_mask field) in the domain and apply 
    42                                  !:       special treatment of freshwater fluxes. 
    43                                  !:  F => suppress closed seas (defined by closea_mask field) from the bathymetry 
    44                                  !:       at runtime. 
    45                                  !:  If there is no closea_mask field in the domain_cfg file or we do not use 
    46                                  !:  a domain_cfg file then this logical does nothing. 
    47                                  !: 
    48    LOGICAL, PUBLIC :: l_sbc_clo  !: T => Closed seas defined, apply special treatment of freshwater fluxes. 
    49                                  !: F => No closed seas defined (closea_mask field not found). 
    50    LOGICAL, PUBLIC :: l_clo_rnf  !: T => Some closed seas output freshwater (RNF or EMPMR) to specified runoff points. 
    51    INTEGER, PUBLIC :: jncs       !: number of closed seas (inferred from closea_mask field) 
    52    INTEGER, PUBLIC :: jncsr      !: number of closed seas rnf mappings (inferred from closea_mask_rnf field) 
    53    INTEGER, PUBLIC :: jncse      !: number of closed seas empmr mappings (inferred from closea_mask_empmr field) 
    54     
    55    INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask       !: mask of integers defining closed seas 
    56    INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_rnf   !: mask of integers defining closed seas rnf mappings 
    57    INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_empmr !: mask of integers defining closed seas empmr mappings 
    58    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surf         !: closed sea surface areas  
    59                                                                   !: (and residual global surface area)  
    60    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfr        !: closed sea target rnf surface areas  
    61    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfe        !: closed sea target empmr surface areas  
    62  
    63    !! * Substitutions 
    64 #  include "vectopt_loop_substitute.h90" 
     34   PUBLIC clo_msk      ! called in domzgr module 
     35 
     36   LOGICAL, PUBLIC :: ln_maskcs        !: logical to mask all closed sea 
     37   LOGICAL, PUBLIC :: ln_mask_csundef  !: logical to mask all undefined closed sea 
     38   LOGICAL, PUBLIC :: ln_clo_rnf       !: closed sea treated as runoff (update rnf mask) 
     39 
     40   LOGICAL, PUBLIC :: l_sbc_clo  !: T => net evap/precip over closed seas spread outover the globe/river mouth 
     41   LOGICAL, PUBLIC :: l_clo_rnf  !: T => Some closed seas output freshwater (RNF) to specified runoff points. 
     42 
     43   INTEGER, PUBLIC :: ncsg      !: number of closed seas global mappings (inferred from closea_mask_glo field) 
     44   INTEGER, PUBLIC :: ncsr      !: number of closed seas rnf    mappings (inferred from closea_mask_rnf field) 
     45   INTEGER, PUBLIC :: ncse      !: number of closed seas empmr  mappings (inferred from closea_mask_emp field) 
     46 
     47   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_opnsea, mask_csundef  !: mask defining the open sea and the undefined closed sea 
     48  
     49   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csglo , mask_csgrpglo !: mask of integers defining closed seas 
     50   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csrnf , mask_csgrprnf !: mask of integers defining closed seas rnf mappings 
     51   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csemp , mask_csgrpemp !: mask of integers defining closed seas empmr mappings 
     52 
    6553   !!---------------------------------------------------------------------- 
    6654   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7664      !! ** Purpose :   Closed sea domain initialization 
    7765      !! 
    78       !! ** Method  :   if a closed sea is located only in a model grid point 
    79       !!                just the thermodynamic processes are applied. 
    80       !! 
    81       !! ** Action  :   Read closea_mask* fields (if they exist) from domain_cfg file and infer 
    82       !!                number of closed seas from closea_mask field. 
    83       !!                closea_mask       : integer values defining closed seas (or groups of closed seas) 
    84       !!                closea_mask_rnf   : integer values defining mappings from closed seas or groups of 
    85       !!                                    closed seas to a runoff area for downwards flux only. 
    86       !!                closea_mask_empmr : integer values defining mappings from closed seas or groups of 
    87       !!                                    closed seas to a runoff area for net fluxes. 
    88       !! 
    89       !!                Python code to generate the closea_masks* fields from the old-style indices 
    90       !!                definitions is available at TOOLS/DOMAINcfg/make_closea_masks.py 
    91       !!---------------------------------------------------------------------- 
    92       INTEGER ::   inum    ! input file identifier 
    93       INTEGER ::   ierr    ! error code 
    94       INTEGER ::   id      ! netcdf variable ID 
    95  
    96       REAL(wp), DIMENSION(jpi,jpj) :: zdata_in ! temporary real array for input 
    97       !!---------------------------------------------------------------------- 
    98       ! 
     66      !! ** Action  :   Read mask_cs* fields (if needed) from domain_cfg file and infer 
     67      !!                number of closed seas for each case (glo, rnf, emp) from mask_cs* field. 
     68      !! 
     69      !! ** Output  :   mask_csglo and mask_csgrpglo  : integer values defining mappings from closed seas and associated groups to the open ocean for net fluxes. 
     70      !!                mask_csrnf and mask_csgrprnf  : integer values defining mappings from closed seas and associated groups to a runoff area for downwards flux only. 
     71      !!                mask_csemp and mask_csgrpemp  : integer values defining mappings from closed seas and associated groups to a runoff area for net fluxes. 
     72      !!---------------------------------------------------------------------- 
     73      INTEGER ::   ios     ! io status 
     74      !! 
     75      NAMELIST/namclo/ ln_maskcs, ln_mask_csundef, ln_clo_rnf 
     76      !!--------------------------------------------------------------------- 
     77      !! 
     78      READ  ( numnam_ref, namclo, IOSTAT = ios, ERR = 901 ) 
     79901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namclo in reference namelist' ) 
     80      READ  ( numnam_cfg, namclo, IOSTAT = ios, ERR = 902 ) 
     81902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namclo in configuration namelist' ) 
     82      IF(lwm) WRITE ( numond, namclo ) 
     83      !! 
    9984      IF(lwp) WRITE(numout,*) 
    10085      IF(lwp) WRITE(numout,*)'dom_clo : read in masks to define closed seas ' 
    10186      IF(lwp) WRITE(numout,*)'~~~~~~~' 
     87      IF(lwp) WRITE(numout,*) 
     88      !! 
     89      !! check option compatibility 
     90      IF( .NOT. ln_read_cfg ) THEN 
     91         CALL ctl_stop('Suppression of closed seas does not work with ln_read_cfg = .true. . Set ln_closea = .false. .') 
     92      ENDIF 
     93      !! 
     94      IF( (.NOT. ln_maskcs) .AND. ln_diurnal_only ) THEN 
     95         CALL ctl_stop('Special handling of freshwater fluxes over closed seas not compatible with ln_diurnal_only.') 
     96      END IF 
    10297      ! 
    10398      ! read the closed seas masks (if they exist) from domain_cfg file (if it exists) 
    10499      ! ------------------------------------------------------------------------------ 
    105100      ! 
    106       IF( ln_read_cfg) THEN 
    107          ! 
    108          CALL iom_open( cn_domcfg, inum ) 
    109          ! 
    110          id = iom_varid(inum, 'closea_mask', ldstop = .false.) 
    111          IF( id > 0 ) THEN  
    112             l_sbc_clo = .true. 
    113             ALLOCATE( closea_mask(jpi,jpj) , STAT=ierr ) 
    114             IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask array') 
    115             zdata_in(:,:) = 0.0 
    116             CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in ) 
    117             closea_mask(:,:) = NINT(zdata_in(:,:)) * tmask(:,:,1) 
    118             ! number of closed seas = global maximum value in closea_mask field 
    119             jncs = maxval(closea_mask(:,:)) 
    120             CALL mpp_max('closea', jncs) 
    121             IF( jncs > 0 ) THEN 
    122                IF( lwp ) WRITE(numout,*) 'Number of closed seas : ',jncs 
    123             ELSE 
    124                CALL ctl_stop( 'Problem with closea_mask field in domain_cfg file. Has no values > 0 so no closed seas defined.') 
    125             ENDIF 
    126          ELSE  
    127             IF( lwp ) WRITE(numout,*) 
    128             IF( lwp ) WRITE(numout,*) '   ==>>>   closea_mask field not found in domain_cfg file.' 
    129             IF( lwp ) WRITE(numout,*) '           No closed seas defined.' 
    130             IF( lwp ) WRITE(numout,*) 
    131             l_sbc_clo = .false. 
    132             jncs = 0  
    133          ENDIF 
    134  
    135          l_clo_rnf = .false. 
    136  
    137          IF( l_sbc_clo ) THEN ! No point reading in closea_mask_rnf or closea_mask_empmr fields if no closed seas defined. 
    138  
    139             id = iom_varid(inum, 'closea_mask_rnf', ldstop = .false.) 
    140             IF( id > 0 ) THEN  
    141                l_clo_rnf = .true.             
    142                ALLOCATE( closea_mask_rnf(jpi,jpj) , STAT=ierr ) 
    143                IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_rnf array') 
    144                CALL iom_get ( inum, jpdom_data, 'closea_mask_rnf', zdata_in ) 
    145                closea_mask_rnf(:,:) = NINT(zdata_in(:,:)) * tmask(:,:,1) 
    146                ! number of closed seas rnf mappings = global maximum in closea_mask_rnf field 
    147                jncsr = maxval(closea_mask_rnf(:,:)) 
    148                CALL mpp_max('closea', jncsr) 
    149                IF( jncsr > 0 ) THEN 
    150                   IF( lwp ) WRITE(numout,*) 'Number of closed seas rnf mappings : ',jncsr 
    151                ELSE 
    152                   CALL ctl_stop( 'Problem with closea_mask_rnf field in domain_cfg file. Has no values > 0 so no closed seas rnf mappings defined.') 
    153                ENDIF 
    154             ELSE  
    155                IF( lwp ) WRITE(numout,*) 'closea_mask_rnf field not found in domain_cfg file. No closed seas rnf mappings defined.' 
    156                jncsr = 0 
    157             ENDIF 
    158   
    159             id = iom_varid(inum, 'closea_mask_empmr', ldstop = .false.) 
    160             IF( id > 0 ) THEN  
    161                l_clo_rnf = .true.             
    162                ALLOCATE( closea_mask_empmr(jpi,jpj) , STAT=ierr ) 
    163                IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_empmr array') 
    164                CALL iom_get ( inum, jpdom_data, 'closea_mask_empmr', zdata_in ) 
    165                closea_mask_empmr(:,:) = NINT(zdata_in(:,:)) * tmask(:,:,1) 
    166                ! number of closed seas empmr mappings = global maximum value in closea_mask_empmr field 
    167                jncse = maxval(closea_mask_empmr(:,:)) 
    168                CALL mpp_max('closea', jncse) 
    169                IF( jncse > 0 ) THEN  
    170                   IF( lwp ) WRITE(numout,*) 'Number of closed seas empmr mappings : ',jncse 
    171                ELSE 
    172                   CALL ctl_stop( 'Problem with closea_mask_empmr field in domain_cfg file. Has no values > 0 so no closed seas empmr mappings defined.') 
    173                ENDIF 
    174             ELSE  
    175                IF( lwp ) WRITE(numout,*) 'closea_mask_empmr field not found in domain_cfg file. No closed seas empmr mappings defined.' 
    176                jncse = 0 
    177             ENDIF 
    178  
    179          ENDIF ! l_sbc_clo 
    180          ! 
    181          CALL iom_close( inum ) 
    182          ! 
    183       ELSE ! ln_read_cfg = .false. so no domain_cfg file 
    184          IF( lwp ) WRITE(numout,*) 'No domain_cfg file so no closed seas defined.' 
    185          l_sbc_clo = .false. 
    186          l_clo_rnf = .false. 
    187       ENDIF 
    188       ! 
     101      ! load mask of open sea 
     102      CALL alloc_csmask( mask_opnsea ) 
     103      CALL read_csmask( cn_domcfg, 'mask_opensea' , mask_opnsea  ) 
     104      ! 
     105      IF ( ln_maskcs ) THEN 
     106         ! closed sea are masked 
     107         IF(lwp) WRITE(numout,*)'          ln_maskcs = T : all closed seas are masked' 
     108         IF(lwp) WRITE(numout,*) 
     109         ! no special treatment of closed sea 
     110         ! no redistribution of emp unbalance over closed sea into river mouth/open ocean 
     111         l_sbc_clo = .false. ; l_clo_rnf = .false. 
     112      ELSE 
     113         ! redistribution of emp unbalance over closed sea into river mouth/open ocean 
     114         IF(lwp) WRITE(numout,*)'          ln_maskcs = F : net emp is corrected over defined closed seas' 
     115         ! 
     116         l_sbc_clo = .true. 
     117         ! 
     118         ! river mouth from lakes added to rnf mask for special treatment 
     119         IF ( ln_clo_rnf) l_clo_rnf = .true. 
     120         ! 
     121         IF ( ln_mask_csundef) THEN 
     122            ! closed sea not defined (ie not in the domcfg namelist used to build the domcfg.nc file) are masked  
     123            IF(lwp) WRITE(numout,*)'          ln_mask_csundef = T : all undefined closed seas are masked' 
     124            ! 
     125            CALL alloc_csmask( mask_csundef ) 
     126            CALL read_csmask( cn_domcfg, 'mask_csundef', mask_csundef ) 
     127            ! revert the mask for masking of undefined closed seas in domzgr  
     128            ! (0 over the undefined closed sea and 1 elsewhere) 
     129            mask_csundef(:,:) = 1 - mask_csundef(:,:) 
     130         END IF 
     131         IF(lwp) WRITE(numout,*) 
     132         ! 
     133         ! allocate source mask for each cases 
     134         CALL alloc_csmask( mask_csglo ) 
     135         CALL alloc_csmask( mask_csrnf ) 
     136         CALL alloc_csmask( mask_csemp ) 
     137         ! 
     138         ! load source mask of cs for each cases 
     139         CALL read_csmask( cn_domcfg, 'mask_csglo', mask_csglo ) 
     140         CALL read_csmask( cn_domcfg, 'mask_csrnf', mask_csrnf ) 
     141         CALL read_csmask( cn_domcfg, 'mask_csemp', mask_csemp ) 
     142         ! 
     143         ! compute number of cs for each cases 
     144         ncsg = MAXVAL( mask_csglo(:,:) ) ; CALL mpp_max( 'closea', ncsg ) 
     145         ncsr = MAXVAL( mask_csrnf(:,:) ) ; CALL mpp_max( 'closea', ncsr ) 
     146         ncse = MAXVAL( mask_csemp(:,:) ) ; CALL mpp_max( 'closea', ncse ) 
     147         ! 
     148         ! allocate closed sea group masks  
     149         !(used to defined the target area in case multiple lakes have the same river mouth (great lakes for example)) 
     150         CALL alloc_csmask( mask_csgrpglo ) 
     151         CALL alloc_csmask( mask_csgrprnf ) 
     152         CALL alloc_csmask( mask_csgrpemp ) 
     153 
     154         ! load mask of cs group for each cases 
     155         CALL read_csmask( cn_domcfg, 'mask_csgrpglo', mask_csgrpglo ) 
     156         CALL read_csmask( cn_domcfg, 'mask_csgrprnf', mask_csgrprnf ) 
     157         CALL read_csmask( cn_domcfg, 'mask_csgrpemp', mask_csgrpemp ) 
     158         ! 
     159      END IF 
    189160   END SUBROUTINE dom_clo 
    190161 
    191  
    192    SUBROUTINE sbc_clo( kt ) 
    193       !!--------------------------------------------------------------------- 
    194       !!                  ***  ROUTINE sbc_clo  *** 
    195       !!                     
    196       !! ** Purpose :   Special handling of closed seas 
    197       !! 
    198       !! ** Method  :   Water flux is forced to zero over closed sea 
    199       !!      Excess is shared between remaining ocean, or 
    200       !!      put as run-off in open ocean. 
    201       !! 
    202       !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt 
    203       !!---------------------------------------------------------------------- 
    204       INTEGER         , INTENT(in   ) ::   kt       ! ocean model time step 
    205       ! 
    206       INTEGER             ::   ierr 
    207       INTEGER             ::   jc, jcr, jce   ! dummy loop indices 
    208       REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon 
    209       REAL(wp)            ::   zfwf_total, zcoef, zcoef1         !  
    210       REAL(wp), DIMENSION(jncs)    ::   zfwf      !: 
    211       REAL(wp), DIMENSION(jncsr+1) ::   zfwfr     !: freshwater fluxes over closed seas 
    212       REAL(wp), DIMENSION(jncse+1) ::   zfwfe     !:  
    213       REAL(wp), DIMENSION(jpi,jpj) ::   ztmp2d   ! 2D workspace 
    214       !!---------------------------------------------------------------------- 
    215       ! 
    216       IF( ln_timing )  CALL timing_start('sbc_clo') 
    217       ! 
    218       !                                                   !------------------!  
    219       IF( kt == nit000 ) THEN                             !  Initialisation  ! 
    220          !                                                !------------------! 
    221          IF(lwp) WRITE(numout,*) 
    222          IF(lwp) WRITE(numout,*)'sbc_clo : closed seas ' 
    223          IF(lwp) WRITE(numout,*)'~~~~~~~' 
    224  
    225          ALLOCATE( surf(jncs+1) , STAT=ierr ) 
    226          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array') 
    227          surf(:) = 0.e0_wp 
    228          ! 
    229          ! jncsr can be zero so add 1 to avoid allocating zero-length array 
    230          ALLOCATE( surfr(jncsr+1) , STAT=ierr ) 
    231          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfr array') 
    232          surfr(:) = 0.e0_wp 
    233          ! 
    234          ! jncse can be zero so add 1 to avoid allocating zero-length array 
    235          ALLOCATE( surfe(jncse+1) , STAT=ierr ) 
    236          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfe array') 
    237          surfe(:) = 0.e0_wp 
    238          ! 
    239          surf(jncs+1) = glob_sum( 'closea', e1e2t(:,:) )   ! surface of the global ocean 
    240          ! 
    241          !                                        ! surface areas of closed seas  
    242          DO jc = 1, jncs 
    243             ztmp2d(:,:) = 0.e0_wp 
    244             WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 
    245             surf(jc) = glob_sum( 'closea', ztmp2d(:,:) ) 
    246          END DO 
    247          ! 
    248          ! jncs+1 : surface area of global ocean, closed seas excluded 
    249          surf(jncs+1) = surf(jncs+1) - SUM(surf(1:jncs)) 
    250          ! 
    251          !                                        ! surface areas of rnf target areas 
    252          IF( jncsr > 0 ) THEN 
    253             DO jcr = 1, jncsr 
    254                ztmp2d(:,:) = 0.e0_wp 
    255                WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 
    256                surfr(jcr) = glob_sum( 'closea', ztmp2d(:,:) ) 
    257             END DO 
    258          ENDIF 
    259          ! 
    260          !                                        ! surface areas of empmr target areas 
    261          IF( jncse > 0 ) THEN 
    262             DO jce = 1, jncse 
    263                ztmp2d(:,:) = 0.e0_wp 
    264                WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:) 
    265                surfe(jce) = glob_sum( 'closea', ztmp2d(:,:) ) 
    266             END DO 
    267          ENDIF 
    268          ! 
    269          IF(lwp) WRITE(numout,*)'     Closed sea surface areas (km2)' 
    270          DO jc = 1, jncs 
    271             IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jc, surf(jc) * 1.0e-6 
    272          END DO 
    273          IF(lwp) WRITE(numout,FMT='(A,ES12.2)') 'Global surface area excluding closed seas (km2): ', surf(jncs+1) * 1.0e-6 
    274          ! 
    275          IF(jncsr > 0) THEN 
    276             IF(lwp) WRITE(numout,*)'     Closed sea target rnf surface areas (km2)' 
    277             DO jcr = 1, jncsr 
    278                IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jcr, surfr(jcr) * 1.0e-6 
    279             END DO 
    280          ENDIF 
    281          ! 
    282          IF(jncse > 0) THEN 
    283             IF(lwp) WRITE(numout,*)'     Closed sea target empmr surface areas (km2)' 
    284             DO jce = 1, jncse 
    285                IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jce, surfe(jce) * 1.0e-6 
    286             END DO 
    287          ENDIF 
    288       ENDIF 
    289       ! 
    290       !                                                      !--------------------! 
    291       !                                                      !  update emp        ! 
    292       !                                                      !--------------------! 
    293  
    294       zfwf_total = 0._wp 
    295  
    296       ! 
    297       ! 1. Work out total freshwater fluxes over closed seas from EMP - RNF. 
    298       ! 
    299       zfwf(:) = 0.e0_wp            
    300       DO jc = 1, jncs 
    301          ztmp2d(:,:) = 0.e0_wp 
    302          WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 
    303          zfwf(jc) = glob_sum( 'closea', ztmp2d(:,:) ) 
    304       END DO 
    305       zfwf_total = SUM(zfwf) 
    306  
    307       zfwfr(:) = 0.e0_wp            
    308       IF( jncsr > 0 ) THEN 
    309          ! 
    310          ! 2. Work out total FW fluxes over rnf source areas and add to rnf target areas.  
    311          !    Where zfwf is negative add flux at specified runoff points and subtract from fluxes for global redistribution. 
    312          !    Where positive leave in global redistribution total. 
    313          ! 
    314          DO jcr = 1, jncsr 
    315             ! 
    316             ztmp2d(:,:) = 0.e0_wp 
    317             WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 
    318             zfwfr(jcr) = glob_sum( 'closea', ztmp2d(:,:) ) 
    319             ! 
    320             ! The following if avoids the redistribution of the round off 
    321             IF ( ABS(zfwfr(jcr) / surf(jncs+1) ) > rsmall) THEN 
    322                ! 
    323                ! Add residuals to target runoff points if negative and subtract from total to be added globally 
    324                IF( zfwfr(jcr) < 0.0 ) THEN  
    325                   zfwf_total = zfwf_total - zfwfr(jcr) 
    326                   zcoef    = zfwfr(jcr) / surfr(jcr) 
    327                   zcoef1   = rcp * zcoef 
    328                   WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0) 
    329                      emp(:,:) = emp(:,:) + zcoef 
    330                      qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    331                   ENDWHERE 
    332                ENDIF 
    333                ! 
    334             ENDIF 
    335          END DO 
    336       ENDIF  ! jncsr > 0     
    337       ! 
    338       zfwfe(:) = 0.e0_wp            
    339       IF( jncse > 0 ) THEN 
    340          ! 
    341          ! 3. Work out total fluxes over empmr source areas and add to empmr target areas.  
    342          ! 
    343          DO jce = 1, jncse 
    344             ! 
    345             ztmp2d(:,:) = 0.e0_wp 
    346             WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) 
    347             zfwfe(jce) = glob_sum( 'closea', ztmp2d(:,:) ) 
    348             ! 
    349             ! The following if avoids the redistribution of the round off 
    350             IF ( ABS( zfwfe(jce) / surf(jncs+1) ) > rsmall ) THEN 
    351                ! 
    352                ! Add residuals to runoff points and subtract from total to be added globally 
    353                zfwf_total = zfwf_total - zfwfe(jce) 
    354                zcoef    = zfwfe(jce) / surfe(jce) 
    355                zcoef1   = rcp * zcoef 
    356                WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0.0) 
    357                   emp(:,:) = emp(:,:) + zcoef 
    358                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    359                ENDWHERE 
    360                ! 
    361             ENDIF 
    362          END DO 
    363       ENDIF ! jncse > 0     
    364  
    365       ! 
    366       ! 4. Spread residual flux over global ocean.  
    367       ! 
    368       ! The following if avoids the redistribution of the round off 
    369       IF ( ABS(zfwf_total / surf(jncs+1) ) > rsmall) THEN 
    370          zcoef    = zfwf_total / surf(jncs+1) 
    371          zcoef1   = rcp * zcoef 
    372          WHERE( closea_mask(:,:) == 0 ) 
    373             emp(:,:) = emp(:,:) + zcoef 
    374             qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    375          ENDWHERE 
    376       ENDIF 
    377  
    378       ! 
    379       ! 5. Subtract area means from emp (and qns) over closed seas to give zero mean FW flux over each sea. 
    380       ! 
    381       DO jc = 1, jncs 
    382          ! The following if avoids the redistribution of the round off 
    383          IF ( ABS(zfwf(jc) / surf(jncs+1) ) > rsmall) THEN 
    384             ! 
    385             ! Subtract residuals from fluxes over closed sea 
    386             zcoef    = zfwf(jc) / surf(jc) 
    387             zcoef1   = rcp * zcoef 
    388             WHERE( closea_mask(:,:) == jc ) 
    389                emp(:,:) = emp(:,:) - zcoef 
    390                qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:) 
    391             ENDWHERE 
    392             ! 
    393          ENDIF 
    394       END DO 
    395       ! 
    396       emp (:,:) = emp (:,:) * tmask(:,:,1) 
    397       ! 
    398       CALL lbc_lnk( 'closea', emp , 'T', 1._wp ) 
    399       ! 
    400    END SUBROUTINE sbc_clo 
    401  
    402162   SUBROUTINE clo_rnf( p_rnfmsk ) 
    403163      !!--------------------------------------------------------------------- 
    404       !!                  ***  ROUTINE sbc_rnf  *** 
     164      !!                  ***  ROUTINE clo_rnf  *** 
    405165      !!                     
    406166      !! ** Purpose :   allow the treatment of closed sea outflow grid-points 
     
    412172      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow) 
    413173      !!---------------------------------------------------------------------- 
     174      !! subroutine parameter 
    414175      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array) 
    415       !!---------------------------------------------------------------------- 
    416       ! 
    417       IF( jncsr > 0 ) THEN 
    418          WHERE( closea_mask_rnf(:,:) > 0 .and. closea_mask(:,:) == 0 ) 
    419             p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp ) 
    420          ENDWHERE 
    421       ENDIF 
    422       ! 
    423       IF( jncse > 0 ) THEN 
    424          WHERE( closea_mask_empmr(:,:) > 0 .and. closea_mask(:,:) == 0 ) 
    425             p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp ) 
    426          ENDWHERE 
    427       ENDIF 
     176      !! 
     177      !! local variables 
     178      REAL(wp), DIMENSION(jpi,jpj) :: zmsk 
     179      !!---------------------------------------------------------------------- 
     180      ! 
     181      ! zmsk > 0 where cs river mouth defined (case rnf and emp) 
     182      zmsk(:,:) = ( mask_csgrprnf (:,:) + mask_csgrpemp(:,:) ) * mask_opnsea(:,:) 
     183      WHERE( zmsk(:,:) > 0 ) 
     184         p_rnfmsk(:,:) = 1.0_wp 
     185      END WHERE 
    428186      ! 
    429187   END SUBROUTINE clo_rnf 
    430     
    431188       
    432    SUBROUTINE clo_bat( k_top, k_bot ) 
    433       !!--------------------------------------------------------------------- 
    434       !!                  ***  ROUTINE clo_bat  *** 
     189   SUBROUTINE clo_msk( k_top, k_bot, k_mask, cd_prt ) 
     190      !!--------------------------------------------------------------------- 
     191      !!                  ***  ROUTINE clo_msk  *** 
    435192      !!                     
    436193      !! ** Purpose :   Suppress closed sea from the domain 
    437194      !! 
    438       !! ** Method  :   Read in closea_mask field (if it exists) from domain_cfg file. 
    439       !!                Where closea_mask > 0 set first and last ocean level to 0 
     195      !! ** Method  :   Where closea_mask > 0 set first and last ocean level to 0 
    440196      !!                (As currently coded you can't define a closea_mask field in  
    441197      !!                usr_def_zgr). 
     
    443199      !! ** Action  :   set k_top=0 and k_bot=0 over closed seas 
    444200      !!---------------------------------------------------------------------- 
     201      !! subroutine parameter 
    445202      INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices 
    446       INTEGER                           :: inum, id 
    447       INTEGER,  DIMENSION(jpi,jpj) :: closea_mask ! closea_mask field 
    448       REAL(wp), DIMENSION(jpi,jpj) :: zdata_in ! temporary real array for input 
    449       !!---------------------------------------------------------------------- 
    450       ! 
    451       IF(lwp) THEN                     ! Control print 
     203      INTEGER, DIMENSION(:,:), INTENT(in   ) ::   k_mask         ! mask used to mask ktop and k_bot 
     204      CHARACTER(LEN=*),        INTENT(in   ) ::   cd_prt         ! text for control print 
     205      !! 
     206      !! local variables 
     207      !!---------------------------------------------------------------------- 
     208      !! 
     209      IF ( lwp ) THEN 
    452210         WRITE(numout,*) 
    453          WRITE(numout,*) 'clo_bat : suppression of closed seas' 
     211         WRITE(numout,*) 'clo_msk : Suppression closed seas based on ',TRIM(cd_prt),' field.' 
    454212         WRITE(numout,*) '~~~~~~~' 
     213         WRITE(numout,*) 
    455214      ENDIF 
    456       ! 
    457       IF( ln_read_cfg ) THEN 
    458          ! 
    459          CALL iom_open( cn_domcfg, inum ) 
    460          ! 
    461          id = iom_varid(inum, 'closea_mask', ldstop = .false.)       
    462          IF( id > 0 ) THEN 
    463             IF( lwp ) WRITE(numout,*) 'Suppressing closed seas in bathymetry based on closea_mask field,' 
    464             CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in ) 
    465             closea_mask(:,:) = NINT(zdata_in(:,:)) 
    466             WHERE( closea_mask(:,:) > 0 ) 
    467                k_top(:,:) = 0    
    468                k_bot(:,:) = 0    
    469             ENDWHERE 
    470          ELSE 
    471             IF( lwp ) WRITE(numout,*) 'No closea_mask field found in domain_cfg file. No suppression of closed seas.' 
    472          ENDIF 
    473          ! 
    474          CALL iom_close(inum) 
    475          ! 
    476       ELSE 
    477          IF( lwp ) WRITE(numout,*) 'No domain_cfg file => no suppression of closed seas.' 
    478       ENDIF 
    479       ! 
    480       ! Initialise l_sbc_clo and l_clo_rnf for this case (ln_closea=.false.) 
    481       l_sbc_clo = .false. 
    482       l_clo_rnf = .false. 
    483       ! 
    484    END SUBROUTINE clo_bat 
    485  
    486    !!====================================================================== 
     215      !! 
     216      k_top(:,:) = k_top(:,:) * k_mask(:,:) 
     217      k_bot(:,:) = k_bot(:,:) * k_mask(:,:) 
     218      !! 
     219   END SUBROUTINE clo_msk 
     220 
     221   SUBROUTINE read_csmask(cd_file, cd_var, k_mskout) 
     222      !!--------------------------------------------------------------------- 
     223      !!                  ***  ROUTINE read_csmask  *** 
     224      !!                     
     225      !! ** Purpose : read mask in cd_filec file 
     226      !!---------------------------------------------------------------------- 
     227      ! subroutine parameter 
     228      CHARACTER(LEN=256),          INTENT(in   ) :: cd_file    ! netcdf file     name 
     229      CHARACTER(LEN= * ),          INTENT(in   ) :: cd_var     ! netcdf variable name 
     230      INTEGER, DIMENSION(:,:), INTENT(  out) :: k_mskout            ! output mask variable 
     231      ! 
     232      ! local variables 
     233      INTEGER :: ics                       ! netcdf id 
     234      REAL(wp), DIMENSION(jpi,jpj) :: zdta ! netcdf data 
     235      !!---------------------------------------------------------------------- 
     236      ! 
     237      CALL iom_open ( cd_file, ics ) 
     238      CALL iom_get  ( ics, jpdom_data, TRIM(cd_var), zdta ) 
     239      CALL iom_close( ics ) 
     240      k_mskout(:,:) = NINT(zdta(:,:)) 
     241      ! 
     242   END SUBROUTINE read_csmask 
     243 
     244   SUBROUTINE alloc_csmask( kmask ) 
     245      !!--------------------------------------------------------------------- 
     246      !!                  ***  ROUTINE alloc_csmask  *** 
     247      !!                     
     248      !! ** Purpose : allocated cs mask 
     249      !!---------------------------------------------------------------------- 
     250      ! subroutine parameter 
     251      INTEGER, ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: kmask 
     252      ! 
     253      ! local variables 
     254      INTEGER :: ierr 
     255      !!---------------------------------------------------------------------- 
     256      ! 
     257      ALLOCATE( kmask(jpi,jpj) , STAT=ierr ) 
     258      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'alloc_csmask: failed to allocate surf array') 
     259      ! 
     260   END SUBROUTINE 
     261 
    487262END MODULE closea 
    488  
  • NEMO/trunk/src/OCE/DOM/daymod.F90

    r12158 r12377  
    5858      !! 
    5959      !! ** Action  : - nyear        : current year 
    60       !!              - nmonth       : current month of the year nyear 
    61       !!              - nday         : current day of the month nmonth 
    62       !!              - nday_year    : current day of the year nyear 
    63       !!              - nsec_year    : current time step counted in second since 00h jan 1st of the current year 
    64       !!              - nsec_month   : current time step counted in second since 00h 1st day of the current month 
    65       !!              - nsec_day     : current time step counted in second since 00h of the current day 
    66       !!              - nsec1jan000  : second since Jan. 1st 00h of nit000 year and Jan. 1st 00h of the current year 
    67       !!              - nmonth_len, nyear_len, nmonth_half, nmonth_end through day_mth 
    68       !!---------------------------------------------------------------------- 
    69       INTEGER  ::   inbday, idweek   ! local integers 
     60      !!              - nmonth       : current month of the current nyear 
     61      !!              - nday         : current   day of the current nmonth 
     62      !!              - nday_year    : current   day of the current nyear 
     63      !!              - nsec_year    : seconds between 00h jan 1st of the current  year and half of the current time step 
     64      !!              - nsec_month   : seconds between 00h 1st day of the current month and half of the current time step 
     65      !!              - nsec_monday  : seconds between 00h         of the   last Monday and half of the current time step 
     66      !!              - nsec_day     : seconds between 00h         of the current   day and half of the current time step 
     67      !!              - nsec1jan000  : seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h of the current year 
     68      !!              - nmonth_len, nyear_len, nmonth_beg through day_mth 
     69      !!---------------------------------------------------------------------- 
     70      INTEGER  ::   inbday, imonday, isecrst   ! local integers 
    7071      REAL(wp) ::   zjul             ! local scalar 
    7172      !!---------------------------------------------------------------------- 
     
    7677            &           'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 
    7778      ENDIF 
    78       nsecd   = NINT( rday      ) 
     79      nsecd   = NINT(       rday ) 
    7980      nsecd05 = NINT( 0.5 * rday ) 
    8081      ndt     = NINT(       rdt  ) 
     
    9091      nhour   =   nn_time0 / 100 
    9192      nminute = ( nn_time0 - nhour * 100 ) 
    92  
    93       CALL ymds2ju( nyear, nmonth, nday, nhour*3600._wp+nminute*60._wp, fjulday )   
     93      isecrst = ( nhour * NINT(rhhmm) + nminute ) * NINT(rmmss) 
     94 
     95      CALL ymds2ju( nyear, nmonth, nday, REAL(isecrst,wp), fjulday )   
    9496      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error 
    95       IF( nhour*3600 + nminute*60 - ndt05 .lt. 0 ) fjulday = fjulday + 1.                    ! move back to the day at nit000 (and not at nit000 - 1) 
    96        
     97      IF( nhour*NINT(rhhmm*rmmss) + nminute*NINT(rmmss) - ndt05 .LT. 0 ) fjulday = fjulday+1.       ! move back to the day at nit000 (and not at nit000 - 1) 
     98 
    9799      nsec1jan000 = 0 
    98100      CALL day_mth 
     
    112114      nday_year = nday + SUM( nmonth_len(1:nmonth - 1) ) 
    113115 
    114       !compute number of days between last monday and today 
    115       CALL ymds2ju( 1900, 01, 01, 0.0, zjul )  ! compute julian day value of 01.01.1900 (our reference that was a Monday) 
    116       inbday = FLOOR(fjulday - zjul)            ! compute nb day between  01.01.1900 and start of current day 
    117       idweek = MOD(inbday, 7)                  ! compute nb day between last monday and current day 
    118       IF (idweek .lt. 0) idweek=idweek+7       ! Avoid negative values for dates before 01.01.1900 
     116      !compute number of days between last Monday and today 
     117      CALL ymds2ju( 1900, 01, 01, 0.0, zjul )     ! compute julian day value of 01.01.1900 (our reference that was a Monday) 
     118      inbday = FLOOR(fjulday - zjul)              ! compute nb day between  01.01.1900 and start of current day 
     119      imonday = MOD(inbday, 7)                    ! compute nb day between last monday and current day 
     120      IF (imonday .LT. 0) imonday = imonday + 7   ! Avoid negative values for dates before 01.01.1900 
    119121 
    120122      ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step 
    121       IF (nhour*3600+nminute*60-ndt05 .gt. 0) THEN 
     123      IF( isecrst - ndt05 .GT. 0 ) THEN 
    122124         ! 1 timestep before current middle of first time step is still the same day 
    123          nsec_year  = (nday_year-1) * nsecd + nhour*3600+nminute*60 - ndt05  
    124          nsec_month = (nday-1)      * nsecd + nhour*3600+nminute*60 - ndt05     
     125         nsec_year  = (nday_year-1) * nsecd + isecrst - ndt05  
     126         nsec_month = (nday-1)      * nsecd + isecrst - ndt05     
    125127      ELSE 
    126128         ! 1 time step before the middle of the first time step is the previous day  
    127          nsec_year  = nday_year * nsecd + nhour*3600+nminute*60 - ndt05  
    128          nsec_month = nday      * nsecd + nhour*3600+nminute*60 - ndt05    
    129       ENDIF 
    130       nsec_week  = idweek    * nsecd + nhour*3600+nminute*60 - ndt05 
    131       nsec_day   =             nhour*3600+nminute*60 - ndt05  
    132       IF( nsec_day .lt. 0 ) nsec_day = nsec_day + nsecd 
    133       IF( nsec_week .lt. 0 ) nsec_week = nsec_week + nsecd*7 
     129         nsec_year  = nday_year     * nsecd + isecrst - ndt05  
     130         nsec_month = nday          * nsecd + isecrst - ndt05    
     131      ENDIF 
     132      nsec_monday   = imonday       * nsecd + isecrst - ndt05 
     133      nsec_day      =                         isecrst - ndt05  
     134      IF( nsec_day    .LT. 0 ) nsec_day    = nsec_day    + nsecd 
     135      IF( nsec_monday .LT. 0 ) nsec_monday = nsec_monday + nsecd*7 
    134136 
    135137      ! control print 
    136138      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8,a,i8,a,i8)')   & 
    137139           &                   ' =======>> 1/2 time step before the start of the run DATE Y/M/D = ',   & 
    138            &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week, '  & 
     140           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_monday:', nsec_monday, '  & 
    139141           &                   nsec_month:', nsec_month , '  nsec_year:' , nsec_year 
    140142 
     143      nsec000_1jan000 = nsec1jan000 + nsec_year + ndt05 
     144      nsecend_1jan000 = nsec000_1jan000 + ndt * ( nitend - nit000 + 1 ) 
     145       
    141146      ! Up to now, calendar parameters are related to the end of previous run (nit000-1) 
    142147      ! call day to set the calendar parameters at the begining of the current simulaton. needed by iom_init 
     
    160165      !! ** Purpose :   calendar values related to the months 
    161166      !! 
    162       !! ** Action  : - nmonth_len    : length in days of the months of the current year 
    163       !!              - nyear_len     : length in days of the previous/current year 
    164       !!              - nmonth_half   : second since the beginning of the year and the halft of the months 
    165       !!              - nmonth_end    : second since the beginning of the year and the end of the months 
    166       !!---------------------------------------------------------------------- 
    167       INTEGER  ::   jm               ! dummy loop indice 
     167      !! ** Action  : - nyear_len     : length in days of the previous/current year 
     168      !!              - nmonth_len    : length in days of the months of the current year 
     169      !!              - nmonth_half   : second since the beginning of the current year and the halft of the months 
     170      !!              - nmonth_end    : second since the beginning of the current year and the end of the months 
     171      !!---------------------------------------------------------------------- 
     172      INTEGER  ::   jm ,jy                   ! dummy loop indice 
     173      INTEGER, DIMENSION(12) ::   idaymt     ! length in days of the 12 months for non-leap year 
    168174      !!---------------------------------------------------------------------- 
    169175 
    170176      ! length of the month of the current year (from nleapy, read in namelist) 
    171177      IF ( nleapy < 2 ) THEN 
    172          nmonth_len(:) = (/ 31, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31 /) 
     178         ! default values 
     179         idaymt(1:12) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) 
     180         nmonth_len(-11: 25) = (/ idaymt(1:12), idaymt(1:12), idaymt(1:12), idaymt(1) /) 
    173181         nyear_len(:) = 365 
     182         ! 
    174183         IF ( nleapy == 1 ) THEN   ! we are using calandar with leap years 
    175             IF ( MOD(nyear-1, 4) == 0 .AND. ( MOD(nyear-1, 400) == 0 .OR. MOD(nyear-1, 100) /= 0 ) ) THEN 
    176                nyear_len(0)  = 366 
    177             ENDIF 
    178             IF ( MOD(nyear  , 4) == 0 .AND. ( MOD(nyear  , 400) == 0 .OR. MOD(nyear  , 100) /= 0 ) ) THEN 
    179                nmonth_len(2) = 29 
    180                nyear_len(1)  = 366 
    181             ENDIF 
    182             IF ( MOD(nyear+1, 4) == 0 .AND. ( MOD(nyear+1, 400) == 0 .OR. MOD(nyear+1, 100) /= 0 ) ) THEN 
    183                nyear_len(2)  = 366 
    184             ENDIF 
     184            DO jy = -1,1 
     185               IF ( MOD(nyear+jy, 4) == 0 .AND. ( MOD(nyear+jy, 400) == 0 .OR. MOD(nyear+jy, 100) /= 0 ) ) THEN 
     186                  nmonth_len(2 + 12*jy) = 29 
     187                  nyear_len( 1 +    jy) = 366 
     188               ENDIF 
     189            ENDDO 
    185190         ENDIF 
    186191      ELSE 
     
    189194      ENDIF 
    190195 
    191       ! half month in second since the begining of the year: 
    192196      ! time since Jan 1st   0     1     2    ...    11    12    13 
    193197      !          ---------*--|--*--|--*--| ... |--*--|--*--|--*--|-------------------------------------- 
    194198      !                 <---> <---> <--->  ...  <---> <---> <---> 
    195199      ! month number      0     1     2    ...    11    12    13 
    196       ! 
    197       ! nmonth_half(jm) = rday * REAL( 0.5 * nmonth_len(jm) + SUM(nmonth_len(1:jm-1)) ) 
    198       nmonth_half(0) = - nsecd05 * nmonth_len(0) 
    199       DO jm = 1, 13 
    200          nmonth_half(jm) = nmonth_half(jm-1) + nsecd05 * ( nmonth_len(jm-1) + nmonth_len(jm) ) 
     200      nmonth_beg(1) = 0 
     201      DO jm = 2, 25 
     202         nmonth_beg(jm) = nmonth_beg(jm-1) + nsecd * nmonth_len(jm-1) 
    201203      END DO 
    202  
    203       nmonth_end(0) = 0 
    204       DO jm = 1, 13 
    205          nmonth_end(jm) = nmonth_end(jm-1) + nsecd * nmonth_len(jm) 
     204      DO jm = 0,-11,-1 
     205         nmonth_beg(jm) = nmonth_beg(jm+1) - nsecd * nmonth_len(jm) 
    206206      END DO 
    207207      ! 
     
    235235      zprec = 0.1 / rday 
    236236      !                                                 ! New time-step 
    237       nsec_year  = nsec_year  + ndt 
    238       nsec_month = nsec_month + ndt 
    239       nsec_week  = nsec_week  + ndt 
     237      nsec_year    = nsec_year    + ndt 
     238      nsec_month   = nsec_month  + ndt 
     239      nsec_monday  = nsec_monday  + ndt 
    240240      nsec_day   = nsec_day   + ndt 
    241241      adatrj  = adatrj  + rdt / rday 
     
    272272              &   '      New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year 
    273273         IF(lwp) WRITE(numout,'(a,i8,a,i7,a,i5)') '         nsec_year = ', nsec_year,   & 
    274               &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day, '   nsec_week = ', nsec_week 
    275       ENDIF 
    276  
    277       IF( nsec_week > 7*nsecd )   nsec_week = ndt05     ! New week 
    278  
    279       IF(ln_ctl) THEN 
     274              &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day, '   nsec_monday = ', nsec_monday 
     275      ENDIF 
     276 
     277      IF( nsec_monday > 7*nsecd )   nsec_monday = ndt05     ! New week 
     278 
     279      IF(sn_cfctl%l_prtctl) THEN 
    280280         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear 
    281281         CALL prt_ctl_info(charout) 
     
    319319      ! 
    320320      REAL(wp) ::   zkt, zndastp, zdayfrac, ksecs, ktime 
    321       INTEGER  ::   ihour, iminute 
     321      INTEGER  ::   ihour, iminute, isecond 
    322322      !!---------------------------------------------------------------------- 
    323323 
     
    349349               CALL iom_get( numror, 'adatrj', adatrj , ldxios = lrxios ) 
    350350          CALL iom_get( numror, 'ntime' , ktime  , ldxios = lrxios ) 
    351           nn_time0=INT(ktime) 
     351               nn_time0 = NINT(ktime) 
    352352               ! calculate start time in hours and minutes 
    353           zdayfrac=adatrj-INT(adatrj) 
    354           ksecs = NINT(zdayfrac*86400)        ! Nearest second to catch rounding errors in adatrj          
    355           ihour = INT(ksecs/3600) 
    356           iminute = ksecs/60-ihour*60 
     353               zdayfrac = adatrj - REAL(INT(adatrj), wp) 
     354          ksecs = NINT(zdayfrac * rday)          ! Nearest second to catch rounding errors in adatrj          
     355               ihour = ksecs / NINT( rhhmm*rmmss ) 
     356          iminute = ksecs / NINT(rmmss) - ihour*NINT(rhhmm) 
    357357            
    358358               ! Add to nn_time0 
    359359               nhour   =   nn_time0 / 100 
    360360               nminute = ( nn_time0 - nhour * 100 ) 
    361           nminute=nminute+iminute 
     361          nminute = nminute + iminute 
    362362           
    363           IF( nminute >= 60 ) THEN 
    364              nminute=nminute-60 
    365         nhour=nhour+1 
     363               IF( nminute >= NINT(rhhmm) ) THEN 
     364             nminute = nminute - NINT(rhhmm) 
     365        nhour = nhour+1 
    366366          ENDIF 
    367367          nhour=nhour+ihour 
    368           IF( nhour >= 24 ) THEN 
    369         nhour=nhour-24 
    370              adatrj=adatrj+1 
     368          IF( nhour >= NINT(rjjhh) ) THEN 
     369        nhour = nhour - NINT(rjjhh) 
     370             adatrj = adatrj + 1. 
    371371          ENDIF           
    372372          nn_time0 = nhour * 100 + nminute 
    373           adatrj = INT(adatrj)                    ! adatrj set to integer as nn_time0 updated           
     373               adatrj = REAL(INT(adatrj), wp)                    ! adatrj set to integer as nn_time0 updated           
    374374            ELSE 
    375375               ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day) 
     
    377377               nhour   =   nn_time0 / 100 
    378378               nminute = ( nn_time0 - nhour * 100 ) 
    379                IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0) 
     379               isecond = ( nhour * NINT(rhhmm) + nminute ) * NINT(rmmss) 
     380               IF( isecond - ndt05 .lt. 0 )   ndastp = ndastp - 1      ! Start hour is specified in the namelist (default 0) 
    380381               adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 
    381382               ! note this is wrong if time step has changed during run 
     
    386387            nhour   =   nn_time0 / 100 
    387388       nminute = ( nn_time0 - nhour * 100 ) 
    388             IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0) 
     389            isecond = ( nhour * NINT(rhhmm) + nminute ) * NINT(rmmss) 
     390            IF( isecond - ndt05 .LT. 0 )   ndastp = ndastp - 1         ! Start hour is specified in the namelist (default 0) 
    389391            adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 
    390392         ENDIF 
  • NEMO/trunk/src/OCE/DOM/depth_e3.F90

    r10069 r12377  
    3636   PUBLIC   e3_to_depth        ! called by domzgr.F90 
    3737       
    38    !! * Substitutions 
    39 #  include "vectopt_loop_substitute.h90" 
    4038   !!---------------------------------------------------------------------- 
    4139   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
  • NEMO/trunk/src/OCE/DOM/dom_oce.F90

    r10068 r12377  
    1212   !!            3.7  ! 2015-11  (G. Madec) introduce surface and scale factor ratio 
    1313   !!             -   ! 2015-11  (G. Madec, A. Coward)  time varying zgr by default 
     14   !!            4.1  ! 2019-08  (A. Coward, D. Storkey) rename prognostic variables in preparation for new time scheme. 
    1415   !!---------------------------------------------------------------------- 
    1516 
     
    3233   LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time 
    3334   LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc) 
    34    REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice 
    3535   REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer 
    3636   REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter 
    3737   INTEGER , PUBLIC ::   nn_euler       !: =0 start with forward time step or not (=1) 
    38    LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet 
    3938   LOGICAL , PUBLIC ::   ln_crs         !: Apply grid coarsening to dynamical model output or online passive tracers 
    4039 
     
    121120   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point 
    122121   ! 
    123    REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   ff_f  , ff_t                    !: Coriolis factor at f- & t-points  [1/s] 
     122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff_f  , ff_t                    !: Coriolis factor at f- & t-points  [1/s] 
    124123   !!---------------------------------------------------------------------- 
    125124   !! vertical coordinate and scale factors 
     
    129128   LOGICAL, PUBLIC ::   ln_sco       !: s-coordinate or hybrid z-s coordinate 
    130129   LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF  
    131    !                                                        !  ref.   ! before  !   now   ! after  ! 
    132    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m] 
    133    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m] 
    134    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m] 
    135    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m] 
    136    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m] 
    137    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m] 
    138    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m] 
    139  
    140    !                                                        !  ref.   ! before  !   now   ! 
    141    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0 , gdept_b , gdept_n   !: t- depth              [m] 
    142    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m] 
    143    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m] 
     130   !                                                        !  reference scale factors 
     131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0   !: t- vert. scale factor [m] 
     132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0   !: u- vert. scale factor [m] 
     133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0   !: v- vert. scale factor [m] 
     134   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0   !: f- vert. scale factor [m] 
     135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0   !: w- vert. scale factor [m] 
     136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0   !: uw-vert. scale factor [m] 
     137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0   !: vw-vert. scale factor [m] 
     138   !                                                        !  time-dependent scale factors 
     139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e3t, e3u, e3v, e3w, e3uw, e3vw  !: vert. scale factor [m] 
     140   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   e3f                             !: F-point vert. scale factor [m] 
     141 
     142   !                                                        !  reference depths of cells 
     143   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0  !: t- depth              [m] 
     144   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0  !: w- depth              [m] 
     145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0  !: w- depth (sum of e3w) [m] 
     146   !                                                        !  time-dependent depths of cells 
     147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  gdept, gdepw   
     148   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  gde3w   
    144149    
    145    !                                                      !  ref. ! before  !   now   !  after  ! 
    146    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0            ,    ht_n             !: t-depth              [m] 
    147    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m] 
    148    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m] 
    149    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m] 
    150    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m] 
     150   !                                                      !  reference heights of water column 
     151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m] 
     152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m] 
     153   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m] 
     154                                                          ! time-dependent heights of water column 
     155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht                     !: height of water column at T-points [m] 
     156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hu, hv, r1_hu, r1_hv   !: height of water column [m] and reciprocal [1/m] 
    151157 
    152158   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
     
    158164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m) 
    159165 
     166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep, bathy 
    160167 
    161168   !!---------------------------------------------------------------------- 
     
    170177   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book) 
    171178 
    172    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level             (ISF) 
    173    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level (ISF) 
    174    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF) 
     179   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level           (ISF) 
    175180 
    176181   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts 
     
    190195   INTEGER , PUBLIC ::   ndastp        !: time step date in yyyymmdd format 
    191196   INTEGER , PUBLIC ::   nday_year     !: current day counted from jan 1st of the current year 
    192    INTEGER , PUBLIC ::   nsec_year     !: current time step counted in second since 00h jan 1st of the current year 
    193    INTEGER , PUBLIC ::   nsec_month    !: current time step counted in second since 00h 1st day of the current month 
    194    INTEGER , PUBLIC ::   nsec_week     !: current time step counted in second since 00h of last monday 
    195    INTEGER , PUBLIC ::   nsec_day      !: current time step counted in second since 00h of the current day 
     197   INTEGER , PUBLIC ::   nsec_year     !: seconds between 00h jan 1st of the current  year and half of the current time step 
     198   INTEGER , PUBLIC ::   nsec_month    !: seconds between 00h 1st day of the current month and half of the current time step 
     199   INTEGER , PUBLIC ::   nsec_monday   !: seconds between 00h         of the last Monday   and half of the current time step 
     200   INTEGER , PUBLIC ::   nsec_day      !: seconds between 00h         of the current   day and half of the current time step 
    196201   REAL(wp), PUBLIC ::   fjulday       !: current julian day  
    197202   REAL(wp), PUBLIC ::   fjulstartyear !: first day of the current year in julian days 
    198203   REAL(wp), PUBLIC ::   adatrj        !: number of elapsed days since the begining of the whole simulation 
    199204   !                                   !: (cumulative duration of previous runs that may have used different time-step size) 
    200    INTEGER , PUBLIC, DIMENSION(0: 2) ::   nyear_len     !: length in days of the previous/current/next year 
    201    INTEGER , PUBLIC, DIMENSION(0:13) ::   nmonth_len    !: length in days of the months of the current year 
    202    INTEGER , PUBLIC, DIMENSION(0:13) ::   nmonth_half   !: second since Jan 1st 0h of the current year and the half of the months 
    203    INTEGER , PUBLIC, DIMENSION(0:13) ::   nmonth_end    !: second since Jan 1st 0h of the current year and the end of the months 
    204    INTEGER , PUBLIC                  ::   nsec1jan000   !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year 
     205   INTEGER , PUBLIC, DIMENSION(  0: 2) ::   nyear_len     !: length in days of the previous/current/next year 
     206   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_len    !: length in days of the months of the current year 
     207   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_beg    !: second since Jan 1st 0h of the current year and the half of the months 
     208   INTEGER , PUBLIC                  ::   nsec1jan000     !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year 
     209   INTEGER , PUBLIC                  ::   nsec000_1jan000   !: second since Jan 1st 0h of nit000 year and nit000 
     210   INTEGER , PUBLIC                  ::   nsecend_1jan000   !: second since Jan 1st 0h of nit000 year and nitend 
    205211 
    206212   !!---------------------------------------------------------------------- 
     
    256262         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) ) 
    257263         ! 
    258       ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,      & 
    259          &      gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) ,                             & 
    260          &      gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 
    261          ! 
    262       ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) ,   & 
    263          &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   &  
    264          &      e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) ,   &  
    265          &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             & 
    266          !                                                          ! 
    267          &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         & 
    268          &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &                
    269          &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                        
    270          ! 
    271       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
    272          &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     & 
    273          &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     & 
    274          &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  ) 
    275          ! 
    276          ! 
    277       ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) ) 
     264      ALLOCATE( gdept_0(jpi,jpj,jpk)     , gdepw_0(jpi,jpj,jpk)     , gde3w_0(jpi,jpj,jpk) ,      & 
     265         &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(4) ) 
     266         ! 
     267      ALLOCATE( e3t_0(jpi,jpj,jpk)     , e3u_0(jpi,jpj,jpk)     , e3v_0(jpi,jpj,jpk)     , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk)     ,   & 
     268         &      e3t  (jpi,jpj,jpk,jpt) , e3u  (jpi,jpj,jpk,jpt) , e3v  (jpi,jpj,jpk,jpt) , e3f  (jpi,jpj,jpk) , e3w  (jpi,jpj,jpk,jpt) ,   &  
     269         &      e3uw_0(jpi,jpj,jpk)     , e3vw_0(jpi,jpj,jpk)     ,         & 
     270         &      e3uw  (jpi,jpj,jpk,jpt) , e3vw  (jpi,jpj,jpk,jpt) ,    STAT=ierr(5) )                        
     271         ! 
     272      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj)    , hv_0(jpi,jpj)     ,                                             & 
     273         &      ht  (jpi,jpj) , hu(  jpi,jpj,jpt), hv(  jpi,jpj,jpt) , r1_hu(jpi,jpj,jpt) , r1_hv(jpi,jpj,jpt) ,   & 
     274         &                      STAT=ierr(6)  ) 
     275         ! 
     276      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  )  
     277         ! 
     278      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) ) 
    278279         ! 
    279280      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &  
     
    281282         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
    282283         ! 
    283       ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     & 
    284          &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) ) 
     284      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) 
    285285         ! 
    286286      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &  
  • NEMO/trunk/src/OCE/DOM/domain.F90

    r11536 r12377  
    3030   USE trc_oce        ! shared ocean & passive tracers variab 
    3131   USE phycst         ! physical constants 
    32    USE closea         ! closed seas 
    3332   USE domhgr         ! domain: set the horizontal mesh 
    3433   USE domzgr         ! domain: set the vertical mesh 
     
    3837   USE c1d            ! 1D configuration 
    3938   USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine) 
    40    USE wet_dry,  ONLY : ll_wd 
     39   USE wet_dry, ONLY : ll_wd 
     40   USE closea , ONLY : dom_clo ! closed seas 
    4141   ! 
    4242   USE in_out_manager ! I/O manager 
     
    5858CONTAINS 
    5959 
    60    SUBROUTINE dom_init(cdstr) 
     60   SUBROUTINE dom_init( Kbb, Kmm, Kaa, cdstr ) 
    6161      !!---------------------------------------------------------------------- 
    6262      !!                  ***  ROUTINE dom_init  *** 
     
    7373      !!              - 1D configuration, move Coriolis, u and v at T-point 
    7474      !!---------------------------------------------------------------------- 
     75      INTEGER          , INTENT(in) :: Kbb, Kmm, Kaa          ! ocean time level indices 
     76      CHARACTER (len=*), INTENT(in) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables 
     77      ! 
    7578      INTEGER ::   ji, jj, jk, ik   ! dummy loop indices 
    7679      INTEGER ::   iconf = 0    ! local integers 
    7780      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"  
    78       CHARACTER (len=*), INTENT(IN) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables 
    7981      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
    8082      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0 
     
    134136      ENDIF 
    135137      ! 
    136       CALL dom_hgr                     ! Horizontal mesh 
    137       CALL dom_zgr( ik_top, ik_bot )   ! Vertical mesh and bathymetry 
    138       CALL dom_msk( ik_top, ik_bot )   ! Masks 
    139       IF( ln_closea )   CALL dom_clo   ! ln_closea=T : closed seas included in the simulation 
    140                                        ! Read in masks to define closed seas and lakes  
    141       ! 
    142       DO jj = 1, jpj                   ! depth of the iceshelves 
    143          DO ji = 1, jpi 
    144             ik = mikt(ji,jj) 
    145             risfdep(ji,jj) = gdepw_0(ji,jj,ik) 
    146          END DO 
    147       END DO 
     138      CALL dom_hgr                      ! Horizontal mesh 
     139 
     140      IF( ln_closea ) CALL dom_clo      ! Read in masks to define closed seas and lakes 
     141 
     142      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry 
     143 
     144      CALL dom_msk( ik_top, ik_bot )    ! Masks 
    148145      ! 
    149146      ht_0(:,:) = 0._wp  ! Reference ocean thickness 
     
    161158      ! 
    162159         !       before        !          now          !       after         ! 
    163             gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
    164             gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          ! 
    165                                    gde3w_n = gde3w_0   !        ---          ! 
     160            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   ;   gdept(:,:,:,Kaa) = gdept_0   ! depth of grid-points 
     161            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   ;   gdepw(:,:,:,Kaa) = gdepw_0   ! 
     162                                   gde3w = gde3w_0   !        ---          ! 
    166163         !                                                                   
    167               e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors 
    168               e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! 
    169               e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
    170                                      e3f_n =   e3f_0   !        ---          ! 
    171               e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          ! 
    172              e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          ! 
    173              e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          ! 
     164              e3t(:,:,:,Kbb) =   e3t_0  ;     e3t(:,:,:,Kmm) =   e3t_0   ;   e3t(:,:,:,Kaa) =  e3t_0    ! scale factors 
     165              e3u(:,:,:,Kbb) =   e3u_0  ;     e3u(:,:,:,Kmm) =   e3u_0   ;   e3u(:,:,:,Kaa) =  e3u_0    ! 
     166              e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    ! 
     167                                     e3f =   e3f_0   !        ---          ! 
     168              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   ;    e3w(:,:,:,Kaa) =   e3w_0   !  
     169             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   ;   e3uw(:,:,:,Kaa) =  e3uw_0   !   
     170             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   ;   e3vw(:,:,:,Kaa) =  e3vw_0   ! 
    174171         ! 
    175172         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF 
     
    177174         ! 
    178175         !        before       !          now          !       after         ! 
    179                                       ht_n =    ht_0   !                     ! water column thickness 
    180                hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  
    181                hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   ! 
    182             r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness 
    183             r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   ! 
     176                                      ht =    ht_0   !                     ! water column thickness 
     177               hu(:,:,Kbb) =    hu_0  ;      hu(:,:,Kmm) =    hu_0   ;    hu(:,:,Kaa) =    hu_0   !  
     178               hv(:,:,Kbb) =    hv_0  ;      hv(:,:,Kmm) =    hv_0   ;    hv(:,:,Kaa) =    hv_0   ! 
     179            r1_hu(:,:,Kbb) = z1_hu_0  ;   r1_hu(:,:,Kmm) = z1_hu_0   ; r1_hu(:,:,Kaa) = z1_hu_0   ! inverse of water column thickness 
     180            r1_hv(:,:,Kbb) = z1_hv_0  ;   r1_hv(:,:,Kmm) = z1_hv_0   ; r1_hv(:,:,Kaa) = z1_hv_0   ! 
    184181         ! 
    185182         ! 
    186183      ELSE                       != time varying : initialize before/now/after variables 
    187184         ! 
    188          IF( .NOT.l_offline )  CALL dom_vvl_init  
     185         IF( .NOT.l_offline )  CALL dom_vvl_init( Kbb, Kmm, Kaa ) 
    189186         ! 
    190187      ENDIF 
     
    192189      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point 
    193190      ! 
    194       IF( ln_meshmask .AND. .NOT.ln_iscpl )                        CALL dom_wri     ! Create a domain file 
    195       IF( ln_meshmask .AND.      ln_iscpl .AND. .NOT.ln_rstart )   CALL dom_wri     ! Create a domain file 
    196       IF(                                       .NOT.ln_rstart )   CALL dom_ctl     ! Domain control 
    197       ! 
    198       IF( ln_write_cfg )   CALL cfg_write         ! create the configuration file 
     191      IF( ln_meshmask    )   CALL dom_wri       ! Create a domain file 
     192      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control 
     193      ! 
     194      IF( ln_write_cfg   )   CALL cfg_write     ! create the configuration file 
    199195      ! 
    200196      IF(lwp) THEN 
     
    292288         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     & 
    293289         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
    294          &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 
    295       NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
     290         &             ln_cfmeta, ln_xios_read, nn_wxios 
     291      NAMELIST/namdom/ ln_linssh, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
    296292#if defined key_netcdf4 
    297293      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
     
    306302      ! 
    307303      ! 
    308       REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run 
    309304      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901) 
    310305901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist' ) 
    311       REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run 
    312306      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 ) 
    313307902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist' ) 
     
    343337         WRITE(numout,*) '      overwrite an existing file      ln_clobber      = ', ln_clobber 
    344338         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz      = ', nn_chunksz 
    345          WRITE(numout,*) '      IS coupling at the restart step ln_iscpl        = ', ln_iscpl 
    346339         IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    347340            WRITE(numout,*) '      READ restart for a single file using XIOS ln_xios_read =', ln_xios_read 
     
    404397#endif 
    405398 
    406       REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 
    407399      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 
    408400903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist' ) 
    409       REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep) 
    410401      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
    411402904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist' ) 
     
    417408         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh 
    418409         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
    419          WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]' 
    420410         WRITE(numout,*) '      ocean time step                         rn_rdt      = ', rn_rdt 
    421411         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp 
     
    436426#if defined key_netcdf4 
    437427      !                             ! NetCDF 4 case   ("key_netcdf4" defined) 
    438       REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF 
    439428      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907) 
    440429907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namnc4 in reference namelist' ) 
    441       REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF 
    442430      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 ) 
    443431908   IF( ios >  0 )   CALL ctl_nam ( ios , 'namnc4 in configuration namelist' ) 
  • NEMO/trunk/src/OCE/DOM/dommsk.F90

    r11536 r12377  
    4444 
    4545   !! * Substitutions 
    46 #  include "vectopt_loop_substitute.h90" 
     46#  include "do_loop_substitute.h90" 
    4747   !!---------------------------------------------------------------------- 
    4848   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    103103      !!--------------------------------------------------------------------- 
    104104      ! 
    105       REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition 
    106105      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 ) 
    107106901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlbc in reference namelist' ) 
    108       REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition 
    109107      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 ) 
    110108902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlbc in configuration namelist' ) 
     
    133131      ! 
    134132      tmask(:,:,:) = 0._wp 
    135       DO jj = 1, jpj 
    136          DO ji = 1, jpi 
    137             iktop = k_top(ji,jj) 
    138             ikbot = k_bot(ji,jj) 
    139             IF( iktop /= 0 ) THEN       ! water in the column 
    140                tmask(ji,jj,iktop:ikbot  ) = 1._wp 
    141             ENDIF 
    142          END DO   
    143       END DO 
     133      DO_2D_11_11 
     134         iktop = k_top(ji,jj) 
     135         ikbot = k_bot(ji,jj) 
     136         IF( iktop /= 0 ) THEN       ! water in the column 
     137            tmask(ji,jj,iktop:ikbot  ) = 1._wp 
     138         ENDIF 
     139      END_2D 
    144140      ! 
    145141      ! the following call is mandatory 
     
    148144 
    149145     ! Mask corrections for bdy (read in mppini2) 
    150       REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries 
    151146      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
    152147903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' ) 
    153       REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
    154148      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
    155149904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' ) 
     
    159153         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 
    160154         CALL iom_close( inum ) 
    161          DO jk = 1, jpkm1 
    162             DO jj = 1, jpj 
    163                DO ji = 1, jpi 
    164                   tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj) 
    165                END DO 
    166             END DO 
    167          END DO 
     155         DO_3D_11_11( 1, jpkm1 ) 
     156            tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj) 
     157         END_3D 
    168158      ENDIF 
    169159          
     
    173163      DO jk = 1, jpk 
    174164         DO jj = 1, jpjm1 
    175             DO ji = 1, fs_jpim1   ! vector loop 
     165            DO ji = 1, jpim1   ! vector loop 
    176166               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk) 
    177167               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk) 
     
    247237         DO jk = 1, jpk 
    248238            zwf(:,:) = fmask(:,:,jk)          
    249             DO jj = 2, jpjm1 
    250                DO ji = fs_2, fs_jpim1   ! vector opt. 
    251                   IF( fmask(ji,jj,jk) == 0._wp ) THEN 
    252                      fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
    253                         &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
    254                   ENDIF 
    255                END DO 
    256             END DO 
     239            DO_2D_00_00 
     240               IF( fmask(ji,jj,jk) == 0._wp ) THEN 
     241                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
     242                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
     243               ENDIF 
     244            END_2D 
    257245            DO jj = 2, jpjm1 
    258246               IF( fmask(1,jj,jk) == 0._wp ) THEN 
  • NEMO/trunk/src/OCE/DOM/domvvl.F90

    r11536 r12377  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
     10   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    1314   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
    1415   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
    15    !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid 
     16   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid 
    1617   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
    1718   !!   dom_vvl_rst      : read/write restart file 
     
    3637 
    3738   PUBLIC  dom_vvl_init       ! called by domain.F90 
     39   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90 
    3840   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    39    PUBLIC  dom_vvl_sf_swp     ! called by step.F90 
     41   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
    4042   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    4143 
     
    6264 
    6365   !! * Substitutions 
    64 #  include "vectopt_loop_substitute.h90" 
     66#  include "do_loop_substitute.h90" 
    6567   !!---------------------------------------------------------------------- 
    6668   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    9395 
    9496 
    95    SUBROUTINE dom_vvl_init 
     97   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 
    9698      !!---------------------------------------------------------------------- 
    9799      !!                ***  ROUTINE dom_vvl_init  *** 
     
    102104      !! ** Method  :  - use restart file and/or initialize 
    103105      !!               - interpolate scale factors 
     106      !! 
     107      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     108      !!              - Regrid: e3[u/v](:,:,:,Kmm) 
     109      !!                        e3[u/v](:,:,:,Kmm)        
     110      !!                        e3w(:,:,:,Kmm)            
     111      !!                        e3[u/v]w_b 
     112      !!                        e3[u/v]w_n       
     113      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 
     114      !!              - h(t/u/v)_0 
     115      !!              - frq_rst_e3t and frq_rst_hdv 
     116      !! 
     117      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     118      !!---------------------------------------------------------------------- 
     119      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     120      ! 
     121      IF(lwp) WRITE(numout,*) 
     122      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
     123      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     124      ! 
     125      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
     126      ! 
     127      !                    ! Allocate module arrays 
     128      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
     129      ! 
     130      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     131      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
     132      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
     133      ! 
     134      CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 
     135      ! 
     136   END SUBROUTINE dom_vvl_init 
     137   ! 
     138   SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 
     139      !!---------------------------------------------------------------------- 
     140      !!                ***  ROUTINE dom_vvl_init  *** 
     141      !!                    
     142      !! ** Purpose :  Interpolation of all scale factors,  
     143      !!               depths and water column heights 
     144      !! 
     145      !! ** Method  :  - interpolate scale factors 
    104146      !! 
    105147      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     
    115157      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    116158      !!---------------------------------------------------------------------- 
     159      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     160      !!---------------------------------------------------------------------- 
    117161      INTEGER ::   ji, jj, jk 
    118162      INTEGER ::   ii0, ii1, ij0, ij1 
     
    120164      !!---------------------------------------------------------------------- 
    121165      ! 
    122       IF(lwp) WRITE(numout,*) 
    123       IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
    124       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    125       ! 
    126       CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
    127       ! 
    128       !                    ! Allocate module arrays 
    129       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    130       ! 
    131       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    132       CALL dom_vvl_rst( nit000, 'READ' ) 
    133       e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    134       ! 
    135166      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    136167      !                                ! Horizontal interpolation of e3t 
    137       CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U 
    138       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    139       CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V  
    140       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    141       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F 
     168      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
     169      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     170      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
     171      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     172      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
    142173      !                                ! Vertical interpolation of e3t,u,v  
    143       CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W 
    144       CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  ) 
    145       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW 
    146       CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    147       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW 
    148       CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     174      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
     175      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
     176      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW 
     177      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     178      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW 
     179      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    149180 
    150181      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 
    151       e3t_a(:,:,:) = e3t_n(:,:,:) 
    152       e3u_a(:,:,:) = e3u_n(:,:,:) 
    153       e3v_a(:,:,:) = e3v_n(:,:,:) 
     182      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 
     183      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 
     184      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 
    154185      ! 
    155186      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep) 
    156       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration) 
    157       gdepw_n(:,:,1) = 0.0_wp 
    158       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg 
    159       gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 
    160       gdepw_b(:,:,1) = 0.0_wp 
    161       DO jk = 2, jpk                               ! vertical sum 
    162          DO jj = 1,jpj 
    163             DO ji = 1,jpi 
    164                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    165                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    166                !                             ! 0.5 where jk = mikt      
    167 !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
    168                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    169                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    170                gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
    171                   &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk))  
    172                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    173                gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 
    174                gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  & 
    175                   &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk))  
    176             END DO 
    177          END DO 
     187      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration) 
     188      gdepw(:,:,1,Kmm) = 0.0_wp 
     189      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg 
     190      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
     191      gdepw(:,:,1,Kbb) = 0.0_wp 
     192      DO_3D_11_11( 2, jpk ) 
     193         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     194         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     195         !                             ! 0.5 where jk = mikt      
     196!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
     197         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
     198         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     199         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     200            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     201         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     202         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     203         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     204            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     205      END_3D 
     206      ! 
     207      !                    !==  thickness of the water column  !!   (ocean portion only) 
     208      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) .... 
     209      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 
     210      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 
     211      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 
     212      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 
     213      DO jk = 2, jpkm1 
     214         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     215         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 
     216         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 
     217         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 
     218         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    178219      END DO 
    179220      ! 
    180       !                    !==  thickness of the water column  !!   (ocean portion only) 
    181       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) .... 
    182       hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 
    183       hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 
    184       hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 
    185       hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 
    186       DO jk = 2, jpkm1 
    187          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    188          hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
    189          hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
    190          hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
    191          hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
    192       END DO 
    193       ! 
    194221      !                    !==  inverse of water column thickness   ==!   (u- and v- points) 
    195       r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF 
    196       r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 
    197       r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
    198       r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     222      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF 
     223      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 
     224      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 
     225      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 
    199226 
    200227      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies) 
     
    211238         ENDIF 
    212239         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    213             DO jj = 1, jpj 
    214                DO ji = 1, jpi 
     240            DO_2D_11_11 
    215241!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    216                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    217                      ! values outside the equatorial band and transition zone (ztilde) 
    218                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    219                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    220                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    221                      ! values inside the equatorial band (ztilde as zstar) 
    222                      frq_rst_e3t(ji,jj) =  0.0_wp 
    223                      frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
    224                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    225                      !                                      ! (linearly transition from z-tilde to z-star) 
    226                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    227                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    228                         &                                          * 180._wp / 3.5_wp ) ) 
    229                      frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
    230                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
    231                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    232                         &                                          * 180._wp / 3.5_wp ) ) 
    233                   ENDIF 
    234                END DO 
    235             END DO 
     242               IF( ABS(gphit(ji,jj)) >= 6.) THEN 
     243                  ! values outside the equatorial band and transition zone (ztilde) 
     244                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
     245                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     246               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
     247                  ! values inside the equatorial band (ztilde as zstar) 
     248                  frq_rst_e3t(ji,jj) =  0.0_wp 
     249                  frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
     250               ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
     251                  !                                      ! (linearly transition from z-tilde to z-star) 
     252                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
     253                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     254                     &                                          * 180._wp / 3.5_wp ) ) 
     255                  frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
     256                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
     257                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     258                     &                                          * 180._wp / 3.5_wp ) ) 
     259               ENDIF 
     260            END_2D 
    236261            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    237262               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
     
    263288      ENDIF 
    264289      ! 
    265    END SUBROUTINE dom_vvl_init 
    266  
    267  
    268    SUBROUTINE dom_vvl_sf_nxt( kt, kcall )  
     290   END SUBROUTINE dom_vvl_zgr 
     291 
     292 
     293   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )  
    269294      !!---------------------------------------------------------------------- 
    270295      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
     
    288313      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling. 
    289314      !!---------------------------------------------------------------------- 
    290       INTEGER, INTENT( in )           ::   kt      ! time step 
    291       INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
     315      INTEGER, INTENT( in )           ::   kt             ! time step 
     316      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step 
     317      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence 
    292318      ! 
    293319      INTEGER                ::   ji, jj, jk            ! dummy loop indices 
     
    321347      !                                                ! --------------------------------------------- ! 
    322348      ! 
    323       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     349      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 
    324350      DO jk = 1, jpkm1 
    325          ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 
    326          e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     351         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 
     352         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 
    327353      END DO 
    328354      ! 
     
    337363         zht(:,:)   = 0._wp 
    338364         DO jk = 1, jpkm1 
    339             zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    340             zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     365            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 
     366            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    341367         END DO 
    342368         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
     
    348374               DO jk = 1, jpkm1 
    349375                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    350                      &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
     376                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 
    351377               END DO 
    352378            ENDIF 
     
    361387         IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    362388            DO jk = 1, jpkm1 
    363                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
     389               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    364390            END DO 
    365391         ELSE                         ! layer case 
    366392            DO jk = 1, jpkm1 
    367                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
     393               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    368394            END DO 
    369395         ENDIF 
     
    381407         zwu(:,:) = 0._wp 
    382408         zwv(:,:) = 0._wp 
    383          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    384             DO jj = 1, jpjm1 
    385                DO ji = 1, fs_jpim1   ! vector opt. 
    386                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    387                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    388                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    389                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    390                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    391                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    392                END DO 
    393             END DO 
    394          END DO 
    395          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    396             DO ji = 1, jpi 
    397                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    398                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    399             END DO 
    400          END DO 
    401          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    402             DO jj = 2, jpjm1 
    403                DO ji = fs_2, fs_jpim1   ! vector opt. 
    404                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    405                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    406                      &                                            ) * r1_e1e2t(ji,jj) 
    407                END DO 
    408             END DO 
    409          END DO 
     409         DO_3D_10_10( 1, jpkm1 ) 
     410            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
     411               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     412            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
     413               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     414            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     415            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     416         END_3D 
     417         DO_2D_11_11 
     418            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     419            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     420         END_2D 
     421         DO_3D_00_00( 1, jpkm1 ) 
     422            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     423               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
     424               &                                            ) * r1_e1e2t(ji,jj) 
     425         END_3D 
    410426         !                       ! d - thickness diffusion transport: boundary conditions 
    411427         !                             (stored for tracer advction and continuity equation) 
     
    476492            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    477493         END DO 
    478          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     494         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 
    479495         DO jk = 1, jpkm1 
    480             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     496            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 
    481497         END DO 
    482498 
     
    486502      !                                           ! ---baroclinic part--------- ! 
    487503         DO jk = 1, jpkm1 
    488             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     504            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    489505         END DO 
    490506      ENDIF 
     
    501517         zht(:,:) = 0.0_wp 
    502518         DO jk = 1, jpkm1 
    503             zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    504          END DO 
    505          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
     519            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     520         END DO 
     521         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 
    506522         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    507          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 
     523         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 
    508524         ! 
    509525         zht(:,:) = 0.0_wp 
    510526         DO jk = 1, jpkm1 
    511             zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 
    512          END DO 
    513          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
     527            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 
     528         END DO 
     529         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 
    514530         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    515          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 
     531         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 
    516532         ! 
    517533         zht(:,:) = 0.0_wp 
    518534         DO jk = 1, jpkm1 
    519             zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 
    520          END DO 
    521          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 
     535            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 
     536         END DO 
     537         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 
    522538         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    523          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 
    524          ! 
    525          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) ) 
     539         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 
     540         ! 
     541         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) ) 
    526542         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    527          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 
    528          ! 
    529          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) ) 
     543         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 
     544         ! 
     545         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) ) 
    530546         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    531          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 
    532          ! 
    533          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) ) 
     547         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 
     548         ! 
     549         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) ) 
    534550         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    535          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
     551         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 
    536552      END IF 
    537553 
     
    540556      ! *********************************** ! 
    541557 
    542       CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 
    543       CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 
     558      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 
     559      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 
    544560 
    545561      ! *********************************** ! 
     
    547563      ! *********************************** ! 
    548564 
    549       hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 
    550       hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 
     565      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 
     566      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 
    551567      DO jk = 2, jpkm1 
    552          hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 
    553          hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 
     568         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 
     569         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 
    554570      END DO 
    555571      !                                        ! Inverse of the local depth 
    556572!!gm BUG ?  don't understand the use of umask_i here ..... 
    557       r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 
    558       r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     573      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 
     574      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 
    559575      ! 
    560576      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt') 
     
    563579 
    564580 
    565    SUBROUTINE dom_vvl_sf_swp( kt ) 
    566       !!---------------------------------------------------------------------- 
    567       !!                ***  ROUTINE dom_vvl_sf_swp  *** 
     581   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 
     582      !!---------------------------------------------------------------------- 
     583      !!                ***  ROUTINE dom_vvl_sf_update  *** 
    568584      !!                    
    569       !! ** Purpose :  compute time filter and swap of scale factors  
     585      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors  
    570586      !!               compute all depths and related variables for next time step 
    571587      !!               write outputs and restart file 
    572588      !! 
    573       !! ** Method  :  - swap of e3t with trick for volume/tracer conservation 
     589      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 
    574590      !!               - reconstruct scale factor at other grid points (interpolate) 
    575591      !!               - recompute depths and water height fields 
    576592      !! 
    577       !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 
     593      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step 
    578594      !!               - Recompute: 
    579595      !!                    e3(u/v)_b        
    580       !!                    e3w_n            
     596      !!                    e3w(:,:,:,Kmm)            
    581597      !!                    e3(u/v)w_b       
    582598      !!                    e3(u/v)w_n       
    583       !!                    gdept_n, gdepw_n  and gde3w_n 
     599      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w 
    584600      !!                    h(u/v) and h(u/v)r 
    585601      !! 
     
    587603      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    588604      !!---------------------------------------------------------------------- 
    589       INTEGER, INTENT( in ) ::   kt   ! time step 
     605      INTEGER, INTENT( in ) ::   kt              ! time step 
     606      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices 
    590607      ! 
    591608      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    595612      IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
    596613      ! 
    597       IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp') 
     614      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update') 
    598615      ! 
    599616      IF( kt == nit000 )   THEN 
    600617         IF(lwp) WRITE(numout,*) 
    601          IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 
    602          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step' 
     618         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 
     619         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 
    603620      ENDIF 
    604621      ! 
     
    615632         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    616633      ENDIF 
    617       gdept_b(:,:,:) = gdept_n(:,:,:) 
    618       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
    619  
    620       e3t_n(:,:,:) = e3t_a(:,:,:) 
    621       e3u_n(:,:,:) = e3u_a(:,:,:) 
    622       e3v_n(:,:,:) = e3v_a(:,:,:) 
    623634 
    624635      ! Compute all missing vertical scale factor and depths 
     
    626637      ! Horizontal scale factor interpolations 
    627638      ! -------------------------------------- 
    628       ! - ML - e3u_b and e3v_b are allready computed in dynnxt 
    629       ! - JC - hu_b, hv_b, hur_b, hvr_b also 
     639      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 
     640      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    630641       
    631       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
     642      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  ) 
    632643       
    633644      ! Vertical scale factor interpolations 
    634       CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  ) 
    635       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    636       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    637       CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  ) 
    638       CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    639       CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     645      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
     646      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     647      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
     648      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  ) 
     649      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     650      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    640651 
    641652      ! t- and w- points depth (set the isf depth as it is in the initial step) 
    642       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    643       gdepw_n(:,:,1) = 0.0_wp 
    644       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    645       DO jk = 2, jpk 
    646          DO jj = 1,jpj 
    647             DO ji = 1,jpi 
    648               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    649                                                                  ! 1 for jk = mikt 
    650                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    651                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    652                gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  & 
    653                    &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) )  
    654                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    655             END DO 
    656          END DO 
    657       END DO 
     653      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     654      gdepw(:,:,1,Kmm) = 0.0_wp 
     655      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
     656      DO_3D_11_11( 2, jpk ) 
     657        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     658                                                           ! 1 for jk = mikt 
     659         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     660         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     661         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     662             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     663         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     664      END_3D 
    658665 
    659666      ! Local depth and Inverse of the local depth of the water 
    660667      ! ------------------------------------------------------- 
    661       hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:) 
    662       hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:) 
    663       ! 
    664       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 
     668      ! 
     669      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 
    665670      DO jk = 2, jpkm1 
    666          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     671         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    667672      END DO 
    668673 
    669674      ! write restart file 
    670675      ! ================== 
    671       IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' ) 
    672       ! 
    673       IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp') 
    674       ! 
    675    END SUBROUTINE dom_vvl_sf_swp 
     676      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 
     677      ! 
     678      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update') 
     679      ! 
     680   END SUBROUTINE dom_vvl_sf_update 
    676681 
    677682 
     
    704709         ! 
    705710      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    706          DO jk = 1, jpk 
    707             DO jj = 1, jpjm1 
    708                DO ji = 1, fs_jpim1   ! vector opt. 
    709                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    710                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    711                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    712                END DO 
    713             END DO 
    714          END DO 
     711         DO_3D_10_10( 1, jpk ) 
     712            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     713               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     714               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     715         END_3D 
    715716         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    716717         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    717718         ! 
    718719      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    719          DO jk = 1, jpk 
    720             DO jj = 1, jpjm1 
    721                DO ji = 1, fs_jpim1   ! vector opt. 
    722                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    723                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    724                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    725                END DO 
    726             END DO 
    727          END DO 
     720         DO_3D_10_10( 1, jpk ) 
     721            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     722               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     723               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     724         END_3D 
    728725         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    729726         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    730727         ! 
    731728      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    732          DO jk = 1, jpk 
    733             DO jj = 1, jpjm1 
    734                DO ji = 1, fs_jpim1   ! vector opt. 
    735                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    736                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    737                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    738                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    739                END DO 
    740             END DO 
    741          END DO 
     729         DO_3D_10_10( 1, jpk ) 
     730            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     731               &                       *    r1_e1e2f(ji,jj)                                                  & 
     732               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     733               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     734         END_3D 
    742735         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    743736         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     
    783776 
    784777 
    785    SUBROUTINE dom_vvl_rst( kt, cdrw ) 
     778   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 
    786779      !!--------------------------------------------------------------------- 
    787780      !!                   ***  ROUTINE dom_vvl_rst  *** 
     
    795788      !!                they are set to 0. 
    796789      !!---------------------------------------------------------------------- 
    797       INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    798       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     790      INTEGER         , INTENT(in) ::   kt        ! ocean time-step 
     791      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices 
     792      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag 
    799793      ! 
    800794      INTEGER ::   ji, jj, jk 
     
    806800         IF( ln_rstart ) THEN                   !* Read the restart file 
    807801            CALL rst_read_open                  !  open the restart file if necessary 
    808             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    ) 
     802            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
    809803            ! 
    810804            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
     
    813807            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    814808            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     809            ! 
    815810            !                             ! --------- ! 
    816811            !                             ! all cases ! 
    817812            !                             ! --------- ! 
     813            ! 
    818814            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    819                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    820                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
     815               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     816               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    821817               ! needed to restart if land processor not computed  
    822                IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 
     818               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 
    823819               WHERE ( tmask(:,:,:) == 0.0_wp )  
    824                   e3t_n(:,:,:) = e3t_0(:,:,:) 
    825                   e3t_b(:,:,:) = e3t_0(:,:,:) 
     820                  e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
     821                  e3t(:,:,:,Kbb) = e3t_0(:,:,:) 
    826822               END WHERE 
    827823               IF( neuler == 0 ) THEN 
    828                   e3t_b(:,:,:) = e3t_n(:,:,:) 
     824                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    829825               ENDIF 
    830826            ELSE IF( id1 > 0 ) THEN 
    831                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 
     827               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    832828               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    833829               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    834                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    835                e3t_n(:,:,:) = e3t_b(:,:,:) 
     830               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     831               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    836832               neuler = 0 
    837833            ELSE IF( id2 > 0 ) THEN 
    838                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 
     834               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    839835               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    840836               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    841                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
    842                e3t_b(:,:,:) = e3t_n(:,:,:) 
     837               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
     838               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    843839               neuler = 0 
    844840            ELSE 
    845                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 
     841               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    846842               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    847843               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    848844               DO jk = 1, jpk 
    849                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
     845                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
    850846                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    851847                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    852848               END DO 
    853                e3t_b(:,:,:) = e3t_n(:,:,:) 
     849               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    854850               neuler = 0 
    855851            ENDIF 
     
    888884               IF( cn_cfg == 'wad' ) THEN 
    889885                  ! Wetting and drying test case 
    890                   CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
    891                   tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    892                   sshn (:,:)     = sshb(:,:) 
    893                   un   (:,:,:)   = ub  (:,:,:) 
    894                   vn   (:,:,:)   = vb  (:,:,:) 
     886                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     887                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     888                  ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
     889                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     890                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
    895891               ELSE 
    896892                  ! if not test case 
    897                   sshn(:,:) = -ssh_ref 
    898                   sshb(:,:) = -ssh_ref 
    899  
    900                   DO jj = 1, jpj 
    901                      DO ji = 1, jpi 
    902                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    903  
    904                            sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    905                            sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    906                            ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    907                         ENDIF 
    908                      ENDDO 
    909                   ENDDO 
     893                  ssh(:,:,Kmm) = -ssh_ref 
     894                  ssh(:,:,Kbb) = -ssh_ref 
     895 
     896                  DO_2D_11_11 
     897                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     898                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     899                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     900                     ENDIF 
     901                  END_2D 
    910902               ENDIF !If test case else 
    911903 
    912904               ! Adjust vertical metrics for all wad 
    913905               DO jk = 1, jpk 
    914                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) & 
     906                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) & 
    915907                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    916908                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    917909               END DO 
    918                e3t_b(:,:,:) = e3t_n(:,:,:) 
     910               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    919911 
    920912               DO ji = 1, jpi 
     
    930922               ! Just to read set ssh in fact, called latter once vertical grid 
    931923               ! is set up: 
    932 !               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  ) 
     924!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    933925!               ! 
    934926!               DO jk=1,jpk 
    935 !                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) & 
     927!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
    936928!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 
    937929!               END DO 
    938 !               e3t_n(:,:,:) = e3t_b(:,:,:) 
    939                 sshn(:,:)=0._wp 
    940                 e3t_n(:,:,:)=e3t_0(:,:,:) 
    941                 e3t_b(:,:,:)=e3t_0(:,:,:) 
     930!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     931                ssh(:,:,Kmm)=0._wp 
     932                e3t(:,:,:,Kmm)=e3t_0(:,:,:) 
     933                e3t(:,:,:,Kbb)=e3t_0(:,:,:) 
    942934               ! 
    943935            END IF           ! end of ll_wd edits 
     
    957949         !                                           ! all cases ! 
    958950         !                                           ! --------- ! 
    959          CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios ) 
    960          CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) 
     951         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 
     952         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 
    961953         !                                           ! ----------------------- ! 
    962954         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases ! 
     
    991983      !!----------------------------------------------------------------------  
    992984      ! 
    993       REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :  
    994985      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 
    995986901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' ) 
    996       REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run 
    997987      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 
    998988902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' ) 
     
    10331023      ! 
    10341024      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    1035       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    10361025      ! 
    10371026      IF(lwp) THEN                   ! Print the choice 
  • NEMO/trunk/src/OCE/DOM/domwri.F90

    r11532 r12377  
    1616   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 
    1717   !!---------------------------------------------------------------------- 
     18   ! 
    1819   USE dom_oce         ! ocean space and time domain 
    1920   USE phycst ,   ONLY :   rsmall 
     
    3233 
    3334   !! * Substitutions 
    34 #  include "vectopt_loop_substitute.h90" 
     35#  include "do_loop_substitute.h90" 
    3536   !!---------------------------------------------------------------------- 
    3637   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    99100       
    100101      CALL dom_uniq( zprw, 'T' ) 
    101       DO jj = 1, jpj 
    102          DO ji = 1, jpi 
    103             zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
    104          END DO 
    105       END DO                             !    ! unique point mask 
     102      DO_2D_11_11 
     103         zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
     104      END_2D 
    106105      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )   
    107106      CALL dom_uniq( zprw, 'U' ) 
    108       DO jj = 1, jpj 
    109          DO ji = 1, jpi 
    110             zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
    111          END DO 
    112       END DO 
     107      DO_2D_11_11 
     108         zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
     109      END_2D 
    113110      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )   
    114111      CALL dom_uniq( zprw, 'V' ) 
    115       DO jj = 1, jpj 
    116          DO ji = 1, jpi 
    117             zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
    118          END DO 
    119       END DO 
     112      DO_2D_11_11 
     113         zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
     114      END_2D 
    120115      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )   
    121116!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil 
     
    155150       
    156151      ! note that mbkt is set to 1 over land ==> use surface tmask 
    157       zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 
     152      zprt(:,:) = REAL( mbkt(:,:) , wp ) 
    158153      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points 
    159       zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 
     154      zprt(:,:) = REAL( mikt(:,:) , wp ) 
    160155      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points 
    161       zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 
    162       CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points 
    163156      !                                                         ! vertical mesh 
    164157      CALL iom_rstput( 0, 0, inum, 'e3t_1d', e3t_1d, ktype = jp_r8  )    !    ! scale factors 
  • NEMO/trunk/src/OCE/DOM/domzgr.F90

    r10425 r12377  
    4444 
    4545  !! * Substitutions 
    46 #  include "vectopt_loop_substitute.h90" 
     46#  include "do_loop_substitute.h90" 
    4747   !!---------------------------------------------------------------------- 
    4848   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7171      INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices 
    7272      ! 
    73       INTEGER  ::   jk                  ! dummy loop index 
     73      INTEGER  ::   ji,jj,jk            ! dummy loop index 
     74      INTEGER  ::   ikt, ikb            ! top/bot index 
    7475      INTEGER  ::   ioptio, ibat, ios   ! local integer 
    7576      REAL(wp) ::   zrefdep             ! depth of the reference level (~10m) 
     
    118119      ! Any closed seas (defined by closea_mask > 0 in domain_cfg file) to be filled  
    119120      ! in at runtime if ln_closea=.false. 
    120       IF( .NOT.ln_closea )   CALL clo_bat( k_top, k_bot ) 
     121      IF( ln_closea ) THEN 
     122         IF ( ln_maskcs ) THEN 
     123            ! mask all the closed sea 
     124            CALL clo_msk( k_top, k_bot, mask_opnsea, 'mask_opensea' ) 
     125         ELSE IF ( ln_mask_csundef ) THEN 
     126            ! defined closed sea are kept 
     127            ! mask all the undefined closed sea 
     128            CALL clo_msk( k_top, k_bot, mask_csundef, 'mask_csundef' ) 
     129         END IF 
     130      END IF 
    121131      ! 
    122132      IF(lwp) THEN                     ! Control print 
     
    138148      !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 
    139149      CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1 
    140        
    141  
     150      ! 
     151      !                                ! ice shelf draft and bathymetry 
     152      DO_2D_11_11 
     153         ikt = mikt(ji,jj) 
     154         ikb = mbkt(ji,jj) 
     155         bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1) 
     156         risfdep(ji,jj) = gdepw_0(ji,jj,ikt  ) 
     157      END_2D 
     158      ! 
    142159      !                                ! deepest/shallowest W level Above/Below ~10m 
    143160!!gm BUG in s-coordinate this does not work! 
     
    296313      !                                    ! N.B.  top     k-index of W-level = mikt 
    297314      !                                    !       bottom  k-index of W-level = mbkt+1 
    298       DO jj = 1, jpjm1 
    299          DO ji = 1, jpim1 
    300             miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
    301             mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
    302             mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
    303             ! 
    304             mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
    305             mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
    306          END DO 
    307       END DO 
     315      DO_2D_10_10 
     316         miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
     317         mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
     318         mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
     319         ! 
     320         mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
     321         mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
     322      END_2D 
    308323      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
    309324      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
  • NEMO/trunk/src/OCE/DOM/dtatsd.F90

    r11536 r12377  
    3535   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
    3636 
     37   !! * Substitutions 
     38#  include "do_loop_substitute.h90" 
    3739   !!---------------------------------------------------------------------- 
    3840   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6567      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
    6668      ! 
    67       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
    6869      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
    6970901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' ) 
    70       REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
    7171      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
    7272902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' ) 
     
    186186         ENDIF 
    187187         ! 
    188          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    189             DO ji = 1, jpi 
    190                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
    191                   zl = gdept_0(ji,jj,jk) 
    192                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    193                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    194                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    195                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    196                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    197                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    198                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    199                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    200                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    201                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    202                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    203                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
    204                         ENDIF 
    205                      END DO 
    206                   ENDIF 
    207                END DO 
    208                DO jk = 1, jpkm1 
    209                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    210                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    211                END DO 
    212                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    213                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
     188         DO_2D_11_11 
     189            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     190               zl = gdept_0(ji,jj,jk) 
     191               IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
     192                  ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
     193                  zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
     194               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
     195                  ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
     196                  zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
     197               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
     198                  DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
     199                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
     200                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
     201                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
     202                        zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     203                     ENDIF 
     204                  END DO 
     205               ENDIF 
    214206            END DO 
    215          END DO 
     207            DO jk = 1, jpkm1 
     208               ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
     209               ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
     210            END DO 
     211            ptsd(ji,jj,jpk,jp_tem) = 0._wp 
     212            ptsd(ji,jj,jpk,jp_sal) = 0._wp 
     213         END_2D 
    216214         !  
    217215      ELSE                                !==   z- or zps- coordinate   ==! 
     
    221219         ! 
    222220         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
    223             DO jj = 1, jpj 
    224                DO ji = 1, jpi 
    225                   ik = mbkt(ji,jj)  
    226                   IF( ik > 1 ) THEN 
    227                      zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    228                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 
    229                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
    230                   ENDIF 
    231                   ik = mikt(ji,jj) 
    232                   IF( ik > 1 ) THEN 
    233                      zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
    234                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
    235                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
    236                   END IF 
    237                END DO 
    238             END DO 
     221            DO_2D_11_11 
     222               ik = mbkt(ji,jj)  
     223               IF( ik > 1 ) THEN 
     224                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     225                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 
     226                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
     227               ENDIF 
     228               ik = mikt(ji,jj) 
     229               IF( ik > 1 ) THEN 
     230                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
     231                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
     232                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
     233               END IF 
     234            END_2D 
    239235         ENDIF 
    240236         ! 
  • NEMO/trunk/src/OCE/DOM/istate.F90

    r10499 r12377  
    2828   USE dtauvd         ! data: U & V current             (dta_uvd routine) 
    2929   USE domvvl          ! varying vertical mesh 
    30    USE iscplrst        ! ice sheet coupling 
    3130   USE wet_dry         ! wetting and drying (needed for wad_istate) 
    3231   USE usrdef_istate   ! User defined initial state 
     
    4342 
    4443   !! * Substitutions 
    45 #  include "vectopt_loop_substitute.h90" 
     44#  include "do_loop_substitute.h90" 
    4645   !!---------------------------------------------------------------------- 
    4746   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    5150CONTAINS 
    5251 
    53    SUBROUTINE istate_init 
     52   SUBROUTINE istate_init( Kbb, Kmm, Kaa ) 
    5453      !!---------------------------------------------------------------------- 
    5554      !!                   ***  ROUTINE istate_init  *** 
     
    5756      !! ** Purpose :   Initialization of the dynamics and tracer fields. 
    5857      !!---------------------------------------------------------------------- 
     58      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices 
     59      ! 
    5960      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    6061!!gm see comment further down 
     
    7677      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    7778      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    78       tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk 
     79      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
    7980      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    8081#if defined key_agrif 
    81       ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization 
    82       va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization     
     82      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization 
     83      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization     
    8384#endif 
    8485 
    8586      IF( ln_rstart ) THEN                    ! Restart from a file 
    8687         !                                    ! ------------------- 
    87          CALL rst_read                        ! Read the restart file 
    88          IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry 
     88         CALL rst_read( Kbb, Kmm )            ! Read the restart file 
    8989         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    9090         ! 
     
    9797         ! 
    9898         IF( ln_tsd_init ) THEN                
    99             CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000 
     99            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000 
    100100            ! 
    101             sshb(:,:)   = 0._wp               ! set the ocean at rest 
     101            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest 
    102102            IF( ll_wd ) THEN 
    103                sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
     103               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
    104104               ! 
    105105               ! Apply minimum wetdepth criterion 
    106106               ! 
    107                DO jj = 1,jpj 
    108                   DO ji = 1,jpi 
    109                      IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN 
    110                         sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
    111                      ENDIF 
    112                   END DO 
    113                END DO  
     107               DO_2D_11_11 
     108                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
     109                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
     110                  ENDIF 
     111               END_2D 
    114112            ENDIF  
    115             ub  (:,:,:) = 0._wp 
    116             vb  (:,:,:) = 0._wp   
     113            uu  (:,:,:,Kbb) = 0._wp 
     114            vv  (:,:,:,Kbb) = 0._wp   
    117115            ! 
    118116         ELSE                                 ! user defined initial T and S 
    119             CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )          
     117            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
    120118         ENDIF 
    121          tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    122          sshn (:,:)     = sshb(:,:)    
    123          un   (:,:,:)   = ub  (:,:,:) 
    124          vn   (:,:,:)   = vb  (:,:,:) 
    125          hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
    126          CALL div_hor( 0 )                    ! compute interior hdivn value   
    127 !!gm                                    hdivn(:,:,:) = 0._wp 
     119         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     120         ssh (:,:,Kmm)     = ssh(:,:,Kbb)    
     121         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     122         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     123         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
     124         CALL div_hor( 0, Kbb, Kmm )         ! compute interior hdiv value   
     125!!gm                                    hdiv(:,:,:) = 0._wp 
    128126 
    129127!!gm POTENTIAL BUG : 
    130 !!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
     128!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
    131129!!             as well as gdept and gdepw....   !!!!!  
    132130!!      ===>>>>   probably a call to domvvl initialisation here.... 
     
    138136!            ALLOCATE( zuvd(jpi,jpj,jpk,2) ) 
    139137!            CALL dta_uvd( nit000, zuvd ) 
    140 !            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:) 
    141 !            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:) 
     138!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb) 
     139!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb) 
    142140!            DEALLOCATE( zuvd ) 
    143141!         ENDIF 
    144142         ! 
    145143!!gm This is to be changed !!!! 
    146 !         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 
     144!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here 
    147145!         IF( .NOT.ln_linssh ) THEN 
    148146!            DO jk = 1, jpk 
    149 !               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
     147!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
    150148!            END DO 
    151149!         ENDIF 
     
    157155      ! Do it whatever the free surface method, these arrays being eventually used 
    158156      ! 
    159       un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp 
    160       ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
     157      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp 
     158      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp 
    161159      ! 
    162 !!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked 
    163       DO jk = 1, jpkm1 
    164          DO jj = 1, jpj 
    165             DO ji = 1, jpi 
    166                un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 
    167                vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
    168                ! 
    169                ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 
    170                vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 
    171             END DO 
    172          END DO 
    173       END DO 
     160!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked 
     161      DO_3D_11_11( 1, jpkm1 ) 
     162         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
     163         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
     164         ! 
     165         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
     166         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
     167      END_3D 
    174168      ! 
    175       un_b(:,:) = un_b(:,:) * r1_hu_n(:,:) 
    176       vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:) 
     169      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 
     170      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 
    177171      ! 
    178       ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 
    179       vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 
     172      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb) 
     173      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb) 
    180174      ! 
    181175   END SUBROUTINE istate_init 
Note: See TracChangeset for help on using the changeset viewer.